BACKGROUND

TOP ABSTRACT BACKGROUND COMPUTATIONAL MODEL AND METHODS RESULTS SUMMARY AND CONCLUSIONS APPENDIX A APPENDIX B REFERENCES

The flow of excess vibrational energy into or out of vibrational modes is important to many reactions in which there is bond breaking or formation. One such reaction is carbon monoxide (CO) dissociation and rebinding in heme proteins. Investigation of CO relaxation in heme proteins can provide vital information concerning the cooperative nature of ligand binding and rebinding, in particular, and protein dynamics, in general.

Vibrational relaxation processes of carbonyl compounds have received considerable attention with the advent of ultra-fast time resolved laser methods (King et al., 1993). The population relaxation times for metal carbonyls vary considerably from molecule to molecule, ranging from the picosecond to nanosecond timescales. This demonstrates that the local environment influences the vibrational relaxation time (T₁) of the carbonyl. To date, there have been several experimental probes into the vibrational relaxation of CO in model heme compounds (Hill et al., 1995) and heme proteins (Hill et al., 1994; Owrutsky et al., 1995). These studies have focused on the population relaxation of the bound CO, which, in the heme protein, myoglobin, is called the A-states.

The vibrational lifetime of bound CO in myoglobin is strongly dependent on the structure of the surrounding protein. Vibrational relaxation of bound CO in heme proteins is very fast compared to other CO-ligated compounds (Hill et al., 1995; Owrutsky et al., 1995). For example, Owrutsky et al. (1995) have found that the CO vibrational lifetime in heme proteins is ~20 ps, while King and coworkers (1993) have found a 150-ps time scale in the relaxation of the CO stretch in first row transition metal carbonyl compounds.

To justify the fast relaxation of bound CO in heme proteins, it has been suggested that the manifold of states within the heme are responsible (Hill et al., 1996a). Vibrational relaxation tends to be more efficient through strong covalent interactions, and this appears to be true in MbCO. Although there is only one such interaction between the CO and the heme, studies indicate that intermolecular transfer of vibrational energy is relatively unimportant in the relaxation process of bound CO (Hill et al., 1996a).

Interestingly, vibrational lifetimes also differ within the same protein but in different conformational substates (Hill et al., 1994). Hill and coworkers (1996b) have shown that changes in the vibrational lifetime of bound CO are correlated with changes in the frequency of the CO molecule that accompany such transitions. Moreover, Fayer and coworkers have found that the T₁ population relaxation time for the carbonyl ligand is only slightly dependent on temperature, decreasing ~10% from 50 to 300 K. They concluded that the low frequency bath modes do not play a significant role in the relaxation of bound CO (Hill et al., 1994).

In contrast to this detailed understanding of CO energy relaxation in metal carbonyls and the bound A-states of carboxy myoglobin, little is known about the relaxation time scale or mechanism in the photolyzed state. After dissociation of the CO from the heme at room temperature, the CO quickly finds itself in one of two B-states. These states are identified by a change in the CO vibrational frequency from ~1960 to ~2131 cm¹ or ~2119 cm¹ for B₁ and B₂, respectively. Anfinrud and coworkers have investigated the B-states and contend that the CO is in one of two orientations within a binding site located in the heme pocket (Lim et al., 1995).

The CO in the B-state has no covalent interactions. Any contribution of the heme or protein to its relaxation may give a picture of the role of the protein solvent in the relaxation of the CO in the bound and photolyzed states. This paper describes a computational study of the CO population relaxation and mechanism of the photolyzed state using molecular dynamics and normal mode analysis of the fully solvated, room-temperature myoglobin system.

Calculation of the vibrational relaxation time (T₁)

The process of vibrational relaxation involves the dissipation of excess vibrational energy into the surroundings. Zwanzig (1973) showed that, when a system can be represented by the generalized Langevin equation, the time decay of the vibrational energy relaxation is a single exponential (the Landau-Teller result)

(1)

where the brackets ··· indicate an average over initial zero time states. Using this formalism along with the assumption that the oscillator is anharmonic and is linearly coupled to the solvent, an expression for T₁ (Oxtoby, 1979; Oxtoby, 1981; Zwanzig, 1961) can be obtained

(2)

where µ is the reduced mass

(3)

and ₀ is the frequency of the oscillator. The time dependent friction, (t), is a measure of the dissipation force along the oscillator's vibrational coordinate.

Time correlation functions of harmonic quantum systems can be directly related to those of the corresponding classical systems. Bader and Berne (1994) have shown this result in the context of vibrational relaxation that, if the system can be described by a set of harmonic normal modes, the relaxation time obtained in a fully quantum system (quantum solute in a quantum bath) is the same as that obtained from a fully classical simulation (classical solute in a classical bath). Furthermore, they suggest that this result may also be applicable when anharmonicities are present. Therefore, it is conceivable that accurate values for vibrational population relaxation times in myoglobin can be obtained from classical molecular dynamics simulations, provided the bath and system-bath coupling are well described as harmonic.

The time-dependent friction is proportional to the time correlation function of the fluctuating random force on the bond. When the bond is of high frequency, it can be considered rigid, and molecular dynamics are run with the oscillator (CO molecule) constrained at its equilibrium distance (Berne et al., 1990). The friction kernel in Eq. 2 can then be easily computed in the rigid bond approximation where (t) can be written,

(4)

F(t) = F(t)   is the fluctuating solvent force on the oscillator bond. The force correlation function includes the effects of the density of states and coupling strength of the surrounding solvent. The scalar force along the bond is determined using

(5)

where F_i is the force felt by atom i of the CO molecule due to the solvent, and _co is the CO bond unit vector.

In the case of CO, the assumptions made above are not unreasonable. CO is ground-state dominated, and the bottom of its bond potential can be accurately modeled as harmonic. Furthermore, in myoglobin, the coupling between the protein and the free CO is believed to be weak. Finally, the CO bond is a high frequency oscillator (₀ ~ 2140 cm¹), for which the rigid bond approximation holds. (This limit has also been found to be reasonable for low frequency oscillators, an indication of the robustness of the approximation (Berne et al., 1990).) The applicability of equilibrium molecular dynamics in the study of relaxation processes is supported by work done by Whitnell (1990, 1992), in which both equilibrium and nonequilibrium techniques were used to determine vibrational relaxation times. With each method, a similar value for T₁ was found, suggesting that the vibrational relaxation is adequately described by the linear response approximation.

Normal mode analysis

To examine the mechanism of energy relaxation, it is essential to be able to identify those protein modes that are most strongly coupled to the CO bond stretching mode. This can be accomplished by a normal mode analysis based on quenched normal modes (QNM) or instantaneous normal modes (INM). The set of QNM is computed by a diagonalization of the Hessian matrix for a set of equilibrium configurations minimized prior to diagonalization. The QNM method offers a straightforward way to separate and examine the coupling between the system and bath modes.

A good deal of attention has been given to the use of INM in the determination of short time dynamical properties of simple solutes in liquids (Cho et al., 1994; Goodyear and Stratt, 1996, 1997; Ladanyi and Stratt 1998). The INM method involves diagonalizing the Hessian matrix for a set of instantaneous configurations taken from a dynamical trajectory (Seeley and Keyes, 1989). Because these configurations are not necessarily at the potential minimum, there are both real frequencies corresponding to stable oscillations of the system and imaginary frequencies corresponding to unstable motion of the system. These imaginary modes tend not to couple strongly to solute molecules (Wann and Stratt, 1994). (As a result, the imaginary modes were excluded in the calculation of the friction kernel presented in this work.) The INM formalism can be used to describe short-time dynamics. In fact, if the time increment t is infinitesimally small, the INM method is exact. The INM method can also be used to calculate the fluctuating friction along a bond as shown by Goodyear and Stratt (1996, 1997). When using a normal mode model to calculate the friction along a vibrational coordinate, we begin with Zwanzig's model of an anharmonic system oscillator r, bilinearly coupled to a bath of harmonic oscillators x_i. In other words,

(6)

where

(7)

Changing to mass-weighted coordinates, the Hamiltonian becomes

(8)

Zwanzig derived the time-dependent friction from the Hamiltonian (Eq. 8) and found

(9)

The sum in Eq. 9 is over the normal modes of the bath oscillators i with mode frequency _i. The key to Eq. 9 is the coupling constants, C_i, between the bath coordinates and the oscillator (CO) stretching coordinate. When considering vibrational relaxation, the couplings should be a measure of how the bath modes affect the force along the CO bond,

(10)

where q_i are the mass-weighted normal mode coordinates. The force along the CO bond is

(11)

where V is the interaction potential between the carbon monoxide molecule and the surroundings. The expression for the coupling constants can then be written with respect to the mass-weighted Hessian matrix used in the determination of the normal modes

(12)

(13)

The x_j/q_i terms are the coefficients of the eigenvector matrix of the normal modes, which will be denoted as U_i,j.

Normal mode calculations can also be used to examine the role of collective motions in the dynamics of the system. We compute the participation ratios defined as

(14)

and

(15)

The sum in Eq. 14 is over all 3N degrees of freedom in the system. In Eq. 15, the summation is divided into two components. The sum over l extends over the M residues in the NM calculation; the second sum over j extends over the 3N_l degrees of freedom of residue l.

The inverses of R_i^I and R_i^II are a measure of the degree of localization of each mode. 1/R_i^I corresponds to the number of atoms participating in the ith mode, while 1/R_i^II corresponds to the number of residues and water molecules participating in the ith mode. If mode i is completely delocalized, 1/R_i^I should be a good estimate of the number of degrees of freedom of the system, 3N, and 1/R_i^II should be a good estimate of the number of residues and water molecules that contribute to that mode. If a mode is completely localized so that only one local mode participates in the mode, only one of the eigenvector coefficients would be nonzero. Because the eigenvectors are normalized, that coefficient and 1/R_i^I would be equal to unity. Conversely, if the mode is fully delocalized, each degree of freedom would be equally involved. Therefore,

(16)

and

(17)

In this work, the classical equilibrium molecular dynamics method is used to calculate the T₁ time of CO in the heme pocket of myoglobin. Second, using normal mode techniques, the coupling of the protein-solvent environment to the CO oscillator has been investigated. From these couplings, a friction kernel has been calculated and compared with the molecular dynamics results. The coupling constants have been used to identify bath modes important to the mechanism of CO vibrational relaxation in myoglobin.

	COMPUTATIONAL MODEL AND METHODS

TOP ABSTRACT BACKGROUND COMPUTATIONAL MODEL AND METHODS RESULTS SUMMARY AND CONCLUSIONS APPENDIX A APPENDIX B REFERENCES

The sperm whale myoglobin molecule (Schlichting et al., 1994) was placed in a 56.70 × 56.70 × 37.712 Å³ box and simulated using periodic boundary conditions. A previously equilibrated box of TIP3P water molecules (Jorgensen et al., 1983) was overlayed and any water molecule whose oxygen atom was within 2.5 Å of a nonhydrogen protein atom was removed. The all-hydrogen parameter set (version 25) within CHARMM (MacKerell et al., 1992) was used. The system, in its entirety, was composed of 11,499 atoms; 2551 myoglobin atoms, 2982 water molecules, and the two atoms of CO. To minimize any bad contacts or unusual strain, ~1500 steepest descent energy minimization were performed.

Molecular dynamics

The temperature was increased slowly to 300 K using the following protocol. From 0 to 150 K, the temperature was increased 10% after every 2-ps of molecular dynamics by randomly sampling the atomic velocities from a Maxwell distribution. The criteria for velocity resampling was based on a 5% temperature window. From 150 to 215 K, the temperature was scaled by 5% and the window became 2.5% of the current temperature. When the temperature reached 215 K, the temperature window became constant at 5 K with an interval of 4 ps between velocity resamplings.

After a temperature of 300 K was reached, 20 ps of constant temperature molecular dynamics was run, in which the temperature was checked every 10 fs. During the last 10 ps, the average temperature remained within the 5 K window and there was no need to resample the atomic velocities. At this point, it was assumed that an equilibrium state had been reached and data could be collected from constant energy dynamics. Microcanonical trajectories of varying lengths (10-20 ps) were run for a total of 350 ps for the -tautomer and 200 ps for -tautomer. The entire 350 ps of the  trajectories were used in determining the T₁ times; however, only 150 ps were used in the decomposition of the force autocorrelation function and its comparison with the -tautomer. The molecular dynamics timestep was 1.0 fs and the intermolecular potential was truncated at 10.5 Å using a group switching function from 9.5 Å.

The three-site model of Straub and Karplus (Ma et al., 1997; Straub and Karplus, 1991), in which charges are placed on the carbon and oxygen atoms as well as the center of mass, was used to represent the CO. This simple charge model, when combined with the Huffaker RRKR potential for the bond (Huffaker, 1976), has been shown to reproduce the experimental dipole and quadrupole moments of the CO molecule, as well as ab initio interaction energies of CO with a variety of molecules, including water and formamide. Interactions between the CO molecule and the protein-solvent environment involve a Lennard-Jones contribution and a three-site electrostatic term,

(18)

where Q_i is the charge on site i and r_ij is the distance between sites i and j. The Lennard-Jones parameters, _ij and _ij, are determined via (Lorentz-Berthelot combining rules)

(19)

The first term in Eq. 18 describes the Coulombic interactions in which the second summation is over the three sites of the CO molecule (the carbon and oxygen atoms and the center of mass). The second term is the van der Waals contribution in which the summations involve interactions between the carbon and oxygen atoms of the CO molecule and the surrounding protein and solvent atoms.

The CO bond was constrained to its equilibrium position (1.128 Å) using the SHAKE algorithm. The fluctuating force autocorrelation function was calculated, and its Fourier transform was carried out using fast Fourier transform (FFT) techniques (Press et al., 1989).

Normal modes

The density of states of a given system can provide insight into possible modes available for vibrational relaxation of the CO molecule, and is given by

(20)

To identify relaxation pathways for the CO molecule in the myoglobin pocket, several quenched and instantaneous normal mode calculations were performed on subsets of the solvated myoglobin system. Several configurations were picked from equilibrium trajectories. For each of these configurations, any residue whose center of mass was outside a 14.0-Å radius from the center of mass of the CO molecule was removed. For the QNM calculations, this system subset was then minimized using the adopted basis Newton-Raphson (ABNR) method (Brooks et al., 1983). To help maintain the system integrity in the vicinity of the CO molecule, a harmonic restoring force was applied to all atoms beyond a 12.0-Å radius from the CO molecule. The above distances were chosen after several tests indicating that a 14.0-Å cutoff was sufficient to achieve converged minimized geometries near the CO as well as the density of states. The reduced system was used to reduce the associated computational overhead. The harmonic restoring force constant was set at 1.0 kcal/(mol Ų). The INM calculations were performed on the system subset without additional constraints. The system subset was diagonalized without potential energy truncation or periodic boundary conditions.

	RESULTS

TOP ABSTRACT BACKGROUND COMPUTATIONAL MODEL AND METHODS RESULTS SUMMARY AND CONCLUSIONS APPENDIX A APPENDIX B REFERENCES

In this section, values of T₁ are estimated, normal mode methods are used to determine important doorway modes in the protein and solvent, and the dominant contributors to the relaxation process are identified.

CO population relaxation times

As described above, relaxation times of high frequency oscillators can be directly related to the Fourier transform of the fluctuating force-force autocorrelation function F(0)F(t) of the force along a rigid bond. In determining the vibrational relaxation time, the value of the friction kernel at the frequency of the oscillator is used. The friction kernel () as a function of frequency calculated from the above simulations for the  and -tautomers is shown in Fig. 1. The arrow points to the frequency of the free CO bond and the inlay depicts the time dependence of F(0)F(t). Using this method and

(21)

the relaxation time of the CO oscillator was estimated.

View larger version (17K): [in this window] [in a new window]   FIGURE 1   The frequency-dependent friction kernels, (), for the carbon monoxide ligand in the heme pocket of sperm whale myoglobin ( and -tautomers) are shown. Displayed in the inlay is the fluctuating force autocorrelation function for the carbon monoxide ligand stretch. These functions were averaged over 350 () and 175 () ps of CO dynamics in the heme pocket at 300 K. The  marks the CO oscillator frequency at 2100 cm1. The () data have been smoothed for clarity.

The spectrum of the time-dependent friction was somewhat noisy. To remove some of the noise, the spectra were smoothed by locally averaging over the data points. This provided an average value, , at the frequency of the oscillator, as well as a standard deviation, , for each point. We related the standard deviation to the uncertainty in our estimate of the relaxation time as

(22)

Table 1 lists the values obtained by averaging over 2, 5, and 10 points. As can be seen, the values obtained are all within the calculated error of each other. We chose to use a value of 10 for averaging over all spectra presented in this work.

View this table: [in this window] [in a new window]   TABLE 1   CO vibrational relaxation times computed for the  and  tautomers of sperm whale myoglobin

In Table 1, we present several vibrational relaxation times for the various B-states. Because interconversion between B₁ and B₂ is a spontaneous process at 300 K, the force along the CO molecule calculated is probably due to both B₁ and B₂ substates. Therefore, care should be taken not to attribute too much significance to the differences in the T₁ times for the different substates. The results reported are an average over the substates and should be viewed as such. Nonetheless, these values for the vibrational relaxation time are in very good agreement with the experimental value of ~600 ps determined by Anfinrud (submitted for publication). It is interesting to note that the T₁ time calculated for the -tautomer is in closer agreement with the experimental result than that for the -tautomer. Because the time scale for tautomerization is much longer than that for the photolyzed CO to enter the B-states, our results suggest that the -tautomer may be present at the time of photolysis. Although this interpretation contradicts the neutron crystal data of Schoenborn (Hanson and Schoenborn, 1981; Cheng and Schoenborn, 1991), recent resonance Raman work indicates the presence of the H in carbonmonoxy myoglobin (Unno et al., 1998).

The use of Eq. 2 and the equilibrium molecular dynamics methods illustrated in this paper have proven to be an effective method in determining population relaxation times for a number of other systems (Gnanakaran and Hochstrasser, 1996; Whitnell et al., 1992), demonstrating the wide applicability of the method and supporting our use of the Landau-Teller approach.

From the data in Table 1, it appears that the T₁ time of the B-state CO population relaxation is considerably slower than the ~17 ps found for CO in its bound state. (Owrutsky et al., 1995). This suggests that the Fe-CO bond and the heme itself contribute to A-state CO relaxation. This comes as no surprise. Owrutsky and co-workers (1995) and Hill et al. (1996a, b) have suggested that -backbonding plays a role in the relaxation process. In studies of vibrational relaxation of model heme compounds, Fayer and coworkers (Hill et al., 1996a, b) go further in saying that anharmonic coupling through the  bond provides the dominant pathway of relaxation.

Normal mode calculations

The density of states determined using the QNM and INM formalisms is shown in Fig. 2. The direct contribution of the restraining force used in the QNM calculation has been removed. The resulting peak appears at ~155 cm¹ and does not change any of the observations mentioned. As expected, the INM spectrum possesses imaginary modes, plotted here in the standard way along the negative frequency axis. These imaginary frequencies make up approximately 5% of D(). Although Lennard-Jones fluids exhibit much higher relative numbers of imaginary modes (Wan and Stratt, 1994; Wu and Loring, 1992; Madan et al., 1990; Madan and Keyes, 1993; Seeley et al., 1991), the result presented here is similar in magnitude to results seen in crystals or ordered liquids such as liquid water (Cho et al., 1994). The most noticeable difference between the two methods is an apparent smearing of states in our INM calculation. This behavior is due in part to the thermal energy of the atoms, a feature lacking in the QNM calculation. Perhaps more importantly, the spreading of states seen in Fig. 2 is an indication of the increasing importance of anharmonicities as trajectories move up the potential energy surface. This is especially evident at lower frequencies, in particular, from 1200 to 1600 cm¹. In both spectra, there is a clear separation of states between 2000 and 2800 cm¹. This separation effectively isolates the CO vibration from the remainder of the system, as a result CO vibrational coupling to the system is weak.

View larger version (19K): [in this window] [in a new window]   FIGURE 2   The density of states as defined by Eq. 20 from QNM and INM calculations. The data have been smoothed for clarity and normalized so that they have unit area.

Figure 3 a and b, displays the inverse participation ratios for the INM calculations. The low frequency modes are quite delocalized with the lowest frequency modes corresponding to translational and rotational motion. With increasing frequency, the modes become more localized. Modes from 2000 to 3200 cm¹ involve only 1-2 residues. At higher frequencies, we see a decrease in localization due to the numerous OH and NH stretching modes. Nonetheless, the overall degree of localization remains considerable.

View larger version (19K): [in this window] [in a new window]   FIGURE 3   The inverse participation ratios as determined from the eigenvectors of the QNM and INM calculations described in the text. (a) The number of degrees of freedom involved in a given mode. (b) The number of residues or water molecules involved in a given mode. The data have been smoothed for clarity.

The distribution of the square of the coupling constants found in Eq. 10, also called the influence spectrum, is shown in Fig. 4 for the -tautomer. The most noticeable peak is at a frequency of ~3515 cm¹ corresponding to the N-H stretch of the distal histidine, His⁶⁴. The second most prominent peak can be seen in the region of 1500-1700 cm¹. The modes found here are a combination of angle bending and bond stretching motions with His⁶⁴ playing a key role. The influence spectrum for the -tautomer exhibited similar couplings to the vibrational modes of the distal histidine (with the exception of the absent N-H stretch).

View larger version (18K): [in this window] [in a new window]   FIGURE 4   The square of the vibrational coupling constants against a backdrop of the INM density of states for myoglobin at 300 K. The results above are averaged over ten configurations taken from separate molecular dynamics trajectories. The spike at ~2100 cm1 corresponds to the CO vibrational frequency.

Although some solvent coupling is evident between 3200 and 3300 cm¹, the effect appears to be minimal. This information, combined with that from the participation ratios shown in Fig. 3 a and b, suggest that the principle modes responsible for CO relaxation in myoglobin are highly localized. Large scale collective motions are relatively unimportant in the relaxation process. This is in agreement with the interpretation of experiments done by Hill et al. (1996a), in which they attributed high frequency, localized modes to the relaxation of the CO ligand in myoglobin. Furthermore, the residues most strongly involved in the relaxation tend to be those closest to the CO molecule. Hard collisions are most efficient and this is done via high frequency localized modes near the CO. This behavior is similar to that found in liquids, in which the fluctuations in the first solvation shell of an oscillator were found to contribute most strongly to vibrational relaxation.

Given that the protein environment is not a liquid, there is no first solvation shell per se. The surrounding protein is more like a flexible cage. The closest residues act as the first solvation shell, and their flexibility enhances the relaxation process. Table 2 lists the residues most strongly coupled to the CO. From Table 2, the role of side chain flexibility is evident. For example, the heme is the closest residue to the CO, however, its coupling to the CO bond is small. This is due to the more or less rigid nature of the heme, which lacks the necessary flexibility for relaxation of the free CO. Interestingly, the coupling of the heme evident near 1090 cm¹ involves a collective of bond stretches, angle bends, torsions, and improper torsions of the heme.

Using the INM approach, we have calculated a friction kernel for our system. It is pictured in Fig. 5 and compared with the direct result from the MD simulation. It should be stressed that the INM calculations were performed for 10 configurations and may not be converged. The density of states for the different configurations showed little variation from one to another. However, the influence spectra and fluctuating force autocorrelation functions varied considerably. This is in agreement with the work of Goodyear and Stratt (1996), who have demonstrated how dramatically INM friction spectra can differ from configuration to configuration. The oscillations in the INM function and its zero time value are a result of the lack of convergence.

View larger version (14K): [in this window] [in a new window]   FIGURE 5   The friction kernels from the INM and molecular dynamics calculations at 300 K for the -tautomer.

The INM method was used in hopes of gaining qualitative insight behind the relaxation and not for the purpose of assigning definitive values relating to it. Considering the underlying approximations, it is impressive that the INM friction kernel approximates the time dependence of the full MD trajectory average so well. For example, there are clearly two time scales involved in the relaxation. The first is a result of nearest neighbor collisions and the second is derived from longer range collective motion. Furthermore, the zero time value, although not correct, is of the right order of magnitude.

Testing the assumption of a harmonic bath

The molecular dynamics results presented here are encouraging in that they agree quite well with experimental measurement. A simple test to probe the validity of the harmonic approximation in treating the bath is given by the distribution of the fluctuating force along the bond. If the harmonic approach is appropriate, this distribution should be Gaussian. The result of this test can be seen in Fig. 6, in which a Gaussian fit to the data has been overlayed. As can be seen, the data, although not strictly Gaussian, do exhibit some Gaussian character and are reasonably approximated by a Gaussian distribution shifted to the left by 0.6 kcal/(mol Å). This is consistent with the good agreement between the Landau-Teller estimate for T₁ and the experimental result.

View larger version (15K): [in this window] [in a new window]   FIGURE 6   Comparison of fluctuating force distribution of the -tautomer with a Gaussian. The Gaussian fit was performed about the center of the original distribution and not about zero.

The 0.6 kcal/(mol Å) shift in Fig. 6 indicates an imbalance between forces that compress the bond and those that act to extend it. Which forces act to stretch the bond and which act to compress it? If the force along the bond is less than zero, the contributing force acts to expand the bond. Conversely, if the force along the bond is greater than zero, the effect is one of compression. Table 3 lists the force along the CO bond for the - and -tautomers. From this data, it is clear that the Coulombic force tends to stretch the CO bond whereas the LJ force tends to compress it.

View this table: [in this window] [in a new window]   TABLE 3   The average Coulombic and Lennard-Jones forces felt along the CO vibrational coordinate for both the - and -tautomers of His64 in myoglobin at 300 K

The decomposition of the fluctuating force autocorrelation function

Further insight into the relaxation mechanism of the CO molecule in myoglobin can be gained via the decomposition of the fluctuating force autocorrelation function. Assuming that the potential can be expressed as pairwise interactions, the friction kernel may be decomposed into Lennard-Jones and Coulombic contributions as

(23)

where F_QQ and F_LJ represent the Coulombic and Lennard-Jones contributions to the force. This should be possible due to the long range, long time scale Coulombic interactions and the short range, short time scale changes in the Lennard-Jones interactions. A similar decomposition may be accomplished by residue or entire segments of the system. For example, we can write

where F_prot, F_heme, and F_solv indicate the force felt along the CO bond due to the protein, the heme and the solvent respectively.

The independent binary collision model (IBC) (Litovitz, 1957) has been used with success to describe vibrational relaxation in solution. This simple model views the events contributing to vibrational relaxation as separate independent collisions. Although there are criticisms of the model (Harris et al., 1990), it is proven to be useful in systems in which many bodied or long range interactions are minimal. This is not always the case in a protein bath. However, the CO molecule is relatively nonpolar molecule in a relatively nonpolar pocket. It may be that the IBC model is appropriate for this system in which the solvent is a protein. A more detailed decomposition by residue allows us to test the IBC model in which it is assumed that the effect of the cross terms is negligible. Although cross terms may contribute significantly to the friction, this may not be the case in the frequency domain. Therefore, the IBC model should be tested after the transform has been done and not before. In short, the IBC model assumes that

(24)

where the sum over i is over the residues of the system and () is the Fourier transform of the friction kernel (t). Such an analysis can also give clues as to the relative importance of specific regions of the protein in accepting excess vibrational energy from the CO molecule.

The fluctuating force distribution for the - and -tautomers is given in Fig. 7. Somewhat surprisingly, the distributions of the two tautomers are almost identical in height, width, and shift. Yet the distributions of the Coulombic and van der Waals components for the two tautomers differ. For example, the -tautomer's electrostatic force distribution is broader than the corresponding -tautomer's distribution. Therefore, although the magnitude of the average electrostatic force along the CO bond is greater for the -tautomer (as seen in Table 3), the fluctuations about this average are smaller. As a result the zero-frequency contribution to the friction kernel should be larger for the -tautomer (Gnanakaran and Hochstrasser, 1996). The Lennard-Jones force distributions of the two tautomers are fairly similar. They differ slightly in shift and height but are roughly equivalent in their widths. The long tail of the force distributions seen in Fig. 7 is due to close-range repulsive interactions felt by the CO molecule in the myoglobin pocket. Interestingly, the differences between the Coulombic and van der Waals distributions in the two tautomers balance each other so that the overall distributions are identical.

View larger version (14K): [in this window] [in a new window]   FIGURE 7   The 300 K distribution of the fluctuating force felt along the CO bond for the - and -tautomers as labeled in the figure.

Although the vibrational relaxation of CO in myoglobin is due to fluctuations in both the electrostatic and Lennard-Jones forces, the major contribution is from relatively weak Lennard-Jones forces. This is demonstrated in Fig. 8 a and b, where the electrostatic and Lennard-Jones contributions to the friction spectrum for the two tautomers are presented. Clearly, the major contribution is due to the van der Waals interactions. This is not necessarily surprising. The CO molecule has a very small permanent dipole moment and moderate permanent quadrupole moment. However, it is likely that the electrostatic effects play a much larger role in the process of pure dephasing.

View larger version (27K): [in this window] [in a new window]   FIGURE 8   The decomposition of the friction spectra for the (a) - and (b) -tautomers into Coulombic and Lennard-Jones contributions. Also shown for comparison are the total spectra. The data were smoothed for visual clarity.

To learn more about the vibrational mechanism of CO in the myoglobin pocket, several decompositions of the fluctuating force autocorrelation function were performed for the -tautomer. The first involved separating the correlation function according to segments of the system. In one such decomposition the system was divided into three segmentsthe protein, the heme, and the solvent. From Fig. 9 a, it is obvious that the self terms of the protein, heme, and solvent very closely reproduce the total friction spectrum. Any cooperative interactions between these groups is negligible.

View larger version (28K): [in this window] [in a new window]   FIGURE 9   A look at the protein, heme, and solvent contributions on the total friction spectrum of the CO bond at 300 K for the -tautomer. (a) The total spectrum and the spectrum obtained by decomposing the force autocorrelation. (b) The friction spectra of each of the components as labeled in the figure. All spectra have been smoothed for visual clarity.

A second decomposition based on individual residues was performed to investigate the applicability of the IBC model. Cross correlations, although contributing to the friction kernel, have little influence on the friction spectrum in the vicinity of the CO vibrational frequency. Fig. 10, a and b, demonstrates this fact. This indicates that the relaxation of the CO molecule takes place via a series of successive, uncorrelated collisions with residues in the myoglobin pocket.

View larger version (22K): [in this window] [in a new window]   FIGURE 10   (a) The fluctuating force autocorrelation function of the protein on the CO bond and the friction kernel calculated from the self terms with respect to the individual residues calculated at 300 K for the -tautomer. (b) The corresponding friction spectra. The spectra have been smoothed for visual clarity.

The pocket residues

The residues that line the myoglobin pocket are Leu²⁹, Leu³², Phe⁴³, Val⁶⁸, His⁶⁴, and Ile¹⁰⁷. Because these residues act as nearest neighbors to the CO molecule, their respective contributions to the friction spectrum were assessed. The resulting spectra are pictured in Fig. 11 a-f. It is clear from Fig. 11 that the tautomeric state of the distal histidine does not affect the friction exerted by Phe⁴³ and Ile¹⁰⁷ on the CO molecule. Although the average distances of the CO molecule from Phe⁴³ in each of the tautomers are comparable, this is not so for Ile¹⁰⁷. The CO tends to be considerably further from Ile¹⁰⁷ in the -tautomer pocket than in the -tautomer pocket. Residues Leu²⁹, Val⁶⁸, and the heme contribute more to CO relaxation in the -tautomer than they do in the -tautomer.

View larger version (25K): [in this window] [in a new window]   FIGURE 11   The individual components of the friction spectrum due to (a) Leu29, (b) Leu32, (c) Phe43, (d) Val68, (e) Ile107, and (f) the heme. The solid and dashed lines represent the - and -tautomers, respectively. The simulations were carried out at 300 K and the data were smoothed for visual clarity.

The dominant contribution to the relaxation for each of these components results from van der Waals contacts. The drop in the heme friction spectrum for the -tautomer shown in Fig. 11 f is also due to a subsequent drop in the van der Waals fluctuations, as well as the electrostatic contributions. This is interesting because, on the average, the CO molecule is closer to the heme in the -tautomer (see Table 4).

The observed distance dependence could be due to different locations of the CO molecule in the pocket or a change in pocket geometry resulting from the tautomeric state. We found that the average distance between the center of mass of Ile¹⁰⁷ and the heme was ~10.3 Å for both tautomers. The distance between the N and the Ile¹⁰⁷ center-of-mass was found to be ~1.5 Å closer in the -tautomer pocket than in the -tautomer pocket. The average distances within the pocket indicate that the CO molecule in the -tautomer has more room in which to move about. To confirm this possibility, the coefficient of diffusion D of the center-of-mass of the CO molecule was calculated. This was done noting that, at intermediate times,

(25)

where D is the diffusion coefficient, R(t) is the position of the center-of-mass of the CO molecule at time t. In the limit of hard sphere collisions, the mean collision frequency  can be extracted easily from the diffusion coefficient of the Enskog model via

(26)

where  = 1/k_bT and m is the mass of the oscillator. Using Eq. 25, diffusion coefficients of 4.53 × 10⁶ cm²/s and 1.43 × 10⁶ cm²/s were calculated for the  and -tautomer, respectively. This translates to mean collision frequencies of 200 and 635 ps¹. These values indicate that the CO molecule in the -pocket does experience more translational restraint and undergoes more collisions per unit time than the CO in the -tautomer pocket. This extra mobility of the CO molecule in the -tautomer pocket may be a cause of the drop in the heme contribution to the friction as is addressed in the Appendix.

The role of the distal histidine 64

The direct effect of the distal histidine on CO dynamics and relaxation was examined. Table 5 lists the average distances of the carbon and oxygen of the CO molecule from the N and, if applicable, the corresponding H. The CO molecule is considerably closer to N in the -tautomer, with the oxygen atom being closer than the carbon atom. Naively, we might find that this is in accord with a simple point charge model fit to the experimental dipole moment showing the oxygen to be the positive pole of the CO molecule. However, ab initio quantum chemical calculations have demonstrated that it is the midbond of the CO molecule that is most strongly attracted to the imidazole's N in the gas phase, in which the minimum distance was found to be 3.46 Å (Straub and Karplus, 1991).

View this table: [in this window] [in a new window]   TABLE 5   The average distance of the carbon and oxygen atoms of the CO molecule from the N and H of the distal histidine, His64, at 300 K

In Fig. 12, the histidine components of the force correlation functions for both tautomers are presented. As expected, the zero-time fluctuation for the -tautomer is larger than for the -tautomer. As is shown in Fig. 12 b and c, the fluctuations of the pure electrostatic components are similar 0.413 versus 0.375 kcal/(mol Å) for the - and -tautomers, respectively. However, the Lennard-Jones components differ considerably 0.075 versus 0.392 kcal/(mol Å). This seems to imply that cross correlations are significant because the overall zero-time fluctuation for the -tautomer is higher. The Coulombic and Lennard-Jones fluctuations are positively correlated in the -tautomer, but are actually anticorrelated in the -tautomer counterpart. The origin of this effect stems from the CO interaction with the N of the histidine.

View larger version (31K): [in this window] [in a new window]   FIGURE 12   (a) The friction kernel contributions of the distal histidine in its two tautomeric states. (b) The electrostatic component of the friction kernel due to the distal histidine in its two tautomeric states. (c) The Lennard-Jones component of the friction kernel due to the distal histidine in its two tautomeric states.

In the -tautomer, the carbon and oxygen atoms each feel a strong electrostatic attraction to the H. The CO tends to align itself with either its carbon or oxygen atom pointing toward the nitrogen. This is demonstrated in Fig. 13 a, in which the N CO angle is plotted as a function of time, indicating that there is a barrier to rotation. This flipping behavior is reminiscent of the proposed two-state orientational behavior of CO in its photolyzed state (Ma et al., 1997; Meller and Elber, 1998).

View larger version (33K): [in this window] [in a new window]   FIGURE 13   The angle formed by N, C, and O atoms as a function of time for the (a) - and (b) -tautomers at 300 K.

The or