Oenin and Quercetin Copigmentation: Highlights From Density Functional Theory

Making use of anthocyanin copigmentation, it is possible to effectively improve color quality and stability of red wines and other foods. This can be done by selecting strong copigments, but a 1-fold experimental screening usually entails a high cost and a low efficiency. The aim of this work is to show how a theoretical model based on density functional theory can be useful for an accurate and rapid prediction of copigmentation ability of a copigment. The present study, concerning the copigmentation between oenin and quercetin under the framework of implicit solvent, indicates that, in these conditions, the intermolecular hydrogen bonds play an important role in the system stabilization. The dispersion interaction slightly affects the structure, energies and UV-Vis spectral properties of the copigmentation complex.

Different from experiments, a quantum mechanical (QM) screening allows to get structures and properties of copigmentation systems (Quartarolo and Russo, 2011;Di Meo et al., 2012;Kalisz et al., 2013;Rustioni et al., 2013;Trouillas et al., 2015Trouillas et al., , 2016, such as binding energies and spectral shifts, in a quicker, time-saving and money-saving way. This implies the use of a reliable theoretical model which consists in a good strategy for searching the most stable conformers of copigmentation complex in the conformational space, in an appropriate choice of the QM methods for geometry optimization, energetic and spectral calculations and in a reliable description of solvation effects (Li et al., 2011a,b;Nave et al., 2012;Rustioni et al., 2013;Trouillas et al., 2015Trouillas et al., , 2016Marpaung et al., 2017). Actually, finding out the most stable conformer can be done by calculations totally based on a QM approach, starting from few most probable orientations, or by a preliminary molecular dynamics (MD) simulation followed by a QM refinement (Di Meo et al., 2012;Trouillas et al., 2015Trouillas et al., , 2016. The former way is expected to be more reliable and more time-saving than the latter one as long as reasonable QM methods are utilized and appropriate initial guess of orientations are adopted (Trouillas et al., 2016).
The nature of the non-covalent interaction between anthocyanin and copigment is not quite clear and the suggested π-π stacking configuration may be driven by dispersion forces, hydrogen bonds (HBs), hydrophobic effects, etc (Dimitrić Marković et al., 2005;Kunsági-Máté et al., 2006;Di Meo et al., 2012;Kalisz et al., 2013;Teixeira et al., 2013;Trouillas et al., 2015Trouillas et al., , 2016Zhang et al., 2016). This complexity calls to consider various possibilities for the initial guess, and a careful selection of computational methods. In some previous reports, anthocyanin glycone was substituted by a methyl group in order to explore the general nature of copigmentation and reduce computational cost (Di Meo et al., 2012;Trouillas et al., 2015Trouillas et al., , 2016. However, it is quite reasonable to take the holo glycoside into account when we screen strong copigments against a specific anthocyanin. The hydroxyl-rich sugar moiety ought to have notable impact. The purpose of this paper is to investigate the copigmentation process between malvidin-3-O-glucoside (oenin), normally the highest content anthocyanin in Vitis vinifera young red wines and a determining factor of wine color (Han et al., 2015), and quercetin (see Scheme 1 in Supporting Information), a representative and intensively studied copigment in red wines (Lambert et al., 2011). Our QM study includes conformer search strategy, methods for geometry optimization, energy evaluation, spectral shift and solvation effects estimation.

METHODS AND COMPUTATIONAL DETAILS
According to previous researches based on density functional theory (DFT) (Grimme et al., 2010;Anouar et al., 2012;Di Meo et al., 2012;Trouillas et al., 2015Trouillas et al., , 2016, the hybrid B3LYP-D3 functional combined with the 6-31+G(d) basis set with polarization and diffusion functions on heavy atoms was adopted for geometry optimization of individual oenin and individual quercetin, as well as the conformer search of oenin/quercetin complex. Vibrational frequency calculations were performed with the same method to verify that each obtained structure is a minimum on the potential surface, and to get corrections for thermodynamic functions. Integral equation formalism polarizable continuum model (IEFPCM) was used to describe water solvent effect in all calculations except spectral evaluation (Li et al., 2011a,b;Trouillas et al., 2015Trouillas et al., , 2016.
The binding energies for selected conformers were obtained according to Equation 1: where i stands for oenin or quercetin. The geometries of the complex conformers, oenin and quercetin were optimized individually. In order to explore the impact of functional and basis set on the binding energy, benchmark single point computations on previously optimized geometries were performed with normal hybrid GGA functionals B3LYP-D3(BJ) and B3LYP-D3, range-separated hybrid functionals ωB97X-D and CAM-B3LYP-D3, and meta-GGA functional M06-2X-D3, each of which was combined with 6-311++G(d,p), aug-cc-pVDZ, aug-cc-pVTZ, and aug-cc-pVQZ basis sets. The basis set superposition error (BSSE) (Simon et al., 1996) of complex energy was assessed by counterpoise method and calculated with CAM-B3LYP-D3/aug-cc-pVTZ. Binding energy is the sum of total distortion energy (Equation 2) of oenin and quercetin and their interaction energy (Equation 3) (Scaranto et al., 2011): where "complexed" means the geometry of oenin or quercetin in the complex. Thus, the binding, distortion and interaction energies for 12 selected conformers are estimated. The Boltzmann distribution of the nine most stable conformers was evaluated according to relative Gibbs free energies, compared with the free energy of the optimal conformer. The dispersion contribution to binding and interaction energies of the 12 conformers were calculated with CAM-B3LYP-D3/aug-cc-pVTZ. To explore the CT character, the electronic population analysis was achieved by CHelpG formalism (Breneman and Wiberg, 1990) with CAM-B3LYP-D3/aug-cc-pVTZ for ground states and TD-ωB97X-D/cc-pVDZ for excited states. Pearson correlation analysis was made by IBM SPSS Statistics to investigate the relation between energies and CT.
Functionals, with different HF exchange component, of B3LYP, PBE0, B3PW91, CAM-B3LYP, ωB97X-D, and M06-2X coupled with cc-pVDZ basis set and state-specific PCM (SS-PCM) were employed to predict the spectral shift of the optimal SCHEME 1 | Chemical structures of (A) oenin and (B) quercetin, of which the backbone atoms and rings are numbered.
The visualization of the non-covalent contributions were obtained by using the NCIPLOT software version 3.0 (Johnson et al., 2010;Contreras-Garcia et al., 2011).

RESULTS AND DISCUSSION
Twelve most probable orientations of complex oenin/quercetin were considered (see Figure S1). Orientations 1∼4 show the backbone of quercetin being parallel or antiparallel to that of oenin, like largely suggested by previous studies (Kunsági-Máté et al., 2006;Di Meo et al., 2012;Kalisz et al., 2013;Teixeira et al., 2013;Trouillas et al., 2015Trouillas et al., , 2016. This interaction mode should favor a great extent of electron delocalization between oenin and quercetin, with a good chance of getting a stable conformer. Intuitively, orientations 1 and 2 are beneficial for dispersion interaction and charge transfer (CT), namely the apparent π-π stacking, while orientations 3 and 4 are also favorable for HB interactions between quercetin and hydroxylrich sugar moiety of oenin. Since HB interactions could play a key role, more orientations offering this possibility were taken into examination. In orientations 5∼8, quercetin takes a triangular-shaped arrangement to simultaneously interact with the sugar part of oenin mainly by HB, and with backbone of oenin mainly by HB or dispersion force. In orientations 9∼12, quercetin exhibits as a sandwiched layer located between sugar and backbone segments of oenin, which allows similar or weaker interactions as in the triangular-shaped arrangement.
More importantly, the HB effect is critical to copigmentation complex stability, which has not been highlighted with substantial evidence before (Di Meo et al., 2012;Kalisz et al., 2013;Zhang et al., 2016). Complex 5 is the most stable one although it does not presents the B. . . B stacking and has a shifted and inclined AC. . . AC stacking. Its stability could be attributed to two strong HBs connecting the 6 ′′ -OH of oenin with 3 ′ -OH of quercetin and 3 ′′ -OH of oenin with 7-OH of quercetin, and another weak HB between 7-OH of oenin and 4-OH of quercetin. In addition, a weak interaction between 3 ′ -OCH 3 of oenin with conjugate π-electrons of B-ring of quercetin is observed. In Figure 2, different kinds of non-covalent interactions in complex 5 are indicated. Other than HBs it is possible to note the presence of many van der Walls contributions between the two moieties. So, also these interactions contribute to stabilize the complex.
Similar situation happens in conformers 1∼3 and 6∼9, in which, however, there is insufficient ring overlap. On the contrary, in spite of a good AC. . . AC stacking, 10 and 12 have higher relative energies (Table 1), which should be partly ascribed to weak HBs formed between oenin and quercetin. Oenin 3 ′′ -OH group is most likely to form a strong HB, and 6 ′′ -OH and 3 ′ -OCH 3 may also have great possibility. For quercetin, 3 ′ -OH has an absolute advantage to form a strong HB. Thus, a substitution of any of these groups shall dramatically affect copigmentation stability. As in our case, another recent study found a great influence of HB on the structure of catechol dimer (Barone et al., 2017). Besides, it's worth noting that the glycone of anthocyanin can have a double-effect. On one hand it may 1 | Main parameters for HB and ring-stacking (RS) and relative Gibbs free energies G relative for the examined conformers.

Conformer
HB Distance in Å, angle in degree and energy in kcal/mol. a Parameters are shown only for the strongest hydrogen bond, while others are shown in Figure 1 and Figure S2. D and A stand for hydrogen donor and acceptor, respectively. The HB type is O-H…O for all conformers, except that conformer 5 has another weak C-H…π. b Distance between H and hydrogen acceptor, or between two rings. c Angle of a hydrogen bond, or dihedral of two rings. d A, B and AC denote the A-ring, B-ring and AC-rings, respectively. e The energy of conformer 5 was taken as the reference. enhance the interaction between pigment and copigment via its hydroxyls to form HBs, on the other hand, its steric effect may impede copigment approaching the skeleton of anthocyanin, as uncovered by experiments (Gras et al., 2018). A great structural distortion takes place when oenin and quercetin associates with each other. Table 2 shows that the total torsion angle between AC-rings and B-ring is 23.0 ± 7.4 • for all conformers, including 7.7 ± 4.7 • distortion from oenin and 15.3 ± 10.7 • from quercetin. This mainly comes from the flexibility of glucoside and B-ring, and results in the formation of strong HBs. Thus, intermolecular HBs, causing significant structural distortion of the two monomers, shall further affect the stability and spectra properties of the copigmentation complex.
It is noteworthy that this study is based on the implicit solvent model (Li et al., 2011a,b). An explicit description of the solvent (Trouillas et al., 2015(Trouillas et al., , 2016 could modify the obtained results in that it may change the hydrogen pattern and the relative stability of the conformers.

Energetic Behavior
Binding energy is a key index for evaluation of copigmentation ability of a copigment. The binding energy values of conformer 5 estimated with different functionals and basis sets are listed in Table 3. All the methods predict a prominent binding energy, in agreement with experimental observations (Lambert et al., 2011) and other calculations (Di Meo et al., 2012). The impact of functional is not significant, even if CAM-B3LYP-D3 seems better according to experimental binding enthalpy. For all tested functionals, a larger basis set provides a lower binding energy, implying the necessity of employing high-level basis set to capture the complicated interactions between oenin and quercetin as comprehensively as possible. Still, 6-311++G(d,p) and aug-cc-pVDZ basis sets are good enough to give reasonable results.
Actually, binding energy consists of two parts: distortion of copigmentation system and interaction between anthocyanin and copigment (Scaranto et al., 2011). It is easy to confuse binding energy and interaction energy, of which the latter one really represents the stability of both the complex and the color. The more negative the interaction energy, the more stable the complex and the color. Thus, interaction energy should be taken as the index for copigmentation ability evaluation of a copigment. CAM-B3LYP-D3/ aug-cc-pVTZ approach was used to evaluate the binding, distortion and interaction energies for 12 selected conformers, and the Boltzmann weights (BW) of the nine most stable conformers (see Table 4). As can be seen from values of BW, conformer 5 has overwhelming possibility to exist in real wine solution, due to strong HBs formed between oenin and quercetin as described before. This indicates that neglecting the sugar moiety of oenin, or miss-capturing the optimal conformer assisted by HBs may lead to a higher binding energy, which is much closer to the experimental value but not rational. Noteworthy distortion energies of 0.41∼5.11 kcal/mol for conformers 1∼9 are obtained, in accordance with the structural distortion. This suggests that the use of interaction energy rather than    a Geometry is based on B3LYP-D3/6-31+G(d). Dispersion energies and q GS were calculated with CAM-B3LYP-D3/aug-cc-pVTZ. q ES were computed by TD-ωB97X-D/cc-pVDZ. Eelectronic population analysis was achieved by CHelpG formalism. b "−" means the electron is transferred from quercetin to oenin, vice versa.
binding energy is more appropriate to describe copigmentation stability. Dispersion contribution to binding and interaction energies, as well as intermolecular CT of ground and excited states for 12 conformers are displayed in Table 5. The dispersive interaction contributes −22.01 ± 0.79 kcal/mol to the binding energy and −22.41 ± 0.85 kcal/mol to the interaction energy. This manifests that the difference of dispersion contribution between conformers or between binding and interaction energies is indistinctive. Since the dispersion contribution is of the same magnitude as the binding or interaction energy itself for each conformer, as previous studies revealed ( Kunsági-Máté et al., 2006;Di Meo et al., 2012;Teixeira et al., 2013), one may conclude that the dispersion plays an uppermost contribution to copigmentation. Nevertheless, a Pearson correlation analysis shows that the correlation coefficient between dispersion energy and binding energy is 0.348 (p = 0.268), while it is 0.440 (p = 0.153) between dispersion contribution and interaction energy (see Table 6). This means that the dispersion contribution is significant to neither binding nor interaction energy, just as a similar finding in catechol dimer (Barone et al., 2017). The most likely scenario is that HBs play a more important role to binding or interaction energy, and thus to stability of copigmentation (Kalisz et al., 2013;Zhang et al., 2016;Barone et al., 2017). This is supported by the former structural analysis. To elucidate this issue, a reliable interaction energy decomposition is necessary in a future study.

Spectral Shift
UV-Vis spectral shift is another key index for giving insight to the copigmentation process. The more red-shift of the spectra, the stronger the copigment. The maximum absorption wavelengths, λ max , of oenin and conformer 5, as well as the experimental and predicted spectral shifts estimated with different functionals are listed in Table 7. All functionals underestimate λ max of oenin and conformer 5 compared with experiments (Lambert et al., 2011), as it happens in previous DFT studies (Di Meo et al., 2012;Trouillas et al., 2015Trouillas et al., , 2016. However, it is found that range-separated hybrid functionals CAM-B3LYP and ωB97X-D, and meta-GGA functional M06-2X seem to better describe the electronic excitation of copigmentation complex (Di Meo et al., 2012;Trouillas et al., 2015Trouillas et al., , 2016, since the bathochromic spectral shifts from 7.2 to 9.4 nm are much closer to the experimental value of 8.5 nm (Lambert et al., 2011). Previous evaluation of the spectral shift of 3-O-methylcyanidin complexed with quercetin by ωB97X-D also gave reasonable results (Di Meo et al., 2012). Thus, TD-ωB97X-D/cc-pVDZ combined with SS-PCM solvent model is trustworthy to describe the electronic transition and corresponding spectral shifts of copigmentation complex of oenin and quercetin.
Vertical excitation energies, maximum absorption wavelengths and corresponding shifts, oscillator strengths and molecular orbital (MO) contributions for 12 conformers are exhibited in Table 8. Compared to the spectrum of oenin, each conformer presents an evident bathochromic shift in the range 3.8∼14.5 nm, except a slight and notable hypochromatic shift of 0.7 nm and 8.7 nm for conformers 8 and 11, respectively. Each conformer also has a high oscillator strength. There is a significant dependency (correlation coefficient = 0.781, p = 0.003) between the extent of shift and the charge transfer degree in the ground state. A CT from quercetin to oenin (shown in Table 5) can cause a bathochromic shift. More fundamentally, the CT of equilibrium ground state is determined by multiple factors, of which the HB effect rather than the dispersive interaction should be the key, as a Superscripts of significance (two-tailed) ** and *stand for highly significant (p < 0.01) and significant (p < 0.05), respectively. discussed above under the circumstance of implicit solvent scheme. For all conformers, the maximum absorptions is due to the transition from the second highest occupied molecular orbital (HOMO-1) to the lowest unoccupied MO (LUMO). MO correlation of oenin, conformer 5 and quercetin (see Figure 3), shows that HOMO, HOMO-1 and LUMO of conformer 5 are close to the HOMO of quercetin, and to the HOMO and LUMO of oenin. The HOMO→ LUMO energy gap of the complex is only 5.77 eV, much lower than the gap of quercetin (7.11 eV), and very similar to the gap of oenin. HOMO, HOMO-1, and LUMO of the complex locate mainly on B-ring and C-ring of quercetin, B-ring and AC-rings of oenin. Thus, the HOMO→ LUMO transition corresponds to an intermolecular CT from quercetin to oenin (Di Meo et al., 2012;Trouillas et al., 2015Trouillas et al., , 2016, while HOMO-1→ LUMO refers to an intramolecular CT of oenin. Considering the structure of conformer 5 (Figure 1), B-ring of quercetin almost has no superposition with oenin, and C-ring of quercetin is aslant and translated parallel to AC-rings of oenin with a dihedral of 15.3 • . This consistently decreases the overlap between HOMO and LUMO, and thus diminishes the HOMO→ LUMO contribution to the maximum absorption of conformer 5, while protrudes the HOMO-1→ LUMO contribution. Such an intramolecular excitation does not have an apparent intermolecular CT (see Table 5). However, for  conformers 1, 2, 3, 6, 7, 8, 10, and 12, HOMO-1→ LUMO transition corresponds to an intermolecular excitation that has a quite big CT (Di Meo et al., 2012;Trouillas et al., 2015Trouillas et al., , 2016.

CONCLUSIONS
In summary, copigmentation of oenin with quercetin was carefully investigated by using density functional theory. From our study the following conclusions can be drawn: -Among the 12 considered initial copigmentation orientations, 10 are characterized by HB interactions between quercetin and glucoside of oenin. Quercetin assumes a (anti)parallel, triangular-shaped, or sandwiched form located between sugar and backbone segments of oenin. The above 10 orientations include the most stable conformer, which is stabilized by two strong HBs and other two weak interactions. Thus, the intermolecular HB can be considered the driving force of copigmentation process between oenin and quercetin when FIGURE 3 | Molecular orbital correlation diagram of oenin, conformer 5 and quercetin. A tube molecular representation was adopted for oenin in red and quercetin in cyan.
implicit solvent effect is considered. But it may be different if specific water molecules are taking into account. -Based on our results it appears that a reliable description of the copigmentation process should comprise the sugar moiety of anthocyanin to allow to take into account the different HB interactions that can occur between hydroxyl-rich glycone of oenin and quercetin. -Grimme dispersion correction are mandatory to describe dispersive interaction between pigment and copigment. -CAM-B3LYP-D3 in conjuction with IEFPCM and 6-311++G(d,p) or aug-cc-pVDZ basis set is good enough to provide reasonable binding and interaction energies. -TD-ωB97X-D/cc-pVDZ partnered with SS-PCM is able to afford trustworthy description of the electronic spectral properties of the copigmentation complex.
The strategy which we have used in this case could further be applied to estimate other copigmentation systems. In this way, it will be possible to screen strong copigments from large-scale samples efficiently, and sequentially improve color quality and stability of red wines and other foods.

AUTHOR CONTRIBUTIONS
All authors listed have made a substantial contribution to the work and approved its publication. YL designed the protocol, made calculations and wrote the paper. MP made calculations. MT and NR designed the protocol and wrote the paper.