Oenin/Syringic Acid Copigmentation: Insights From a Theoretical Study

On the basis of the dispersion-corrected density functional theory, a computational model is proposed to describe the oenin/syringic acid copigmentation and to explore the non-covalent interaction between the anthocyanin and the copigment in the framework of implicit solvent approach. The predicted binding free energy and visible spectrum shift of this copigmentation complex are in accordance with the experimental observations. The used model provides a good structural description of oenin/syringic acid complex and suggests that the intermolecular hydrogen bonding, in which the hydroxyl-rich sugar moiety in oenin plays a key role, may be the determinant for the formation and nature of the copigmentation complex.

Quantum mechanics (QM) screening is able to provide microscopic interactive conformation, spectrum and binding free energy of any copigmentation system (Quartarolo and Russo, 2011;Di Meo et al., 2012;Kalisz et al., 2013;Rustioni et al., 2013;Trouillas et al., 2016;Khalifa et al., 2018;Bayach et al., 2019;He et al., 2019), and is generally less moneyand time-consuming compared with experimental approaches. The QM screening calls for a robust theoretical model composed of an efficient conformer-scanning strategy in search of the copigmentation conformers with the lowest energy in conformational space, appropriate algorithms for structural optimization, energetic and spectral evaluation, and solvent effect description (Li et al., 2011a,b;Nave et al., 2012;Rustioni et al., 2013;Trouillas et al., 2016;Marpaung et al., 2017). The seek of the most stable copigmentation conformers could be achieved by sequential molecular dynamics (MD) simulation and QM refinement (Di Meo et al., 2012;Trouillas et al., 2016), or totally by QM calculations with several preferential orientations as the starting point. The latter way may be less timeconsuming than the former one if suitable QM approaches and reasonable initial guess of preferential orientations are adopted (Trouillas et al., 2016).
The non-covalent nature of the interactions between pigments (anthocyanins) and copigments is still controversial and lacks of solid evidences. Many studies indicate that the π-π stacking interaction is driven mainly by dispersion forces, and strengthened by 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;Zhang et al., 2015;Trouillas et al., 2016). In some previous theoretical studies the glycosyl group in anthocyanins is replaced by a methyl group so as to uncover the chemical nature for copigmentation with a lower cost (Di Meo et al., 2012;Trouillas et al., 2016). However, you can't rule out that polyhydroxyl sugar moiety in anthocyanins, which could have appreciable impact upon copigmentation.
Following our computational strategy used in another recent work on oenin/quercetin copigmentation (Li et al., 2018a), we thought it interesting to explore the copigmentation of syringic acid (Zhang et al., 2015) (see Figure 1) and malvidin-3-Oglucoside (also known as oenin) that is the most concentrated anthocyanin in young red wines obtained from Vitis vinifera grapevine varieties.

METHODS AND COMPUTATIONAL DETAILS
The first step to study the oenin/syringic acid copigmentation is to acquire the optimal conformer with the lowest energy of the complex, which is the footstone of the following energetic and spectral calculations. For this purpose, nine most preferable interaction orientations for the copigmentation complex of oenin/syringic acid were taken into account ( Figure S1) (Li et al., 2018a). In orientations 1∼4, just as previous studies suggested ( Kunsági-Máté et al., 2006;Di Meo et al., 2012;Kalisz et al., 2013;Teixeira et al., 2013;Trouillas et al., 2016), the syringic acid backbone is parallel to the oenin backbone. Among these orientations some seem to be suitable for charge transfer (CT) and dispersion interactions (1∼2), whilst others (3∼4) should facilitate the formation of HBs between polyhydroxyl sugar moiety in oenin and syringic acid. Orientations 5∼9, in which the embedded syringic acid is able to interact simultaneously with the polyhydroxyl sugar segment in oenin by HB and with the oenin backbone by dispersion force or HB, were also examined.
For each orientation, a potential energy curve was built up performing a relaxed scan along Z-direction (varying the distance between the planes of the rings involved in the stacking from 2.6 to 4.6 Å with an interval of 0.2 Å) to get a conformational minimum. The conformational minima thus acquired were subjected to further relaxed scan in the XY-plane (varying the angle of reciprocal orientation of the rings involved in stacking from 0 to 360 • with an interval of 10 • ). In such a way, hundreds of conformers were examined to obtain the most stable conformers.
As in previous studies performed at density functional level of theory (DFT) (Anouar et al., 2012;Di Meo et al., 2012;Trouillas et al., 2016), geometry optimizations for individual syringic acid (both protonated and deprotonated form), oenin and their complexes were carried out by employing the hybrid functional with Grimme dispersion correction B3LYP-D3 (Grimme et al., 2010) in connection with the 6-31+G(d) basis set, with diffusion and polarization functions on the heavy atoms (Li et al., 2018a). The followed vibration frequency analysis was accomplished by the same approach to confirm that the optimal structures obtained are potential surface minima and to acquire corrections for zero point energy, thermal energy, enthalpy, and Gibbs free energy thermodynamic functions. Polarizable continuum model with integral equation formalism (IEFPCM) was adopted to account for the solvent effect except than in the spectra evaluation (Li et al., 2011a,b;Trouillas et al., 2016).
The binding thermal energy E (and similarly the binding enthalpy H and binding Gibbs free energy G) (Li et al., 2018b) for the selected conformers were obtained on the basis of Equation 1: where i denotes syringic acid (both protonated and deprotonated form) or oenin. The binding entropy S is determined by G= H-T S at 293K. The geometries of oenin, syringic acid and their complex were individually optimized. For the examination of the influence of basis sets and functionals, single point computations on previously optimized geometries were performed using B3LYP-D3, ωB97X-D, B3PW91-D3, CAM-B3LYP-D3, M06-2X-D3, and PBE0-D3 functionals, combined with cc-pVDZ, cc-pVTZ, aug-cc-pVDZ, 6-311++G(2d,2p), 6-311++G(d,p), 6-311G(d,p), 6-311++G and 6-311G basis sets. The basis set superposition error (BSSE) (Simon et al., 1996) for interaction energy was evaluated by the counterpoise method with CAM-B3LYP-D3/aug-cc-pVDZ. The distortion energy E distortion for the most representative conformers was obtained applying the Equation 2: where "complexed" stands for the geometry of syringic acid or oenin in copigmentation complex. The Boltzmann weights for conformers were assessed by relative binding Gibbs free energies. The contribution arising from dispersion interaction to binding energy was calculated at B3PW91-D3/aug-cc-pVDZ level. The electronic population analysis was implemented by the CHelpG formalism (Breneman and Wiberg, 1990) to capture the CT characteristic of the complexes, using B3PW91-D3/augcc-pVDZ and TD-CAM-B3LYP-D3/cc-pVDZ for ground and excited states, respectively. Normal hybrid functionals B3LYP-D3, PBE0-D3 and B3PW91-D3, range-separated hybrid functionals CAM-B3LYP-D3 and ωB97X-D, and meta-GGA functional M06-2X-D3, combined with cc-pVDZ basis set and state-specific polarizable continuum model (SS-PCM) were utilized to evaluate the spectral shift for the optimal conformer (Simon et al., 1996;Yanai et al., 2004;Chai and Head-Gordon, 2008;Di Meo et al., 2012;Trouillas et al., 2016). Since CAM-B3LYP-D3 turned out to be one of the best performing functional for the optimal conformer, it was further employed to estimate the spectral shift for other ten conformers. All the calculations were accomplished by Gaussian09 packages (Frisch et al., 2009).

Structure Feature
Starting from the orientations depicted in Figure S1, eleven optimal conformers were achieved by exploring the conformational space of the copigmentation complex. A threeview drawing of the most stable conformer 1 is reported in Figure 2. Figure S2 shows all other conformers. All interesting structural parameters are reported in Table 1. The stacking form includes translated-parallel (1, 3, 4, 8, 9), also called "parallel-displaced" (Di Meo et al., 2012), aslant-parallel (2, 5, 6, 7, 11) having a dihedral over 10 • . In the complex (10) there is no stacking although aslant-parallel form for this conformer was reported before (Li et al., 2018a). The distance between the backbone planes of the pigment and copigment molecules falls in the range 3.21∼3.51 Å for the translated-parallel stacking forms and 3.11∼5.59 Å for the aslant-parallel ones. The dihedral angle of the backbone planes of the pigment and copigment molecules is 1.57∼9.38 • and 11.12∼43.74 • for the translated-parallel stacking and for the aslant-parallel stacking forms, respectively. The presence of aslant-parallel (2, 5, 6, 7, 11) ring stacking entails the decrease of complex stability.
The HB interaction appears to be of great importance for the copigmentation complex stability in spite of that previous studies has not been able to provide solid support (Di Meo et al., 2012;Kalisz et al., 2013;Zhang et al., 2015;Li et al., 2018a). Conformer 1 is the most stable one although it presents a weak stacking between the C-ring of oenin and S-ring of syringic acid, as in the case of the oenin/quercetin copigmentation complex (Li et al., 2018a). Its stability could be attributed to the presence of intermolecular HBs (one O-H. . . O and other two C-H. . . O), as shown in Figure 2 and  (Johnson et al., 2010) in complex 1 can be found in Figure 3. The three strong intermolecular HBs, as well as several others of lesser importance, are recognizable in the blue isosurfaces with spike value of sign(λ 2 )ρ between−0.024 and−0.005 (Johnson et al., 2010). The strongest spike at ca.−0.024 is connected with the strongest O-H. . . O HB. In the same figure, some weak van der Walls interactions between the oenin and syringic acid partners that contribute to the complex stabilization, appear in the green isosurfaces with spike value of sign(λ 2 )ρ from−0.005 to 0.005 (Johnson et al., 2010).
HBs also happen in conformers 2∼7 and 10∼11 where the stacking is scarce. Conformer 8 and 9 present relative binding FIGURE 2 | Front (A), side (B), and top (C) views of the most stable conformer 1 with a tube molecular representation. Carbon atoms are colored in green for oenin and in cyan for syringic acid. Oxygen and hydrogen (involved in hydrogen bonds) atoms are depicted in red and blue, respectively. HBs, presented in blue dashed lines and numbered, are exhibited with key parameters. The stacking distance and dihedral between oenin rings and syringic acid ring are also included. A fogging depth-cueing is used to improve perception.

Conformer
HB (D-H…A) a RS (Oenin…syringic acid)  Figure 2 and Figure S2. D and A stand for hydrogen donor and acceptor, respectively. b Distance between H and A, or between two rings. c Angle of a hydrogen bond, or dihedral of two rings. d B, C, AC, and S denote the B-ring, C-ring and AC-rings of oenin and S-ring of syringic acid, respectively.
free energy values that are positive (see Table 3), although they possesses a good stacking of B. . . S or AC. . . S and HBs. Since, in oenin, the 3 ′′ -OH, 5-OH and 7-OH hydroxyl groups can form strong HBs and in syringic acid, 1-COOH and 4-OH groups can do the same, we suppose that these interaction can contribute significantly in determine the stability of copigmentation complexes. This assumption is supported in a recent work in which the impact of HBs on the formation and stability of a catechol dimer (Barone et al., 2017) was investigated. Other contributions to the stability of copigmentation complexes can derive from glycosyl group in oenin despite it is well-known that it can play a double role. In fact, with its many hydroxyl groups it may strengthen the interactions between oenin and the copigment through the formation of HBs, but its steric hindrance can hinder the electronic conjugation between the copigment and oenin (Gras et al., 2018).
As showed in Table S1, an obvious structure distortion arises when syringic acid and oenin conjugate with each other (Li et al., 2018a). These distortions are mainly due to structural flexibility of B-ring and glucoside in oenin   Zhang et al. (2015).
and of 1-COOH group in syringic acid, leading to the generation of HBs. Thus, not only the intermolecular HBs cause distortion of the partners' structure, it may further make important contributions to the stability and spectral behavior of the complex. It is worth mentioning that the copigmentation complex is embedded in aqueous solution, whose influence in this work (Li et al., 2011a,b) is accounted by implicit solvent model approach. An explicit description of the solvent (Trouillas et al., 2016) may affect the HBs' strength and the conformers' stability.

Energetic Feature
Thermodynamic aspects of the binding process are important to evaluate the copigmentation ability of a copigment. The binding energy, binding enthalpy, binding Gibbs free energy, and binding entropy for conformer 1 computed using various functionals and basis sets are presented in Table 2 and Table S2, together with available experimental values. As can be observed, all methods we have used to make calculations overestimate strongly both binding enthalpy and entropy with respect to the measured values. Instead, the binding free energy seems to be always quite well reproduced. This problem concerns also the theoretical determination of Zhang et al. (2015) obtained through a MD simulation to sample conformers of oenin/syringic acid in explicit water solvent followed by a QM estimation of the energies.
B3PW91-D3/aug-cc-pVDZ protocol was employed to determine the Boltzmann weights (BW) for the 11 selected conformers. Results are listed in Table 3. The binding free energy of conformer 1 is 0.95∼3.09 kcal/mol lower than the other conformers, thus, it turns out to be the one with the largest possibility of existing (BW = 61%). The evaluation of Pearson correlation coefficients indicate that the difference of the binding free energy between conformers mainly comes from the binding enthalpy rather than from binding entropy [0.817 (p = 0.002) vs. 0.414 (p = 0.206)]. As already mentioned, the dispersion forces and HB interactions may give important contributions to the stability of the copigmentation complexes. As shown in Table 4, the dispersion energy contribution is−21.72±2.68 kcal/mol. For the HB contribution, it is not easy to get an accurate evaluation. However, as recommended by Barone et al. (2017), the influence  As demonstrated in Figure S3, the rotation starts from the geometry of conformer 1 (ψ = 195 • ), fixing all atoms except the hydrogen atom in the strongest O-H···O intermolecular HB. The rotation is also designed to bring about a minimal perturbation to other interactions, such as dispersion and repulsion interactions. Ten geometries, from ψ = 195 • to ψ = 120 • , were sampled along with the rotation trajectory. The variations of the binding energy, binding enthalpy, binding free energy and dispersion contribution values are listed in Table 5, while the change of the non-covalent interactions are illustrated in Figure S3. Plotting the reduced density gradient vs. the electron density ρ multiplied by the sign of the second Hessian eigenvalue sign (λ 2 ) facilitates to uncover the types and strength of non-covalent interactions, as detailedly introduced by Johnson and coworkers (Johnson et al., 2010). Very low density values (i.e.,−0.005 a.u. < sign (λ 2 )ρ < 0 a.u.) generally map to weak dispersion interactions; while higher density values (i.e.,−0.05 a.u. < sign (λ 2 )ρ < −0.005 a.u.) map to stronger HB interactions. In our case, the spike value of sign(λ 2 )ρ changes from−0.027 to−0.012 and the isosurface color of the HB varies from blue to light green when ψ changes from 180 • to 120 • , suggesting its gradual breaking of HB interaction. Along with the HB weakening, the binding energy, binding enthalpy and binding free energy rise rapidly from−18.32 to 14.67 kcal/mol, from−18.91 to 14.08 kcal/mol and from−2.73 to 30.25 kcal/mol, respectively, while the dispersion contribution almost remain unchanged, from−24.77 to−23.84 kcal/mol. Although an exact value of the HB contribution cannot be given one can still see that the formation of this strong intermolecular interaction influences the binding free energy of conformer 1, similarly to what occurs for catechol dimer (Barone et al., 2017) and in the copigmentation process between epicatechin and pelargonidin-3-O-glucoside (Zou et al., 2018).
The HBs and dispersion interactions cause structural distortions in both oenin and syringic acid, which are accompanied by significant CT. As shown in Table 4, the distortion energy is−2.59 ± 3.20 kcal/mol, and the CT for ground and excited states is−0.05 ± 0.05 |e| and−0.16 ± 0.18 |e|, respectively. The distortion and CT are related to the sugar moiety of oenin, which is key for HBs assisting the stability of the copigmentation complex. 6 | Vertical excitation energy (E max ), maximum absorption wavelength (λ max ), spectral shift ( λ max ), oscillator strength f and MO contribution (%) for the 11 conformers (energy and wavelength are in eV and nm, respectively) a . a Geometry is based on B3LYP-D3/6-31+G(d) computational scheme. The spectra were calculated with SS-PCM, TD-CAM-B3LYP-D3/cc-pVDZ computational protocol. b Compared to λmax(oenin).

Conformer
FIGURE 4 | Molecular orbital correlation diagram of oenin, conformer 1 and syringic acid. A tube molecular representation was adopted for oenin (red) and syringic acid (cyan).

Spectral Feature
Visible spectrum shift along with the copigmentation process is another important index to evaluate the copigmentation ability of a copigment. The wavelength of maximum absorption (λ max ) for conformer 1 and oenin, and the corresponding shifts λ max evaluated by various functionals are showed in Table S3. It seems that, compared with experiments (Zhang et al., 2015), each functional underestimates λ max both for conformer 1 and oenin, likewise to other DFT studies (Di Meo et al., 2012;Trouillas et al., 2016). In spite that the functionals of B3LYP, B3PW91 and B3P86 predict relatively closer λ max of oenin to the experimental value, they failed to describe the spectral shifts reasonably (Oliveira et al., 2017). In previous reports, although the functionals of B3LYP and B3P86 were found to perform particularly well at predicting λ max for a series of flavonoids (Trouillas et al., 2015) and anthraquinones (Anouar et al., 2014), ωB97X-D seems better to account for the spectral shift when quercetin forms copigmentation complex with 3-O-methylcyanidin (Di Meo et al., 2012). In our case, RSH functionals of CAM-B3LYP-D3 and ωB97X-D, as well as M06-2X-D3 appear to be better than other functionals for determine the electronic excitation of the copigmentation complex. Thus, it was considered that TD-DFT (ωB97X-D or CAM-B3LYP-D3)/cc-pVDZ integrated with the implicit solvent model SS-PCM could be reliable enough to characterize the ultrafast electron transition and relevant spectral shift for the copigmentation of syringic acid and oenin. The vertical excitation energy E max , the relevant λ max and λ max , the corresponding oscillator strength f and molecular orbital (MO) description for each conformer are listed in Table 6. With respect to the λ max in oenin, the conformers exhibit notable bathochromic shift of 2.0∼14.8 nm, excluding a little hypsochromic shift for conformer 7 or 10. A weighted mean of 3.9 nm of the bathochromic shift is basically in line with the experimental value of 14.0 nm (Zhang et al., 2015). Strong oscillator strengths along with the electron transition were obtained for the conformers.
The λ max could be attributed to the electron transition from the highest occupied MO, i.e., HOMO, to the lowest unoccupied TABLE 7 | Binding energies ( E) and its dispersion contribution ( E disp ), binding enthalpies ( H), binding Gibbs free energies ( G), binding entropies ( S), relative binding Gibbs free energies ( G) and Boltzmann weights (BW) for selected conformers of the complex of oenin/deprotonated syringic acid (energies in kcal/mol, entropies in J/(K•mol) and BW in %). MO, i.e., LUMO, for all conformers except 10. As illustrated in Figure 4, the MO correlation analysis reveals that the frontier MOs in conformer 1 are closer to the corresponding orbitals in oenin rather than to those in syringic acid. The energy gap between LUMO and HOMO in conformer 1 is 5.01 eV, obviously lower than the gap of 7.31 eV in syringic acid but almost the same as the gap in oenin. More specifically, the LUMO and HOMO in conformer 1 mainly distribute on ACrings and B-ring of the oenin segment, respectively. In other words, the transition from HOMO to LUMO should have an intramolecular CT in oenin (Di Meo et al., 2012;Trouillas et al., 2016), accompanied with different degrees of intermolecular CT (-0.05±0.05 |e| for the ground states and−0.16±0.18 |e| for the excited states) (see Table 4). Actually, the bathochromic shift is highly correlated (correlation coefficient = −0.850, p = 0.004) with the intermolecular CT in excited states.

Deprotonated Syringic Acid
Syringic acid is able to exist either in a neutral protonated state, i.e., carboxyl form, or in an anionic deprotonated state, i.e., carboxylate form. The reported dissociation constant pKa 1 is 4.0 (Almeida et al., 2014), 4.34 (Wang et al., 2006), 4.33 (Hyder and Jönsson, 2012) or 4.30∼4.47 in 14%vol wine (Erdemgil et al., 2007). Thus, the dominant form shall be protonated state for pH < 2, deprotonated state for pH > 6, otherwise, the two states tend to coexist. Since a normal pH of red wines is around 3.5, syringic acid shall mainly exhibit in the protonated state under such a condition, which is just what we studied and discussed above. However, in other circumstances, syringic acid may possess a deprotonated form to associate with anthocyanins. Therefore, conformers 1, 2, 5, 9 and 10 were selected and their syringic acid carboxyls were mutated into carboxylate anions. Geometry optimization and frequency analysis were carried out to obtain corresponding conformers of oenin/deprotonated syringic acid, marked as 1 ′ , 2 ′ , 5 ′ , 9 ′ , and 10 ′ , respectively. Afterwards, thermodynamic energies and spectral properties of these conformers were calculated employing the methods mentioned in section Methods and Computational Details.
The HB parameters, energetic and spectral features of 1 ′ , 2 ′ , 5 ′ , 9 ′ , and 10 ′ were collected in Table S4, Tables 7, 8, respectively. On the whole, the deprotonation has not caused major changes in structure but in some specific interactions between oenin and syringic acid, e.g., the HB interactions. From the HB parameters standpoint, the deprotonation of syringic acid appears to strengthen the HBs in 1 ′ , 2 ′ , and 5 ′ , compared with those in 1, 2, and 5. This is also confirmed by the binding free energies in Table 7, which exihibits a much higher affinity between oenin and deprotonated syringic acid. The anionic carboxylate plays a major role for the enhancement of the HBs, especially in conformer 2 ′ , which makes it 3.37 kcal/mol of the binding free energy more stable than conformer 1 ′ . For the optical properties, the deprotonation of syringic acid is inclined to hypochromatic shift the maximum absorption of the spectra.

CONCLUSIONS
To sum up, a study was carried out to clarify the copigmentation process between syringic acid and oenin, on the basis of dispersion-corrected DFT computations and using implicit solvent model. Conclusions can be made as follows: -Eleven preferable conformers for the copigmentation complex oenin/syringic acid were selected after an accurate sampling of the conformational space. Among them, nine have HB interactions between syringic acid and the glycosyl group in oenin. The stabilization of the lowest lying energy conformer was proved to be due in large part to the intermolecular HBs interactions. -It was found that the dispersion forces represent an important contribution to the intermolecular interactions which govern the formation of copigmentation complex oenin/syringic acid, thus, a Grimme dispersion correction is compulsory to characterize copigmentation process. On the basis of our results we think we can suggest using B3PW91-D3, M06-2X-D3, or ωB97X-D functionals in conjuction with IEFPCM and aug-cc-pVDZ for the calculation of the thermodynamic properties, and the combination of TD-DFT (ωB97X-D or CAM-B3LYP-D3)/cc-pVDZ with SS-PCM the for the characterization of the electron transition and spectral feature. -The polyhydroxyl sugar moiety in oenin appeared to be significantly involved in the formation of intermolecular HBs, therefore, it is necessary to take the holo glycosyl group in anthocyanin into account to obtain an accurate description of copigmentation.
The present findings are obtained under the scheme of implicit solvent but it is important to remember that the description with an explicit solvent model could present significative differences. On the basis of obtained results we think that the computational approach used in this work may be quite suitable to a good description of copigmentation processes. It may further be applied to estimate other copigmentation systems. In such a way, it will be possible to screen strong copigments from large-scale samples and, sequentially, to improve the stability of natural food pigments, which is an important subject for food quality and safety.

DATA AVAILABILITY
The datasets generated for this study are available on request to the corresponding author.