Harmonic and anharmonic dynamics of Fe-CO and Fe-O2 in heme models

ABSTRACT We present density-functional molecular dynamics simulations of FeP(Im)(AB) heme models (AB = CO, O2, Im = imidazole) as a way of sketching the dynamic motion of the axial ligands at room temperature. The FeP(Im)(CO) model is characterized by an essentially upright FeCO unit, undergoing small deviations with respect to its linear equilibrium structure (bending and tilting up to 10 deg and 7 deg, often occur). The motion of the carbon monoxide ligand is found to be quite complex and fast, its projection on the porphyrin plane sampling all the porphyrin quadrants in a short time (~0.5 ps). Simultaneously, the imidazole ligand rotates slowly around the Fe-N^sub E^. bond. In contrast to carbon monoxide, the oxygen ligand in FeP(Im)(O2) prefers a conformation where the projection of the O-O axis on the porphyrin plane bisects one of the porphyrin quadrants. A transition to other quadrants takes place through an O-O/Fe-N^sub P^ overlapping conformation, within 4-6 ps. Further details of these mechanisms and their implications are discussed.

INTRODUCTION

The reaction of carbon monoxide and 02 with the heme iron of myoglobin (Mb) and hemoglobin (Hb) has long been a topic of discussion and debate (Springer et al., 1994; Olson and Phillips, 1996; Sage and Champion, 1996; Ormos et al., 1998). Of special interest are the properties of the Fe-CO and Fe-O2 units, due to their possible relation with the protein function mechanisms. In this respect, the structure and vibrational properties of the Fe-CO and Fe-O2 moieties have often been used as a sensitive probe for heme-ligand binding and for electrostatic interactions in the distal pocket (Hirota et al., 1996; Unno et al., 1998). The structure of the Fe-CO unit, in particular, was thought from an early stage to be relevant for carbon monoxide binding control (Stryer, 1995), although the quite large distortions reported in x-ray studies (Cheng and Schoenborn. 1991; Yang and Phillips, 1996) have been challenged (see, for instance, Ray et al., 1994; Slebodnick and Ibers, 1997). Spectroscopic studies (Lim et al., 1995; Sage and Jee, 1997) give a much smaller distortion (~7 deg deviation from linearity) and theoretical calculations have demonstrated that the energetic cost for such a small deformation is marginal (Vangberg et al., 1997; Rovira et al., 1997; Spiro and Kozlowski, 1998; Havlin et al., 1998). In addition, experiments on isolated heme models have been undertaken with success (Collman, 1997, and references therein). However, little is known on dynamic aspects like the distribution of the carbon monoxide orientation at room temperature, which could help in the modeling of the experimental data.

The Fe-O bond exhibits a bent end-on geometry both in MbO^sub 2^ and HbO^sub 2^ (Shaanan, 1982), as it is also found in isolated hernes (Jameson et al., 1978). Neutron and x-ray diffraction measurements have found that the 02 is hydrogen-bonded with the N^sub epsilon^-H of His-64 in both MbO^sub 2^ and alpha-Hb (Phillips and Schoenborn, 1981; Shaanan 1983). However, there is no evidence of a hydrogen bond in Beta-HbO^sub 2^, where free rotation of the ligand around its equilibrium position is expected (Shaanan, 1993). Recent electron paramagnetic resonance measurements in cobalt-substituted Hb (in both alpha and Beta subunits) have found evidence of O2 rotation (Bowen et al., 1997). Evidence of O2 rotation is also provided by investigations of synthetic models. In particular, the fourfold disorder found in the crystal structure of picket-fence oxygen systems (Jameson et al., 1978, Collman, 1997) has been interpreted as a dynamic 02 Motion by both M6ssbauer and nuclear magnetic resonance experiments (Spartalian et al., 1975; Mispelter et al., 1983; Oldfield et al., 1991).

The thermal motion of imidazole axial ligands has also been investigated by electron paramagnetic and nuclear magnetic resonance spectra in proteins and synthetic models in homogeneous solution (Nakamura et al., 1996; Momot and Walker, 1997). In the heme proteins, the imidazole is constrained by covalent bonding to the protein backbone. However, in synthetic heme models, it appears to rotate freely and adopt a variety of orientations.

To understand the relation of these motions to the properties of the protein, it is necessary to transcend a purely static point of view and fully examine the influence of thermal fluctuations. This is the objective of the present study, in which first-principles molecular dynamics (MD) will be used to extract information on Fe-CO, Fe-O^sub 2^ and Fe-imidazole motion in FeP(Im)(O^sub 2^) and FeP(Im)(CO) models at room temperature. This study is part of our investigation into ligand binding properties of heme models (Rovira et al., 1997; Rovira and Parrinello 1999). To the best of our knowledge, it is the first study in which the heme-ligand dynamics is examined from first principles.

THE COMPUTATIONAL METHOD

A detailed description of the Car-Parrinello method can be found in several publications (Car and Parrinello, 1995-, Galli and Parrinello, 1991). The calculations were performed with periodic boundary conditions, using an orthorhombic supercell of dimensions a = b = 15 Angstrom, c = 12 Angstrom [FeP(lm)(CO)J and a = b = 15 Angstrom, c = 10 Angstrom (FeP(Im)(O^sub 2^)]- We used the generalized gradient-corrected approximation of spin-dependent density functional theory, following the prescription of Becke (1986) and Perdew (1986). The electronic wave functions were expanded in plane waves up to a kinetic energy cutoff of 70 Ry. Martins-Troullier norm-conserving pseudopotentials were used (Troullier and Martins, 1991), supplemented with nonlinear core corrections for the iron atom (Louie et al., 1982). This is the same setup as we used for the study of structural and equilibrium properties of related heme models (Rovira et al., 1997; Rovira and Parrinello, 1999). The deuterium mass for the hydrogen atoms was used in the dynamic simulations, which allows the use of a longer time step (0.121 fs). The fictitious mass of the Car-Parrinello Lagrangian was set to 700 au. The four meso hydrogens of the porphyrin were kept fixed, because almost all experimental heme models contain bulky substituents in this position. The initial configurations were taken as an overlapping O-O/Fe-N^sub P^ conformation in the case of the O2 complex and the minimum energy structure in the case of the carbon monoxide complex. The systems were allowed to evolve during 2 ps to achieve vibrational equilibration and to lose information of the initial configuration. Further methodological details can be found in Rovira et al. (1997).

RESULTS AND DISCUSSION

Static properties

Detailed information on the structure and electronic properties of FeP(Im)(XY) complexes (XY = CO, O^sub 2^) can be found in our previous work (Rovira et al., 1997; Rovira and Parrinello, 1999) and in several other theoretical studies of heme models (for instance, Dedieu et al. 1983; Jewsbury et al., 1994; Nagatsuji et al., 1996; Vangberg et al., 1997; Havlin et al., 1998; Maseras, 1998; Spiro and Kozlowski, 1998; Sigfridson and Ryde, 1999). Here, we just briefly recall the main stereochemical properties of these complexes. Both systems have a quite similar structure, except for the Fe-XY unit. As depicted in Scheme 1, carbon monoxide binds linearly to the iron, whereas O^sub 2^ prefers a bent end-on type of bond (<Fe-O-O = 121), with the O-O projection on the porphyrin plane bisecting one of the porphyrin quadrants (Rovira et al., 1997). The energetic cost for just a small distortion of the Fe-CO bond is very small: a 7 deg deviation of the Fe-CO angle involves less than one kcal/mol; however larger deviations (>= 15 deg) are prohibitive (Vangberg et al., 1997). Small changes in the Fe-O-O angle can also be made at a small energetic cost, although the change to an upright Fe-O-O is costly in terms of energy (15 kcal/mol) and involves a change of spin state (from singlet to triplet) (Rovira and Parrinello, 1999).

Additional calculations were performed to explore in more detail the shape of the potential energy surface of the FeP(Im)AB systems (AB = CO, O^sub 2^) with respect to the axial ligand position. The structure of the FeP is kept fixed in a planar configuration. The orientation of the axial ligands with respect to the porphyrin was kept fixed at the orientations of interest (see Fig. 1), while the remaining degrees of freedom were optimized. The imidazole ligand was oriented either along the bisector of the NP-Fe-N^sub p^ angle (Fig. 1 a) or overlapping a Fe-N^sub p^ bond (Fig. 1 b). Concerning the diatomic ligands, only the orientational configuration corresponding to the potential energy surface minimum was considered for carbon monoxide (Scheme 1). In the case Of O^sub 2^, which at equilibrium is bent (Scheme 1), the orientational configurations of Fig. 1 were calculated.

From these calculations, we can infer that, for both carbon monoxide and O^sub 2^, the energy is insensitive to the imidazole orientation. In the case Of O^sub 2^, there is an energy barrier of 1.3 kcal/mol (independent of the imidazole rotation) when the O2 is oriented along a Fe-N^sub p^ bond (Fig. 1, a1, b1, b2). Additional calculations relaxing all degrees of freedom gave a value of 1.2 kcal/mol for this barrier, which indicates that also porphyrin distortions will not affect the torsional motion of the ligand. This splits the system into four essentially equivalent porphyrin quadrants.

Given the small energy barrier for O2 rotation around Fe-O, and considering that the torsional mode associated with this motion should have a very low frequency, jumping of the O^sub 2^ from one porphyrin quadrant to another at room temperature is expected to occur in a timescale of a few picoseconds. We should finally point out that the above results do not change significantly regardless of whether the system is in an open-shell or closed-shell electronic state (O-O = 1.29 Angstrom, <Fe-O-O = 122 deg, Fe-O = 1.74 Angstrom at the minimum structure). Thus, for computational reasons, we decided to perform the simulations on the closed-shell surface of this system. It is also worth noting that the above structure/energy relations and the main equilibrium properties of both carbon monoxide and O^sub 2^ models are quite similar to those found for picket-fence systems (Rovira and Parrinello, 1999). Thus, we expect that the conclusions drawn here can also be applied quantitatively to cases that are as yet unmanageable with first-principles MD.

Heme-CO dynamics

An MD simulation of FeP(Im)(CO) (Fig. 2) was performed for a total period of 18 ps, with an average temperature of 300 K. The Fe-N^sub epsilon^ was constrained to its equilibrium value (2.08 Angstrom) during an initial time interval of 2 ps after equilibration of the system.

The most salient features that can be observed from the ligand dynamics is the fast complex motion of the carbon monoxide and the much slower motion of the imidazole. As a way to display the motion of the carbon monoxide ligand, we monitored the projection of the carbon and oxygen atoms on the average plane defined by the four porphyrin nitrogens (hereafter referred to as the N^sub p^ plane). Figure 3 shows the trajectory of both carbon and oxygen projections on this plane. The iron atom is located at the center of each plot, with the x and y axes aligned with the Fe-N^sub p^ bonds (see Fig. 2). As shown in Fig. 3 A, the trajectory of the carbon atom appears to be rather complex and it is concentrated around the iron atom. The trajectory of the oxygen atom (Fig. 3 B) has similar features but is more spread (~0.4 Angstrom from the center). Nevertheless, it should be noted that, with respect to the size of the porphyrin (Fe-N^sub p^ = 2.02 Angstrom), the whole spread of values shown in Fig. 3 B corresponds to just a very small area over the iron atom. Further analysis of the time evolution of the carbon monoxide orientation (data not shown here) reveals that the projection of the C-O axis on the porphyrin plane visits all the porphyrin quadrants in a very short time (~0.5 ps). Therefore, the global picture that can be inferred from our simulation is that of an essentially upright Fe-CO unit, with the carbon monoxide ligand undergoing a fast complex motion within a very small region around its equilibrium position.

Deviations of the Fe-CO unit from its linear equilibrium structure are commonly described in terms of the Fe-C-0 tilt (delta^sub c^) and bend (Theta^sub o^,) angles (Collman, 1997; Schlichting et al., 1994; Yang and Phillips, 1996), depicted in Fig. 2, which have been traditionally related to the protein discrimination for carbon monoxide. The frequency distribution of the delta^sub c^ and theta^sub o^, angles obtained from our dynamics is shown in Fig. 4 A. Small fluctuations of these angles during the MD simulation are quite frequent: delta^sub c^ <= 5 deg, theta^sub o^ <= 8 deg, but larger deformations are unlikely to occur. This is consistent with the results of static calculations (Vangberg et al., 1997; Rovira et al., 1997; Havlin et al., 1998), which have predicted that small delta^sub c^ - theta^sub o^ variations raise the energy of the system by less than 1.5 kcal/mol, approximately the thermal energy available to the Fe-CO moiety.

However, we should point out that, given the complex motion of the ligand described above (Fig. 3), the instantaneous carbon monoxide position cannot be easily defined only in terms of the 5c and 0,, angles. To further illustrate this point, we characterized the trajectory of the carbon monoxide using the 8, and 8. angles defined in Fig. 2. Note that the case 80 > bc corresponds to an in-phase combination of the two angles (i.e., outward bending of the carbon monoxide). As shown in Fig. 4 B, these conformations are the most frequently observed. However, those characterized with an inward bending (50 < 5,) also contribute significantly to the carbon monoxide motion. Therefore, care should be taken when labeling the structure of the Fe-CO unit with the conventional in-phase combination of 8, and 00 angles. The problem should be best regarded as that of a highly dynamic Fe-CO moiety, sampling many different conformations with different probability in a short time.

Concerning the movement of the imidazole, it is dominated by the rotation around the Fe-NE bond and is much slower: while Im rotates 1800 from its initial position (which takes place in one ps) the carbon monoxide samples all porphyrin quadrants. No correlation was observed between the motion of the two axial ligands, as to be expected from the static calculations reported in the previous section. Both clockwise and counterclockwise types of rotation are found, with frequent changes from one to the other mode of rotation and without a clear preference for either an overlapping or bisecting conformation with respect to the Fe-NP bonds (Fig. 1, a and b, respectively). Simultaneously, the porphyrin ring undergoes large amplitude out-of-plane displacements (up to 0.9-A displacements of the pyrrolic carbon atoms are observed). Many peculiar modes of vibration are excited, like a boat distortion and ruffling of the porphyrin ring.

Heme-02 dynamics

An MD simulation for the oxygen analogue, FeP(Im)(02), was performed for a total time of 15.5 ps. As for the carbon monoxide complex, our main interest was to characterize the dynamic motion of the 02 and imidazole axial ligands. To define the orientation of the 02 with respect to the porphyrin plane, we have used one NP-Fe-01-02 torsional angle (phi) and the projection of the O^sub 2^ center of mass on the average porphyrin plane.

The evolution of the phi angle during the MD simulation is presented in Fig. 5. During the first period of the simulation, the 0-0 axis projection on the porphyrin plane lies on the first porphyrin quadrant (1) but undergoes large oscillations between the Fe-N, and Fe-N^sub 2^ bonds. After 2.2 ps, the 02 jumps over the Fe-N^sub 2^ bond toward the second quadrant. Apparently, the energy accumulated in the Fe-O^sub 2^ rotational mode is high enough for the ligand to skip the second and third quadrants and end up in the fourth quadrant (I --> IV, counterclockwise). Two more transitions take place at 8 ps (TV - III) and 13.5 ps (111 ---> IV). All transitions take place via rotation of 02 around the Fe-O axis and involve a conformation with a more open Fe-0-0 angle (124 deg) and the Fe-O bond slightly tilted (3-5 deg ) with respect to the z axis. The Fe-02 tilt (delta) and Fe-0-0 angle (theta) show oscillations similar to the ones found for carbon monoxide (delta <= 10 and 0 = 117-130, as shown in Fig. 6). It is at one of these oscillation maxima that the transition takes place. Therefore, our results provide evidence for the Fe-02 dynamic motion proposed to explain the fourfold disorder found in the crystal structure of picket-fence oxygen systems (Jameson et al., 1978, Collman, 1997). They also confirm the fact that the O-O/Fe-N overlapping configuration is the transition state for the dynamic motion of 02 between the porphyrin quadrants (Spartalian et al., 1975; Mispelter et al., 1983; Oldfield et al., 1991).

Concerning the motion of both axial ligands, we do not observe any correlation between the movement of the 02 and the conformation of the imidazole (for instance, complete fixing of the Im position for 0.5 ps during the second picosecond of the simulation did not change the evolution of the 02 motion), as it was found in the case of the carbon monoxide analogue. The conformation of Im at each jump (i.e., when the 0-0 axis overlaps one Fe-N^subp^ bond) is very similar to al, b2, al, and al for the first four jumps, respectively (Fig. 1). The last transition (111 ---> IV) is characterized by an almost eclipsing conformation of the imidazole similar to b L

It is important to note that the limited time sampled in the simulation precludes a rigorous statistical analysis of the rotational motion of O^sub 2^ about the Fe-O bond and the one of Im with respect to Fe-N.. However, it is expected that, for longer times, the axial ligands would sample all porphyrin quadrants with equal probability. Thus, as a way to enhance our signal, we averaged our data over the porphyrin symmetry operations. Figure 7 A shows the probability distribution corresponding to the 02 center-of-mass projection on the porphyrin plane. For the sake of comparison, the same type of distribution for the carbon monoxide ligand is shown (Fig. 7 B). In the case of carbon monoxide, the distribution is characterized by having a probability maximum for a small region around the iron. In contrast, the distribution for the 02 is qualitatively different from that of carbon monoxide. Four probability maxima are found in this case, which are concentrated on regions along the bisector of each NP-Fe-NP quadrant. This conformation corresponds to the equilibrium structure of the Fe-02 bond. The small barrier among the minima (1.2 kcal/mol, as reported in the section, Static properties) is the reason for the lack of density in the direction along the Fe-NP bonds. In contrast, the high energy of a linear Fe-0-0, 15 kcal/mol (Rovira and Parrinello, 1999), precludes the sampling of this configuration at room temperature. The shape of the density pockets of Fig. 7 A reflects the frequent oscillations of the ligand within each quadrant, changing the Np-Fe01-02 torsional angle. This is consistent with the results of static calculations that show that the torsion around the Fe-O bond is a soft degree of freedom. To extract additional information on the 02 libration from the MD trajectory, the temporal Fourier transform of the velocity autocorrelation function for the iron and oxygen atoms was computed. The torsional mode corresponding to the 02 rotation around Fe-O is indeed characterized by a very low frequency (97 cm-'). The Fe-0-0 bending is found at 380 cm-'. To the best of our knowledge, this mode has not yet been assigned for oxymyoglobin, but it has been observed at 273 cm-' in heme models (Nakamoto, 1990). The 02 stretching frequency obtained (1180 cm- 1) is very close to the values reported for Mb02 (1103 cm-'; see Momenteau and Reed, 1994). The closest experimental realization of our model calculation is the picket-fence oxygen synthetic model (Jameson et al., 1978). It is therefore important to compare our results with the reported x-ray structure of the Fe-02 moiety: 0-0 = 1.15, 1.17 A, Fe-O = 1.75 A, <Fe-0-0 = 133, 129. These data are in apparent disagreement with our results: 0-0 = 1.30 A, Fe-O = 1.96 A, <Fe-0-0 = 1230, which are obtained by averaging these properties over our MD run. This could be attributed to the effect of the picket substituents, but our previous study on FeP(Im)02 and oxypicket fence (Rovira and Parrinello, 1999) showed these effects to be rather small, which rules out this possibility. However, we can reconcile theory and experiment if we recall that x-ray experiments do not measure distances and angles directly, but rather the positions of maximum probability, from which the other structural properties are deduced. We can prove this point from our trajectories, and, in fact, if we extract distances and angles from the position of maximum probability, we find 0-0 = 1. 19 A, Fe-O = 1.72 A, <Fe-0-0 = 1390, which is in much better agreement with experiment. Note that the experimentally reported 0-0 distances are very surprising because they are smaller than the reported gas phase value, which, for 02, is 1.23 A. Similar considerations hold for the reported data on the proteins, in which comparison with our results is further complicated by other affects, such as the presence of residues close to the heme pocket. In particular, the Fe-0-0 angles reported by x-ray and neutron structures of oxymyoglobin (1150) (Phillips, 1980; Phillips and Schoenborn, 1981) and oxyhemoglobin (1530 for the a subunit and 1590 for the 0 subunit) (Shaanan, 1983) are quite far from our average value (1230) and are not sampled in the simulation (Fig. 6).

In summary, our simulations of FeP(Im)(AB) (AB = CO, 02) reveal a picture of highly dynamic axial ligands. In the case of carbon monoxide, an essentially harmonic dynamics is found. The carbon monoxide ligand undergoes small, albeit extremely complex, displacements around its equilibrium position and its dynamics is not influenced by the imidazole ligand. The distributions obtained (Figs. 3, 4, and 7) represent the expected carbon monoxide dynamics in myoglobin in the absence of interaction with the polypeptide framework. This can be useful in the interpretation of the x-ray and spectroscopic measurements to determine the average Fe-CO structure, although, as mentioned before, care should be taken when associating a rigid structure with the highly dynamic Fe-CO moiety. In particular, the propensity for the oxygen to bend out from the porphyrin axial direction should be noted. Further interactions within the heme pocket could perturb the distributions reported here. For instance, a weak hydrogen bond (Unno et al., 1998) or an electrostatic interaction with His-64 would probably lead to the nonequivalence of the distribution of Fig. 7 A over the four porphyrin quadrants.

In the case of FeP(Im)(02), our simulation reveals a highly anharmonic dynamics for the 02 ligand, which undergoes large amplitude oscillations within one porphyrin quadrant and jumps from one to the other every 4-6 ps. This is consistent with the highly dynamic nature of 02 bound to heme proposed by several experiments in proteins and synthetic models (Jameson et al., 1978; Bowen et al., 1997; Spartalian et al., 1975; Mispelter et al., 1983; Oldfield et al., 1991) especially those that lack a hydrogen bond at the terminal oxygen. Ligand rotation in these models has been shown to occur by nuclear magnetic resonance experiments (Mispelter et al., 1983; Oldfield et al., 1991), on the basis of the equivalence of the pyrrole proton resonances. Our results suggest that, for a nonhydrogen bonded 02, precise determination of the rate of rotation would require picosecond time resolution.

In contrast, a possible hydrogen bond with His-64, which has been suggested several times as the origin of the discrimination between carbon monoxide and 0, would alter this picture. In this respect, restriction of the librational motion of the bound ligand caused by hydrogen bonding has been discussed in relation with the loss of entropy upon 02 binding to myoglobin (Filiaci and Nienhaus, 1997). Given the low frequency of the torsional mode of 02 rotation (97 cm^sup -1^), a weak hydrogen bond to His-64 (Phillips and Schoenborn, 1981; Shaanan, 1983) is indeed very likely to slow down the rotational motion of the ligand.

CONCLUSIONS

The room temperature motion of ligand binding to heme has been depicted by means of Car-Parrinello MD simulations on FeP(Im)(CO) and FeP(Im)(02) models. Our results illustrate the fluxionality of the Fe-CO and FC-02 bonds in the absence of the protein environment. From a technical point of view, they provide a reference that can be used to interpret future QM/MM calculations aimed at quantifying the role of the protein framework. The simulations reported here reflect a quite rigid Fe-CO bond, which undergoes little distortion around its equilibrium position. Rotation of the imidazole axial ligand around the Fe-N^sub epsilon^ bond is also observed, in agreement with experiments on heme models with imidazole-based axial ligands. Several features are shared with the analogous 02 model, but the bent Fe-02 bond shows a preference for a bisecting conformation, jumping over one Fe-NP bond toward a different porphyrin quadrant within 4-6 ps. Oscillations of the Fe-02 tilt and bend angles are quite large and, because of the anharmonicity of the Fe-02 motion, results obtained by an averaging value analysis might differ from those obtained from the maximum probability position. This reconciles the values reported in x-ray analysis with the average values we obtain in the dynamics. The rotational motion of the ligand is characterized by a frequency mode at 97 cm^sup -1^. Given this low value, the axial rotation of O2 would probably be slowed down by a weak hydrogen bond with His-64. Thus, in the case of Mb02 and a-Hb02 (Phillips and Schoenborn, 1981; Shaanan, 1983), only one of the quadrants would probably be selected.

We thank the Garching Computer Center (Garching, Germany) for computing support. C. Rovira acknowledges the financial support of the Training and Mobility of Researchers programme of the European Union under contract No. ERBFMBICT96-0951. We also thank S. Raugei, G. Lippert, R. Rousseau and E. Canadell for useful discussions.

[Reference]

REFERENCES

[Reference]

Becke, A. D. 1986. Density functional calculations of molecular bond energies. J. Chem. Phys. 84:4524-4529.

Bowen, J. H., N. V. Shokhirev, A. M. Raitsimring, D. H. Buttlaire, and F. A. Walker. 1997. EPR studies of the dynamics of rotation of dioxygen in model cobalt(II) hemes and cobalt-containing hybrid hemoglobins. J. Am. Chem. Soc. 101:8683-8691.

Car, R., and M. Parrinello. 1985. Unified approach for molecular dynamics and density-functional theory. Phys. Rev. Lett. 55:2471-2474.

Cheng, X., and B. P. Schoenborm. 1991. Neutron diffraction study of carbonmonoxymyoglobin. J. MoL BioL 220:381-399.

Collman, I P. 1997. Functional analogs of heme protein active sites. Inorg. Chem. 36:5145-5155.

[Reference]

Dedieu, A., M.-M. Rohmer, and A. Willard. 1983. Structure and properties of a model of deoxyheme. An ab initio SCF calculation. Chem. Phys. 77:449-462.

Filiaci. M., and G. U. Nienhaus. 1997. The role of entropy in the discrimination between CO and 02 in myoglobin. Eur. Biophys. 1. 26:209-214. Galli, G., and M. Parrinello. 1991. Computer Simulation in Materials

Science. V. Pontikis and M. Meyer, editors. Kluwer, Dordrecht, The Netherlands.

Havlin. R. H., N. Godbout, R. Salzmann, M. Wojdelski, W. Arnold. C. E. Schulz, and E. Oldfield. 1998. An experimental and density functional investigation of iron-57 Mossbauer quadrupole splittings in organometallic and heme-model compounds: applications to carbonmonoxy-heme protein structure. J. Am. Chem. Soc. 120:3144-3151.

Hirota, S., Tiansheng, L., Phillips, G. N., Olson, J. S., Mukai, M., Kitagawa, T. 1996. Perturbation of the Fe-02 bond by nearby residues in heme pocket: observation of VF-02 raman bands for oxymyoglobin mutants. J. Am. Chem. Soc. 118:7845-7846.

[Reference]

Jameson, G. B., G. A. Rodley, W. T. Robinson, R. R. Gagne, C. A. Reed, and J. P. Collman. 1978. Structure of a dyoxygen adduct of (1methylimidazole)-meso-tetrakis(a,a,a,a-o-pivalamidophenyl) porphinatoiron(II). An iron dioxygen model for the heme component of oxymyoglobin. Inorg. Chem. 17:850-857.

Jewsbury, P., S. Yamamoto, T. Minato, M. Saito, and T. Kitagawa. 1994. The proximal residue largely determines the CO distortion in carbonmonoxy globin proteins. An ab initio study of a heme prosthetic group. J. Am. Chem. Soc. 116:11586-11587.

Kroon, J., and J. A. Kanters. 1974. Non-linearity of hydrogen bonds in molecular crystals. Nature. 248:667-669.

Lim, M., T. A. Jackson, and P. A. Anfinrud. 1995. Binding of CO to myoglobin from a heme pocket docking site to form nearly linear Fe-C-0. Science. 269:962-966.

Louie, S. G., S. Froyen, and M. L. Cohen. 1982. Nonlinear ionic pseudopotentials in spin-density-functional calculations. Phys. Rev. B. 26: 1738-1742.

Maseras, F. 1998. Binding of dioxygen in a picket-fence porphyrin complex of iron. A theoretical QM/MM study. New. I Chem. 22:327-332. Mispelter, J., M. Momenteau, D. Lavalette, and J.-M. Lhoste. 1983. Fly

drogen-bond stabilization of oxygen in hemoprotein models. I Am. Chem. Soc. 105:5165-5166.

[Reference]

Momenteau, M., and C. A. Reed. 1994. Synthetic heme dioxygen complexes. Chem. Rev. 94:659-698.

Momot, K. I., and F. A. Walker. 1997. Investigations of rotation of axial ligands in six-coordinate low-spin iron(III) tetraphenylporphyrinates: measurement of rate constants from saturation transfer experiments and comparison with molecular mechanics calculations. J. Phys. Chem. A. 101:2787-2795.

Nakamoto, K. 1990. Vibrational spectra of dioxygen adducts of metal chelate compounds. Coord. Chem. Rev. 100:363-402.

Nakamura, M., T. I\keue, S. Neya, N. Funasaki, and N. Nakamura. 1996. Fixation of the 2-methylimidzaole ligand and anomalous pyrrole chemical shifts in bis(2-methylimidazole)(meso-tetraalkylporphyrinato) iron(III) chloride caused by the nonplanar porphyrin ring. Inorg. Chem. 35:3731-3732.

Nakatsuji, H., J. Hasegawa, H. Ueda, and M. Hada. 1996. Ground and excited states of oxyheme: SAC/SAC-Cl study. Chem. Phys. Lett. 250: 379-386.

[Reference]

Oldfield, E.. H. C. Lee, C. Coretsopoulos, F. Adebodum, K. D. Park, S. Yang, J. Chung, and B. Phillips. 1991. Solid state 170 nuclear-magneticresonance spectroscopic studies of 102- 1701 picket-fence porphyrin, myoglobin, and hemoglobin. J Am. Chem. Soc. 113:8680-8685.

Olson, I S., and G. N. Phillips. 1996. Kinetic pathways and barriers for ligand binding to myoglobin. J. BW Chem. 271:17593-17596. Ormos, P., S. Szaraz, A. Cupanes, and U. N. Nienhaus. 1998. Structural

factors controlling ligand binding to myoglobin: a kinetic hole-burning study. Proc. Nad. Acad Sci. USA. 95:6762-6767.

Perdew, J. P. 1986. Density-functional approximation for the correlation energy of the inhomogeneous electron gas. Phys. Rev. B. 33:8822- 8824. Phillips, S. E. V. 1980. Structure and refinement of oxymyoglobin at 1.6

resolution. I Mol. Biol. 142:531-554.

[Reference]

Phillips, S. E. V., and B. P. Schoenborn. 1981. Neutron diffraction reveals oxygen-histidine hydrogen bond in oxymyoglobin. Nature. 292:81-82.

Ray, G. B., X.-Y. Li, J. A. Ibers, J. L. Sessler, and T. G. Spiro. 1994. How far proteins bend the FeCO unit? Distal polar and steric effects in heme proteins and models. J. Am. Chem. Soc. t 16-.162-176.

Rovira, C., K. Kunc, J. Hotter, P. Ballone, and M. Parrinello. 1997. Equilibrium geometries and electronic structure of iron-porphyrin complexes: a density functional study. J. Phys. Chem. A. 101: 8914-8925.

Rovira, C., and M. Parrinello. 1999. Factors influencing ligand binding properties of heme models@ a first principles study of picket-fence and protoheme complexes. Chem. Eur. J. 5:250-262.

[Reference]

Sage, J. T., and P. M. Champion. 1996. Small substrate recognition in heme proteins. In Comprehensive Supramolecular Chemistry. K. S. Suslick, editor. vol. 5, ch. 8. Pergamon, Oxford. 171-213,

Sage, I T., and W. Jee. 1997. Structural characterization of the myoglobin active site using infrared crystallography. J. Mol. Biol. 274:21-26. Schlichting, I., J. Berendzen, G. N. Phillips, and R. M. Sweet. 1994.

Crystal structure of photolysed carbonmonoxy-myoglobin. Nature. 371: 808-812.

Shaanan, B. 1982. The iron-oxygen bond in human oxyhaemoglobin. Nature. 296:683-684.

Shaanan, B. 1983. Structure of human oxyhaemoglobin at 2.1 A resolution. I MoL Biol. 171:31-59.

[Reference]

Sigfridson, E., and U. Ryde. 1999. On the significance of hydrogen bonds for the discrimination between CO and 02 by myoglobin. I BioL Inorg. Chem. 4:99 -110.

Slebodnick, C., and J. A. Ibers. 1997. Myoglobin models and steric origins of the discrimination between 02 and CO. J. Biol. Inorg. Chem. 2:521-525.

Spartalian, K., G. Lang, J. P. Collman, R. R. Gagne, and C. A. Reed. 1975. Mossbauer spectroscopy of hemoglobin model compounds: evidence for conformational excitation. J Chem. Phys. 63:5375-5382.

Spiro, T. G., and P. M. Kozlowski. 1998. Discordant results on FeCO deformability in heme proteins reconciled by density functional theory. J. Am. Chem. Soc. 120:4524-4525.

[Reference]

Springer, B. A., S. G. Sligar, J. S. Olson, and G. N. Phillips. 1994. Mechanisms of ligand recognition in myoglobin. Chem. Rev. 94: 699-714.

Stryer, L. 1995. Portrait of an allosteric protein. In Biochemistry. 4th ed. ch. 7. Freeman, New York.

Troullier, N., and J. L. Martins. 1991. Efficient pseudopotentials for plane-wave calculations. Phys. Rev. B. 43:1993.

Vangberg, T., D. F. Bocian, and A. Ghosh. 1997. Deformability of Fe(ll) CO and Fe(III)CN groups in heme protein models: nonlocal density functional theory calculations. J. Biol. Inorg. Chem. 2:526-530.

Unno, M., J. F. Christian, J. S. Olson, J. T. Sage, and P. M. Champion. 1998. Evidence for hydrogen bonding effects in the iron ligand vibrations of carbonmonoxy myoglobin. J. Am. Chem. Soc. 120:2670 -2671.

Yang, F., and G. N. Phillips. t996. Crystal structures of CO- deoxy- and met-myoglobins at various pH values. J. Mol. Biol. 256:762-774.

[Author Affiliation]

Carme Rovira and Michele Parrinello

Max-Planck Institut fur Festkorperforschung, Heisenbergstrasse 1, 70569 Stuttgart, Germany

[Author Affiliation]

Received for publication 20 April 1999 and in final form 20 October 1999.

Address reprint requests to Michele Parrinello, Max-Planck Institut fur Festkorperforschung, Heisenbergstrasse 1, 70569 Stuttgart, Germany. Tel.: +49-711-689-1700; Fax: +49-711-689-1702; E-mail: parrinello@prr.mpi. stuttgart.mpg.de.

[Author Affiliation]

Dr. Rovira's present address is Departament de Quimica Fisica, Universitat de Barcelona, Marti' i Franques 1, 08028 Barcelona, Spain.

(c)2000 by the Biophysical Society