DFT study of the stabilization effect on anthocyanins via secondary interactions.

Anthocyanins, which are the labile flavonoid pigments in botanical food, are attracting intensive attention because they can reduce the risk of noncommunicable diseases. Thus, many dietary molecules have been explored to minimize anthocyanin degradation. This study developed a novel model based on the density functional theory (DFT) and conceptual density functional theory (CDFT) to screen small dietary compounds that can stabilize aqueous anthocyanins. The progression of anthocyanin degradation, which was modeled as an aqueous food system, was illustrated using thermodynamic computation and relaxed scanning. The nucleophilic index and dipole moment were applied to quantify van der Waals interaction between anthocyanins and stabilizers. Two equations based on first-order kinetics were established to demonstrate that the equilibrium constant and free energy of the binding reaction between anthocyanins and stabilizers were theoretically important. The change in binding free energy change (ΔG) may be the best indicator of the protection offered by dietary stabilizers on anthocyanins, which was demonstrated by comparisons of computational ΔG with the thermal half time from the previous study on the effects of gallic/ferulic/caffeic acids on anthocyanin stability. Based on established forecasting methods, trans-resveratrol (ΔG = -35.63 kJ/mol) was found to be the best stabilizer among dietary compounds.

Introduction

Anthocyanins (ACNs) are widely distributed pigments, responsible for the red, purple, blue, and deep-almost black colors of many fruits, vegetables, cereal grains, and flowers. Six ACN aglycones (pelargonidin, cyanidin, peonidin, malvidin, delphinidin, and petunidin) frequently occur in nature (Giusti & Jing, 2007). The number and position of the hydroxy group in aglycones and the identity, number, and position of sugar moieties and acyl groups provide distinct ACNs with potential structural and functional variants (Fig. 1).

Fig. 1

Energy barriers of C-O bond newly formed between C2 position of C3G and O of nucleophile in presence/absence of anthraquinone via relaxed scans. a) Hydrolysis reaction of C3G and the nucleophile was water. b) Hydrolysis reaction of C3G with anthraquinone and the nucleophile was water. c) Hydrolysis reaction of C3G and the nucleophile was hydroxyl ion. d) Hydrolysis reaction of C3G with anthraquinone and the nucleophile was hydroxyl ion. e) six major anthocyanidins in nature.

ACNs are receiving intensive attention for their potential health benefits (Gamel, Wright, Pickard, & Abdel-Aal, 2020; Jing Li et al., 2019a, Tsuda, 2012, Wu et al., 2018, Zhang et al., 2021). However, ACNs are susceptible to pH, heat, light, and metal due to the pyran ring, which limits their commercial use. For example, the ACN half-time life (t1/2) decreases with increasing pH values (Kırca, Özkan, & Cemeroğlu, 2007), and under this circumstance, hydroxyl ions are the major nucleophilic reagents that attack the pyran ring. When the temperature increases during heat treatment, water molecules became more aggressive and cause the thermal instability of ACNs even at low pH (Jie Li et al., 2013). The rich π-electrons from ligands can associate with electron-poor ACN cations via secondary interactions (Castañeda-Ovando, Pacheco-Hernández, Páez-Hernández, Rodríguez, & Galán-Vidal, 2009), which help stabilize ACN. Bąkowska, Kucharska, and Oszmiański (2003) suggested that flavones, tannic acid, chlorogenic acid, etc. could prevent ACN color from fading during heating. Chung, Rojanasasithara, Mutilangi, and McClements (2016) found that Vanillin, quillaja saponins, epigallocatechin gallate, and green tea extract could also stabilize ACNs (Table 1).

Table 1

Comparison between computational data at 368.15 K, p and experimental data of cyanidin-3-O-sophoroside-5-O-glucoside complexes.

Stabilizers	Experimental data*	Computational data
	t1/2,95℃ (hr)	Nucleophilic index (eV)	Dipole moment (Debye)	ΔG (kJ/mol)
Caffeic acid	3.39	2.954	5.643	−17.53
Ferulic acid	3.50	3.020	6.762	−19.01
Gallic acid	9.64	2.578	6.914	−21.82

*Data from Qian, B. J., Liu, J. H., Zhao, S. J., Cai, J. X., & Jing, P. (2017). The effects of gallic/ferulic/caffeic acids on colour intensification and anthocyanin stability. Food Chemistry, 228, 526–532.

The hyperchromicity or hypochromicity of ACN was attributable to the typical π-π complex or charge transfer complex in the presence of these ligands (Li et al., 2019a, Li et al., 2019b, Wu et al., 2018). The secondary interactions are sometimes be so strong that they generate complicated complexes with distinguished properties such as commelinin (Shiono, Matsugaki, & Takeda, 2008). Notably, the interactions are intermolecular and intramolecular, especially in acylated anthocyanins (Bąkowska-Barczak, 2005). Acyl groups, especially aroyl groups, may provide anthocyanins with additional stability (Hoshino, Matsumoto, & Goto, 1980). Furthermore, some acylated anthocyanins with sandwich-type conformations were reported, and no hydration was detected (Figueiredo, George, Tatsuzawa, Toki, Saito, & Brouillard, 1999).

Given the ACN complex richness, we aimed to develop a DFT computational method to quantify the stability of anthocyanin complexes with ligands. First, the mechanisms of ACN degradation and its thermodynamics were analyzed based on intermediate/final products to screen natural stabilizers. Then, the theoretic screening model was established and verified by comparison with experimental data from the previous literature.

Materials and methods

Materials

Anthocyanin was selected as our target molecule, and several dietary molecules with biological activity were selected as potential stabilizers: 1) cis-/trans-resveratrol, which is a naturally occurring phytoalexin produced by some spermatophytes such as grapevines (Frémont, 2000; Salehi et al., 2018); 2) genipin, which is an active constituent isolated from the fruit of Gardenia jasminoides Ellis and some fruits, is widely used in traditional medicine, and even has potential anticancer effects (Shanmugam et al., 2018); 3) carnosic-8-lactone, which is a derivative from carnosic acid and carnosol and highly abundant in rosemary (Rosmarinus officinalis)(Loussouarn, Krieger-Liszkay, Svilar, Bily, Birtić, & Havaux, 2017; Temerdashev, Milevskaya, Shpigun, Prasad, Vinitskaya, & Ryaboko, 2020); 4) myrtenol, which is a monoterpene with multiple pharmacological activities extracted from myrtle and was committeed by WHO expert and recognized as a safe favoring ingredient (Britto et al., 2018; Selvaraj, Valliammai, Sivasankar, Suba, Sakthivel, & Pandian, 2020); 5) R-4-acetyl-2-carene, which is a derivative from carene (Ho, 1983; Macaev et al., 2007).

Computational methods

Computations were performed using the Gaussian09 (revision D.01) suite (Frisch et al., 2013) of the quantum chemical program. All structures were optimized in an implicit solvent model using the M062X hybrid meta exchange–correlation functional (Zhao & Truhlar, 2008) with D3 correlation (Sure & Grimme, 2013) in conjunction with Pople’s valence triple-ζ 6-311G(d,p) basis set (Krishnan, Binkley, Seeger, & Pople, 1980) for all atoms. All structures were in a local minimum potential energy surface with zero imaginary frequency, or at the first-order saddle point on the potential surface with one imaginary frequency. Intrinsic reaction coordinate (IRC) calculations were additionally performed to further characterize the true nature of the TSs (Gonzalez & Schlegel, 1989). Harmonic vibrational frequencies, thermal corrections, and entropic corrections at specific temperatures and p were obtained from frequency calculations. A global multiplicative harmonic frequency scaling factor for M062X/6-311G(d,p) (0.9460) was used to correct the calculated harmonic frequencies and thermal data (Kashinski et al., 2017). The quasi-harmonic approximation proposed by Grimme (Grimme, 2012) considers entropic terms for all normal-mode frequencies below the cutoff of 100.0 cm−1 using the free-rotor approximation; above this threshold, the RRHO approximation was retained. All frequency and thermal corrections were computed by using Shermo (Lu & Chen, 2020). The single point energy and dipole moment were computed using the M062X-D3 functional with the def2-TZVP basis set (Weigend & Ahlrichs, 2005). The SMD continuum universal solvation model (Marenich, Cramer, & Truhlar, 2009) was used to consider the water solvation effects. Because many structures are possible when discussing the interactions between two molecules, the identification of the most stable structure of the complex demanded detailed conformational sampling or scanning. Therefore, 1000 initial structures were generated by using Molclus (Kırca et al., 2007); after the pre-optimization by using MOPAC2016 (Stewart, 1990, Stewart, 2016) and PM7 (Stewart, 2013), a semi-experimental method, only the structure with the lowest energy was used for further optimization and calculation. The wave function analysis and various quantities from concept density functional theory (CDFT) were computed using Multiwfn (Version 3.7) (Lu & Chen, 2012) at the same calculation level as the optimization.

Results and discussion

Mechanism and thermodynamic computation of anthocyanin degradation

The mechanism of anthocyanin (i.e., cyanidin-3-O-glucoside, C3G) degradation in Scheme 1 was deduced based on the intermediate/final products reported in previous studies (Adams, 1973, Patras et al., 2010, Seeram et al., 2001, Sinela et al., 2017). The Gibbs free energy profile was analyzed in Scheme 2. To investigate the possible transition state in the C3G degradation, we designed relaxed scans of newly-formed C-O bonds of anthocyanin and nucleophile at 0.005 Å or 0.010 Å intervals (depending on the structure), by fixing the bond length and optimizing the remainder to a local minimum point at the same level as the geometry optimization.

Scheme 1

Mechanism of cyanidin-3-O-glucoside degradation. (i) Nucleophilic attacking; (ii) Acid-base reaction; (iii) Intramolecular elimination reaction; (iv) 1,4-addition reaction; (v) Keto-enol tautomerism; (vi) Inverse aldol reaction; (vii) Oxidation.

Scheme 2

Gibbs free energy profile of cyanidin-3-O-glucoside degradation in presence/absence of anthraquinone at 298.15 K, p. Structure 6 and 7 are shown using Pymol due to their complex structures. Units of Gibbs free energies were kJ/mol.

First, a water molecule undertook a nucleophilic attack on the C2 position of the pyran ring (step i). Then the pyran ring lost a proton and an intramolecular elimination reaction occurred to produce an α,β-unsaturated ketone (steps ii and iii). After an addition reaction with a water molecule, an inverse aldol condensation reaction occurred within the molecule, and then the pyrolysis products were finally formed (steps iv-vi). Theoretically, glycoside may be removed from C3G at any steps in the degradation reaction. Since there were no confirmable experiments and the computed Gibbs free energy of the hydrolysis of the glucosidic bond was 12.05 kJ/mol, we were inclined to keep glucoside at the first few steps in the degradation progress.

After identifying all steps of the C3G degradation reaction (Scheme 1) and the associated energy features, the overall Gibbs free energy profile was analyzed and shown in Scheme 2. The reactions in protons had to be fitted due to computational difficulty on the dissolution free energy of the proton. The fitting method (Kheirjou, Imanzadeh, Rezaei, & Safari, 2019) was applied while six molecules with experimental pKa values were adapted from references (Freitas et al., 2018, Mazza and Brouillard, 1987, Vidot et al., 2016) (see supplementary document: Table S-1, page S43). In Scheme 2, the C3G degradation process appeared to be thermodynamically and kinetically favorable. However, in the presence of anthraquinone, the C3G degradation reaction unlikely occurred because the free energy change of the reaction increased by ∼ 40 kJ/mo and the relaxed scanning result of the reaction process in Fig. 2b indicated an impossible transition state or intermediate.

Fig. 2

Surface fitting of ΔH, nucleophilic index and dipole moment, respectively. Function of the fitting result would be z = -54.44 + 4.59x-2.85y-2.12x2 + 0.041y2 + 0.80xy, R2 = 0.5050. Surface fitting result of z = z0 + ax + by + cx2 + dy2 + fxy, where z, x and y were ΔH, nucleophilic index and dipole moment, respectively. Function of the fitting result would be z = -54.44 + 4.59x-2.85y-2.12x2 + 0.041y2 + 0.80xy, R2 = 0.5050.

Fig. 2a shows the energy barriers during the hydrolysis reaction of C3G alone when the nucleophile was water. However, Fig. 2b suggests that there was no local minimum potential energy surface, i.e., a stable structure for potential intermediates while water molecular approached C3G in the presence of anthraquinone. Fig. 2c and 2d indicate that the hydrolysis of C3G can be rapid regardless of the anthraquinone presence in a system where hydroxyl ions prevailed. This result may attribute to the lack of a potential barrier when hydroxyl ions approach C3G.

Quantify secondary interactions and their energy

Major secondary interactions are generally hydrogen bonds, π-π stacking, charge transfer interactions, and electrostatic interactions. Hydrogen bonds have been sufficiently discussed, and the method to estimate the energy of intramolecular and intermolecular hydrogen bonds was already given by Grabowski, 2001, Wendler et al., 2010, respectively. Hydrogen bonds will not be discussed here since they were nothing special for anthocyanins.

Most studies used phenolic acids as a stabilizer of anthocyanin (Fan, Wang, Xie, Zhang, Li, & Zhou, 2019), and anthocyanin was believed to be a good electron acceptor (Ferreira da Silva, Lima, Freitas, Shimizu, Maçanita, & Quina, 2005). In this study, cyanidin with benzene/benzene derivatives was chosen to prove and quantify this characteristic of ACN. As an index of polarity, dipole moment was positively correlated with the interaction energy (Becke & Johnson, 2005). The nucleophilic index, which is defined as the difference in HOMO (highest occupied molecular orbital) energy of a molecule and tetrachloroethylene, was a quantified index for the nucleophilicity of a molecule. Thus, the dipole moment and nucleophilic index were selected as indicators that may relate to the dispersion interaction energy between anthocyanin (i.e., cyanidin) and benzene derivatives. The data are shown in Table S-3 (Part D in the supplementary document on page S45) and were fitted to a surface that monotonically decreased in both × and y directions (Fig. 3), although the R2 was only 0.5050, which made this fitting only an exhibition of tendency. Nevertheless, molecules with high nucleophilic index and high dipole moment were supposed to have more enthalpy changes when binding to cyanidin. This characteristic can also make phenolic acids potential stabilizers of ACN due to their high nucleophilic index since their hydroxy groups act as electron-donating groups.

Kinetics: Binding free energy change as a gold standard indicator for stabilizers

ACN degradation regardless of stabilizers has been proven to be a first-order kinetic reaction in previous studies (Adams, 1973, Qian et al., 2017). If the system is highly acidic, water is the attacking molecule. If k1, k2, and k3 are the rate constants of the degradation of anthocyanin cations, anthocyanin-stabilizer complexes, and measured data, respectively, then the following first-order rate equation can be established (Eq. (1)).

In the expression, [A], [AM], and [A]sum are the concentrations of anthocyanin, anthocyanin-stabilizer complex, and anthocyanin in total, respectively. There was an equilibrium to eliminate variables, where K was the equilibrium constant of the binding reaction between anthocyanin and stabilizer, and [M] was the stabilizer concentration (Eq. (2)).

In Section 3.1, we had proved that anthocyanin with a stabilizer could not be attacked by a water molecule in an acidic system, thus k2 can be set as 0. After the above expressions were fully simplified, the measured rate constant was expressed in Eq. (3).

Similarly, in a system with enough hydroxide anion, the degradation of the anthocyanin-stabilizer complex cannot be neglected, which leads to Eq. (4).

The coefficient of the variable in Eq. (4) was positive since with a stabilizer, k1 was usually larger than k2. Therefore, both Eqs. (3), (4) show that a larger K leads to lower rate constants of degradation. The binding free energy (ΔG) of ACN and stabilizer were calculated in Table 2 to verify whether the computational data were consistent with experimental data from the reference. The ΔG was −21.82 kJ/mol < -19.01 kJ/mol < -17.53 kJ/mol for cyanidin-3-O-sophoroside-5-O-glucoside with gallic acids, ferulic acids, and caffeic acids, respectively. The binding free energy changes of cyanidin-3-O-sophoroside-5-O-glucoside with stabilizers were consistent with the sequence of half-life values of the complexes (Qian et al., 2017), which suggests that ΔG predicted the stability of ACN with stabilizers.

Table 2

Prediction about some potential natural micromolecular stabilizers for cyanidin-3-O-glucoside at 298.15 K, p.

Stabilizer	Number of hydrogen bonds with C3G	Nucleophilic index (eV)	Dipole moment (Debye)	ΔG (kJ/mol)
trans-Resveratrol	1	2.171	4.783	−35.63
Genipin	2	2.204	2.785	−19.71
Carnosic-8-lactone	1	2.464	5.181	−18.22
cis-Resveratrol	1	2.390	5.317	−14.73
Myrtenol	1	3.184	1.853	9.33
R-4-acetyl-2-carene	1	3.412	3.239	11.41

Prediction of stabilizers

With the established and clarified forecasting method, a few natural compounds with the potential ability to stabilize ACN were chosen for comparison, as shown in Table 3. Among all selected molecules. Trans-resveratrol was a protential natural compound to stabilize ACN with a binding free energy of −35.62 kJ/mol. In addition, genipin, carnosic-8-lactone, and cis-resveratrol may be good stabilizers.

Conclusion

The thermodynamic computation of C3G degradation alone or with anthraquinone demonstrated that C3G with the stabilizer appeared unlikely to degrade in a highly acidic system. Stabilizers such as anthraquinone could also increase the free energy of the hydrolysis reaction of C3G, which made the degradation thermodynamically much harder. In the analysis of cyanidin with benzene and benzene derivatives, the nucleophilic index and dipole moment were demonstrated to indicate the intensity of van der Waals interactions. Equations based on first-order kinetics demonstrated that the equilibrium constant and free energy of the binding reaction between anthocyanin and stabilizer theoretically mark the capability of stabilizers to stabilize anthocyanins. Based on the established forecasting methods, molecules such as trans/cis- resveratrol, genipin, and carnosic-8-lactone were predicted as potential anthocyanin stabilizers. Among them, trans-resveratrol was screened as the best potential anthocyanin stabilizer, which suggests its potential in the storage of wine, juice, or other liquid foods. Additionally, the established model can be applied to screen other more potential stabilizers.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.