The SPASIBA Force Field for Studying Iron-Tannins Interactions: Application to Fe3+/Fe2+ Catechol Complexes

The SPASIBA force field parameters have been obtained for Fe3+/Fe2+-Oxygen interactions occuring between non-heminic iron and hydroxyl groups of polyphenols found in tannins. These parameters were derived from normal modes analyses based on quantum chemical calculations using the Density Functional Theory (DFT). Four models involving complexation of iron with water ([Fe(H2O)6]3+, [Fe(H2O)6]2+) and with cathechol molecules ([Fe(cat)2(H2O)2]−1, [Fe(cat)2(H2O)2]−2) were studied using the Density Functional Theory and the B3LYP hybrid functional under high spin states of iron.


Introduction
Iron and Tannins (polyphenolic compounds) present in vegetables are part of very important nutrient components for human organisms because of the participation of iron in a large variety of metabolisms and the benefical effects of tannins in the treatment of cardiovasular diseases [1][2][3][4].
Non heminic iron (Fe 3+ \Fe 2+ ) ions present in food display a weak intestinal absorption, worsened by a strong decrease in its availability originating from the appearance of very stable and irreversible complexes with hydroxyl groups present in tannins (polyphenols).
Ealier conformational analyses and complexation investigations on Tannins were first carried out using X-Rays crystallographic studies , NMR experiments using Molecular Mechanics to extract low energy-minimized structures reproducing coupling constants and NOE effects [5,6]. Molecular Modelling methods such as Docking-Scoring , 3D QSAR or Ligand Base Virtual Screening have shown their utility to approach the conformational structure of Tannins in interactions with specific receptors [7]. These methods lead only to a static point of view describing the ability of polyphenols to interact with a large variety of biomolecules and inorganic ions. Empirical methods and, to a less extent, quantum chemical studies are able to give a dynamical point of view on the way that polyphenols can interact with their ligands , both considering their electronic properties and their intrinsic potential energy properties .
A comprehensive description of the geometrical and electronic properties implied in such complex appears as a first task to gain some knowledge about the mechanisms of complex formation. For this purpose, the Density Functional Theory (DFT) was used, in one hand, to identify various possible types of complexations and, in a second hand, to obtain the empirical SPASIBA force field parameters from vibrational spectra in order to perform further molecular dynamical studies.

Density Functional Theory
DFT optimizations of the different complexes were performed using the B3LP hybrid density functional. Iron-water and iron-catechol complexes are well known to appear under a high spin state for iron [8][9] and different basis sets have been shown to be adequate when using hybrid functionals [10].
An extensive review on the structural and spectroscopic (UV absorption, NMR) characteristics of catecholato iron III complexes in catechol dioxygenases observed in various solvents or under tri or tetradentate ligands forms was reported by Yamahara et al. [9]. These last authors reported cis form for catechol rings bound to iron (III) as observed by X-rays studies on protocatechuate 3,4 dioxygenase (3,4-PCD). This form was found to be less stable than the trans form which was chosen in this study.
Iron III adopts a 3d 5 electronic configuration and in the presence of a strong field environment displays a repartition of its d orbitals in two kinds of Molecular Orbitals i.e ; e g (d z2 ,d x2-y2 ) and t 2g (d xy , d xz , d yz ).
For hexaaqua and catechol complexes, a LAVCP** basis set was used for heavy atoms and a 6-31g (tm) ** one for the other atoms.  2 ] -2 . These basis sets exhibit polarization and diffuse functions well adapted to inorganic ions and are implemented in the Jaguar program used in the present work [11].

Empirical force field
The empirical SPASIBA (Spectroscopic Potential Algorithm for Simulating Biomolecular Conformational Adaptability) potential energy function includes classical and Urey-Bradley-Shimanouchi terms adapted to Infrared and Raman spectroscopic studies. It includes redundancies terms originating from tetrahedral constraints and additional trans-gauche terms. The particular form of this potential energy has been largely described elsewhere [12][13][14] and the related parameters (force constants) have been obtained for a large variety of chemical groups and biomolecular families (peptides, lipids, saccharides) from Raman and Infrared vibrational spectroscopies. 6 6 ] 3+ complexes displays an octahedral ferric/ferrous combination with six water molecules adopting a pratically pure T h molecular point group symmetry leading to a spherical symmetry (Figure 1     The final optimized geometry of the Fe(III)-catechol complex (Fe(Cat) 2 (H 2 O) 2 ] -1 displays a C 2h molecular point group symmetry with the Fe atom located at the inversion center. As can be observed in Figure 2, the two ring planes adopt quasi parallel orientations displaying an mean interplanar distance of 0.6 Ǻ and alternative angular values of 84.7° and 95.3° for the two direct O hyd -Fe-O wat angles belonging to a same ring. By examination of Table 2b, one can observe a reasonable agreement between the values of the optimized internal parameters and those extracted from X-Rays studies [15], ESR spectroscopy [16] or other experimental and calculated values reported by various auhors [17][18][19][20] . According to Funabiki et al. [19] a mean value of 2.Å for the Fe-O bond length must correspond to a minimum in the structural energy.  The SPASIBA parameters (Table 3) associated to the two Fe(III) complexes were obtained from normal coordinate analyses using experimental and DFT vibrational frequencies partly obtained from the work of Jarzecki et al. [8].

Fe(III) complexes [Fe(H 2 O)
The present DFT vibrational frequencies and the associated potential energy distribution (P.E.D) among internal coordinates were derived via the Redong program (Allouche and Pourcin [21]) using a general scaling factor of 0.962 applied to each internal force constants to fit the theoretical and the experimental data. Such a factor is generally used to give correct agreements with experiments [22].  (Table 3).
Since no exprimental vibrational data are available, comparizons were directly done with the DFT derived vibrational frequencies in one hand, and, in a second hand, with the characteristic values observed for the various chemical groups [18][19][20]. The empirical values associated with the internal coordinates of the catechol complex after minimization are reported in Table 2b and display a general good agreement with the experimental and predicted DFT values with a very weak distance of 0.13-0.20 Å between the two catechol molecular planes .

DFT calculations
In the [Fe(H 2 O) 6 ] 2+ complex, the Fe 2+ ion displays a d 6 electronic configuration and forms an octahedral complex subject to Jahn-Teller distorsions constraining the molecular point group to adopt a C i symmetry to form a stable conformation [8]. Table 5   The DFT optimized geometry related to the[Fe(Cat) 2 (H 2 O) 2 )] -2 complex displays a similar C2 h molecular point group symmetry with an interplanar distance of about 0.5 angströms between the two catechol molecular planes with O hyd -Fe-O wat valence angles (belonging to a same ring) adopting values of 72.5° and 107.5°.

SPASIBA parameters for the Ferrous complexes
After vibrational normal modes analyses and vibrational frequencies refinements on the ferrous complexes, the final derived empirical Fe-O stretching force constant appears as being the most affected with an associated force constant decreasing from K=117.0 (Fe 3+ ) to 65.0 (Fe 2+ ) kcalmol -1 Å -2 while in plane valence angle bending and torsional parameters remain essentially unchanged (Table 6).   Torsional C-O motions are generally coupled with C-C torsional contributions and have been calculated at 261.7(SPASIBA), 237.6(DFT) cm -1 in the present woork. These values can be correlated to the experimental value (309cm -1 ).
The Other stretching modes such as C-O have vibrational wavenumbers in the 1149-1284 (SPASIBA) and 1151-1309 (DFT) cm -1 ranges. These two regions include the characteristic experimental band assigned to this mode observed at 1275 cm -1 .
High frequency vibrational ranges related to C-H symmetric and antisymmetric stretching motions (3114-3121 /SPASIBA and 3099-3160 cm -1 /DFT) can be fitted with the 3027,3049 and 3070-3075 cm -1 ) observed ones. O-H stretching modes have comparable predicted DFT/SPASIBA wavenumbers in the 3647-3766 cm -1 ) range. Both DFT and the SPASIBA empirical force field under use reproduce satifactorily the mean vibrational features related to the phenolic groups.

Conclusion
A preliminary determination of the empirical SPASIBA force field devoted to aqueous and cathecolate complexes with Fe(III) and Fe(II) was carried out. Geometrical parameters and vibrational frequencies obtained from normal mode analyses on hexaaqua and dicatechol-water (Fe 3+ /Fe 2+ ) complexes have been compared to experimental and theoretical quantum DFT values . The same set of force constants could be used for both FeIII/FeII complexes apart for that related to the Fe-O stretching which has to be chosen according to the electronic state of spin under consideration A generally good agreement between all methods was observed. Further calculations are in progress to apply this force field for Molecular Dynamics which would permit to follow the complex formation and its evolution along physicochemical constraints (pH, temperature, counter-ions or external competitive ligands) keeping in mind the biochemical importance of iron releasing in organisms.