A theoretical antioxidant pharmacophore for natural hydroxycinnamic acids

The structure-activity relationship analysis has been performed for transand cis-hydroxycinnamic acids, to determine their theoretical antioxidant pharmacophore. Based on the detailed conformational studies, the most stable rotamers have been selected. We have analyzed the descriptors of four antioxidant mechanisms important in free radical scavenging: hydrogen atom transfer, sequential proton loss electron transfer, single electron transfer – proton transfer and transition metal chelation, based on the B3LYP/6-311++G(2d,2p) calculations in vacuum and polar media. The results explain the activity difference between cinnamic acid and its derivatives. The descending order of antioxidant potential is as follows: caffeic > sinapinic ~ ferulic > p-coumaric > o-coumaric > m-coumaric ~ phenol. The results have shown that transisomers indicate higher reactivity than cisand may be considered as good antioxidants. It has been determined that the highest antioxidant ability is related to the hydroxyl group in para position, supported by planar structure and stability of radical forms. π-Type delocalization of unpaired electron on aromatic ring, double bond and para O-atom is the key to radical stabilization. The ortho-dihydroxy substitution in benzene ring positively influences the ability to neutralize free radicals and makes caffeic acid the antioxidant pharmacophore of hydroxycinnamic acids.

Phenolic antioxidants were proved to interact with free radicals according to the four mechanisms [17-28]:

TMC (Transition Metals Chelation)
The mechanism is connected to the Fenton Reaction 1.6. According to the TMC, each molecule that can dissociate to anion and proton has the ability to chelate heavy metals, especially polyphenols. The ability of the compound to donate a proton is calculated due to the fact that metal chelation often occurs with the presence of deprotonated hydroxyls in the polyphenols (Reaction 1.7). HOMO and LUMO distributions, SOMO distribution and spin density (SD) of molecules are also related to the ability of scavenging free radicals. Low HOMO energy values are related with weak abilities to donate protons and provide information about the parts of molecule, which are very prone to free radical attack. The difference between HOMO and LUMO energy indicates the chemical activity of the compound. Low ∆E(LUMO-HOMO) values are related with higher chemical activity [29]. SOMO is related to the distribution of non-paired electrons. Together with SD values, it determines the stability of the radical form. Numerous theoretical studies dedicated to some CiA derivatives have been performed [17,[29][30][31][32][33][34][35][36][37][38][39]. Table 1 summarizes which structural and functional parameters have been calculated. The calculated values of these parameters are also provided (Table 1). Till now, the properties of hydroxycinnamic acids have been studied using different functionals (AM1, HF, B3LYP and M05-2X) and basis sets, e.g. 6-31+G(d); 6-311+G(d,p); 6-311++G(3df,2p), which makes the comparative analysis very difficult. Moreover, not all antioxidant descriptors are available for cinnamic acid derivatives in the literature. In this paper, we provide comprehensive analysis of structure -antioxidant activity relationship for both isomers of CiA, as well as for its derivatives such as: o-, m-, p-coumaric, caffeic, ferulic and sinapinic acids. The DFT method is used to calculate bond dissociation enthalpy (BDE), proton affinity (PA), adiabatic ionization potential (AIP), highest occupied molecular orbital, (HOMO), lowest unoccupied molecular orbital (LUMO), singly occupied molecular orbital (SOMO), spin density (SD), stabilization energy (ΔE ISO ), deprotonation energy (∆H acidity , ∆G acidity ), second order perturbation energy E(2) and NMR proton shift. The goal of our study is to describe molecular features of hydroxycinnamic acids, which are necessary to obtain the highest antioxidant potential. Determining the pharmacophore will be beneficial for selection and design of novel and more potent antioxidants.

Computational details
The calculations were performed with the Gaussian 09W computational package [40]. All structures of trans-and cis-isomers were submitted to HF/6-31G(d,p) geometry conformational search. Dihedral angles α = C 5 -C 6 -C 7 -C 8 (describing the rotation of benzene ring around C 6 -C 7 single bond) and β = C 7 -C 8 -C 9 -O 9` (describing the rotation of carboxyl group around C 8 -C 9 single bond) were scanned to construct the potential energy surface scans (PES Scans). Moreover, the rotation of hydroxyl and methoxyl groups around corresponding C Y -O Y` bonds was also analyzed (Fig. 1). The conformers with the lowest HF/6-31G(d,p) energies were fully optimized in their ground states in vacuum using the Berny algorithm and tight convergence in the SCF, employing the DFT method with B3LYP hybrid functional [41-43] and 6-311++G(2d,2p) basis set. The force constants were computed at the first point using the current method. For neutral and monoanionic forms the restricted level of theory was applied while the unrestricted to radicals and radical cations. The obtained electronic geometries indicated the approximate location of the energy minimum structures, confirmed by the absence of imaginary harmonic vibrational frequencies. Calculations conducted in polar media: water, ethanol and DMSO were performed according to the same methodology employing C-PCM solvation model (conductor-like polarizable continuum model) [44], which allows to produce good quality results for polar media, when calculations of antioxidant parameters are concerned. The solute cavity was built up using radii from the UFF force field, and the standard value of the dielectric constant (ε) for water (ε = 78.3553), ethanol (ε = 24.852) and DMSO (ε = 46.826). Natural Bond Orbital (NBO) analysis implementation to the Gaussian 09W computational package has been used for calculation of second order perturbation energies E(2) [45]. Spin-spin coupling constants [46], using Gauge-Independent Atomic Orbital (GIAO) method [47], were calculated with B3LYP/6-311++G(2d,2p) to obtain NMR proton shifts.

Geometry optimization
To understand the relationship between the molecular structure and antioxidant activity of cinnamic acids, the conformational and electronic features have been studied. A complete geometry optimization was carried out for cis-and trans-isomers of cinnamic (CiA)  Kinetic studies on radical scavenging activity M05-2X/6-311+G(d,p); water and pentylethanoate [39] with previously published data of geometry optimization [17,[29][30][31][32]. The following structural parameters were investigated (according to [36]), to study the effect of conformational and rotational changes on the overall stability of the compounds: (1) rotation around C 6 -C 7 and C 8 -C 9 bonds (for all compounds); (2) syn and anti geometries represented by C 7 -C 8 -C 9 -O 9`` dihedral angle equal to 0º or 180º, respectively (for all compounds); (3) orientation of H 9` atom relative to the carbonyl group, resulting in s-cis or s-trans conformation (for all compounds); (4) orientation of the hydroxyl and methoxyl ring substituents, characteristic for the compound, determined by rotation of the respective C-O bonds. Fig. 1 shows the graphical representation of the parameters studied. For the further calculations, geometries of the most stable conformers were used. For all the compounds analyzed, the most stable structures indicate vinyl bond and carboxylic moiety coplanar to the aromatic ring ( Fig. 2), which is confirmed by 0º value of dihedral angles corresponding to rotation around C 6 -C 7 and C 8 -C 9 bonds. Moreover, these bonds are partially double (data not shown), which is an effect of π-electron delocalization between the aromatic ring, vinyl bond and carboxylic moiety. The lowest energy conformers, both trans-and cis-, have syn geometries and s-cis orientation of H 9` atom, with favored values of 0º for C 7 -C 8 -C 9 -O 9`` and O 9``-C 9 -O 9`-H 9` dihedral angles.
Regarding the orientation of the hydroxyl and methoxyl ring substituents, the most favorable rotamers are coplanar to the aromatic ring ( Fig. 2), which is reflected by the 0º value of respective dihedral angles -C X -C Y -O Y`-H Y` and C X -C Y -O Y`-C z . In case of KA, FA and SA it allows the formation of a stabilizing intramolecular interactions: for FA and SA ( Fig. 1). However, the contribution of these intramolecular hydrogen bonds into the molecule stabilization is small. To estimate more precisely the nature of proper hydrogen bonds, we have performed detailed NMR proton shifts and NBO analysis with B3LYP/6-311++G(2d,2p). Table 3 compiles donor-acceptor interactions, proton shifts and second order perturbation energies E(2) [45]. The values of the NMR proton shifts in vacuum and all media studied reveal minimal influence of electrostatic effect on the molecular stability and suggest that intramolecular hydrogen bonds are very weak. The interactions between the first lone pair of the oxygen and the O-H antibonding (σ*) orbital (Table 3) are also small. Based on the E(2) values, we can estimate that corresponding hydrogen  The linear regression analysis between the most stable optimized conformers (obtained with B3LYP/6-311++G(2d,2p)) and experimental data was performed, to determine the reliability of calculations. Most of the crystal structures correspond to the less stable conformers because of rotation of hydroxyl bonds and carboxylic moiety, with exception of t-m-CA and t-FA. Hence, less stable conformers were also included in the validation process. The geometries were superimposed on the crystal structures (with no changes in the bond lengths and angle values). The root-mean square deviation (RMSD) of the overlay between equivalent atoms was calculated (Table 4). Comparing results for the compounds in vacuum with experimental data, we do not observe significant differences because the RMSD values are very low. The increase in RMSD for the most stable conformers is connected to the bond`s rotation.

Bond dissociation enthalpy (BDE)
Bond dissociation enthalpy of hydroxyl groups (BDE OH ) is related to the hydrogen atom transfer. HAT is one of the most commonly studied antioxidant mechanisms. BDE describes the thermodynamic stability of the OH bond and, at the same time, the susceptibility to homolytic fission. The compounds with the lowest BDE OH are more active. The results of BDE OH calculations are shown in Table 5 and Supplementary Table 1. Phenol was used as a reference system, to determine the reliability of our calculations. Comparing gas and polar phase BDE values, no dramatic differences were found (all were less than 10 kJ mol -1 ). However, the BDEs were lower in vacuum than  Table 1). Polar solvents make charge separation easier and it increases with the higher solvent polarity. Thus in gas phase, homolytic cleavage will be preferred while heterolytic cleavage in polar phase. The values obtained are in good agreement with the experiment [61-67].
As can be seen from Table 5, molecules with the reactive OH in ortho or para position have gas-phase BDE OH values decreased in relation to phenol, in the following order: p-CA > o-CA > FA ~ SA > KA. In polar media, a similar indication can be found, but some changes in the order can be noticed: p-CA ~o-CA > FA > SA > KA. Compounds with p-OH have t-isomers more active than c-isomers according to HAT mechanism. The unpaired electron is distributed over whole structures. The presence of an additional functional group in ortho position is advantageous for antioxidant properties and decreases BDE OH in vacuum, as well as in polar media. The o-hydroxylation (in KA) is more favorable than o-methoxylation (in FA and SA). Moreover, homolytic fission of p-OH bonds gives compounds with more resonance structures, which is in opposition to the hydrogen abstraction in meta position. It explains the results for m-CA and KA (with reactive meta-OH). Both isomers of m-CA are less active than phenol, which suggests that this is not a good antioxidant. The meta-OH group in KA is not preferred in HAT mechanism, because 4`-O-KA radical is less stable than 3`-O-KA radical (ΔE=25 kJ mol -1 in vacuum) due to resonance concentration only on the benzene ring.
Our results are in accordance with the work of 25], where trans-caffeic acid and other hydroxycinnamic acids were mentioned as excellent free radical scavengers by H atom donation, alongside with hydroxytyrosol, gallic acid, epicatechin

Proton Affinity (PA)
Proton affinity of hydroxyl groups (PA OH ) is a numerical descriptor related with two-step antioxidant mechanism named sequential proton loss electron transfer (SPLET). PA describes the thermodynamic stability of OH bond, as well as the susceptibility to heterolytic fission that is a first step of SPLET. The compounds with the lowest PA OH values are the most active. The results of PA OH calculations are shown in Table 5 and Supplementary Table 2. As expected, the PAs are much higher in vacuum than in polar media, as much as ten times (Supplementary Table 2) because gas-phase is less favorable for proton donation. Similarly to BDE OH calculations, phenol was used as a reference system to determine the reliability of our calculations (Supplementary Table 2). The values obtained are in good agreement with the experiment [69,70]. As can be seen from Table 5, compounds have tisomers more active than c-isomers according to SPLET mechanism. In gas-phase PA OH values of molecules with the reactive OH are decreased in relation to phenol in the following order: 4`-OH-KA > m-CA > FA ~ SA > o-CA > p-CA >3`-OH-KA. In polar media, major changes in the order can be noticed: 4`-OH-KA > m-CA > p-CA ~ o-CA > FA > SA > 3`-OH-KA. Based on the results obtained, the m-hydroxylation of cinnamic acid is not preferred according to SPLET mechanism in vacuum, as well as in polar media. In gas-phase, the highest antioxidant potential have compounds only with OH groups in ortho and para position, with 3`-OH-KA as the most active one. The presence of methoxyl group (in FA and SA) decreases antioxidant potential in vacuum. In polar media, 3`-OH-KA is still the most active compound according to SPLET mechanism. It can be explained by the presence of another OH group in 4` position in the aromatic ring (ortho position in relation to 3`-OH), which acts as a strong electron-donating group and contributes to the stabilization of 3`-O-KA monoanion in polar media, also by formation of intramolecular hydrogen bond. The presence of additional functional group in ortho position is advantageous for antioxidant properties and decreases PA OH in vacuum, as well as in polar media.
The o-hydroxylation (in KA) is more favorable than o-methoxylation (in FA and SA). In polar media, FA and SA are more active than o-CA and p-CA, opposed to vacuum. Methoxyl group in ortho position to 3`-OH acts as moderately activating substituent with electron donating abilities. However, its contribution to the stabilization of 3`-O monoanion is weaker than for KA. SA has two OCH 3 groups, which causes the synergistic effect and suggests why lower PA OH value in polar media than FA with only one methoxyl group.

Adiabatic Ionization Potential (AIP)
Adiabatic ionization potential (AIP) is a numerical descriptor of the first reaction of two-step single electron transfer -proton transfer mechanism (SET-PT). It describes the ability of compounds to electron donation. The ease of an unpaired electron abstraction is very important in reactions with free radicals. Molecules with lower AIP are more active according to SET-PT mechanism. This parameter has been calculated for CiA and its derivatives in comparison to phenol. Table 6 shows ionization potential values related to phenol. The absolute AIP are presented in Supplementary Table 3. Our results showed that in vacuum AIP values are significantly higher than in polar media, which suggests that single electron abstraction is preferred in the presence of solvent. Moreover, these values are higher than BDE OH (in vacuum and polar media) and PA OH (in polar media), but lower than gas-phase PA OH . Thus, SET is less preferred in gas-phase than HAT and SET-PT is less preferred in polar media than SPLET. The isomers of CiA and hydroxycinnamic acids have comparable ionization potentials in polar media, but, in vacuum, cis are more active than trans. The planar and bent cis conformation of phenylopropanoic carbon chain in isolated structures promotes close distance between benzene ring and carboxylic group. As an effect, we observe pseudo-diaromatic structure that might additionally stabilize cation radical in vacuum.
Except from CiA (in vacuum and polar media), o-CA and m-CA (in polar media), other compounds have lower AIP than phenol. Thus they are more active according to SET-PT. p-Hydroxylation increases the electron donating capacity, as follows: p-CA < KA < FA < SA in vacuum and p-CA < KA < FA ~ SA in polar media. The lowest AIP values of compounds with p-OH group are connected to the stabilization of radical cation by resonance distributed on whole phenylopropanoic structure. m-OHcation radicals are less stable because of the minor stabilization by resonance. Additional o-hydroxylation or o-methoxylation has inductive effect, which also contributes in the stabilization of cation radical (for KA, FA and SA). The o-hydroxylation (in KA) is more favorable than o-methoxylation (in FA and SA) and it decreases AIP in vacuum, as well as in polar media.

Deprotonation energies
The formation of free radicals is also connected to the Fenton reaction, catalyzed by transition metals in they reduced state e.g. Cu + , Fe 2+ [18]. The presence of O-atoms in phenolic compounds offers chelating sites, like hydroxyl, carbonyl and methoxyl groups. The ability of compound to form metal-ligand complexes is described by transition metals chelation mechanism (TMC). The most preferred chelation sites are connected to deprotonated hydroxyls. Metal chelating properties can be measured quantitatively based on the acidity values (∆H acidity in vacuum and ∆G acidity in polar media). The results obtained are shown in Table 7 (in relation to phenol) and in Supplementary Table 4 (the absolute values).
Acidity of hydroxycinnamic acids is lower than for phenol. The results for c-isomers are higher than for t-in vacuum and comparable in polar media. Acidities of mono-hydroxylated o-CA and p-CA are comparable with the values for FA and SA. At the same time, they are higher than acidity of KA deprotonated in 3`-O position, which suggests that only KA can act as effective chelator. Moreover, deprotonating in meta position for m-CA and 4'-OH-KA is connected with the highest ∆H acidity and ∆G acidity. Thus the meta position may be less active in TMC mechanism.
Our results correspond to the experimental and theoretical studies [70] describing Al(III) complexes with caffeic acid and caffeate. The preferred chelating site for Al(III) was cateholate moiety (o-dihydroxy substitution of aromatic ring), with deprotonated O-atoms contributing in a complexation mechanism.

Stabilization energies (ΔE ISO )
Determination of the stabilization energy is a simple method, which allows one to predict the scavenging ability of phenolic antioxidants [68]. It is connected to trapping free radicals via HAT mechanism. ΔE ISO has been calculated according to the reactions with DPPH and phenol radicals (Reactions 3.1 and 3.2).
The results are shown in Table 9. ΔE ISO values of the compounds studied are compared with phenol. Based on the data obtained, it is possible to estimate the stability for the involved groups in pharmacophore prediction.
The ΔE ISO results for CA suggest that the position of OH group strongly influences the antioxidant`s stability. Hydroxylation in ortho or para position gives better activity than in meta. The presence of 3`-OH group and o-hydroxylation or o-methoxylation decreases ΔE ISO value and enhances the compound`s stability -3`-OH-FA > 3`-OH-SA > 3`-OH-KA (in the decreasing ΔE ISO order). In the case of 3`-O-radicals, an unpaired electron is delocalized on the aromatic ring, vinyl bond and carboxyl group. The π-electron distribution and contribution from o-Oatom (especially strong in KA) stabilizes the free radical of the antioxidant, in vacuum as well as in polar media. The relatively lowest stability is connected to the hydroxyl homolytic fission in meta position for m-CA and KA. It can be explained by the fact that the unpaired electron of m-Oatom radical can only be delocalized in benzene ring. phenol -6.4 -0.6 5.8 -6.5 -0.6 5.9 -6.5 -0.6 5.9 -6.5 -0.6 5.9 c-CiA -6.8 -2.3 4.5 -6.9 -2.3 4.6 -6.9 -2.3 4.6 -6.9 -2.3 4.6 t-CiA -6.

HOMO and LUMO -distribution and energy
Highest occupied molecular orbital (HOMO) and lowest unoccupied molecular orbital (LUMO) are important parameters of molecular electron structure. HOMO represents the ability to donate electrons, whereas LUMO acts as an electron acceptor. Higher frontier orbital energy -E HOMO corresponds to good electron donating ability. The molecule which has lower E HOMO is considered to be a weak electron donor. Moreover, distribution of the frontier molecular orbital is very important in the analysis of antioxidant potential and can indicate the active site prone to free radical attack. The HOMO and LUMO distributions in vacuum for cinnamic acid and its derivatives are shown in Fig. 3 (results in polar media are included in Supplementary Fig. 1). Additionally to the graphical representation, the energy gaps between HOMO and LUMO have been calculated (Table 8).
For all the compounds studied (both c-and tisomers), the electron density of HOMO is concentrated on the benzene ring as well as vinyl bond (C 7 -C 8 ). Mono-hydroxylation in ortho, meta or para position in the aromatic ring, as well as additional hydroxy-/ methoxylation contributes to the HOMO plot (Fig. 3). O-atoms of hydroxyl and methoxyl groups are strongly involved in electron distribution. The electron density of LUMO is mostly located on the C 6 -C 7 and C 8 -C 9 bonds and carboxylic moiety of the phenylopropanoic structure. Mono-hydroxylation in ortho and para positions in the aromatic ring also contributes to the LUMO plot (Fig. 3). The effect from hydroxylation in meta position (for m-CA and KA), as well as methoxylation (for FA and SA) is not observed, opposite to the HOMO plot. Thus, in the transition from ground state to the excited state the whole molecule is involved. Delocalization of the electron density stabilizes the molecular structure. Moreover, the transition event suggests that O-atom in para position is the main active center for hydroxycinnamic acids.
∆E(E LUMO -E HOMO ) informs how easy the transition from ground state to excited state may occur. It is an additional parameter informing about the compound`s activity. The results are shown in Table 8. According to the ∆E values the differences between c-and t-isomers are very small in vacuum and in polar media. The following activity order might be established: phenol < CiA < o-CA ~ m-CA ~ p-CA < KA ~ FA < SA. It can be proposed that among cinnamic acid derivatives KA, SA and FA have very good electron donating abilities with p-OH group prone to the attack of radicals and other electrophilic agents.

Spin density (SD) and SOMO distribution of radicals
Hydroxyl radicals of the compounds studied are the final products of reactions with free radicals, according to the following antioxidant mechanisms -HAT, SPLET and SET-PT [48]. Distribution of canonical singly occupied molecular orbitals (SOMO) represents the resonance structures` formation of the free radicals studied, obtained by hydrogen abstraction of the OH groups ( . For these two compounds the SD values from C 7 -C 8 double bond and carboxylic moiety are 10 to 100-times smaller, which is related to the results of SOMO distribution.

Structure -antioxidant properties of cinnamic acid derivatives
Hydroxycinnamic acids form a group of antioxidants that can act according to the common free radical scavenging mechanisms, namely HAT, SPLET, SET-PT and TMC. The structure-activity relationship analysis indicates the most important features of potent antioxidants: OH groups attached to benzene ring; planar phenylopropanoic structure stabilized by resonance (mainly in radical and radical cation forms) and conjugation effects; additional electron donating functional groups (hydroxyl and methoxyl).

Identification of antioxidant pharmacophore for small cinnamic acids
The SAR analysis of the hydroxycinnamic acids points the t-KA as the most reactive compound with the biggest antioxidant potential. Correlating results of BDE, PA, AIP, ΔE ISO , SD and distribution of different molecular orbitals indicate that 3`-O • -t-KA is the most stable free radical, with an unpaired electron distributed on p-O-atom, benzene ring, vinyl bond and carboxylic moiety. As a consequence, π-type electron delocalization leads to six resonance structures -the key to enhanced free radical stabilization. Based on these results, trans-caffeic acid can be used as pharmacophore antioxidant in search of more effective derivatives. We suggest that using antioxidants as leadcompounds in drug development may be an excellent link between drug development and simultaneous protection against dangerous free radicals.