A comparative study of O2, CO and CN binding to heme-IX protein models

Parametrization of a molecular-mechanics program to include terms specific for five- and six-coordinate transition metal complexes results in computer-simulated structures of heme complexes. The principal new feature peculiar to five and six coordination is a term that measures the effect of electron-pair repulsion modified by the ligand electronegativity and takes into account the different structural possibilities. The model system takes into account the structural differences of the fixing centre in the haemoglobin subunits. The customary proximal histidine is added. The prosthetic group heme IX is wholly considered in our model. The calculations show clearly that certain conformations are much more favourable that others for fixing O2. From the O2 binding in haemoglobin, myoglobin and simple Fe porphyrin models it is concluded that the bent O2 ligand is best viewed as bound superoxide O2-. Axial ligands are practically free-rotating. A small modification of the model in both crystal and protein matrix affects the orientation of the ligands in experimental systems.


Introduction
The heme prosthetic group is in the active centre of a number of relevant proteins as the oxygen (O 2 ) transport proteins haemoglobin (Hb) and myoglobin (Mb) [1,2], as well as enzymes involved in catabolism as peroxidases [3], catalases, oxidases [4] and cytochromes [5].The replacement of Fe by Mg in heme produces chlorophyll [6], and the replacement of Fe by other transition metals coupled with modifications in the aromatic ring produces species as vitamin B 12 [7] and cofactor F-430 [8,9].
Rohmer's group characterized the electronic state of Fe(P) ( P = porphyrin) complexes [10][11][12], and predicted an electronic structure, which was later proven by experiment [13].Other theoretical studies are the real position of the CO group in Fe(P)(Im)(CO) (Im = imidazole) [14,15], the position of the CN group in Fe(mdi) 2 (py)(CN) (mdi = malondialdiminate) complexes [16], the role of distal and proximal histidines (His) on the binding of O 2 in Hb [17][18][19][20], and structural aspects of the binding of O 2 and other ligands to heme [21][22][23][24][25][26][27].The amount of information obtained from the calculations is seriously limited by the size of the heme group, which has only recently allowed the appearance of theoretical studies on reactivity [28][29][30][31][32][33].The coordination of O 2 to Fe(P)(Im) 5-coordinate species leads to 6-coordinate species with octahedral geometry: the biomimetic forms of Mb-O 2 and Hb-O 2 .X-ray data have been reported only for two complexes: Fe(T piv PP) Im](O 2 ) [34] and Fe(T piv PP)[2-(Me)Im](O 2 ) [35].Both are quite similar, sharing the same porphyrin T piv PP = meso-tetrakis(α,α,α,α-o-pivalamidophenyl)porphyrin.Maseras' group optimized the geometry of Fe(T piv PP) Im](O 2 ) with hybrid quantum mechanics/molecular mechanics (QM/MM) IMOMM(DFT-B3LYP:MM3) [36] and pure QM DFT-B3LYP [37].Ghosh and Bocian optimized the geometry of Fe(P)(Im)(CO) with density functional theory (DFT) [19].Salzmann et al. optimized under constraint the geometry of a Fe(T piv PP) Im](CO) model with DFT-B3LYP [38].Han et al. calculated a heme model-CO system employing the ab initio pseudopotential method with local density approximation (LDA) exchange correlation [39].The reversible binding of O 2 and carbon monoxide (CO) plays a central role in studies of hemeprotein structure and function.The affinity of Hb for CO is 200 times that for O 2 , which accounts for the risk for CO poisoning.Numerous encumbered Fe II porphyrin models were synthesized in an effort to elucidate the structural details of small ligand binding.The steric bulk of certain axial ligands bonded to synthetic Fe II porphyrins provided model compounds of reduced O 2 and CO affinity, and models of the tense (T) state of hemoproteins.There is only one single-crystal X-ray structure determination on such a complex [40].There was much discussion on the mechanistic basis of the variation of affinity values in heme proteins and model compounds.This focused on the nature of the axial ligand, distal steric effects, distal polar effects, and enforced doming and ruffling of the porphyrin skeleton.Johansson et al. showed by QM calculations on a heme a model that, upon reduction, the spin pairing at Fe is accompanied by effective delocalization of electrons from the Fe towards the periphery of the porphyrin ring, including its substituents [41,42].In previous works, both non-interacting (NID) and interacting (ID) induced-dipole polarization models were implemented in the program molecular mechanics (MM2) [43].The polarizing force field (MMID2) [44] was applied to For-Gly-NH 2 [45].In the present report, terms specific for five-and six-coordination have been included in MMID2.The new program is called MMIDX.The next section presents the computational method.Following that the results are discussed.Finally, the conclusions are summarized.

Computational method
The molecular polarizability, α ab mol , is defined as the linear response to an external electric field, where µ a ind is the induced molecular dipole moment.Considering a set of N interacting atomic polarizabilities, the atomic µ ind has a contribution also from the other atoms, where T (2)  pq,bc is the interaction tensor as modified by Thole [46] T pq ,ab where v pq = r pq /s pq if r pq < s pq ; otherwise v pq = 1, with s = 1.662(α p α q ) 1/6 ; δ represents the Kronecker delta function: δ ab = 1 if a = b, and δ ab = 0 if a ≠ b.Molecular polarizability can be written as where B is the relay matrix defined as (in a supermatrix notation) 1,3-Interactions between atoms bonded to a common atom are not specifically included in Equation ( 6) in the MM2+polarization approach, because they are effectively already included in the bondlength and bond-angle strainless parameters.
For 5-coordinate structures, there is a need to consider the effect of 1,3-interactions because an energy term is needed for ligand-ligand repulsion to account for (1) the stability difference of various possible geometries and (2) structural effects due to the variation of ligand electronegativity.The geometries of the methylfluorophosphoranes (CH 3 ) n PF 5-n (n = 0→3) were qualitatively correlated with Gillespie's VSEPR theory, which was adopted as a model for the present approach.Bond electron-pair repulsion (EPR) terms were introduced via the non-bonded term E nb of Equation (6) modified to express EPR for atoms bonded to 5-or 6-coordinated atoms.The unmodified term is where P = [r VDW(A) +r VDW(B) ]/r AB , r VDW is the van der Waals radius of the specific atom and r AB is the non-bonded distance between A and B; ε = (ε A ε B ) 1/2 , where ε A and ε B are parameters specific to atoms A and B, and are related to the hardness of the atoms.Hill evaluated the constants in Equation (7), which give the energy E nb(AB) in units of kilocalories per mole [47].The modification of Equation (7) to express EPR terms is [48] E 1, 3 The addition of a scaling factor D, to obtain a suitable balance between this term and the other terms in Equation ( 6), and the replacement of P with P * (where P * = [r VDW(A) +r VDW(B) ]/r * AB and r * AB is the distance between atoms A and B calculated from modified bond lengths, d * CA and d * CB , between the central atom C and either atom A or B) provide the necessary adjustments to quantitatively reproduce the (CH 3 ) n PF 5-n structures [49].The variation in ligand electronegativity is introduced by a distance factor R A in the relation The magnitude of R is inversely related to the electronegativity difference between atoms C and A. The R factors are the means of including EPR between atoms A and B. If the electronegativity difference ∆X CA is large, the bonding electron pair can be considered to move away from atom C, thus decreasing EPR between C-A and C-B bonds.When ∆X CA > ∆X CA' , EPR term E (1,3)AB < E (1,3)A'B even when the actual bond lengths are equal.A set of distance factors R may be obtained from the bond ionic character I, and using the relation ) where r A and r C are covalent radii of atoms A and C.

Calculation Results and Discussion
The structure of the heme(-His)-O 2 model is shown in Figure 1.Heme is the prosthetic group of oxy-Hb.The molecular mechanics calculations use the MM2/MMX+polarization force fields.Van der Waals parameters for the Fe atom have been taken from the force field UFF [50].Torsional contributions involving dihedral angles with the metal atom and the bending terms involving Fe in central position have been set to zero.Molecular mechanics dipole moments µ are given in Table 1.In general, µ increases with the oxidation state of Fe.In particular, µ heme(FeIII)(-His)-CN results the greatest due to the highly polar Fe δ+ -C-N δcomplex.The µ heme(FeIII)(-His)-O2 is relatively large due to the extremely polar Fe δ+ -O-O δ-.The inclusion of polarization in MM2 corrects, in general, µ in the correct direction when compared with MMX+ID.In particular, µ heme(-His) remains almost constant (MM2+polarization).The binding of His in heme increases µ by a factor of 5.Moreover, the subsequent binding of CN doubles µ and the binding of O 2 or CO increases µ 7-16%.

The 4-coordinate Fe(P) system
The heme group is a natural starting point for both experimental and theoretical studies.Crystal structures were reported for a number of heme derivatives with different substituents in the ring.In particular, Fe(TPP) (TPP = meso-tetraphenylporphyrin) has been chosen.Its electronic state is well known experimentally to correspond to a low spin triplet (S = 1).Selected structural parameters are listed in Table 2.The agreement in bond angles between both MM2/MMX+polarization geometries and the X-ray structure is good, with discrepancies smaller than 5°.Differences in bond distances are larger in a number of cases as for the C-C bridge , C-C' and C'-C'' distances.These have values of 1.395Å, 1.439Å and 1.365Å, respectively, in X-ray, and ≈1.336Å, 1.340Å and 1.332Å, correspondingly, in MM2/MMX+polarization.  Agreement between experiment and MM2/MMX+polarization is good (≈0.06Å).All these atoms are, in part, purely described with MM2.The optimized MM2/MMX+polarization values are close to the optimal bond distance for these types of atoms in the applied force field, which is 1.337Å.Another discrepancy in the geometries appears in the Fe-N distance.This is more puzzling, because the calculated distance 1.877-1.881Å is smaller than the experimental value of 1.966Å.The calculated displacement of the Fe atom to the mean plane of the porphyrin N-atoms ring (0.1Å) is not observed experimentally.
The 5-coordinate Fe(P)(Im) system  Coordination of an Im ligand to the heme group leads to a 5-coordinate species with a square pyramidal geometry.These compounds are good biomimetic models of Mb and Hb, with Im replacing the proximal His of the biological systems.The need to avoid dimerization and formation of 6-coordinate species with two axial ligands poses serious restrictions on the nature of the porphyrins able to give this kind of complexes.For this study, the species Fe(Piv 2 C 8 )[1-(Me)Im] {Piv 2 C 8 = α,α,5,15-[2,2'-(octanediamido)diphenyl]-α,α,10,20-bis(o-pivalamidophenyl)porphyrin} has been chosen.This species has the advantage of having 1-methylimidazole as axial ligand, in constrast with the more common 2-methylimidazole, which is more sterically demanding.Neither for Fe(Piv 2 C 8 )[1-(Me)Im] nor for other 5-coordinate derivatives of heme is the electronic state experimentally known.Electronic spectroscopy, magnetic susceptibility and Mössbauer measurements identify it as high spin (S = 2).Selected parameters are collected in Table 3.The Fe-N porph distances are longer than those in the 4-coordinate system by approximately 0.02Å (MMX+ID).This trend is in agreement with the reference values.The result is fully consistent with the shift from low spin to high spin in the metal.Overall agreement in the geometric parameters is correct.Moreover, one has to view with suspicion the X-ray parameters of Im, which would make the N=C double bond N Im -C ε in Im longer than the N-C single bond N Im -C δ .Furthermore, the MM calculations are, in general, in agreement with the QM/MM reference, which provides the expected result.
The sharper discrepancy concerns the N porph -Fe-N Im -C ε dihedral angle.This angle measures the rotation around the Fe-N Im single bond, and rules the placement of the Im plane with respect to the porphyrin ring.Its sign is arbitrary, because the x and y directions are equivalent in absence of an axial ligand.In this work, a positive sign has been chosen for consistency with data on the 6-coordinate complexes presented below.An angle of 90° (as in the pure QM reference) means that the Im plane is eclipsing one of the Fe-N porph bonds, while an angle of 135° (≈135.2° in MMX) indicates a staggered orientation of Im with respect to the Fe-N porph bonds.Therefore both pure QM and MMX values are just opposite with the experimental value, 126.0°, lying in between, although closer to MMX results.The structural minima lead to structures where the Im ligand is located about the bisector of an angle N porph -Fe-N porph .The MMX+polarization results lie, in general, in the range 90-133° of the references.The importance of the large discrepancy between different values is, however, arguable, because there is also a large dispersion in different experimental 5-coordinate derivatives of heme, as well as in experimental reports of Mb and Hb.The corresponding interpretation is that the rotation around this single bond has a rather low barrier.
A final number to mention is the displacement of the Fe atom out of the porphyrin N-atoms plane.This is one of the most important parameters, because it seems to play a critical role in the cooperativity of Hb involved in the allosteric transition.In this case, MMX+NID calculation of 0.274Å lies near the experimental and reference values ca.0.3Å.The positive sign corresponds to the fact that the displacement is towards Im, and away from the empty site.
The geometry agrees with Pauling's prediction of a bent FeO 2 bond, and the O-O distance is close to that of 1.27Å, predicted by him [54].It is slightly shorter than that of 1.34Å in O 2 -close to the electron spin resonance results, which show that no more than 2/3 of the density of one electron is transferred from the metal to the antibonding π * orbitals of O 2 .
The greatest discrepancy concerns the N porph -Fe-O-O dihedral angle.This angle measures the rotation around the Fe-O single bond, and rules the placement of the Fe-O-O plane with respect to the porphyrin ring.An angle of 0° (approximately the -3.3° in MM2) means that the Fe-O-O plane is eclipsing one of the Fe-N porph bonds, while an angle of -45° (approximately the -44.6° in pure QM) indicates a staggered orientation of the O 2 with respect to the Fe-N porph bonds.The structural minima lead to structures where the O 2 ligand is located about the bisector of an angle N porph -Fe-N porph .The MMX/MMX+NID results lie near the reference results.The structural minimum correspond to a trans isomer, where the ending O is placed above the opposite quadrant where the N δ is located.Finally, the Fe atom is computed by MMX+ID to lie in the energy minimum rather close to the porphyrin N-atoms plane, with a deviation of 0.044Å, always towards the O atom.This is again in agreement with the experiment and reference calculations, and rather different from the behaviour of the 5-coordinate system.
Structural parameters of Fe(P)(Im)(CO) are summarized in Table 5.The geometric parameters concerning the coordination of CO, which are probably the most critical for the biochemical activity of Hb, are well reproduced.The computed values for the Fe-C distance, 1.97-2.01Å,are close to the experimental value of 1.770Å.The calculated values for the C-O distance, 1.11-1.14Å,are also close to the experimental report of 1.120Å.The C-O distance is similar in the free molecule, 1.171Å, calculated with AM1 [55] and in Fe(P)(Im)(CO), 1.143Å (MMX).The interpretation is that electronic charge is not transferred from FeP to CO.This is in agreement with the experimental observation that the Fe-CO bond can be formally described as Fe II -CO [9].The most significant feature of the present structure is the linear Fe-C-O bond.While such a linear bond is to be expected based on the extensive literature on transition metal CO complexes, the result is nonetheless highly significant since bent Fe-C-O bonds with bond angles of 135-145º appear to be the rule in various CO complexes of hemoproteins [56][57][58].The global energy minimum is linear for Fe-C-O (MM2/MMX+polarization).These calculations are in agreement with the experiment and calculation references.The Im ring is rotated so that the N atom is directed toward the Fe atom, and the rotation angle N porph -Fe-N Im -C ε reaches approximately 176° (MM2+polarization) in agreement with the LDA reference (174.2°).The Fe atom is computed by MMX to lie in the energy minimum rather close to the porphyrin N-atoms plane, with a deviation of 0.037Å, always towards the C atom.This is again in agreement with the experiment and reference calculations.Selected bond lengths and angles of Fe(P)(Im)(CN) are indicated in Table 6.The geometric parameters concerning the coordination of CN are well reproduced.The computed values for the Fe-C distance, 1.9-2.0Å,are in agreement with the experimental value of 1.934Å.The calculated values for the C-N distance, 1.16-1.20Å,are close to the experimental report of 1.150Å.The C-N distance increases from 1.147Å in the free molecule (AM1) to 1.202Å in Fe(P)(Im)(CN) (MMX+NID).The interpretation is that electronic charge is transferred from FeP to CN.This is in agreement with the experimental observation that the Fe-CN bond can be formally described as Fe III -(CN) - [9].The most significant feature of the present structure is the linear Fe-C-N bond.The global energy minimum is linear for Fe-C-N (MM2/MMX+polarization).These calculations are in agreement with the experiment and calculation references.The Fe atom is computed by MMX to lie in the energy minimum rather close to the porphyrin N-atoms plane, with a deviation of 0.188Å, always towards the C atom.This distance is hardly observed in the experiment.The MMX+ID displacement of the Fe atom to the porphyrin N-atoms plane is shown in Table 7 for the heme-IX adducts as a function of different scaling factors D. The average unsigned error for D = 0.05 (0.094Å) is smaller than for D = 0.1 (0.260Å) and D = 0.2 (0.113Å).The results for D = 0.05 lie in or are closer to the range of the three references.

Conclusions
From the preceding results the following conclusions can be drawn: 1.For the heme-IX adducts, NID or ID polarization energy represents 73% of the total energy MM2+polarization.EPR corresponds to 68% of the total energy MMX+polarization.
2. The optimized structure of the 5-coordinated complexes is a square-base pyramid with the Fe atom 0.3-0.4Åout of the porphyrin plane.This result is in good agreement with the structures of the biomimetic systems that have a sterically unvoluminous axial ligand, but differs from the position of Fe in the α-subunit of deoxy-Hb (0.6Å).
3. Three different Fe-binding models are proposed for O 2 , CO and CN: bent superoxide Fe III -O 2 -, linear Fe II -CO and linear Fe III -CN -.The nature of O 2 binding in Hb, Mb and simple Feporphyrin models is becoming clear.When O 2 is bound as a bent, rather than a triangular, ligand, it is best described as bound superoxide.This bent geometry may be critical to biological functioning because it permits discrimination between O 2 and CO.
4. The distal His protects the heme Fe II acting as a proton trap.The distal His has a pK a of ca.5.5; at neutral pH it is protonated only at N δ , which faces the solvent.Any proton entering the heme pocket of Hb would be bound by N ε , and simultaneously N δ would release its proton to the solvent.When the His side chain swings out of the hemo pocket, the protons would interchange, restoring the previous state.No other amino acid side chain could function in this way.
5. The fact that in a close analogue of the CO complexes of the hemoproteins the Fe-C-O linkage is linear strongly suggests that another interpretation of the results of the protein studies is in order.The allegedly bent Fe-C-O linkage in these proteins was derived from Fourier maps on which the C and O atoms remained unresolved.These maps were interpreted on the assumption that the Fe-C vector was perpendicular to the porphyrin plane.It is much more reasonable to expect that owing to the fixed nature of the globin pocket, bending will occur at the Fe atom, leading to a linear Fe-C-O bond, which is not perpendicular to the porphyrin plane.
6. Unlike Fe-CO, the Fe-O 2 complex is highly polar and the bound O 2 is selectively stabilized by H-bonding to distal His in Mb and Hb.The geometry optimizations performed on the experimental structures are in good agreement with X-ray results, in spite of the fact that only certain residues of the fixing centre have been taken into account.In Hb-O 2 , the distal His H-bonds with the O directly linked with Fe.The H-bond in the fixing centres of Hb differs from that customarily observed in the biomimetic systems and Mb, in which the distal His H-bonds the ending O. The electronic density of the binding O is lower than that of the ending O, suggesting that the H-bond between O and the distal His is weaker than in Mb.

Figure 1 .
Figure 1.Stereo view of the heme(-His)-O 2 model.Heme is the prosthetic group of oxyhaemoglobin.
c ID: polarization by interacting induced dipoles.d Scaling factor D = 0.05.e Average values.f Distance of the Fe atom to the mean plane of the porphyrin N atoms.
b NID: polarization by non-interacting induced dipoles.c ID: polarization by interacting induced dipoles.d Scaling factor D = 0.05.e Average values.f Distance of the Fe atom to the mean plane of the porphyrin N atoms.g Frozen in calculation.h Corresponds to N-H in this calculation.
b NID: polarization by non-interacting induced dipoles.c ID: polarization by interacting induced dipoles.d Scaling factor D = 0.05.e Average values.f Distance of the Fe atom to the mean plane of the porphyrin N atoms.g Frozen in calculation.The state of this system is a low spin open-shell singlet (S = 1) resulting in a Fe III--O 2 -charge distribution.Selected parameters are reported in on Fe(T piv PP)Im](CO) system are provided for comparison[38].b NID: polarization by non-interacting induced dipoles.c ID: polarization by interacting induced dipoles.d Scaling factor D = 0.05.e Average values.f Distance of the Fe atom to the mean plane of the porphyrin N atoms.g Frozen in calculation.
b NID: polarization by non-interacting induced dipoles.c ID: polarization by interacting induced dipoles.d Scaling factor D = 0.05.e Average values.f Distance of the Fe atom to the mean plane of the porphyrin N atoms.
a NID: polarization by non-interacting induced dipoles.b ID: polarization by interacting induced dipoles.

Table 4
[34], but the experimental value, even shorter than the 1.21Å for free O 2 , is suspect, because of the disorder on the placement of the ending O atom within the crystal, as admitted by the authors of the X-ray experiment[34].The O-O interatomic distance increases from 1.21Å in free O 2 to 1.261Å (MMX+NID), suggesting that electronic charge is transferred from FeP to O 2 , in agreement with the experimental observation that the Fe-O 2 bond can be formally described as FeIII-O 2 -[9].The most significant feature of the present structure is the bent Fe-O-O bond.The Fe-O-O bond angles are indicative of a bent η 1 coordination mode, where only one O atom is directly attached to the metal.These calculations are in agreement with the experiment and calculation references.Atomic net charges have been calculated with our program POLAR . The parameters concerning the coordination of O 2 , which are probably the most critical for the biochemical activity of Hb, are well reproduced.The computed values for the Fe-O distance, 1.8-1.9Å,are close to the experimental value of 1.746Å.The calculated values for the O-O distance, 1.21-1.26Å,on the other hand, are far from the experimental value of 1.

Table 5 .
Selected geometric parameters (Å and degrees) from the complete geometry

Table 7 .
Molecular mechanics (MMX+ID) results for heme-IX adducts: distance of the Fe atom to the porphyrin N-atoms plane.