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

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.


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.
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 −1 12 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, tuneoptimization 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.

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: Parr and Yang (1984) Global Hardness Parr and Yang (1984) Electrophilicty Parr et al. (1999) Electrodonating Power ω − = (3I+A) 2 16(I−A) ≈ (3ǫ H +ǫ L ) 2 16η Gázquez et al. (2007) Electroaccepting Power ω + = (I+3A) 2 16(I−A) ≈ (ǫ H +3ǫ L ) 2 16η Gázquez et al. (2007) Net electrophilicity et al. (2009) 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: Parr and Yang (1984) Electrophilic Fukui Function f − (r) = ρ N (r) − ρ N−1 (r) Parr and Yang (1984) Dual Descriptor Morell et al. (2005Morell et al. ( , 2006 Nucleophilic Parr Function P − (r) = ρ rc s (r) Domingo et al. (2013) Electrophilic Parr Function P + (r) = ρ ra s (r) Domingo et al. (2013) 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 ρ rc s (r) and ρ ra s (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).

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: Long-range corrected density functional Chai and Head-Gordon (2008b) ωB97XD ωB97X version including empirical dispersion Chai and Head-Gordon (2008a) 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).

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(Halgren, ,b,c, 1999Halgren 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 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.
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: 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.
As Tables 1-4 provide, the KID procedure applies accurately from MN12SX and N12SX density functionals that are rangeseparated 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 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.   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.
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  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 | 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. 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 , 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. 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.

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 rangeseparated 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.

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).