Density Functional Study of Aflatoxin B1, B2, G1 and G2: NBO and AIM Analyses

Computational Details

To achieve the optimal structure of the aflatoxin composition, it is crucial to obtain the most stable structure with a global minimum. While Monte Carlo methods and molecular dynamics are commonly employed to search for global minima, there is no guarantee that these methods will always lead to the global minimum. Therefore, it is advantageous to utilize various algorithms and pre- optimization techniques with other software [10, 11, 12].

The initial structure of the aflatoxins under investigation was obtained from x-ray crystallography data. Subsequently, the most stable conformers were optimized using the B3LYP (Becke–Lee–Yang–Parr) version of DFT method with the 6-311G (d, p) basis set. All calculations were performed using the Gaussian 03 program [13]. The chosen basis set contains an appropriate number of basis set functions to accurately reproduce experimental data [14, 15]. The absence of negative frequencies in the normal mode analysis confirms that the optimized structures correspond to real minimum stationary points on the potential energy surfaces.

Visualizations of the optimized structures were generated using Gauss View 5.0 [16].

For further analysis, the Natural Bond Orbital (NBO) version 3.1 [17] was employed to perform NBO analysis. Additionally, the AIM2000 software [18, 19] was utilized for Quantum Theory of Atoms in Molecules (QTAIM) analysis. These analyses provide valuable insights into the electronic and intramolecular bonding properties of the studied aflatoxins.

Energies and Geometries

The optimization of the aflatoxin structures provides valuable information regarding their three-dimensional arrangement. Figure 2 illustrates that the optimized structures of the studied aflatoxins are non-planar, and the orientation of the pentagon ring varies among them.

Table 1 presents the distances between selected atoms and the bond lengths of specific bonds in the aflatoxin molecules. It is evident that in AFB (1, 2), the distance between the oxygen atoms of the carbonyl groups (O1 and O2) is greater compared to AFG (1, 2). Additionally, the bond length between O1 and C1 in AFB (1, 2) is longer than that in AFG (1, 2). This suggests that O1 is more strongly bound to its pentagonal ring, allowing the two oxygen atoms (O1, O2) to be more spatially separated and interact with other molecules more effectively. On the other hand, the bond length between O2 and C2 is less influenced by the adjacent ring.

Furthermore, there is a possibility of forming an intermolecular bond between the O6 and C5 atoms, which will be further discussed in the AIM analysis. Interestingly, the distance between O6 and C5 in AFB (1, 2) is greater than in AFG (1, 2), indicating a weaker intermolecular bond between these atoms in AFB (1, 2).

Thermodynamic Properties

The thermodynamic properties of the studied aflatoxins were computed using the B3LYP/6-311G (d, p) level of theory at the standard conditions of 298.15 K and 1 atmospheric pressure. The calculated thermodynamic properties are summarized in Table 2. Based on the reported results, it is observed that aflatoxin B1, which is known to be more toxic according to previous studies [3], exhibits higher instability compared to the other aflatoxins. Additionally, it is noted that the dipole moment of AFB2 is weaker than that of the other aflatoxins.

*Atomic unit: 1 au = 627.5095 Kcal. mol−1 Table 2: Corrected electron energy values and thermodynamic properties of the studied species at B3LYP/6-311G (d, p) level of theory in the gas phase.

Natural Bond Orbital (NBO) Analysis

The NBO analysis is a method that involves converting the wave function into its constituent parts. This analysis allows for the determination of natural atomic orbitals, natural hybrid orbitals, natural bond orbitals, and natural local molecular orbitals, which are then used for population analysis and NBO energy analysis. Natural bond orbitals are special vectors derived from the density matrix, formed by the first and second blocks. The central blocks describe internal and nonbonding orbitals, while the two-block consists of bonding and nonbonding orbitals [20, 21, 22].

Table 3 presents the electronic properties of the studied aflatoxins, including the lowest unoccupied molecular orbital (LUMO) and its energy (ε LUMO), as well as its occupation number. It also includes the highest occupied molecular orbital (HOMO) and its energy (ε HOMO), along with its occupation number. Other properties such as ionization energy (I), electron affinity (A), chemical hardness( ) η , electronic chemical potential( ) µ , electrophilicity index( ) ω , and softness (S) are also provided.

The electrophilicity index( ) ω , introduced by Parr, et al. [23, 24, 25, 26], and is a descriptor used to quantify the global electrophilic nature of a molecule on a relative scale. Chemical hardness, electronic chemical potential, and softness are universal reactivity descriptors [27, 28]. The stability of a chemical species can be related to its hardness, where soft molecules have a small energy band gap between the LUMO and HOMO, while hard molecules have a larger gap [29, 30].

Based on the results in Table 3, it can be observed that despite the structural and reactivity differences among the aflatoxins, their HOMO and LUMO orbitals have the same characteristics but with different energy levels. Additionally, the band gap for AFG1 and AFG2 (6.5) is smaller than that of AFB1 and AFB2 (6.6). This suggests that AFB1 and AFB2 are harder species compared to AFG1 and AFG2, as indicated by the higher energy band gap.

BD*, 2-Center Antibonding Orbital LP is Valence Lone Pair Table 3: Calculated electronic properties for the studied species at the B3LYP/6-311G (d, p) level of theory in the gas phase.

Figures 3 & 4 depict the HOMO and LUMO orbitals of the studied aflatoxins, illustrating that the HOMO and LUMO orbitals are the same across all aflatoxins.

Table 4 presents the natural charges of the important atoms in aflatoxins B1, B2, G1 and G2, obtained using the B3LYP/6-311G (d, p) level of theory. From Table 4, it can be observed that the partial charge of O1 is greater than that of O2. Additionally, as indicated in Table 3, the HOMO orbital is primarily located on the O1 atom. Therefore, the O1 atom exhibits a higher tendency to react with other molecules based on its higher partial charge and its involvement in the HOMO orbital.

Click to enlarge

Click to enlarge

Atoms in Molecules (AIM) Analysis

In AIM theory, critical points (CPs) are points where the gradient of the electron density becomes zero. There are four categories of critical points, which are determined by the values of the hessian matrix of the electron density. The hessian matrix contains the second derivatives of the electron density with respect to various coordinate axes, and its eigenvalues and eigenvectors are denoted as follows:

Here, i

λ represents the eigenvalues and ui represents

the eigenvectors for the ith critical point. Each critical point

is characterized by its degree and symbol (r, s), where r is

the number of non-zero eigenvalues and s is the sum of the

algebraic signs of the eigenvalues. For example, a bond critical

point (BCP) typically exists between adjacent nuclei and is

represented as (1, 3). This critical point has two negative

$$ \lambda 1 \leq \lambda 2 \leq 0 a $$

and one positive eigenvalue 3

0 λ >

. Table 5 provides a breakdown of critical points based on the

values of r and s.

The trajectory of electron density starting from a BCP and ending at each nucleus is referred to as the pathway. The molecular graph consists of a network of bond paths, which are closely related to the chemical bonding network in the molecule’s equilibrium geometry. The electron density along the bond path at the BCP reaches its lowest value, indicating a strong chemical bond.

The electron density ( ) r ρ is often used to describe the strength of chemical bonds. Higher values of ( ) r ρ correspond to stronger bonding. The Laplacian of the electron density $$ \left(\nabla^ {2} \rho (r)\right) $$

is the sum of the second derivatives of the

electron density along the coordinate axes at a given point

r. The Laplacian values provide information about bond

characteristics. A negative value of

( ) ² r ρ ∇

and high values

of ( ) r ρ

indicate covalent bonding, while a positive value of

( ) ² r ρ ∇

and low values of

( ) r ρ

are associated with closed-

shell interactions, including ionic, coordinate, hydrogen, and

van der Waals bonds. In systems involving hydrogen bond interactions, rare gas dimers, or ion systems, the electron density is approximately one hundredth of an atomic unit, and the Laplacian is positive.

Figure 5 shows the critical points of electron density distribution and the gravitational pathways, where an intermolecular bond between the O6 and C5 atoms is indicated by an arrow. Table 5 summarizes the values of electron density ( ) ( ) r ρ , electron density Laplacian ( ) 2 ( ) r ρ ∇ , total energy ( ) ( ), 2 H r V G ratio, and ellipticity (ε ) at the BCP for the studied aflatoxins.

Click to enlarge

In Table 6, the small values of ρ(r) and the positive sign of the Laplacian indicate that the electron density is concentrated in the interaction region between the O6 and C5 nuclei, suggesting a closed-shell interaction for the O6- C5 molecular bond. Furthermore, by comparing the values of ρ(r) for this bond in different aflatoxin types, it can be concluded that the O6-C5 bond is stronger in aflatoxins AFG (1, 2) compared to the others.

The electron energy density, represented by the Hamiltonian, is negative for collective interactions, with the absolute value indicating the dominance of potential energy in the covalent character of the interaction. The quantities

V(r) and G(r) represent potential energy and kinetic energy,

respectively, and the ratio

2 V

G is always positive. If this ratio

is greater than 1, it indicates a negative Laplacian

$$ \left(\nabla^ {2} \rho (r)\right) $$

and a covalent interaction. If the ratio falls between 0.5 and

1, it implies a positive Laplacian and a negative Hamiltonian

(H(r)), indicating a partial covalence interaction. When the

ratio is less than 0.5, both the Laplacian and Hamiltonian are

positive, suggesting a non-covalent interaction. Therefore,

for the O6-C5 interaction, the ratio falls within the range of

0.5-1, indicating a partial covalence bond.

Regarding the C1-O1, C2-O2, and C3-C4 bonds, they exhibit covalent couplings, with the C3-C4 bond being weaker due to its lower ρ(r) value. The ellipticity values (ε ) indicate the proximity of curves perpendicular to the bond path, with a higher ε value suggesting a higher percentage of π-type bond participation. In the case of the O6-C5 interaction, the ε value is greater than for other bonds, indicating a higher degree of π-type bond character and participation [31, 32].

Aflatoxin BCP ( )* r ρ ( )* ² r ρ ∇ ( )* H r ε AFB1

AFB2

AFG1

AFG2

2 V G C1-O1 0.4135 -0.0585 0.6958 0.0541 4.0429 C2-O2 0.4313 -0.1215 0.743 0.1203 4.0853 C3-C4 0.3499 -1.0442 0.4131 0.4254 7.4337 O6-C5 0.0123 0.0549 -0.0021 1.8102 1.6429 C1-O1 0.4123 -0.0672 0.6931 0.0539 4.0497 C2-O2 0.4308 -0.1205 0.7416 0.1194 4.0847 C3-C4 0.2473 -0.5687 0.1998 0.0279 8.9293 O6-C5 0.0125 0.0559 -0.0021 1.9248 1.6441 C1-O1 0.4258 -0.1422 0.7297 0.1133 4.1024 C2-O2 0.4311 -0.1193 0.7429 0.1202 4.0837 C3-C4 0.3479 -1.0311 0.4082 0.4257 7.4262 O6-C5 0.0188 0.0759 -0.0027 0.3701 1.6674 C1-O1 0.4257 -0.1432 0.7295 0.1132 4.1032 C2-O2 0.4309 -0.1213 0.7423 0.1199 4.0852 C3-C4 0.2473 -0.5686 0.1998 0.0279 8.9297 O6-C5 0.0189 0.0763 -0.0027 0.3706 1.6673