Computational insights into substituent effects on the stability and reactivity of flavylium cation analogs of anthocyanins

Synthetic flavylium salts and natural anthocyanins undergo hydration by water at physiological pHs, fading out the color. Equilibrium constants of hydration ( K h ), tautomerization ( K t ) and isomerization ( K i ) of 21 substituted flavylium salts were rationalized based on molecular properties calculated by density functional theory (DFT)


Introduction
3][4][5] In aqueous media, anthocyanins and flavylium salts can undergo a complex network of chemical reactions, as shown in Scheme 1. Below pH 2, the colored flavylium cation, AH + , is the most stable form in solution.8][9][10][11]  Much progress has been made in the last years in understanding the factors affecting the reactivity of flavylium salts, especially regarding their proton transfer [6][7][8][9][10][11][12] and charge transfer 4,13 reactions in water and at micelle-water interfaces.However, measurements of the thermodynamic equilibrium constants can be demanding because of the competitive formation of B and the chalcone isomers, as well as their deprotonated forms. 11The one-electron reduction of AH + in aqueous solutions is irreversible, leading to large uncertainties in the redox potentials. 13,14Consequently, the rationalization of the reactivity from a physical organic chemical point of view is difficult because there are relatively few compounds for which all of the experimental constants (either kinetic or thermodynamic) have been determined.For these reasons, procedures based on computational methods are an attractive alternative for estimating the reactivity of flavylium salts and anthocyanins.

OH
Previously we applied Density Functional Theory (DFT) with success to calculate the pKas of a series of nine substituted flavylium salts, as well as the one-electron reduction potentials of five flavylium salts. 14We also developed 15 Linear Free Energy Relationships (LSER) between experimental pKh, Kt and Ki values for flavylium salts and Hammett substituent constants for substituents on both activated and non-activated positions, resulting in a model for the prediction of the "apparent pKa" (Kap = Ka + Kh + KhKt + KhKtKi) from structure. 1 Empirical Hammett relationships have been used for decades for the rationalization of organic reactions, 16 their mechanisms and other types of chemical phenomena, namely NMR shifts, 17 frontier orbital energies 18 or core-electron binding energy shifts. 19Our objective in this study was to identify ab initio quantum molecular quantities that can be employed as surrogate Hammett parameters, that is, as oneparameter regression descriptors.Such descriptors should facilitate a more basic understanding of the underlying aspects responsible for the observed variations in the reactivity of flavylium cations and anthocyanins.Thus, going beyond earlier work that employed DFT to estimate pKa values, 14,20 we have now expanded the methodology to compute the hydration equilibrium constants pKh.

Assessment of the hydration equilibrium constants of flavylium salts
The flavylium salts employed in this study and the corresponding experimental equilibrium constants are listed in Table 1.The first step was the evaluation of calculated hydration equilibrium constants in order to ensure that the level of theory employed was adequate.Since the hydration equilibrium constant, Kh, can be written as  ℎ =  ℎ ′   /[ 2 ], the alternative thermodynamic cycle shown in Scheme 2 was employed in the theoretical calculations of  ℎ ′ because it required only the calculation of the free energy of OH -in the gas phase, rather than of the free energies of both water on the left side and H + on the right side of the equation (thus minimizing potential sources of error).The value of Kh was then obtained by using the experimental autoionization equilibrium constant of water, Kw.The calculated free energies are presented in the Supplementary Material (Table S1).
Scheme 2. Thermodynamic cycle for the calculation of the term K'h.
Because the compounds studied are relatively large, with at least 16 non-hydrogen atoms, the four solvation models examined provided good results over a range of ca. 10 pK units (Figure 1).The solvation free energies calculated by the universal solvation model based on solute electron density (SMD), the Cramer-Truhlar solvation model (SM5.4P) and both polarized continuum models (PCM and PPM2) were particularly good, resulting in predicted pKh values with a mean absolute deviation (MAD) of 0.55-0.68relative to the experimental values.
The level of theory utilized in the gas phase calculations was consistent with that employed elsewhere. 14,21The addition of diffuse functions to a double-ζ basis set considerably reduces errors in the calculation of reaction energies by DFT. 22The use of a triple-ζ basis set should therefore have little effect on the calculated electronic energies because the valence basis set is already saturated. 23Nonetheless, we tested the computation of pKh employing single point electronic energies obtained with the triple-ζ basis set 6-

AUTHOR(S)
311++G(3df,3pd).However, in only one case, specifically compound 5, did the calculated value show an improvement relative to the double-ζ basis set result.Compound 5 has a OCH3 group at position C4', indicating that the description of the electron density needs to be enhanced when methoxy groups are present.

Linear correlations of calculated molecular properties with substituent effects
Previously, we demonstrated 15 that the Hammett resonance R and meta-like m parameters were the best descriptors for the effects of substituents at activated and non-activated positions, respectively, on the hydration equilibrium, pKh, of flavylium salts.The electronic effects of multiple substituents at activated and non-activated positions attached to the flavylium moiety were found to be additive, 15 resulting in linear correlations with the total X of the molecule (see Table S2 in Supplementary Material).
Initially, local descriptors such as interatomic distances, bond angles, dihedral angle between benzopyryl and phenyl rings, atomic charges and bond energies were examined.However, none of the geometrical parameters tested correlated with either X or pKh, as perhaps might be expected.Also, parameters such as the natural atomic charges of the O1 and C2 atoms of the AH + form and the corresponding bond energy performed only at a qualitative level.In conclusion, although the hydration reaction of flavylium salts involves large electrostatic interactions, the local charge values did not quantitatively describe the reactivity of flavylium cations.However, a qualitative analysis of the correlations indicated that electronreleasing groups do indeed enhance the charge density at the reaction center of AH + , thereby increasing pKh.
From the concepts of MO theory, the orbital energies -εHOMO and -εLUMO are related to the ionization potential (IP) and electron affinity (EA), respectively. 25Within the Kohn-Sham (KS) framework of DFT, the eigenvalue of the highest-occupied KS orbital corresponds to the negative of the IP and the

AUTHOR(S)
eigenvalue of the lowest-unoccupied KS orbital to the negative of the EA. 26 Because the longest-wavelength absorption maxima of flavylium cations and quinonoidal bases are largely determined by the HOMO-LUMO electronic transition 14,27 and because IP -EA is simply the difference between εLUMO and εHOMO, we shall focus on the energies of these frontier orbitals in order to seek relationships linking molecular properties to the reactivity of flavylium salts.Furthermore, the HOMO-LUMO gap is related to the chemical hardness η 28 and the HOMO and LUMO energies are also associated with the absolute electronegativity χ, because χ = (IP + EA)/2.For cations, the EA effectively becomes the electronegativity, 25 since electron transfer from the cation to another molecule is unlikely.Because the exact functionals for exchange and correlation are unknown, the vertical ionization potential (VIP) and vertical electron affinity (VEA) of AH + were calculated based on the finite difference approximation employing the energies E(N), E(N-1) and E(N+1), where N is the number of electrons in the molecule (vide Supplementary Material, Table S3).
We confirmed that there is indeed a good correlation linking VEAAH+ and X (X = 4.228VEAAH+ , R 2 = 0.975, sd = 0.04, F = 779, n = 21) and that electron-releasing groups (X < 0) decrease the value of VEAAH+, as shown in Figure 2a.The fair correlation for VIPAH+ (X = 1.501VIPAH+ , R 2 = 0.844, sd = 0.25, F = 108, n = 21) indicates that it is much less sensitive to substituents given the cationic nature of the AH + form.The relationship between VEAAH+ and X allows, for instance, the estimation of group contributions to the total X, the evaluation of whether to use the through-conjugation parameter R + (rather than R) for specific cases, or the estimation of theoretical -values for substituents with unknown Hammett parameters.
The good correlation found between σX and VEAAH+ led us to wonder how this approach might perform for the other multiequilibria of the flavylium salts in aqueous solution.According to Pearson's HSAB principle, 25,27 reactions with hard solvents such as water should become less favorable with the decrease of the absolute hardness of the molecule, which implies a decrease in the difference between the ionization potential and the electron affinity.Indeed, Lietz et al. 27 found two distinct equations correlating the apparent pKap and the chemical hardness  of 4,7-substituted and 7-OMe substituted flavylium cations.With regards to pKh, however, there is no reasonable correlation with , reflecting the fact that VIPAH+ is associated with the HOMO orbital energy, which in turn correlates poorly with pKh (R 2 = 0.281).Because pKh depends only weakly on VIPAH+, one might expect that the hydration equilibrium is dependent on VEAAH+ alone.Indeed, a reasonable correlation was found between VEAAH+ and pKh (Figure 3a, represented as the difference ∆VEA between the VEA of AH + and the corresponding value for compound 1), pointing to the stabilization of AH + by electron-donor substituents, i.e., the pKh of flavylium salts increases as the electron affinity of the flavylium cation decreases.This strongly supports the explanation given earlier by us 13 for the charge-transfer mediated stabilization of flavylium ions against hydration upon complexation with electron rich copigments.The weak correlation between ∆VEAAH+ and the negative logarithm of the tautomerization equilibrium -log(Kt) (R 2 = 0.598, n = 10) improves substantially (to R 2 = 0.935) if the outlier 11 is excluded from the correlation (Figure 3b).No relationship was found between ∆VEA of the Z-or E-chalcone forms and -log(Ki) (Figure 3c), in accordance with the absence of correlation between the redox potentials of substituted chalcones and Hammett constants. 29Unfortunately, the set of experimental Ki data is too small to permit an independent treatment of ring effects in chalcones.Likewise, the pKa values of flavylium compounds were found to be insensitive to VEA (or σX) of the AH + or A forms (Figure 3d), reminiscent of the lack of correlation of the ground-state acidities of pyranoflavylium cations with Hammett sigma parameters. 30Pyranoflavylium cations are synthetic analogs of the pyranoanthocyanins, which are formed during the maturation of red wine via condensation reactions between anthocyanins and copigments or yeast metabolic products.In order to gain more insight into substituent effects on flavylium salt multiequilibria we focused on core-electron binding energies (CEBE), because CEBE shifts correlate nicely with Hammett constants. 19A simple way to calculate the binding energy is provided by Koopman's theorem, which states that the negative of the orbital energy is equal to the ionization energy.This crude approximation does not take into account orbital relaxation effects, however.To make sure that the level of theory proposed was adequate one calculated the 1s binding energies for O, N and C atoms in a series of organic compounds with diverse functional groups, totalizing 33 data points, and compared with the experimental data compiled from NIST Xray Photoelectron Spectroscopy Database. 31The results presented in Table S4 and Figure S1 (Supplementary Material) are quite good, with R 2 coefficients greater than 0.980 in all three cases.The dispersion in binding energies of O atom mostly reflect sensibility to the instrument setup, such as the kind of solid substrate were the compound is deposited, and experimental settings.
As seen earlier, the deprotonation equilibrium constant, pKa, of flavylium compounds is insensitive to Hammett constants or VEA values calculated with frontier orbitals, leading us to focus our attention on core electrons in the search for possible correlations.In polyhydroxy flavylium compounds, the acidity of OH groups in the benzopyrylium part of the molecule (A and C rings) is higher than that of OH groups on the phenyl ring (B-ring) by ca. 2 pKa units (vide compounds 2 and 3 in Table 1; refs.32 and 33).For convenience, it was assumed that multiple OH groups on the benzopyryl moiety can be effectively treated as an average hydroxyl group and that OH groups on the phenyl ring will contribute to the pKa only in the absence of an ionizable group on benzopyrylium moiety.The good correlation (R 2 = 0.956, excluding compound 10) shown in Figure 4d again points to the stabilization of the AH + form by electron-donating substituents.The overall picture is that core-electron binding energies are indeed useful for rationalizing the reactivity of this class of compounds, being more sensitive to substituent-induced changes in reactivity than descriptors based on frontier orbitals such as Electron Affinities.

Implications for the prediction of the apparent pKap of flavylium salts
The stability of the AH + form can be regarded as an "apparent" acidity constant, Kap, for the equilibrium between the flavylium cation and a set of "conjugate bases", defined as

[CB] = [A]+[B]+[CZ]+[CE]. If the isomerization is slow and CE
has not yet been formed or if E-chalcone is not present, known as pseudoequilibrium conditions, the "apparent pKa" is given by pKap1 = -log(Ka + Kh + KhKt).However, when significant amounts of CE are present and the system has reached the final thermodynamic equilibrium composition, pKap2 = -log(Ka + Kh + KhKt + KhKtKi).These two definitions of the "apparent pKa" coexisted in the past, generating some confusion in reported pKap values in the literature. 24ased on the correlations found with ∆CEBE, values of pKap1 and pKap2 were calculated and compared with experimental values (Figure 5a).The mixture of definitions of the "apparent pKa" is evident here because part of the experimental pKap data lies close to the calculated values of pKap1, while the remainder falls in the domain of pKap2 values.Using Hammett substituent constants as descriptors, a non-linear correlation with pKap can be expected given the complex nature of the mixture of conjugate bases, CB.Indeed, separating the experimental values into pKap1 and pKap2 values and plotting them as a function of X (Figure 5b) indicated a change in reactivity, with a plateau near the limiting values of 4 and 3 for pKap1 and pKap2, respectively, when X < -1.Both equilibria follow a similar trend with X and the deviation of compound 11 from the general trend can be attributed to the unusual pH dependence of its reported experimental hydration constant, Kh, presumably reflecting significant cross-conjugative stabilization of the positive charge by the dimethylamino substituent.

Behavior of anthocyanin pigments in nature
Table 2 lists experimental data for 10 of the most common anthocyanins found in nature and Figure S2 (Supplementary Material) presents the equilibrium constants as a function of the calculated CEBE differences (Table S6 in Supplementary Material).Based on the findings for the equilibrium constants and apparent pKap of synthetic flavylium salts one can conjecture that the chemical behavior of natural anthocyanin pigments should follow the same trends observed before, only displaced from the original correlations by an offset since these natural pigments are subjected to steric effects at position C3.This hypothesis is also supported by linear free energy relationships (LSER) based on Hammett correlations with synthetic flavylium salts and natural anthocyanins evaluated together. 15Thus, using the correlations previously found and applying the corresponding offset for each process, the estimated pKh, pKa and pKap2 from CEBE calculations gave mean absolute deviations of 0.39, 0.11 and 0.37 pK units, respectively.
The comparison of the CEBE values of compounds in Tables 1 and 2, indicated that the bulky sugar units at C3, characteristic of anthocyanins, displaces the isomerization equilibrium towards the CZ form to such an extent that the concentration of CE species is very small to be detected in most cases. 34The tautomerization equilibrium in anthocyanins, on the other hand, is displaced to the hemiketal form.Also, pKa and pKh decrease ca.1.1 and 3.5 pK units, respectively, when comparing flavylium salts and anthocyanins with similar CEBE values, even observing that the corresponding anthocyanins have lower X values (Table S7 in Supplementary Material).

Conclusions
An interesting panorama emerges from the correlations of the core-electron binding energies with equilibrium constants.One strategy employed to protect the colored AH + species against the nucleophilic attack of water at carbon C2 is the inclusion of electron-donating substituents on the flavylium cation moiety.However, at the same time, this approach also leads to stabilization of the colorless Z-chalcone form.This intricate effect of substituents on the different pH-dependent multiequilibria of flavylium salts results in a non-linear relationship between descriptors of substituent effects and the apparent pKap, which levels off at values of ca. 4 and 3 for pKap1 and pKap2, respectively, when X < -1.In nature, the chemical structures of the flavylium cation chromophore of the most abundant anthocyanin types (malvidin, delphinidin, pelargonin, cyanin, petunidin and peonin) differ essentially by the presence or absence of OH or methoxy groups at positions 3´and 5' (or in some case due to an additional glycosylation of the OH group at the 5-position) and have pKaps between 1.7 and 3.0. 24,35Their X values fall between -1.2 and -2.1, independent of the substitution pattern if the inductive/resonance effects of the glucosyl (small, similar to t-butyl 15 ) or sambubiosyl units are not taken into account.The sugar units attached to the oxygen at the 3-position increase the solubility of anthocyanins, but their steric effect, which promotes hydration of the AH + form, is compensated by smaller X values in nature.Rather than attempting to promote stability via the use of substituents, strategies such as the incorporation of flavylium cations in porous clays 37 may be a more useful way to increase the thermal and photochemical stability of the color of this class of compounds.

Geometry optimization
Geometry optimizations of the flavylium salt structures were performed without any geometry constraint with the hybrid functional 38 mPW1PW91 in vacuum using the 6-31+G(d,p) basis set.The selection of this functional was based on our previous study 14 showing that the electronic transitions for flavylium cations and quinonoidal bases calculated with this functional were in better agreement with experimental results than those calculated with B3LYP.Harmonic frequency calculations indicated that all stationary points found were minima on the electronic potential energy surface, i.e., no imaginary frequencies were found.In the condensed phase, the fully relaxed geometries were obtained at the mPW1PW91/6-31G(d) level and the Integral Equation Formalism for the Polarizable Continuum Model (IEFPCM) 39 described the implicit solvent.
The united atom topological model, 40 UA0, was used to build the molecular cavity.

Natural population analysis
The natural population analysis (NPA) was carried out with the optimized geometries in condensed phase at the mPW1PW91/6-31G(d) level in order to obtain the natural charges, Qn, and orbital energies of all species (AH + , B, Z-and E-chalcones).NPA is a partitioning procedure that transforms a molecular wave function into a localized and orthogonal atomic form satisfying the Pauli exclusion principle.Atomic charge is not physical observable and its value depends on the choice of the basis set and the scheme adopted for the partitioning of the electron density matrix (schemes that do not come from first principles).The advantage of the NPA method is that it solves some shortcomings associated with the Mulliken population analysis procedure such as the basis set dependency problem. 41

Evaluation of the vertical electron affinities (VEA)
The vertical electron affinities were evaluated as VEA(AH+, B) = E(AH+, B) -E(AH+, B) -, where E(AH+, B) is the electronic energy of the flavylium cation (or hemiketal) optimized at the mPW1PW91/6-31G(d) level in the condensed phase and E(AH+, B) -is the electronic energy of the corresponding reduced form calculated at the same geometry.

Calculation of the hydration equilibrium constant pKh
The hydration equilibrium constant Kh can be written as Kh = K'hKw/[H2O], where the ionic product of water is In these equations, ∆   ( 2 ) is the free energy of dissociation of water in the gas phase,    () is the standard free energy of species "i" in the gas phase, ∆  * () is the solvation free energy of "i" and ∆  * ( + ) is the free energy change for the proton transfer reaction in the aqueous phase.The open circle © AUTHOR(S) superscript ( o ) represents free energies that use a standard-state gas-phase pressure of 1 atm and the superscript asterisk ( * ) denotes free energies that use an aqueous phase standard-state of 1 mol•L -1 .The term ∆  → * = −(24.46)was added in the ∆   ( + ) expression to take into account the change of concentration units from the gas phase reference state to the liquid phase reference state (1 atm to 1 mol•L -1 ).At 298.15 K ∆  → * is −1.89 kcal•mol -1 .Note that this term cancels out in reactions where the same number of moles of reactants and products are transferred from the gas phase to solution. 42,43he thermal contributions to the gas phase free energies of the molecules were carried out within the framework of statistical thermodynamics at the mPW1PW91/6-31+G(d,p) level.The electronic contribution to the gas phase free energy was acquired by single-point calculations with the 6-311+G(2d,2p) basis set.The Gibbs free energy of OH -in the gas phase was calculated as follows: first, the H2O molecule was optimized in vacuum at the mPW1PW91/6-31+G(d,p) level and the thermal contributions were carried out at the mPW1PW91/6-311+G(2d,2p) level to calculate the term    ( 2 ).The gas phase free energy for OH -was then calculated from the equation ∆   ( 2 ) =    ( + ) +    ( − ) −    ( 2 ), taking    ( + ) = -6.28kcal•mol - 1 from the Sackur−Tetrode equation 21 and the experimental 44 value for ∆   ( 2 ) = 383.7 kcal•mol -1 .We employed the experimental value for the free energy of OH -in the gas phase for the purpose of minimizing potential sources of error.Since the solvation model directly affects the calculation of pKh, the solvation free energies were evaluated through four different approaches: by IEFPCM single point at the HF/6-31G(d) level with UAHF radii, 45 by IEFPCM single point reparameterized with Bondi 46 radii (PCM2) at the HF/6-31G(d) level, by solvation model SM5.4P at the PM3 level 47 and by the SMD model. 48The experimental value 49 of ∆  * ( − ) = -104.6kcal•mol -1 was utilized in the calculation of the pKh with the PCM2 method, given that this model is not parameterized for ions.The calculations of the terms in the thermodynamic cycle, natural population analysis, vertical electron affinities (VEA) and solvation energies were performed with the Gaussian 03 and 09 packages 50,51 and AMSOL. 47

Figure 1 .
Figure 1.Calculated pKh values for the flavylium salts studied and the respective mean absolute deviations (MAD) for the solvation models utilized.The straight line corresponds to pKh calculated = pKh experimental and is included to guide the eye.

Figure 2 .
Figure 2. Correlation of the sum of the Hammett constants of substituents of flavylium cations with the: (a) vertical electron affinity and (b) core-electron biding energy of C2(1s) of the AH + form.

Figure 3 .
Figure 3. Representations of (a) pKh and (b) -log(Kt) as a function of the vertical electron affinity of AH + , (c)log(Ki) as a function of the VEA of the CE form and (d) pKa as a function of the VEA of the AH + form.The VEAs of the flavylium salt species were calculated at the mPW1PW91/6-31G(d) level of theory.

Figure 4 .
Figure 4. Representations of (a) pKh and (b) -log(Kt) as a function of the difference between the core-electron binding energy of C2(1s) of AH + and the corresponding value for the unsubstituted compound 1, CEBE; (c)log(Ki) as a function of the average CEBE of C2'(1s) and C6'(1s) in the CZ form; and (d) pKa as a function of CEBE of the O(1s) of the oxygen of the hydroxyl groups of the AH + species.The CEBE values of the flavylium salt species were determined from Natural Population analysis at the mPW1PW91/6-31G(d) level of theory.

Figure 5 .
Figure 5. (a) Comparison between experimental and calculated apparent pKaps of the compounds studied and (b) plot of pKap1 and pKap2 as a function of the sum of Hammett substituent constants for the AH + form.

Table 1 .
Substituents and experimental thermodynamic constants of flavylium salts at 25 °C (see the structures in Scheme 1); data compiled from reference 24

Table 2 .
Substituents and experimental thermodynamic constants of common anthocyanins.