Conceptual DFT Study of the Local Chemical Reactivity of the Colored BISARG Melanoidin and Its Protonated Derivative

Frau, Juan; Glossman-Mitnik, Daniel

doi:10.3389/fchem.2018.00136

ORIGINAL RESEARCH article

Front. Chem., 01 May 2018

Sec. Theoretical and Computational Chemistry

Volume 6 - 2018 | https://doi.org/10.3389/fchem.2018.00136

This article is part of the Research Topic Quantum Mechanical/Molecular Mechanical Approaches for the Investigation of Chemical Systems – Recent Developments and Advanced Applications View all 15 articles

Conceptual DFT Study of the Local Chemical Reactivity of the Colored BISARG Melanoidin and Its Protonated Derivative

$\r\nJuan Frau$ Juan Frau¹

Daniel Glossman-Mitnik^1,2^*

¹Departament de Química, Universitat de les Illes Balears, Palma de Mallorca, Spain
²Laboratorio Virtual NANOCOSMOS, Centro de Investigación en Materiales Avanzados, Departamento de Medio Ambiente y Energía, Chihuahua, Mexico

This computational study assessed eight fixed RSH (range-separated hybrid) density functionals that include CAM-B3LYP, LC-ωPBE, M11, MN12SX, N12SX, ωB97, ωB97X, and ωB97XD related to the Def2TZVP basis sets together with the SMD solvation model in the calculation the molecular structure and reactivity properties of the BISARG intermediate melanoidin pigment (5-(2-(E)-(Z)-5-[(2-furyl)methylidene]-3-(4-acetylamino-4-carboxybutyl)-2-imino-1,3-dihydroimidazol-4-ylideneamino(E)-4-[(2-furyl)methylidene]-5-oxo-1H-imidazol-1-yl)-2-acetylaminovaleric acid) and its protonated derivative, BISARG(p). The chemical reactivity descriptors for the systems were calculated via the Conceptual Density Functional Theory. The choice of active sites applicable to nucleophilic, electrophilic as well as radical attacks were made by linking them with Fukui functions indices, electrophilic and nucleophilic Parr functions, and the condensed Dual Descriptor Δf(r). The study found the MN12SX and N12SX density functionals to be the most appropriate in predicting the chemical reactivity of the molecular systems under study starting from the knowledge of the HOMO, LUMO, and HOMO-LUMO gap energies.

1. Introduction

Visual color in processed foods is largely due to colored products of Maillard or nonenzymic browning reactions. In spite of the longstanding aesthetic and practical interest in Maillard derived food coloring materials, relatively little is known about the chemical structures responsible for visual color (Rizzi, 1997). These chemical structures are known as Colored Maillard Reaction Products and can be isolated at intermediate stages during the melanoidin formation process.

Besides their interest as dye molecules which may be useful as food additives, but also as dyes for dye-sensitized solar cells (DSSC), these compounds have also antioxidant capabilities. Thus, they are amenable to be studied by analyzing their molecular reactivity properties.

One of these isolated molecules is named by the acronym BISARG and together with its protonated derivative, BISARG(p) have been experimentally studied as a part of a work related to the formation of melanidins (Hofmann, 1998) and we believe that it could be of interest to study their molecular reactivity by using the ideas of Conceptual DFT, in the same way of our previous works (Alvarado-González et al., 2013; Cervantes-Navarro and Glossman-Mitnik, 2013; Glossman-Mitnik, 2013a,b; Martínez-Araya et al., 2013a,b; Salgado-Morán et al., 2013; Glossman-Mitnik, 2014a,b,c,d; Martínez-Araya and Glossman-Mitnik, 2015; Martínez-Araya et al., 2015; Soto-Rojo et al., 2015; Frau et al., 2016a,b,c; Mendoza-Huízar et al., 2016; Frau et al., 2017a,b,c,d,e; Frau and Glossman-Mitnik, 2017a,b,c,d,e,f,g; Sastre et al., 2017).

The interest in using range-separated (RS) exchange correlation functionals in KS DFT is on the rise (Gledhill et al., 2016). The functionals tend to partition the $r_{12}^{- 1}$ operator and exchange them into long- and short-range parts, whose range separation parameter ω controls the rate of attaining the long-range behavior. It is possible to fix the value of ω. The value can also be nonempirically “tuned” through a system-by-system mechanism that minimizes some tuning norms. The basis of the optimal tuning approach is the knowledge that the energy that the HOMO should have, ϵ_H(N), in exact KS as well as generalized KS theory for an N electron system, ought to be exactly −IP(N). Hence, IP represents the vertical ionization potential that is calculated by considering a particular functional energy difference E(N-1) − E(N). If approximate functionals are used, it is possible to have considerable differences between ϵ_H(N), and −IP(N). Optimal tuning constitutes determining a system-specific range-separation parameter ω non-empirically in an RSE functional. Optionally, it also implied that several other parameters including ϵ_H(N) = −IP(N) are satisfied optimally (Jacquemin et al., 2014). Even though no equivalency exists to match this prescription of electron affinity (EA) coupled with LUMO in the case of neutral species, it is possible to say that ϵ_H(N+1) = −EA(N), that is, the electron affinity of the neutral system is equal to minus the HOMO energy of the anion (SOMO), which facilitates the finding of an optimized value of ω, and is then optimized to establish both properties simultaneously. Some concerns have been raised during the preparation of this paper regarding the validity of the ionization potential theorem (IP) within the context of Generalized Kohn-Sham (GKS) theory. However, it must stressed that Baer et al. (2010) and more recently Baerends et al. (2013) and Karolewski et al. (2013) have given arguments that the same criterion applies in GKS theories and with with hybrid and range-separated hybrid functionals. This will make it easy to predict the Conceptual DFT descriptors. In the past, the simultaneous prescription has been referred to as the “KID procedure” (for Koopmans in DFT), courtesy of the analogy it shares with the Koopmans' theorem within the Hartree-Fock theory. This SOMO energy will not be, in general, equal to the LUMO of the neutral, but if the difference between them, which we have called ΔSL, is small enough to be considered negligible for predictions of the Conceptual DFT descriptors, then the practical KID procedure will have a computational support.

This implies that the appropriateness of a particular density functional in making predictions of the Conceptual DFT descriptors directly by relying on the properties that the neutral molecule can be easily estimated. It only requires one to check the way that it has followed the KID procedure. Nevertheless, tune-optimization depends on the system and must be performed for each molecule one at a time. Therefore, examining the various density functionals exhibiting significant accuracy across various types of databases in physics, chemistry, and where the ω value is fixed will determine how they perform the practical technique.

Thus, in this computational study we will assess eight density functionals in calculating the molecular properties and structure of the BISARG intermediate melanoidin pigment and its protonated derivative, BISARG(p). Following the same ideas of previous works, we will consider fixed RSH functional instead of the optimally-tuned RSH density functionals that have attained great success and have also supported the validity of the IP theorem in the context of the GKS theory (Stein et al., 2009a,b; Karolewski et al., 2011; Kuritz et al., 2011; Refaely-Abramson et al., 2011; Foster and Wong, 2012; Koppen et al., 2012; Kronik et al., 2012; Phillips et al., 2012a,b; Karolewski et al., 2013; Moore and Autschbach, 2013; Egger et al., 2014; Foster et al., 2014; Jacquemin et al., 2014; Niskanen and Hukka, 2014; Sun and Autschbach, 2014; Manna et al., 2015; Lima et al., 2016; Pereira et al., 2017).

2. Theoretical Background

The theoretical background of this study is similar to the previous conducted research presented complete purposes, because this research is a component of a major project that it is in progress. If we consider the KID procedure mentioned in the Introduction together with a finite difference approximation, then the global reactivity descriptors can be written as:

where ϵ_H and ϵ_L are the energies of the highest occupied and the lowest unoccupied molecular orbitals (HOMO and LUMO), respectively.

Applying the same ideas, the definitions for the local reactivity descriptors are:

where ρ_N+1(r), ρ_N(r), and ρ_N−1(r) are the electronic densities at point r for the system with N + 1, N, and N − 1 electrons, respectively, and $ρ_{s}^{r c} (r)$ and $ρ_{s}^{r a}$ (r) are related to the atomic spin density (ASD) at the r atom of the radical cation or anion of a given molecule, respectively (Domingo et al., 2016).

3. Settings and Computational Methods

Following the lines of our previous works, the computational studies were performed with the Gaussian 09 (Frisch et al., 2018) series of programs with density functional methods as implemented in the computational package. The basis set used in this work was Def2SVP for geometry optimization and frequencies, while Def2TZVP was considered for the calculation of the electronic properties (Weigend and Ahlrichs, 2005; Weigend, 2006). All the calculations were performed in the presence of water as the solvent by doing Integral Equation Formalism-Polarized Continuum Model (IEF-PCM) computations according to the SMD solvation model (Marenich et al., 2009).

For the calculation of the molecular structure and properties of the studied systems, we have chosen eight density functionals which are known to consistently provide satisfactory results for several structural and thermodynamic properties:

In these functionals, GGA stands for generalized gradient approximation (in which the density functional depends on the up and down spin densities and their reduced gradient) and NGA stands for nonseparable gradient approximation (in which the density functional depends on the up/down spin densities and their reduced gradient, and also adopts a nonseparable form).

4. Results and Discussion

The three-dimentional molecular structure of the BISARG system was built with the aid of molecular graphics program starting from structure presented in the original article (Hofmann, 1998). Starting from this, the molecular structure of its protonated derivative, BISARG(p) was built with the aid of a chemical visualization software. The pre-optimization of the systems was done using random sampling that involved molecular mechanics techniques and inclusion of the various torsional angles via the general MMFF94 force field (Halgren, 1996a,b,c, 1999; Halgren and Nachbar, 1996) through the Marvin View 17.15 program that constitutes an advanced chemical viewer suited to multiple and single chemical queries, structures and reactions (https://www.chemaxon.com). Afterwards, the structures that the resultant lower-energy conformers assumed for both molecules were reoptimized using the eight density functionals mentioned in the previous section together with the Def2SVP basis set as well as the SMD solvation model using water as the solvent.

The analysis of the results obtained in the study aimed at verifying that the KID procedure was fulfilled. On doing it previously, several descriptors associated with the results that HOMO and LUMO calculations obtained are related with results obtained using the vertical I and A following the ΔSCF procedure. A link exists between the three main descriptors and the simplest conformity to the Koopmans' theorem by linking ϵ_H with -I, ϵ_L with -A, and their behavior in describing the HOMO-LUMO gap as J_I = |ϵ_H+E_gs(N − 1)−E_gs(N)|, J_A = |ϵ_L+E_gs(N)−E_gs(N + 1)|, and $J_{H L} = \sqrt{{J_{I}}^{2} + {J_{A}}^{2}}$ . Notably, the J_A descriptor consists of an approximation that remains valid only when the HOMO that a radical anion has (the SOMO) shares similarity with the LUMO that the neutral system has. Consequently, we decided to design another descriptor ΔSL (the difference between the SOMO and LUMO energies), to guide in verifying how the approximation is accurate.

The results of the calculation of the electronic energies of the neutral, positive and negative molecular systems (in au) of BISARG and BISARG(p), the HOMO, LUMO, and SOMO orbital energies (in eV), J_I, J_A, J_HL, and ΔSL descriptors calculated with the eight density functionals and the Def2TZVP basis set using water as solvent simulated with the SMD parametrization of the IEF-PCM model are presented in Tables 1, 2.

TABLE 1

Table 1. Electronic energies of the neutral, positive and negative molecular systems (in au) of the BISARG molecule, the HOMO, LUMO, and SOMO orbital energies (in eV), J_I, J_A, J_HL, and ΔSL descriptors (also in eV) calculated with the eight RSH density functionals and the Def2TZVP basis set using water as solvent simulated with the SMD parametrization of the IEF-PCM model.

TABLE 2

Table 2. Electronic energies of the neutral, positive and negative molecular systems (in au) of the protonated BISARG(p) molecule, the HOMO, LUMO, and SOMO orbital energies (in eV), J_I, J_A, J_HL, and ΔSL descriptors (also in eV) calculated with the eight RSH density functionals and the Def2TZVP basis set using water as solvent simulated with the SMD parametrization of the IEF-PCM model.

As presented in previous works, we considered four other descriptors that analyze how well the studied density functionals are useful for the prediction of the electronegativity χ, the global hardness η, and the global electrophilicity ω, and for a combination of these Conceptual DFT descriptors, considering only the energies of the HOMO and LUMO or the vertical I and A: J_χ = |χ − χ_K|, J_η = |η − η_K|, J_ω = |ω − ω_K|, and $J_{C D F T} = \sqrt{J_{χ}^{2} + J_{η}^{2} + J_{ω}^{2}}$ , where CDFT stands for Conceptual DFT. The underscript K stands for the descriptor calculated by applying the KID procedure.The results of the calculations of J_χ, J_η, J_ω, and J_CDFT for the low-energy conformers of BISARG and BISARG(p) in water are displayed in Tables 3, 4, respectively.

TABLE 3

Table 3. J_χ, J_η, J_ω, and J_CDFT (in eV) of the BISARG intermediate melanoidin pigment.

TABLE 4

Table 4. J_χ, J_η, J_ω, and J_CDFT (in eV) of the protonated BISARG(p) intermediate melanoidin pigment.

As Tables 1–4 provide, the KID procedure applies accurately from MN12SX and N12SX density functionals that are range-separated hybrid meta-NGA as well as range-separated hybrid NGA density functionals respectively. In fact, the values of J_I, J_A, and J_HL are actually not zero. Nevertheless, the results tend to be impressive especially for the MN12SX density functional. As well, the ΔSL descriptor reaches the minimum values when MN12SX and N12SX density functionals are used in the calculations. This implies that there are sufficient justifications to assume that the LUMO of the neutral approximates the electron affinity. The same density functionals follow the KID procedure in the rest of the descriptors such as J_χ, J_η, J_ω, and J_CDFT.

Having verified that the MN12SX/Def2TZVP model chemistry is a good choice for the calculation of the global reactivity descriptors, we now present the optimized molecular structures of BISARG and BISARG(p) in water in Supplementary Figures 1, 2. Meanwhile, the calculated bond lengths and bond angles for both cases are shown in Supplementary Tables 1–4.

As a summary of the previous results, the global reactivity descriptors for the BISARG and BISARG(p) molecules calculated with the MN12SX/Def2TZVP model chemistry in water are presented in Table 5.

TABLE 5

Table 5. Global reactivity descriptors for the BISARG intermediate melanoidin pigment and its protonated derivative BISARGd(p) calculated with the MN12SX density functional.

The calculations of the condensed Fukui functions and dual descriptor are done by using the Chemcraft molecular analysis program to extract the Mulliken and NPA atomic charges (Zhurko and Zhurko, 2012) beginning with single-point energy calculations involving the MN12SX density functional that uses the Def2TZVP basis set in line with the SMD solvation model, and water utilized as the solvent.

Considering the potential application the studied molecules as antioxidants, it is of interest to get insight into the active sites for radical attack. Graphical representations of the radical Fukui function f⁰ calculated with the MN12SX/Def2TZVP model chemistry for both systems in water are presented in Figures 1, 2.

FIGURE 1

Figure 1. A graphical schematic representation of the radical Fukui function f⁰ over the atomic sites of the BISARG intermediate melanoidin pigment.

FIGURE 2

Figure 2. A graphical schematic representation of the radical Fukui function f⁰ over the atomic sites of the protonated BISARG(p) intermediate melanoidin pigment.

The condensed electrophilic and nucleophilic Parr functions $P_{k}^{+}$ and $P_{k}^{-}$ over the atoms of the BISARG and BISARG(p) molecules in water have been calculated by extracting the Mulliken and Hirshfeld (or CM5) atomic charges using the Chemcraft molecular analysis program (Zhurko and Zhurko, 2012) starting from single-point energy calculations of the ionic species with the MN12SX density functional using the Def2TZVP basis set in the presence of the solvents according to the SMD solvation model.

The results for the condensed dual descriptor calculated with Mulliken atomic charges Δf_k (M), with NPA atomic charges Δf_k (N), the electrophilic and nucleophilic Parr functions with Mulliken atomic spin densities $P_{k}^{+}$ (M) and $P_{k}^{-}$ (M), and the electrophilic and nucleophilic Parr functions with Hirshfeld (or CM5) atomic spin densities $P_{k}^{+}$ (H) and $P_{k}^{-}$ (H) are displayed in Tables 6, 7 for the BISARG and BISARG(p) molecules in water, respectively, while Figures 3, 4 show schematic representations of the molecules with the numbering of the most important reactive sites according to the results in Tables 6, 7.

TABLE 6

Table 6. The condensed dual descriptor calculated with Mulliken atomic charges Δf_k (M), and with NPA atomic charges Δf_k (N), the electrophilic and nucleophilic Parr functions with Mulliken atomic spin densities $P_{k}^{+}$ (M) and $P_{k}^{-}$ (M), and the electrophilic and nucleophilic Parr functions with Hirshfeld (or CM5) atomic spin densities $P_{k}^{+}$ (H) and $P_{k}^{-}$ (H) for the BISARG melanoidin molecule. Hydrogens and atomic sites where the absolute value of the dual descriptor is lower than 1 are not shown.

TABLE 7

Table 7. The condensed dual descriptor calculated with Mulliken atomic charges Δf_k (M), and with NPA atomic charges Δf_k (N), the electrophilic and nucleophilic Parr functions with Mulliken atomic spin densities $P_{k}^{+}$ (M) and $P_{k}^{-}$ (M), and the electrophilic and nucleophilic Parr functions with Hirshfeld (or CM5) atomic spin densities $P_{k}^{+}$ (H) and $P_{k}^{-}$ (H) for the BISARG(p) melanoidin molecule. Hydrogens and atomic sites where the absolute value of the dual descriptor is lower than 1 are not shown.

FIGURE 3

Figure 3. A graphical schematic representation of the BISARG intermediate melanoidin pigment with the numbering of the most important reactive sites.

FIGURE 4

Figure 4. A graphical schematic representation of the protonated BISARG(d) intermediate melanoidin pigment with the numbering of the most important reactive sites.

From the results for the local reactivity descriptors in Table 6, it can be concluded that C2 and C23 will be the preferred sites for a nucleophilic attack and that these atoms will act as electrophilic species in a chemical reaction. In turn, it can be appreciated that N8 will be prone to electrophilic attacks and that this atomic site will act as a nucleophilic species in chemical reactions that involve the BISARG molecule in water. In turn, for the case of the BISARG(p) molecule in water, C2 and C17 will be the preferred sites for a nucleophilic attack while C15 and C29 will be the sites for electrophilic reactions.

5. Conclusions

Eight fixed RSH density functionals, including CAM-B3LYP, LC-ωPBE, M11, MN12SX, N12SX, ωB97, ωB97X, and ωB97XD, were examined to determine whether they fulfill the empirical KID procedure so as to provide computational support for this common practice. The assessment was conducted by comparing the values from HOMO and LUMO calculations to those generated by the ΔSCF technique for the BISARG molecule and its protonated derivative, BISARG(p). BISARG and BISARG(p) are intermediate melanoidin pigments that are of academic and industrial interest. The study has observed that the range-separated and hybrid meta-NGA density functionals tend to be the most suited in meeting this goal. Thus, they can be suitable alternatives to density functionals where the behavior of them is optimally tuned using a gap-fitting procedure. They also exhibit the desirable prospect of benefiting future studies aimed at understanding the chemical reactivity of colored melanoidins with larger molecular weights when reducing sugars react with proteins and peptides.

It is not the goal of Computational Chemistry to perform studies to reproduce known experimental results except in the case that they can be used for the calibration of a particular technique. Instead, it can be useful to predict in advance the structural and chemical reactivity characteristics of new or unknown molecular systems whose properties have not been reported and as guide for future research. As far as we know, there are no reports in the literature about the chemical reactivity properties for the molecular systems considered in this work and it is not possible to perform any kind of comparisons. However, the present study shows that with an adequate choice of the model chemistry we have been able to predict the sites of interaction of the BISARG and BISARG(p) molecules with impressive accuracy starting from the knowledge of the HOMO, LUMO, and HOMO-LUMO gap energies of the studied systems. This involves having DFT-based reactivity descriptors, including Fukui functions, Parr functions, and Dual Descriptor calculations. In conclusion, the Conceptual DFT descriptors are useful in characterizing and describing the preferred reactive sites and in comprehensively explaining the reactivity of the molecules.

Author Contributions

DG-M conceived and designed the research and headed, wrote and revised the manuscript, while JF contributed to the writing and the revision of the article.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

This work has been partially supported by CIMAV, SC, and Consejo Nacional de Ciencia y Tecnología (CONACYT, Mexico) through Grant 219566-2014 for Basic Science Research. DG-M conducted this work while a Visiting Lecturer at the University of the Balearic Islands from which support is gratefully acknowledged. This work was cofunded by the Ministerio de Economía y Competitividad (MINECO) and the European Fund for Regional Development (FEDER) (CTQ2014-55835-R).

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fchem.2018.00136/full#supplementary-material

References

Alvarado-González, M., Flores-Holguín, N., and Glossman-Mitnik, D. (2013). Computational nanochemistry study of the molecular structure and properties of the chlorophyll a molecule. Int. J. Photoener. 2013:8. doi: 10.1155/2013/424620