Molecular Reactivity and Absorption Properties of Melanoidin Blue-G1 through Conceptual DFT

This computational study presents the assessment of eleven density functionals that include CAM-B3LYP, LC-wPBE, M11, M11L, MN12L, MN12SX, N12, N12SX, wB97, wB97X and wB97XD related to the Def2TZVP basis sets together with the Solvation Model Density (SMD) solvation model in calculating the molecular properties and structure of the Blue-G1 intermediate melanoidin pigment. The chemical reactivity descriptors for the system are calculated via the conceptual Density Functional Theory (DFT). The choice of the active sites related to the nucleophilic, electrophilic, as well as radical attacks is made by linking them with the Fukui function indices, the electrophilic Parr functions and the condensed dual descriptor Δf(r). The prediction of the maximum absorption wavelength tends to be considerably accurate relative to its experimental value. The study found the MN12SX and N12SX density functionals to be the most appropriate density functionals in predicting the chemical reactivity of the studied molecule.


Introduction
In food science, the Maillard reaction is well known (named after the French Chemist who first described it): it consists of the reaction produced between a reducing sugar and an amino acid (the link of a small number of amino acids forms peptides, and at a higher number, they form proteins). When a reducing sugar reacts with a protein, a compound called a Schiff base is formed, and the accumulation of several of these compounds degrades proteins until they create Advanced Glycation End-products (AGEs) that originate fibrils that accumulate in the brain. It is hypothesized that when brain proteins degrade, they may cause degenerative diseases such as Alzheimer's, Parkinson's or diabetes.
The completion of the Maillard or nonenzymatic browning reaction leads to melanoidin formation. If reducing sugars are left to react with biological molecules under physiological conditions, a similar reaction occurs in a process referred to as glycation. Notably, glycation has a considerable influence on the aging process among living organisms. It also has an impact on the pathology of several diseases [1]. In recent times, our emphasis has been on understanding how glycation takes place and the chemical reactivity that reducing carbohydrates have with amino acids, as well as the peptides participating in the process, which are often linked to some diseases, such as Parkinson's, diabetes and Alzheimer's .
The main feature notable in the Maillard reaction includes color formation. However, little is known about the colored moieties that cause the coloration. The light-weight intermediate colored products are referred to as Colored Maillard Reaction Products (CMRP), and it is possible to measure molecule can be easily estimated. It only requires one to check the way it has followed the KID procedure. Nevertheless, tune optimization depends on the system and had to 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, as well as where the ω value is fixed is done to determine how to perform the practical technique.
This study seeks to undertake a comparative study of the way a number of recent density functionals perform in reproducing the chemical reactivity descriptors that the Blue-G1 pigment has in the KID formalism to have sufficient insights about their molecular attributes that future studies use on the chemical reactivity that colored melanoidins of larger molecular weights that form from the reaction that reducing sugars have with proteins and peptides.

Theoretical Background
The theoretical background of this study is similar to the previous conducted research presented  and will be shown here to be complete, because this research is a component of a major project that is in progress.
If we consider the KID procedure presented in our previous works  together with a finite difference approximation, then the global reactivity descriptors can be written as [36,[52][53][54][55] : Electron-accepting power Net electrophilicity 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 [37,52,[56][57][58][59][60][61][62][63] Nucleophilic Fukui function f 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 [40].

Settings and Computational Methods
Consistent with the work presented earlier , this study involved performing computational studies with the Gaussian 09 series of programs [64] whose density functional methods are implemented in the computational package. The gradient technique was relied upon to determine the equilibrium geometries of the molecules; whereas the determination of the force constants and vibrational frequencies involves computing analytical frequencies on various stationary points that are obtained as optimization is completed to check whether the minima were real. The basis set that this work used included Def2SVP for both geometry optimization, as well as frequencies; whereas Def2TZVP was involved in calculating the electronic characteristics [65,66].
Eleven density functionals were selected to calculate the molecular structure and properties that the studied systems have. The selected functionals offer satisfactory results consistently in relation to a number of thermodynamic and structural attributes. Among the functionals are included CAM-B3LYP, which entails Handy and co-workers' long-range-corrected version of B3LYP done through the Coulomb-attenuating method [67]. Another functional relates to LC-wPBEconstituting the long-range-corrected wPBE density functional [68]. In addition, M11entails a range-separated hybrid meta-Generalized Gradient Approximation (GGA) [69]. M11L, on the other hand, is a dual-range local meta-GGA [70], while MN12L consists of a non-separable local meta-Non-separable Gradient Approximation (NGA) [71]. MN12SX entails a range-separated hybrid non-separable meta-NGA [72], while N12 is an NGA [73]. N12SX comprise of a range separated hybrid NGA [72] as the wB97 and wB97X long-range corrected density functionals [74] together with the wB97XD version include empirical dispersion [75]. GGA in these functionals represents the generalized gradient approximation whereby the density functional is dependent on the up and down spin densities together with their reduced gradient; whereas NGA denotes the non-separable gradient approximation in which the density functional tends to rely on the up and down spin densities, as well as their lowered gradient, while also assuming a non-separable form. The various calculations were done with water being the solvent. Integral Equation Formalism-Polarized Continuum Model (IEF-PCM) computations were involved as the Solvation Model Density (SMD) provided [76].

Results and Discussion
This study took the molecular structure of the Blue-G1 intermediate melanoidin pigment from PubChem (https://pubchem.ncbi.nlm.nih.gov), a website that acts as the public repository for information pertaining to chemical substances together with the biological activities with which they are associated. The molecular structure with the IUPAC (International Union of Pure and Applied Chemistry) name includes [34]. The preoptimization of the resultant system involved selecting the most stable conformers. The selection was done using random sampling, which involved molecular mechanics techniques and inclusion of the various torsional angles via the general MMFF94force field [77][78][79][80][81], which involves the Marvin View 17.15 program, which constitutes an advanced chemical viewer suited to multiple and single chemical queries, structures and reactions (https://www.chemaxon.com). Afterwards, the structure that the resultant lower-energy conformer assumes was reoptimized using the eleven density functionals previously mentioned in the previous section together with the Def2SVP basis set, as well as the SMD solvation model in which the solvent was water.
The analysis of the results obtained in the study aimed at verifying that the KID procedure was fulfilled. Upon doing this previously, several descriptors associated with the results that HOMO and LUMO calculations obtained are related to the results obtained using the vertical I and A following the ∆SCFprocedure. A link exists between the three main descriptors and the simplest conformity to 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 is an approximation that is only valid if the HOMO of the radical anion (the SOMO) shares similarity with the LUMO of the neutral system. Consequently, we decided to design another descriptor ∆SL, to guide in verifying how the approximation is accurate. Table 1 illustrates the electronic energies of neutral, positive and negative molecular systems of Blue-G1, HOMO, LUMO and SOMO orbital energies (all in au), while the calculation of J I , J A , J HL and ∆SL descriptors involves using the eleven density functionals, as well as the Def2TZVP basis set that would use water as a solvent that is simulated through the SMD parametrization of the IEF-PCM model.  Afterwards, the study focuses on the other four descriptors analyzing the effectiveness of the density functional in predicting the electronegativity χ, the global electrophilicity ω, the global hardness η and for a combination of the conceptual DFT descriptors, taking into account the energies that the LUMO and HOMO or the vertical I and A are related; where CDFT stands for Conceptual DFT. Next, we consider four other descriptors that analyze how useful the studied density functionals are for the prediction of the electronegativity χ, the global hardness η, 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 results of the calculations of J χ , J η , J ω and J CDFT for the Blue-G1 intermediate melanoidin pigment are displayed in Table 2. As Table 1 provides, the KID procedure applies accurately for 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 MN12SXdensity functional. In addition, the ∆SL descriptor reaches the minimum values when the 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.
It can be seen that the same density functional follows the KID procedure in the rest of the descriptors such as J χ , J η , J ω and J D1 . The significance of these results is attributable to their illustration that reliance on J I , J A and J HL would not be sufficient. For instance, if J χ were considered on its own to apply to each density functional covered in this study, the values would considerably be near zero.
In the case of the other descriptors, only the MN12SX and N12SX density functionals exhibit this behavior. This implies that the results that J χ record are likely to be because of the cancellation of errors by chance.
The considerable success of the approach is however undermined by the issue of tuning being system dependent. Therefore, focus should be on establishing the effectiveness of the behaviors of the fixed RSH density functionals in describing the excitation characteristics. In his works, Becke has recently mentioned that the adiabatic connection and the ideas of Hohenberg, Kohn and Sham apply only to electronic ground states comprise a common misconception [101]. Furthermore, consistent with Baerends et al., the KS model is not appreciated for being superior because of its lowest excitation energy in molecules. Physically, it amounts to an excitation of the KS system rather than electron addition, as would be the case in Hartree-Fock. Thus, it can effectively be used as a measure of the optical gap and is an effective approximation to the gap (in molecules) [102]. In their conclusion, van Meer et al. advance that the HOMO-LUMO gap associated with the KS model tends to be an approximation of the lowest excitation energy, a desirable characteristic with no concerns regarding it [103].
Therefore, calculation of the maximum wavelength absorption of the Blue-G1 pigment involved conducting TDDFT calculations with the aforementioned eleven density functionals at the same level of model chemistry and theory as the calculations for establishing the ground state of the molecule. Table S1 of the Supplementary Information presents the results by comparing the values involved in the ground-state approximation derived from the HOMO-LUMO gap, as well as the TDDFT ones together with the experimental value of 629 nm [34]. Moreover, Figure 1 provides an illustration that compares the results graphically. Notably, the presented results suggest that the differences with the experimental value for λ max tend to have the same order in the various functionals that the current study considers, apart from the N12 density functional. If the λ max values that the HOMO-LUMO gap generates were the ones considered, MN12SX and N12SX would appear to be accurate especially in predicting this value. This does not apply in the rest of the density functionals that the study considers. In fact, the results that the TDDFT calculations give can be improved by expanding the basis set and integrating vibronic corrections; despite the level of accuracy in predicting the excited state property being remarkable if calculations are made at the ground state with the suitable choice of density functional.
Upon verifying that the MN12SX density functional is the best suited for the calculation of the global reactivity density descriptors and in predicting the λ max in agreement with the experimental, Figure 2 presents the optimized structure o the Blue-G1 pigment graphically as calculated based on the theory; whereas Tables S2 and S3 of the Supplementary Materials illustrate the bond length and the bond angles.  Table 3 illustrates the results obtained after calculating the electronegativity χ, chemical hardness η, global electrophilicity ω, electron-accepting (ω + ) and electron-donating (ω − ) powers, as well as net electrophilicity powers with the MN12SX density. The Def2TZVP basis set is used with water acting as a solvent in line with the SMD solvation model. The calculations of the condensed Fukui functions and dual descriptor are done by using the Chemcraft molecular analysis program to extract the Mulliken and NPAatomic charges [104] beginning with single-point energy calculations involving the MN12SX density functional that uses the Def2TZVP basis set. In line with the SMD solvation model, water is utilized as a solvent.
Considering the potential application of the Blue-G1 molecule as an antioxidant, it is of interest to get insight into the active sites for radical attack. A graphical representation of the radical Fukui function f 0 is presented in Figure 3. The condensed electrophilic and nucleophilic Parr functions P + k and P − k over the atoms of the Blue-G1 pigment (excluding the H atoms) have been calculated by extracting the Mulliken and Hirshfeld (or CM5) atomic charges using the Chemcraft molecular analysis program [104] starting from single-point energy calculations of the ionic species with the MN12SX density functional using the Def2TZVP basis set in the presence of water as a solvent 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 charges P + k (M) and P − k (M) and the electrophilic and nucleophilic Parr functions with Hirshfeld (or CM5) atomic charges P + k (H) and P − k (H) are displayed in Table 4. From the results for the local descriptors in Table 4, it can be concluded that C1 will be the preferred site for a nucleophilic attack and that this atom will act as an electrophilic species in a chemical reaction. In turn, it can be appreciated that C14 and C15 will be prone to electrophilic attacks and that this atomic sites will act as nucleophilic species in chemical reactions where the Blue-G1 molecule is involved.
It has been already pointed out that although condensed Fukui functions give interesting results, they are not conclusive. In particular, it has been found that when studying metallic clusters, the condensed Fukui functions can predict the results of nucleophilic and electrophilic interactions with poor reliability [105]. However, from the results obtained in our work, we can present four important considerations: (i) in the first place, it is not the same to obtain conclusions about the reliability of the condensed Fukui functions when studying metallic clusters or even solid systems than when considering pure organic molecules, as happens in our case; (ii) the reliability of the results that we have obtained is impressive because the condensed Fukui functions (and thus the dual descriptor) have been calculated using two different schemes for the partition of the electronic density (i.e., the atomic charges), and the same has been done for the Parr functions: the conclusions about the reaction sites are exactly the same; (iii) in our work, we considered and presented the calculation of the dual descriptor rather than the condensed Fukui functions, and it has been shown by Martínez-Araya that the dual descriptor is more reliable for predicting the electrophilic and nucleophilic sites than the condensed Fukui functions [106]; (iv) as we have shown in several previous works , the reliability of the conceptual DFT descriptors for predicting the reactive sites of a given molecular system is heavily dependent on the goodness of the model chemistry employed for the calculations where we understand for goodness the ability to fulfill the KID procedure mentioned in the Introduction section. Table 4. 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 charges P + k (M) and P − k (M) and the electrophilic and nucleophilic Parr functions with Hirshfeld (or CM5) atomic charges P + k (H) and P − k (H) for the Blue-G1 molecule. Atom

Conclusions
The eleven fixed RSH density functionals, which include CAM-B3LYP, LC-wPBE, M11, N12, M11L, MN12L, N12SX, MN12SX, wB97, wB97X and wB97XD, are examined to establish whether they fulfill the empirical KID procedure. The assessment is done by comparing the values from HOMO and LUMO calculations to those that the ∆SCF technique for the Blue-G1 molecule generates. This is an intermediate melanoidin pigment that is of both academic and industrial interest. The study has observed that the range-separated and hybrid meta-NGA density functionals tend to be the most suited to meeting this goal. In this case, they emerge as suitable alternative to the density functionals once it is established that the behavior of the functionals are tuned using a gap-fitting procedure. They also exhibit a desirable prospect of benefiting future studies in understanding the chemical reactivity that colored melanoidins with larger molecular weights have when reducing sugars react with proteins and peptides.
From the results of this work, it becomes evident that it is easy to predict the sites of interaction of the Blue-G1 pigment under study. This would involve having DFT-based reactivity descriptors including Parr functions and dual descriptor calculations. Evidently, the descriptors were useful in characterizing and describing the preferred reactive sites. They were also useful in comprehensively explaining the reactivity of the molecules. Furthermore, it is also possible to predict the maximum absorption wavelength for the Blue-G1 with considerable accuracy. The prediction would involve the MN12SX density functional beginning with the HOMO-LUMO gap instead of TDDFT calculations. Such a finding is particularly crucial considering the likelihood of it being used to inform the alternative determination method on the color of larger systems such as prosthetic chromophore groups. Such becomes necessary in circumstance where it would not be possible to afford the TDDFT calculations.
Thus, it can be concluded that the model chemistry MN12SX/Def2TZVP/SMD(water) is the best combination of density functional, basis set and solvation model for the prediction of maximum absorption wavelength (starting from the KS HOMO-LUMO gap) and for the prediction of the chemical reactivity sites for the Blue-G1 melanoidin pigment, and it is the recommended methodology to follow in future studies of analog compounds.