INTRODUCTION

Top Abstract Introduction Methods Results Discussion Conclusions References

The fundamental role of myoglobin and hemoglobin in the understanding of protein structure and function follows from their well-defined structural features, their conformational changes, and the connection of the latter with their physiological functions (Perutz, 1979).

Myoglobin is one of the proteins that were studied most intensively in the past (Antonini and Brunori, 1971). Its simple reaction (binding with a small ligand) and the wealth of experimental information available on it make it a good system for computational studies. Detailed mechanistic pictures can be obtained from computer simulations, which (in contrast to many other biological molecules) can be tested directly against experimental data.

Recent experimental advances allowed a detailed description of short time ligand dynamics. Short times are more accessible to atomically detailed computer simulations and therefore open the way for more meaningful comparison to experiments. The low-temperature x-ray diffraction experiments or Laue spectroscopy at room temperature (Teng et al., 1994; Schlichting et al., 1994; Hartmann et al., 1996; Srajer et al., 1996) and room temperature infrared time-resolved studies (Lim et al., 1995a,b, 1997) of the photolyzed state of MbCO suggest that the CO is confined to a small region in the heme pocket, approximately parallel to the heme plane.

A number of molecular dynamics simulations, describing ligand photodissociation, protein relaxation, and geminate recombination, has been reported in the past (for a review, see Kuczera, 1997). In particular and more in the direction of the present work, the rotationally constrained environment for the CO motions in the heme pocket was investigated in several studies (Straub and Karplus, 1991; Vitkup et al., 1997; Ma et al., 1997). These studies made a number of intriguing suggestions that, in addition to the specific investigation of ligand diffusion dynamics, can be further employed in the refinement of the computational models.

In previous studies (Straub and Karplus, 1991; Ma et al., 1997) it was suggested that the spectral width of the CO vibrational excitation is dominated by electrostatic interactions, a suggestion that is further examined in the present study. Clearly, if electrostatic interactions are important, it is crucial to obtain an accurate representation of the electric field at the pocket. A widely used approximation in molecular dynamics simulations that is examined here is the use of a cutoff distance for the evaluation of the electrostatic interactions. A computational protocol that overcomes the cutoff approximation is the Ewald sum (Ewald, 1921).

In this work we computed 270 trajectories to examine the dynamics of the photodissociated CO in myoglobin. Myoglobin is embedded in a water box and periodic boundary conditions are used. Long-range electrostatic interactions are accounted for by using Ewald sum and PME (particle mesh Ewald) (Darden et al., 1993; Essmann et al., 1995). Ewald is advocated as one of the methods of choice for treating long-range electrostatic interactions. It is therefore of considerable interest to apply it to the study of CO orientation and translational diffusion in myoglobin. Furthermore, the electric field (E) at the CO is an essential quantity in the computation of spectroscopic observables. E as obtained from calculations with cutoffs is approximate, and alternative computational protocols are very much to be desired.

Despite numerous successes of the Ewald approach, it is important to emphasize that the Ewald sum imposes translational symmetry on the model. Because the experiment is usually done in solution, in which the symmetry is broken, the Ewald approach is another (although very different) approximation. Therefore, rather than seeking a "better" modeling of electrostatic interactions, a more modest (and realistic) goal would be the examination of the same phenomena with different computational and theoretical tools. Similar results obtained with different computational protocols are more trustworthy.

To further investigate the sensitivity of our results to the model used and to obtain more meaningful comparison to experiment, we also computed diffusion in a mutant of myoglobin. We consider ligand diffusion in a mutant in which Leu²⁹ was replaced by phenylalanine (Carver et al., 1992). Both residues have significant interactions with the photodissociated ligand, interactions that are primarily nonbonded.

For each of the mutants we consider two different models of the CO. In one of the models the CO molecule is presented by two point masses. Each point mass is also assigned a van der Waals radius, a well depth, and a charge. The second model uses similar van der Waals interactions, but three point chargesone additional charge is placed at the center of mass of this molecule to reproduce the quite large quadrupole moment (Straub and Karplus, 1991).

To further examine other possible factors that may influence the simulation, two different rectangular boxes are used. The boxes differ in the orientation of myoglobin in the box and in the number of water molecules. The duplicate simulations provide a useful guide to the effects of different simulation set-ups and to the convergence of the computational protocols.

The spatial and orientational distributions of the CO molecule are extracted from multiple, short (10 ps) time trajectories. We also computed the electric field and the van der Waals forces that are felt by the ligand in the pocket. The latter are the necessary ingredients for the computations of the vibrational lineshape, as discussed in the Methods section.

This manuscript is organized as follows. In the next section we describe the techniques used in the computations. In Results the computational data are provided and described. Comparison to experiments and to other computations is provided in the Discussion.

	METHODS

Top Abstract Introduction Methods Results Discussion Conclusions References

The calculations in this paper were performed using the program MOIL (Elber et al., 1994). Most of the potential energy function of MOIL is the combination of AMBER (Weiner et al., 1984) and OPLS (Jorgensen and Tirado-Rives, 1988), and the heme model is from Kuczera et al. (1989). The TIP3P water model (Jorgensen et al., 1983) is used.

Treatment of electrostatic and van der Waals interactions

Except for the reference calculations with electrostatic cutoff, the electrostatic interactions were not truncated and the Ewald sum (Ewald, 1921) is calculated by the PME method (Darden et al., 1993). The following PME setup was used. The direct-sum Coulomb cutoff that was used to generate the atom neighbor list was 9.1 Å, and the Ewald  parameter was 0.35 (in one set of calculations another choice was made, namely the direct-sum cutoff of 12.1 Å and a  of 0.26). Third-order Euler spline interpolation and a grid of 32 × 32 × 32 points were used (Essmann et al., 1995).

The van der Waals interactions that decay like 1/r⁶ (r is the distance between the two particles) were truncated at 9 Å in most of the computations. The list of neighbors was generated according to a cutoff distance larger by 2 Å, providing a buffer of neighbors to enhance the stability of the calculations. The lists were updated (recomputed) every 10 steps. To further explore the effect of the van der Waals cutoff, we also computed a set of 30 trajectories with a van der Waals cutoff distance of 11.5 Å. The complete details of all of the computations are provided in Table 1.

Models of solvation boxes and protein structures

The crystal structures of the wild type and of the Phe²⁹ mutant (Carver at al., 1992) were employed to construct two rectangular, fully solvated boxes. The first unit cell (which we denote by Box A) has dimensions 40 Å × 54 Å × 56 Å and contains one myoglobin, one C-O molecule bound to the heme iron, and 2627 water molecules (including the 334 crystallographic water molecules). SO₄ ion, which is present in the crystallographic structure, and two sodium ions (to ensure neutrality of the box) are also included.

The second box (Box B) has dimensions 50 Å × 45 Å × 50 Å. The slightly smaller volume is a consequence of the different orientation of the same protein structure in the alternative box. The number of water molecules is now 2332 (including, as before, the crystallographic water molecules).

The crystal structure of the mutant (Carver et al., 1992) and the solvation box A were employed in the study of the Phe²⁹ mutant. Water molecules (2709 molecules), one sodium ion, and one SO₄ ion are included in the Phe²⁹ simulations. For each case, the two different models of the carbon monoxide molecule and its interaction with the protein were used.

Models of the photodissociated ligand

The photodissociation process was modeled in the past as a sudden transition to the unliganded state (for a review see Kuczera, 1997). Here we employ the same approach.

The Fe-C bond was replaced in some of the computations by an exponential repulsive interaction, as was used in studies of geminate recombination (Li et al., 1993), and in others by the usual nonbonded interactions (electrostatic and van der Waals) between the heme atoms and CO molecule (Vitkup et al., 1997 and Ma et al., 1997).

In the first dissociation protocol (photodissociation model I), the usual nonbonded interactions (van der Waals) between the heme and the CO molecule are activated upon dissociation, whereas a three-site "quadrupolar" CO model is employed to model the electrostatic forces (Straub and Karplus, 1991). The electrostatic model consists of three relatively large point charges located at the carbon and oxygen nuclei (0.75 and 0.85 electron units, respectively) and at the center of mass (1.6 electron units). It was tested against ab initio and experimental data (Straub and Karplus, 1991).

In the second dissociation protocol (referred to as photodissociation model II), with the exponential potential, a model with two point charges (one on the carbon and one on the oxygen) was employed to model the electrostatic interactions (Elber and Karplus, 1990).

The two dissociation protocols are used in the present work to explore their abilities to reproduce experimental data as well as to check whether the conclusions reached are sensitive to the model that was used. It is of course easier to accept results that are less sensitive to the exact computational model or protocol (having different and reasonable models tested before against experiment; Straub and Karplus, 1991; Li et al., 1993).

Molecular dynamics protocols

Starting from the x-ray structures embedded in water boxes as discussed above, 15-ps trajectories in which the protein was fixed and only the water molecules were allowed to move were performed. The final structures were further minimized by using 1500 steps of conjugate gradient minimization with the parameters of the bound heme. Next, the structures were slowly heated from 0 to 300 K, during 60 ps of linear heating. The final structures were further equilibrated at 300 K over 20 or 40 ps. From the final equilibrium runs, we picked liganded structures to initiate trajectories of a photodissociated ligand. The time separation between different sampling points was at least 0.5 ps. Each trajectory was initiated with a different set of velocities sampled from Maxwell distribution at 300 K.

The initiation of the unbound potential (by changing the heme parameters to that of the dissociated state) was assumed to be instantaneous. In all of the simulations, the SHAKE algorithm has been employed to constrain the bonds with hydrogen atoms (Ryckaert et al., 1977; van Gunsteren and Berendsen, 1977). The equations of motion were propagated with the velocity Verlet integrator, employing a time step of 1 fs.

In the simulation, the kinetic temperature was monitored (note that the van der Waals interaction is still discontinuous because of the truncation). If the computed temperature deviated from 300° by more than 20°, the velocities were scaled. Such scaling was performed only very rarely (one scaling in 5-10 trajectories). However, in calculations that employed electrostatic cutoff distance, scaling was performed approximately once in two picoseconds. This is one of the advantages of using the PME approach (instead of cutoff), in which the energy function is continuous and differentiable. The RMS deviation from the x-ray structures was ~1.3-1.5 Å for C atoms and 1.7-2.0 Å for all of the atoms, excluding the C and the N terminal. This is higher than the RMS deviation reported in the simulations with solvation shells (Loncharich and Brooks, 1989). However, it is still within acceptable and typical ranges of molecular dynamics simulations. It is plausible that the solvation shells in which the water cannot exchange with bulk solvent molecules are more rigid as compared to the liquid phase.

For the native structure, 30 photodissociated trajectories, each 10 ps long, have been computed for Box A and Box B as well as for photodissociation model I and photodissociation model II. This gives 120 trajectories (each 10 ps long). Additional 30 trajectories were computed with an electrostatic cutoff distance of 12.1 Å to probe the effect of this approximation. We further explore the effect of the cutoff distance on van der Waals interactions by computing the next 30 trajectories in box A with a van der Waals cutoff of 11.5 Å instead of 9 Å. All of the conditions of the different computations are summarized in Table 1.

Only one box (Box A) is used for the Phe²⁹ mutant and 90 10-ps-long trajectories were used in the investigation of this mutant. Sixty trajectories were employed to study photodissociation models I and II. An additional 30 trajectories were employed in photodissociation simulations in which phenylalanines 29 and 43 of the heme pocket were modeled using explicit hydrogens. The explicit hydrogens represent the quadrupole moment of the aromatic rings (Jorgensen and Severance, 1990) and its effect on the ligand dynamics was tested as well.

Representative coordinate sets were saved each 20 fs. The coordinates were analyzed to characterize structural and dynamical features of the CO molecule and of the protein degrees of freedom that are coupled to the early diffusion. The calculated structures were further employed in computations of quantities that were extracted from experiments, such as the CO average position, orientation with respect to the heme plane, spectra and rotational dynamics. Comparison to experiment (Lim et al., 1997) and to recent computations (Ma et al., 1997; Vitkup et al., 1997) will be described in details in the Discussion.

Computations of vibrational spectra

We have attempted to compute the vibrational lineshape of the photodissociated CO, using the approach outlined by Ma et al. (1997). We extended the same model to also include van der Waals interactions and computed the linewidth as discussed below.

For each coordinate set that was saved during the trajectories, we computed the electric field and the van der Waals energy felt by the CO molecule. Of course, the motion of the CO molecule was computed by classical mechanics. We assume, however, that the time-dependent electric field and the van der Waals energy, as obtained from classical trajectories, are perturbations to the quantum mechanical state of the CO molecule. If we denote the vibrating free rotor Hamiltonian of the CO by H_o (the gas phase Hamiltonian) and the other two components by H_elec and H_vdW, the CO Hamiltonian in our representation is

(1)

where µ is the dipole moment of the CO molecule and E is the electric field. The perturbation is time dependent, because it is sampled as a function of time along the trajectory. In Fig. 1 we show the time dependence of the electric field and the CO van der Waals energy in one of the trajectories. To compute infrared (IR) spectra, an oscillating electric field D = D₀ cos(t) is also needed to induce transitions between different energy levels. The frequency of the oscillating field is scanned during the experiment.

View larger version (38K): [in this window] [in a new window]   FIGURE 1   The electric field (in atomic units), the computed spectral shift (in wavenumbers), and the van der Waals energy of the carbon atom (kcal/mol) along a single sample trajectory of the CO (the protocol includes PME, Box A of solvation, a van der Waals cutoff of 11.5 Å, and photodissociation model I). The first three panels (A-C) are for 10 ps, and the next three panels (D-F) are for 1 ps. Note the rare and sharp spikes of the van der Waals energy and the overall slow variation of the perturbations on the time scale of CO vibration (~15 fs). If the electrostatic field is computed with a 12-Å cutoff (not shown), considerably larger and more rapid fluctuations (similar in pace to the CO vibrations) are obtained.

The perturbation on the gas phase CO molecule is induced by the fluctuating protein and by the solvent environment and is time dependent. It may seem necessary to compute a time-dependent quantum evolution of the CO ro-vibrational state.

However, most of the transitions occur at the resonance frequency dictated by D. This frequency corresponds to very rapid oscillations of ~15 fs. The rapid change in the radiation field should be contrasted with the considerably slower changes in the electric field and van der Waals interactions at the CO (Fig. 1). The latter oscillate on a time scale that is slower by a factor of ~10 compared to the radiation field.

Keeping in mind the separation of time scales, it is possible to introduce an adiabatic approximation in which the electric field and the van der Waals interactions are effectively frozen during the transition. It is therefore possible to consider the perturbation as time independent during the short period of the excitation. It is also possible to consider the spectroscopy as sampling static Hamiltonians that differ in their van der Waals and electrostatic potentials.

Similar to the Born-Oppenheimer approximation, the transitions are assumed to be vertical for fixed values of van der Waals and electrostatic interactions. Hence at each specific time, say t₁, we compute the energy difference between the lowest and the first excited states of the Hamiltonian:

(2)

where the dipole moment is a function of the distance r between the carbon and the oxygen. To account for the inhomogeneous linewidth, the solution of the time-independent energy levels of the set of Hamiltonians sampled during the trajectories is required.

The empirical force field and the classical dynamics we employed are already significant approximations. It therefore seems inappropriate to seek a highly accurate quantum mechanical calculation for the spectra, and a simple model is used. H(t₁) is separated into two parts: 1) H_0V = H₀ + H_dW(t₁) and 2) µ(r)E(t₁). Part 1 is considered first.

To estimate the energy levels of H_0V, we computed the value of the corresponding potential at three points and fitted the results by a parabola. We first used the coordinates as obtained from the classical trajectory, R(t₁). In the second configuration, only the coordinates of the carbon and the oxygen were modified. They were shifted along the bond vector from the actual position of r₁ = r(t₁) to r₂ = r₀. The shift was made so that the center of mass of the diatomic molecule remained unchanged, while the relative distance of the C and the O was the exact equilibrium distance of the unperturbed CO. In the third point a similar shift was made along the bond, but this time in the direction opposite the first position, i.e. r₃ = 2r₂  r₁.

The three points are fitted to a harmonic curve, 1/2k(r  r₀)² + V₀, and the frequency appropriate for this curve is computed by  = , where M is the reduced mass of CO.

In addition to the van der Waals interactions, the electric field (contribution (b)) influences the ground and the excited states differently. The excited state wavefunction is wider, and therefore the expectation value of the dipole moment is different. In the first-order perturbation theory, the lowest energy level is corrected by µ₀E(t₁), and the first excited state is corrected by µ₁E(t), where  _i denotes the expectation value of the dipole operator at the ith quantum state of I. The energy separation between the two states is  =   (µ₁  µ₀)E(t₁) =  + µ · E(t₁).

The probability of probing a specific excitation frequency  is therefore

(3)

where the joint probability density, P_dW,E(, µ · E(t₁)/), is computed from the structures sampled in the trajectory. In practice, we found that the contributions of the van der Waals interactions are very small. Excluding very rare strong collisions, the frequency shifts that were obtained due to van der Waals interactions were on the order of 1 cm¹ and were distributed equally about the unperturbed frequency. We therefore decided not to include them in the calculations and to focus instead on the contribution of the electric fields only.

The final procedure used to compute the distribution of frequencies is the same as that of Ma et al. (1997) and is based on empirical adjustment to ab initio data. The only technical difference is that we computed the frequency shift for each structure and did not extract it from the distribution of the electric field. The direct computation added to the stability of the data.

It is intriguing to note that the fluctuations of the electric field in simulations that employ a cutoff distance are significant and rapid. The time scale of the fluctuations is comparable to the period of the oscillation of the radiation. Thus the adiabatic approximation that we used to calculate the frequency distribution is not valid for computations with cutoff. However, if we decide to ignore the test, and to consider configurations and distributions of electric fields obtained with PME (which passed the test) or with cutoff, the results are remarkably similar. This will be demonstrated in the Results.

	RESULTS

Top Abstract Introduction Methods Results Discussion Conclusions References

Below we summarize the results of the simulations with a particular emphasis on computations that are relevant to experiment.

Orientation of the CO with respect to the normal to the heme plane

After photodissociation the CO loses its excess translational and rotational energy and equilibrates very quickly. A clear demonstration to that effect is the trajectory of the CO orientational angle with respect to the heme normal, . We define the heme normal by the vector product of R(Fe_p-N_c) and R(Fe_p-N_b). N_b and N_c are two of the pyrrole nitrogens, and Fe_p is the projection of the iron coordinate vector onto the average plane defined by the four pyrrole nitrogens. The R's are the vectors connecting their positions (Fig. 2).

View larger version (17K): [in this window] [in a new window]   FIGURE 2   A schematic drawing of the coordinate system that is used to define the orientation of the carbon monoxide molecule with respect to the heme plane. n is the vector normal to the heme plane, which is defined as follows. We span the plane by the two vectors R(Fep-Nb) and R(Fep-Nc), where Nb and Nc are two of the pyrrole nitrogens and Fep is the projection of the iron position onto the plane defined by the average of the four pyrrole nitrogens.  is the angle between the vector n and the CO bond.  is the angle between the projection of the CO bond vector onto the heme plane and the vector R(Fep-Nc). Note that the carbon monoxide molecule is displaced in the figure, so that the carbon will be along the vector n. In reality, this is usually not the case. Therefore we also need the last parameter, R, which we used to characterize the position of the CO molecule. It is the distance of the geometric center of the CO molecule to the vector normal to the heme, n.

In Fig. 3, we show the average  as a function of time (the average is computed using 30 trajectories for each of the curves). We compare the statistics for three different cases: The native protein in Box A, the native protein in Box B, and the Phe²⁹ mutant. The relaxation to the equilibrium value is complete well before 1 ps to a value of ~100°, in accord with polarization experiments (Lim et al., 1997). It is probably safe to say that all of the plots are similar, and the effects of the mutation or the different models for the ligands, if any, are small.

View larger version (23K): [in this window] [in a new window]   FIGURE 3   The changes in the averaged value of  (the angle between the CO bond and the normal to the heme) as a function of time from the 10-ps trajectories. The averages are over 30 trajectories for each of the curves. Three panels are shown. In each panel we show two curves: the solid line corresponds to model I of photodissociation (three-charge models for the CO; see text for more details) and the dotted line to model II of photodissociation (two-charge model for CO; see text for more details). (A) The results for box A of solvation; (B) the results for box B of solvation; (C) computations for the phenylalanine 29 mutant, which were performed in Box A of solvation. Thirty mutant trajectories were computed with model I and with an all-hydrogen potential for phenylalanine 29 and 43. Another set of 30 mutant trajectories was computed for model II and with united atom potential. Note that all of the curves are similar and oscillate mildly between 80° and 120°. Note also the very rapid relaxation to the equilibrium position, which is much faster than a picosecond, in accord with spectroscopic measurements (Lim et al., 1997).

The adventurers may further comment that the orientation for the Phe²⁹ mutant is closer to 90°, a suggestion that can be tested, in principle, by experiment. However, it is important to emphasize that the difference between the two boxes (boxes A and B) is at least as large as the difference between the mutant and the native protein (Fig. 3).

The curves in Fig. 3 are averages over 30 trajectories for each specific time. It is useful to examine the complete distributions in addition to the averages to have an appreciation for the range of fluctuations. Joint distributions of the angle  and the distance R are shown in Fig. 4, where R is the distance of the CO geometric center from the vector n. The vector n originates from the heme center and is perpendicular to the heme plane (see also Fig. 2).

View larger version (16K): [in this window] [in a new window]   FIGURE 4   Joint distribution of the orientation angle  (the angle between the CO bond and the normal to the heme n) and of R (the distance of the geometric center of the CO molecule from heme normal n) (see also Fig. 2). Three panels are presented that correspond to different sets of simulations. We also note that similar results were obtained at other simulation conditions (see Table 1) and therefore were not reported in detail. (A) Statistics of 30 trajectories computed for box A and photodissociation model I. Note the single maximum at ~110°. The contour line with the lowest value is 50. The step between the lines is 125 and the highest contour is 675. (B) a distribution similar to that of A, except that the distribution is for the mutant Phe29. The rest of the simulation conditions are the same. The lowest contour line is 50, the highest is 650, and the step between the contours is 100. Note that the ligand distribution in R is significantly less homogeneous for Phe29 as compared to Leu29 (native). Up to 2 Å, the Phe29 distribution is narrower than the distribution of the native protein and is pressed against the normal. However, beyond 2 Å, the distribution is considerably broader than the native and reaches longer distances. This is in accord with the two very different time scales for NO recombination in this mutant (Gibson et al., 1992; Li et al., 1993). (C) The results for Leu29 in Box A, from photodissociation model II. Note the two maxima, one at a distance of ~3 Å from the heme normal and one at a distance of ~1 Å. The second maximum is more similar to that in A and B. Note also that the distribution is considerably broader. The lowest contour is 100 and the highest is 400. The step is 100. Although it is possible to associate the second maximum at 3 Å with the results of low-temperature x-ray data (Schlichting et al., 1994), we prefer not to do this at present, because the state at 3 Å does not have a preferred orientation (the CO rotates freely in the heme plane) and therefore is unlikely to be the B state reported experimentally. It is also possible that because of the weak signal, the B1 and B2 populations are hard to detect in these cases.

We show in Fig. 4 the time-averaged distributions and not snapshots in times as in Fig. 3. As is clear from Fig. 3, the value of  relaxes very quickly to equilibrium. Therefore the use of equilibrium (time-averaged) statistics is a reasonable idea to follow. The distributions have pronounced and sharp peaks at ~110 and 100°, in accord with the average values presented in Fig. 3. Nevertheless, the spread of the population is quite large, and fluctuations of more than 60° are also found.

The orientation of the CO molecule in the plane parallel to the heme ()

This analysis of the angle  (see also Fig. 2) is of significant interest, because it suggests a clear structural picture of the so-called B states. The B states were detected spectroscopically at low temperature (Alben et al., 1982) and were suggested as different docking sites in the heme pocket by Lim et al. (1997). However, detailed structural interpretation at room temperature was not available. Using our extensive sampling, we are able to identify two distinct orientations of the CO molecule in the plane parallel to the heme, which are suggested as the B states based on the available data.

The straightforward identification of the B states from room-temperature trajectories is in accord with the qualitative proposition made by Lim et al. (1995a,b, 1997). They proposed that the B states correspond to similar locations of the CO molecule in space, in which the positions of the carbon and the oxygen are exchanged. This is exactly the picture that we have obtained.

In Fig. 5 we show a one-dimensional histogram of the probability density as a function of , the orientational angle in the plane. We note that the doubly peaked distribution is obtained by using a variety of conditions, including different van der Waals cutoffs, different models of the CO, and even cutoff for electrostatic interactions. However, in the simulations that employed Box B, we were able to detect only one of the states. This perhaps is another warning, reemphasizing the difficulties in obtaining proper statistics in condensed phase simulations. A two-dimensional contour plot (Fig. 6) that includes the angle  and the distance from the heme normal places the two B states just above the heme.

View larger version (19K): [in this window] [in a new window]   FIGURE 5   A one-dimensional histogram of the probability density of the angle , the angle of the CO in the heme plane (see also Fig. 2). Each of the curves is constructed from 30 trajectories of 10 ps. Six different curves are shown, four in A and two in B. Except for the dotted line in A (curve 2), all of the other curves clearly show two maxima, corresponding to two B states. The curves are (for more details see Table 1 and text): curve 1 (solid line, A), Leu29, Box A, PME, photodissociation model I, van der Waals cutoff 11.5; curve 2 (dotted line, A), Leu29, Box B, PME, photodissociation model I, van der Waals cutoff 9.0; curve 3 (dashed line, A) Leu29, Box A, PME, photodissociation model I, van der Waals cutoff 9.0; curve 4 (dashed-dotted line, A) Leu29, Box A, electrostatic cutoff 12.1, photodissociation model I, van der Waals cutoff 9.0; curve 5 (solid line, B) Phe29, Box A, PME, photodissociation model I, van der Waals cutoff 9.0; curve 6 (dotted line, B), Leu29, Box A, PME, photodissociation model II, van der Waals cutoff 9.0. Note that the peaks of the distribution are significantly less pronounced in B (a different scale). Note also that individual trajectories can be found in different sets that occupy only state B1 or only state B2. Experiments (Alben et al., 1992; Lim et al., 1997) suggest that the barrier for rotation is ~1 kcal/mol. The ratio of the populations at the barrier and at the minima should therefore be ~1:5. A is closer to experiment than B and has a population ratio of ~1:3.

View larger version (18K): [in this window] [in a new window]   FIGURE 6   A joint distribution for the angle of the CO in the heme plane, , and the distance from the normal to the heme, R. Three panels are shown. (A) A distribution extracted from 30 trajectories computed in Box A, using PME, photodissociation model I, and 11.5 Å for the van der Waals cutoff. The lowest contour represents 100 events and the step between the contours is 50. The highest level is 300. (B) Similar to A, except that Box B is used and the lowest contour is of 50 events. (C) Also similar to A, except that photodissociation model II is used and the van der Waals cutoff is 9 Å. The lowest contour in C is of 50 events.

A probability density in the two Cartesian dimensions that define the heme plane (see Fig. 2) is shown in Fig. 7, demonstrating again that the area just above the heme iron is mostly populated. This is excluding Fig. 7 C, which corresponds to box B simulations that did not populate the B2 state. The lack of the B2 state helps to identify it as being closer to the heme normal.

View larger version (13K): [in this window] [in a new window]   FIGURE 7   The distribution of the geometric center of the CO in the heme plane (see Fig. 2 for the definition of the plane). Note the drift of the distribution, which is peaked above the heme iron between Nb and Nc. (A) Distribution extracted from 30 trajectories computed in Box A using PME, a van der Waals cutoff of 11.5 Å, and photodissociation model I. The lowest contour is of 100 events and the step is 400. (B) The same as in A, except that an electrostatic cutoff of 12.1 Å is used. Note the significant similarity between the two distributions of A and B. (C) The same as A, except that Box B is used and the van der Waals cutoff was 9 Å. (D) The same as C, except that photodissociation model II is employed and the contour step is of 150 only. The maximum in D is of 550 events, significantly smaller than the 1300 events observed in all of the other panels. Note also that the scale for C and D is enlarged compared to A and B, reflecting the broader distribution in these cases.

A new and intriguing feature is the "leak" of the density between N_b and N_c, suggesting, perhaps, a preferred direction for an escape from the heme pocket.

Computations of IR spectra

In Fig. 8 we show the IR spectra computed according to the procedure outlined in Methods. A series of lineshapes is presented. Each spectrum was extracted from a set of 30 trajectories, using different simulation parameters (see Table 1 for more details). As is clear from the figure, the spectra are broad and similar in shape. The width is comparable to the experimental measurements (Fig. 9). However, in contrast to experiment, two distinct experimental peaks for the two B states are hard to detect.

View larger version (19K): [in this window] [in a new window]   FIGURE 8   Computed frequency distribution for the ro-vibrational excitation of the CO molecule. Our computed spectrum has an arbitrary origin and was set to fit the center of the experimental curve. See text for the details of the computations. Three curves are shown, and all of the computations were performed for the native protein: curve 1 (solid line) was extracted from a set of 30 trajectories using PME and a van der Waals cutoff of 11.5 Å; curve 2 (dotted line) was computed with an electrostatic cutoff distance of 12.1 Å and a van der Waals cutoff of 9 Å; and curve 3 (dashed line) is similar to curve 1, except that a van der Waals cutoff of 9 Å was employed.

View larger version (18K): [in this window] [in a new window]   FIGURE 9   Comparison of the computed vibrational line shape (dotted line, curve 1 of Fig. 8) to experimental data (solid line, Lim et al., 1997). Our computed spectrum has an arbitrary origin and was set to fit the center of the experimental curve. Although the widths of the two lines are similar, we do not observe separate features of B1 and B2. If the two populations of the two states are separated and the spectrum is calculated for each one of them, the spectrum in the lower panel is obtained. The relative shift, 5 cm1, is computed as the difference between the average frequencies of the B1 and the B2 peaks and is too small compared to the experimental number (~10 cm1). The computations are approximate, and such an error is definitely possible and is hard to reduce, considering the approximations involved. We therefore did not attempt to refine the spectrum further.

In an attempt to separate the contributions of the two states to the line width, conformations of the B1 and the B2 states were collected separately and their corresponding frequency distribution was accumulated (Fig. 9). The lines are broad but are clearly shifted with respect to each other. The separation between the average frequencies of the B1 and B2 states is quite small, however ~5 cm¹, which is smaller than the experimental separation of ~10 cm¹.

The smaller separation and the significant width of the frequency distributions are summed to one broad peak in which the B1 and B2 features are difficult to separate, in contrast to experiment. Nevertheless, the correct overall width, the observation of some shift, and the structural interpretation, which is consistent with x-ray data (see next section), makes, in our view, a suggestive case for the elusive B states of the photodissociated CO. The level of the theory that we employed to compute the frequency distribution can certainly result in errors of a few wavenumbers.

It is of considerable interest to find the location of the two states in the heme pocket. As we have already noted, the distribution of structures within each state is quite wide. Therefore, to present a clear picture, we show in Fig. 10 CO positions that correspond to the most probable configurations. The heme and the additional residues were taken from equilibrated bound structures that were used to initiate the photodissociation runs. As a result (for example), the iron is not yet displaced from the heme plane.

View larger version (20K): [in this window] [in a new window]   FIGURE 10   Structural interpretation of the B states. From the set of 30 Leu29 trajectories that employed PME, photodissociation model I, Box A, and a van der Waals cutoff of 11.5 Å, we extracted the most probable configurations in , , and R. In A.1 and A.2 we considered the B1 state. A.1 shows a projection of the CO molecule onto the heme plane and A.2 a side view. Note that the coordinates of the heme and the other residues are taken from the bound form. From the side view it is clear that the oxygen is pointing to the heme iron, which is consistent with the suggestion made by IR spectroscopy (Lim et al., 1997). The spatial location of the molecule between Leu29, Ile107, and Val68 (Val68 is not shown) is consistent with the room-temperature x-ray data (Srajer et al., 1996). The B2 state (B.1 and B.2) is more parallel to the heme, and the CO is even closer to the heme normal.

Despite the fact that the iron is not displaced in the displayed figure, it does reach the equilibrium position very rapidly during the simulation. In Fig. 11 we show the relaxation of the iron, out of the plane for the two boxes and for the Phe²⁹ mutant. The plots are averaged over 30 trajectories for each point in time.

View larger version (20K): [in this window] [in a new window]   FIGURE 11   The relaxation of the iron out of the heme plane (defined by the average of the four pyrrole nitrogens) immediately after dissociation. The curves are averaged (as usual) over 30 trajectories each, and are remarkably similar, independent of the type of the box or of the mutant.

There are four panels in Fig. 10. The first two correspond to the state B1 and the next two to the state B2. In B1 the oxygen is closer to the iron and is somewhat tilted toward it with an angle  of ~120°. The CO molecule is displaced between the pyrrole nitrogens N_b and N_c in a direction consistent with the time-resolved x-ray data (Srajer et al., 1996).

In the B2 state, a nonnegligible but smaller spatial displacement of the ligand (compared to B1) in the N_b and N_c directions is apparent. The B2 state further places the CO molecule in essentially parallel orientation with respect to the heme plane. The largest difference compared to the B1 state is the 180° change in the  angle (see also Fig. 2), a rotation that results in the carbon being closer to the iron atom, as opposed to the oxygen atom in B1.

Our structures are not that far from the n vector. Nevertheless, it is possible that our "states" are not final and that within the 10 ps of the simulations we were only able to probe the ligands as they diffused toward their final destinationthe "true" B states, if the true states are indeed significantly more displaced from the n vector (see discussion in the next section). Of course, the present B states already capture many of the characteristics expected from the B states and discussed by Lim et al. (1997).

Spatial distributions of the CO in the heme pocket

Recent publications of low-temperature or time-resolved x-ray data of the photodissociated state of CO in myoglobin (Schlichting et al., 1994; Teng et al., 1994; Srajer et al., 1996; Hartmann et al., 1996) made it possible to probe the ligand as an intermediate before its final destination. There are low-temperature (Schlichting et al., 1994; Teng et al., 1994; Hartmann et al., 1996) and room-temperature diffraction experiments (Srajer et al., 1996), and comparison to all of them is, in principle, warranted.

However, the experiments that were done at different conditions do not agree on the exact position of the ligand (see discussion by Hartmann et al., 1996). A useful comparison between experiment and theory without reproducing the experimental conditions is therefore difficult. Furthermore, it is not obvious how the low-temperature data or the nanosecond data from room temperature are related to picosecond trajectories at room temperature, which do not correspond to any of the experimental set-ups.

We therefore decided to focus on qualitative and general features that are shared by the different experiments rather than on a detailed comparison of the structure to a particular experiment. In all of the experiments, the ligand was located above the heme and displaced between N_b and N_c, in accord with our structures at Fig. 10. The low-temperature structures suggest that the ligand is displaced even further from the vector nresults that are consistent with photodissociation model II (Fig. 4 C). Nevertheless, it is still useful to examine the results of photodissociation model I (see Discussion). Given the short time of the simulation, in which the ligand did not diffuse very far, it is suggestive to probe the tails of the distribution for plausible positions that the ligand may adopt at later times. The tails are examined in Fig. 7, which demonstrates the general direction of the diffusion toward the back of the heme, similar to what has been observed in crystallography. It should be noted that the spectroscopy finds that the B states are established extremely quickly; however, this does not exclude some small further diffusion within the pocket somewhat later.

In Fig. 12 we show the distance L of the geometric center of the ligand from the heme plane. The distribution is quite wide and is almost independent of the existence (or nonexistence) of an electrostatic truncation scheme (cutoff or PME). However, it does depend on the model of CO. Photodissociation model II is associated with a significantly broader distribution than that of photodissociation model I, an observation that holds under considerable variation in simulation parameters. Note also that the average and most probable positions extracted for the mutant Phe²⁹ are similar to those of the native protein (Leu²⁹).

View larger version (17K): [in this window] [in a new window]   FIGURE 12   The distribution of the distances from the geometric center of the CO molecule to the heme plane. Only the results for the native protein (Leu29) in box A are shown. Each curve is a histogram generated from a sample of 30 trajectories of 10 ps. Curve 1 employed PME and photodissociation model I, curve 2 used PME and photodissociation model II, and curve 3 used an electrostatic cutoff distance of 12.1 Å and photodissociation model I. Clearly the most profound effect was induced by changing the model of photodissociation. Additional calculations (not shown) showed a small effect of the mutation Leu29 to Phe29.

	DISCUSSION

Top Abstract Introduction Methods Results Discussion Conclusions References

Here, molecular dynamics simulations of the short time dynamics of a photodissociated ligand were performed. Modern computer technology enables us to obtain detailed and apparently converging statistics for numerous properties of the protein and the ligand. The present calculations can therefore serve as comprehensive tests of model potentials, of computational protocols, and of previous approximate theories or other numerical studies. They can further serve as an attempt to correlate the results of different experiments, such as time-resolved diffraction and spectroscopy, information that is not always easy to integrate.

X-ray diffraction and the infrared spectroscopy experiments are, in a sense, complementary. Whereas the latter can provide the distribution of the orientations of the CO molecule, the x-ray data provide the relative position of the CO with respect to the heme. On the other hand, the time resolution of the diffraction experiments lags behind that of the spectroscopic measurements. So far, three different positions of CO have been reported (Teng et al., 1994; Schlichting et al., 1994; Srajer et al., 1996). The first two experiments were done at low temperature. It is not obvious whether they are connected to the fast spectroscopy at room temperature. The room-temperature spectroscopy data have femtosecond time resolution.

A significant advance in resolving room-temperature questions and coming closer to the conditions of the IR spectroscopy was made by Srajer et al. (1996), who solved a structure at room temperature at nanosecond time resolution. However, the time scale is quite long compared to the initial rearrangement and docking of the ligand in the pocket.

The advantage of computer simulations that is exploited here is the comprehensive picture that one obtains for different phenomena. We should bear in mind (of course) that the simulations are approximate. However, subject to experimental verification, a unified model that takes into account many experimental results can be constructed.

We suggest that the docking site of the ligand is (on the few picoseconds time scale) almost parallel to the heme plane, and slightly displaced in the direction of the N_b and N_c pyrrole nitrogens. The two clearly separated features, which we suggested to be the B1 and B2 states, are the head-to-tail switching arrangements of the CO molecule in the heme plane. This is in agreement with the qualitative picture of Lim et al. (1997). The CO position that we observed is similar to the "R" state of the Laue x-ray diffraction experiment (Srajer et al., 1996) that was measured at room temperature. The CO is locked between Leu²⁹, Val⁶⁸, and Ile¹⁰⁷ and is in the direction of the C ring.

We emphasize that in Fig. 10 a single structure is shown. The single conformation reflects the most probable distance from the vector n, and the most probable orientation angles  and . The actual distribution can be quite wide. We used the most probable value (instead of an average) because the x-ray experiment is more likely to reveal the most probable position if the distribution is broad and most of it is not above the noise level.

Other recent studies of the CO dynamics in the heme pocket were published: the work by Vitkup et al. (1997) and Ma et al. (1997). The simulation of Vitkup et al. includes 28 trajectories to collect statistics employing a solvation shell of 492 molecules. It is not clear whether this solvation is sufficient to reproduce the electric field that the CO experiences in the pocket. The computations were successful in reproducing a variety of features of the low-temperature x-ray data but did not address spectroscopy. This is in contrast to the computations of Ma et al. (1997), which included full, atomically detailed solvation and focused on spectroscopic measurements.

The study of Ma et al. (1997) is more similar to ours in the sense that the solvation effects were considered in detail. They differ from our computations (for better or worse) because the Ewald sum was not employed and the electrostatic interactions were either truncated or modeled by continuum. The study of Ma et al. is the foundation of many of the analyses made here, and poses (and answers) many interesting questions. However, because of limited statistics (of four trajectories only), it was difficult to extract the elusive properties of the B states in that study.

Another new aspect of the present manuscript is the comparison of many simulation protocols. We have tried different treatments of electrostatic and different models of the ligand in the hope of being able to say clearly (using the known experimental results) which of the models performed better. To our surprise, it seems that the use of cutoff for electrostatic interactions did not have a profound effect on the structural features of the ligand in the pocket. The PME approach has numerous advantages, such as continuous and differentiable energy, which is conserved significantly better than the electrostatic energy computed with a cutoff distance of 12 Å. However, the statistical properties of the structures and even of the electric field are not too different (if the cutoff is maintained at 12 Å!). The electric field fluctuates more violently and rapidly between updates of the nonbonded lists, but the distribution is not that much different. This is, perhaps, a sign of relief for simulations of proteins that employ cutoff.

We also examined the difference between the three-site model of the CO molecule (Straub and Karplus, 1991) and the earlier model of CO that did not take into account the quadrupole moment. The results of the two models studying the ligand diffusion are quite different, but not sufficiently different to completely exclude one of them. The distribution of the center of mass and the orientational angle  of the two-site model are considerably broader than the corresponding distribution of photodissociation model I. The distribution of the angle in the heme plane, , shows two maxima for the two B states. However, the distribution is considerably broader and flatter in photodissociation model II compared to photodissociation model I.

It is of interest to mention that within the broad distribution of the center of mass of the two-site model of CO, one can find maxima that are consistent with all of the x-ray data. This is similar to the results of Vitkup et al. (1997). Perhaps ignoring the quadrupole moment of the CO (and setting the electrostatic interactions almost to zero) has an effect similar to that of the removal of the solvent and eliminating long-range electrostatic forces, leading to a more rapid diffusion of the small ligand inside the protein. A similar distribution was also observed in the study of Ma et al. (1997).

Whereas the wider distribution is advantageous from the perspective of alternative x-ray sites, it is less favorable spectroscopically. The spectroscopy suggests well-defined states, and the three-site model yields more focused and peaked distributions. Because the simulations of short time dynamics are performed in conditions more similar to the conditions of the spectroscopy experiment, we tend to favor the three-site model (photodissociation model I) and the application of an extensive solvation shell.

	CONCLUSIONS

Top Abstract Introduction Methods Results Discussion Conclusions References

In this paper we presented detailed and extensive calculations of the early events that follow the photodissociation of CO in myoglobin. The results were tested for convergence by a variety of means. Different computational models were also examined.

We have shown that the simulations are able to bridge the results of different experiments and to provide a detailed explanation of experimental features.

Our results, although definitely not in perfect agreement with experiment, are suggestive, and are similar to experimental results or have similar experimental trends for a number of independent observations, such as the orientation of the CO with respect to the heme plane, the existence of the B states with tail-to-head exchange, and the location of the B states above the heme and between Leu²⁹, Val⁶⁸, and Ile¹⁰⁷. It is the collective qualitative agreement between the different approaches that (so we hope) gives the present study more support and makes the emerging picture more consistent.