J. Phys. Chem. B 2007, 111, 3243-3250
3243
Molecular Dynamics Study on the Solvent Dependent Heme Cooling Following Ligand Photolysis in Carbonmonoxy Myoglobin Yong Zhang,† Hiroshi Fujisaki,‡ and John E. Straub* Department of Chemistry, Boston UniVersity, Boston, Massachusetts 02215 ReceiVed: September 8, 2006; In Final Form: December 6, 2006
The time scale and mechanism of vibrational energy relaxation of the heme moiety in myoglobin was studied using molecular dynamics simulation. Five different solvent models, including normal water, heavy water, normal glycerol, deuterated glycerol and a nonpolar solvent, and two forms of the heme, one native and one lacking acidic side chains, were studied. Structural alteration of the protein was observed in native myoglobin glycerol solution and native myoglobin water solution. The single-exponential decay of the excess kinetic energy of the heme following ligand photolysis was observed in all systems studied. The relaxation rate depends on the solvent used. However, this dependence cannot be explained using bulk transport properties of the solvent including macroscopic thermal diffusion. The rate and mechanism of heme cooling depends upon the detailed microscopic interaction between the heme and solvent. Three intermolecular energy transfer mechanisms were considered: (i) energy transfer mediated by hydrogen bonds, (ii) direct vibration-vibration energy transfer via resonant interaction, and (iii) energy transfer via vibration-translation or vibrationrotation interaction, or in other words, thermal collision. The hydrogen bond interaction and vibrationvibration interaction between the heme and solvent molecules dominates the energy transfer in native myoglobin aqueous solution and native myoglobin glycerol solutions. For modified myoglobin, the vibration-vibration interaction is also effective in glycerol solution, different from aqueous solution. Thermal collisions form the dominant energy transfer pathway for modified myoglobin in water solution, and for both native myoglobin and modified myoglobin in a nonpolar environment. For native myoglobin in a nonpolar solvent solution, hydrogen bonds between heme isopropionate side chains and nearby protein residues, absent in the modified myoglobin nonpolar solvent solution, are key interactions influencing the relaxation pathways.
1. Introduction In myoglobin (Mb) or hemoglobin (Hb), ligand dissociation can occur when the ligand-heme complex absorbs a visible or UV photon, which can cause vibrational excitation of the ligand, heme, and the surrounding residues.1,2 Understanding the time scales and mechanisms of vibrational energy relaxation (VER) that follow such events is an essential component of an understanding of the ultrafast conformational changes and the reorganization of protein structures. Much experimental and theoretical work has been done on this subject.3-14 Time-resolved Raman spectra of the heme protein myoglobin were measured by Martin and co-workers following a 500 nm pulse excitation.15 It was found that the excess energy deposited in the heme is transferred to protein modes in roughly 5 ps. Within 15 ps, the excess energy of the heme was fully dissipated. Li, Sage, and Champion identified a heme cooling time constant of 4 ps on the basis of Raman scattering experiments.16 Lim and co-workers studied heme cooling following photoexcitation of heme in myoglobin, and observed a single-exponential decay with a time constant of 6.2 ( 0.5 ps.17 Kitagawa and co-workers monitored the mode specific behavior of heme vibrational relaxation using time-resolved resonance Raman (TR3) spectroscopy.12,13,18,19 The decay time constants were found to be 1.1 ( 0.6 ps for the ν4 band and 1.9 ( 0.6 ps for the ν7 band, implying a thermal decay of the heme within 2 ps.18,19 *
[email protected]. †
[email protected]. ‡
[email protected].
Hochstrasser and co-workers measured the D2O stretching mode shift following photoexcitation of heme in myoglobin solution.20 It was determined that 60% of the deposited energy was transferred on a time scale of 7.5 ( 1.5 ps, whereas 40% of the energy relaxed on a longer time scale of 20 ps. Miller and co-workers studied the energy flow between the protein and solvent using transient phase grating spectroscopy,21 demonstrating that energy transfer to the surrounding solvent occurred within 20 ps. Computer simulations of relaxation in heme proteins began with the pioneering study of myoglobin and cytochrome c in vacuum by Henry and co-workers.22 They observed that the heme cooled within 20-40 ps; approximately 50% of the excess energy relaxed within 4 ps. The heme cooling was subsequently studied using molecular dynamics simulation by Straub and coworkers.23-27 Decay of the heme kinetic energy following ligand dissociation was found to be single exponential with a time constant of 5.9 ps in both native type myoglobin and H93G mutated myoglobin.23,25 When the two isopropionate side chains are replaced by aliphatic hydrogen, the relaxation time constant increased to 8.8 ps. The strong electrostatic interaction between the isopropionate side chains and the solvating water molecules was conjectured to be the dominant pathway for dissipation of the heme’s excess kinetic energy.23 This conjecture was supported by subsequent experimental studies on the H93G mutant and similar modified heme variants of myoglobin.28
10.1021/jp065877k CCC: $37.00 © 2007 American Chemical Society Published on Web 03/06/2007
3244 J. Phys. Chem. B, Vol. 111, No. 12, 2007
Figure 1. Structure of carbonmonoxy myoglobin (MbCO). The eight helices and several residues forming the heme pocket were labeled.
For cytochrome c, the kinetic energy relaxation was found to be a biphasic exponential decay process with relaxation time constants of 1.5and 10.1 ps for the fast process and slow process, respectively. Energy flow from the heme to protein was found to occur by “through bond” and “through space” mechanisms.24 As a model of heme dynamics in proteins, metalloporphyrins have been widely studied.29-33 Intriguing mode specific properties were reported by Kitagawa and co-workers.34,35 The authors found that the ν4 and ν7 bands of nickel octaethylporphyrin (NiOEP) relax according to a double exponential with the fast time constant being 2.6 ( 0.5 ps for several nonpolar solvents. A solvent dependent relaxation rate was observed for iron(III) meso-tetrakis(4-sulfonatophenyl)porphyrin (FeTSPP), which is soluble in polar solvents. In water solution, the anti-Stokes ν4 band showed a single-exponential decay with a time constant of 1.9 ( 0.4 ps. In methanol/benzene mixtures, an additional slow decay phase was observed, with a time constant dependent on the methanol/benzene ratio. Overall, these studies suggest that vibrational energy relaxation in a heme or porphyrin complex is sensitive to both the heme side chain topology and the local solvation environment. In this work, we present the results of a computational study of heme cooling in a wild type heme myoglobin and a modified heme myoglobin, in a variety of solvent environments including normal water, heavy water, normal glycerol, deuterated glycerol, and a nonpolar solvent. The observed solvent dependence of the kinetic energy relaxation rate is explained by considering three microscopic intermolecular energy transfer mechanisms: (i) energy transfer via hydrogen bonds interaction, (ii) direct vibration-vibration energy transfer via resonant interaction, and (iii) energy transfer via vibration-translation, or vibrationrotation interaction. The relative importance of these pathways is found to depend sensitively on the heme-protein topology and nature of the heme-solvent interactions. 2. Computational Model and Methods 2.1. Molecular Dynamics. Sperm whale myoglobin36 (Protein Data Bank entry 1MBC; see Figure 1) was simulated with periodic boundary condition using the CHARMM program.37 A previously equilibrated solvent box was laid over the protein and any solvent molecule whose heavy atom was within 2.5 Å of a non-hydrogen protein atom was removed. The excess potential energy due to bad contacts and strain was then reduced using the steepest-decent energy-minimization method. The allhydrogen protein (version 22) parameter set38 of the CHARMM simulation program was used. Five solvent models, normal water (aqueous solution), heavy water, normal glycerol (glycerol-h8), deuterated glycerol (glycerold8) and a nonpolar solvent, were used in this work. The CHARMM version of the TIP3P water model39 was used for both normal water and heavy water simulation with the masses
Zhang et al.
Figure 2. (a) Structure of heme for wild-type myoglobin, in which the positions of two isopropionate side chains can be seen. (b) Structure of heme for the modified myoglobin, in which the two isopropionate side chains are amputated and replaced by hydrogen atoms (light blue, Fe; blue, N; red, O; gray, C; white, H).
of hydrogen atoms doubled in heavy water. Normal glycerol parameters consistent with the all-atom CHARMM force field were used.40 The same parameters were used for deuterated glycerol, with an increase in mass for all H atoms. An extended atom model41 was used as the nonpolar solvent, in which the atoms have no charge or dipole moment and so interact only through the van der Waals term in the interaction potential. The molecular dynamics was simulated using the Verlet algorithm42 with a time step of 1.0 fs. The nonbonded potential was truncated using a group switching function extending from 9.5 to 11.5 Å. For simulations in nonpolar solvents, weak harmonic constraints (about 6 kcal/mol/Å2) were added to all protein CR atoms to maintain the protein’s native structure. The temperature was increased slowly to 300 K. A molecular dynamics trajectory was run at constant pressure and constant temperature for normal water solution and normal glycerol solution, respectively, after which the volume of the box was found to fluctuate about a constant mean value. The water box was found to be 71.19 × 71.19 × 54.96 Å3 containing 8985 solvent molecules, and the glycerol box was 81.41 × 81.41 × 59.24 Å3 containing 2789 solvent molecules. The final parameters of the normal water box were used for heavy water and the nonpolar solvent box; those of the normal glycerol box were used for deuterated glycerol box. A modified heme myoglobin was used in the studies, in which the heme’s two isopropionate side chains were amputated and replaced by hydrogen atoms (see Figure 2). A molecular dynamics trajectory was run at constant energy until the new equilibrium for this modified myoglobin was reached. For each system, 20 ps of constant temperature molecular dynamics was run in which the temperature was checked every 0.2 ps for the first 10 ps and every 2.0 ps for the other 10 ps. 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. The system was allowed to evolve for an additional 200 ps. Coordinates were saved every 20 ps for a total of 10 configurations. Each of these configurations was then run for 20 ps, and resampling the velocities occurred every 0.2 ps for the first 10 ps and every 2.0 ps for the other 10 ps. Ten 30 ps microcanonical trajectories were generated for both the nonequilibrium (photolyzed) and the equilibrium (non-photolyzed) states. The initial configuration for each photolyzed trajectory was identical to the corresponding nonphotolyzed trajectory at the moment of photolysis with the exception of the velocities of the heme atoms. 2.2. Simulation of Photon Absorption by the Heme. To simulate the absorption of a photon, approximately 88 kcal/ mol of excess kinetic energy, equivalent to two 650 nm photons, was deposited in the heme. It was found that the intramolecular
Energy Relaxation of the Myoglobin Heme Moiety
J. Phys. Chem. B, Vol. 111, No. 12, 2007 3245
vibrational energy relaxation within the heme was sufficiently rapid so that the results were insensitive to the exact means of distributing the excess energy within the heme.25 In this work, the excess kinetic energy was distributed uniformly among all the heme atoms, following which, the temperature increased by 400 K for native heme and 500 K for modified heme. In carbonmonoxy myoglobin, the heme chromophore iron atom is hexacoordinate with the His93 on the proximal side and ligand CO on the distal side. Following photolysis of CO, the heme becomes pentacoordinate and the iron atom moves out of the porphyrin plane by roughly 0.3 Å.43,44 The recently refined pentacoordinate heme parameters45 were used to describe the heme following photolysis. 2.3. Measure of the Rate of Local Kinetic Energy Relaxation. The ergodic measure is a useful means of determining the time scale for the self-averaging of a given property in a many-body system.26,46-48 A particular form of the ergodic measure, the fluctuation metric, is defined as N
Ω(t) )
∑j [fj(t) - hf (t)]2
(1)
where fj(t) is the time average for atom j of property F(t) of the system and hf(t) is the average of property F(t) over all N atoms at time t, hf(t) ) (1/N)∑Nj fj(t). Here we take Fj(t) to be the kinetic energy of the jth atom. If the total system is selfaveraging, then the function ΩKE(t) will decay to zero as ΩKE(t)/ΩKE(0) ≈ 1/DKEt f 0. The slope of ΩKE(0)/ΩKE(t) is proportional to the generalized diffusion constant DKE for the kinetic energy. The mean square of the fluctuations in the kinetic energy metric is related to the generalized diffusion constant describing the rate of exploration within the kinetic energy state space.46,47 If we assume that the microscopic motion is well described by a Langevin dynamics model, it is possible to relate the rate of self-averaging of the atomic kinetic energy, as measured by the ergodic measure, to the static friction of a Langevin model as DKE ) γ0, where γ0 is the static friction.26 It follows that at long times the slope of the ratio ΩKE(0)/ ΩKE(t) is approximately the static friction constant γ0. The computation of the kinetic energy metric is an effective means of assigning time scales for atomic kinetic energy relaxation in the molecular systems and numerically more accurate than integrating over the corresponding velocity autocorrelation function.26 2.4. Normal Modes Analysis. Normal-mode analysis (NMA) can be used to analyze the detailed nature of the coupling between protein and solvent bath modes based on quenched normal modes (QNM) or instantaneous normal modes (INM). In each case, the normal modes are computed by a diagonalization of the Hessian matrix for a set of configurations. The density of states of a given system was calculated using
D(ω) )
1 3N
N
〈
∑i δ[ω - ωi]〉
(2)
as previously described.49 For the QNM, the configuration is optimized to the nearest local minimum of the potential energy, and the normal-mode analysis is performed for the quenched state. The INM analysis is carried out on the snapshot configuration itself. Much attention has been given to the use of INM in the determination of the short time dynamical properties of solutes in liquids49-55 and proteins.56 Because the configurations in INM calculation are not necessarily at the potential minimum, there are both real
Figure 3. Exponential decay of the excess kinetic energy in the heme following ligand photolysis. The simulation data are represented by points; the best single-exponential fit is represented by a line. (a) Native heme and (b) modified heme as shown in Figure 1. The solvents are, from top to bottom, normal water, heavy water, normal glycerol, deuterated glycerol, and a nonpolar atomic solvent model. The decay time constants are summarized in Table 1.
frequencies and imaginary frequencies. The imaginary modes were found not to couple strongly to solute molecules.57 The INM were calculated for each solvent in this work. The imaginary modes were ignored when the heme cooling properties were analyzed. 3. Results and Analysis 3.1. Solvent Dependent Kinetic Energy Relaxation Rate of the Excited Heme. Simulated nonequilibrium dynamics trajectories were used to determine the time scales of kinetic energy relaxation or cooling of the heme in myoglobin. The average kinetic energy deposited over the heme atoms during the excitation process was 88 kcal/mol, which led to an increase of over 400 K for native heme and an increase of over 500 K for modified heme. The average time dependence of the excess kinetic energy decay in the heme for each case is shown in Figure 3. The simulation data are well fitted by a singleexponential decay function (see Table 1). Using the refined pentacoordinate heme parameters,45 the aqueous solution trajectories were regenerated. The kinetic energy relaxation time constants are 5.1 ( 0.2and 8.3 ( 0.3 ps for the native heme and the modified heme, respectively. A small change was observed from our previous simulation results of 5.9 ( 0.2 ps for native heme25 and 8.8 ( 0.3 ps for modified heme.23 For each solvent system, the kinetic energy relaxation time for the native heme is faster than for the modified heme. For
3246 J. Phys. Chem. B, Vol. 111, No. 12, 2007
Zhang et al.
TABLE 1: Excess Kinetic Energy Relaxation Time Constants in Heme Following Ligand Photolysis Simulated at 300 Ka time constant (ps) solvent
native heme Mb
modified heme Mb
normal water heavy water glycerol-h8 glycerol-d8 nonpolar
5.1 ( 0.2 6.5 ( 0.3 4.7 ( 0.1 4.3 ( 0.2 6.7 ( 0.3
8.3 ( 0.3 8.5 ( 0.4 7.3 ( 0.2 6.7 ( 0.3 7.7 ( 0.3
a The simulation data are well fitted by a single-exponential function in each case.
TABLE 2: Kinetic Energy Diffusion Constants of Solvents Derived from Kinetic Energy Metric Calculation Together with Experimental Thermal Diffusivities and Viscosities. solvent
DKE, ps-1 (300 K)
normal water heavy water glycerol-h8 glycerol-d8 nonpolar
9.8 9.4 2.6 2.7 3.6
diffusivity,a
thermal 10-8 (J/m3)/K 14.3 12.7c 9.5
viscosity,b mPa s (298 K) 0.890 1.132d 934
a The listed values of thermal diffusivity κ were calculated from the thermal conductivity λ, specific heat cp, and density F using formula κ ) λ/cpF. The experimental values for λ, cp, and F were taken from the CRC Handbook of Chemistry and Physics, 86th ed. (CRC Press: Boca Raton, FL, 2005-2006). b Experimental data was taken from the CRC Handbook of Chemistry and Physics, 86th ed. (CRC Press: Boca Raton, FL, 2005-2006). c The specific heat value of normal water was used in this calculation. d This value was taken from http://physchem.kfunigraz.ac.at/sm/Service/Water/D2Ovisc.htm#lit1.
the modified heme, it was found that the rate of kinetic energy relaxation decreased more than 50% for aqueous and glycerol solutions (approximately 63% for aqueous solution and 55% for both glycerol solutions); the observed decrease was 31% and 15% for heavy water and nonpolar solvent solutions, respectively. The solvent dependent relaxation rate of the excited heme was observed for both the native myoglobin and the modified myoglobin systems. The kinetic energy relaxation of the heme in glycerol solutions was observed to be faster than in aqueous solutions (by 0.4 and 1.0 ps, respectively). The native heme relaxed slower in nonpolar solvent than in aqueous solution or glycerol, with a time scale similar to that of the heavy water solution. For modified myoglobin, the kinetic energy relaxation time of the heme in the nonpolar solvent is longer than in glycerol but shorter than in water. In water solutions, the heme cooling in the native myoglobin slows by 1.4 ps upon deuteration but no change was observed for the modified myoglobin. In contrast, the heme cooling is more rapid in the deuterated glycerol solutions than in corresponding normal glycerol solutions by 0.4 ps for the native heme and 0.6 ps for the modified heme. Glycerol is known to stabilize the structure of proteins because of its high viscosity58(see Table 2). In the present study, the heme was observed to be more deeply seated in the heme pocket in glycerol solution than in water solution, with the protein in glycerol solution remaining closer to the crystal structure. A larger fluctuation of CR was observed for myoglobin in aqueous solution than in glycerol solution (see Figure 4). Hydrogen bonds (HBs) between the O atoms of one heme isopropionate side chain and H atoms of Arg45 in its crystal structure are preserved in addition to the HBs formed between isopropionate groups and solvent in the glycerol solution. In
Figure 4. Fluctuation of myoglobin CR atoms in aqueous solution and in normal glycerol solution. The results were averaged over ten 30 ps trajectories. The residues were labeled by the corresponding helices. A larger fluctuation was observed for Mb in aqueous solution than in glycerol solution especially for the residues between helices or at terminals.
water solution, the heme isopropionate side chains extend into solvent and form hydrogen bonds with solvent molecules. When the protein is immersed in the nonpolar solvent, the negatively charged isopropionate side chains turn inward and form hydrogen bonds with the charged or polar side chains of protein residues Arg45, Lys63, Ser92, and Lys96. The structural alteration is shown in Figure 5. The glycerol-induced alteration of the heme pocket structure relative to that observed in pure aqueous solution was reported by Sage and co-workers on the basis of resonance Raman spectroscopy studies.59 In myoglobin, the heme is mainly buried in the protein matrix with roughly 90 van der Waals interactions with heme pocket residues and a covalent bond between its iron atom and His93 imidazole.6 The two charged isopropionate side chains of the heme are accessible to solvent. It appears that the observed solvent dependent energy relaxation rate of the excited heme results from different heme-solvent interactions in different solutions and/or from the different heme-protein interactions induced by solvents. Further analysis will demonstrate that the heme-solvent interactions are determinant. 3.2. Energy Redistribution in Myoglobin at 10 K Following Ligand Photolysis. To test whether protein structure alteration is responsible for the difference in heme cooling rate, we analyzed the excess kinetic energy redistribution in the protein. For normal water solutions and normal glycerol solutions, the 10 configurations prepared at 300 K were quenched and re-equilibrated at 10 K. The structural differences observed in differing solvation environments at 300 K were retained. Following the procedure described above, heme cooling following ligand photolysis was simulated. Coordinates and velocities were saved every 10 fs. The excited heme was observed to relax within 10 ps in all cases (see Table 3). The kinetic energy increase of the whole protein, residues 36-57 (helices C to D) and residues 58-95 (helices E to F) as a function of time for the first 5 ps is shown in Figure 6. We chose these residues because they are close to the heme (see Figure 1). In addition, helices E and F were suggested to be important for the protein function.43,60,61 Kinetic energy transfer from the heme to the surrounding residues was found to occur through three channels: “through projectile” transfer by collision between freed ligand (CO here) and heme pocket residues, “through bond” transfer through covalent link between the heme and His93, and “through space”
Energy Relaxation of the Myoglobin Heme Moiety
J. Phys. Chem. B, Vol. 111, No. 12, 2007 3247
Figure 5. Hydrogen bonds formed by heme isopropionate side chains in different solvent environments. (a) The two isopropionate side chains of heme extend into solvent in water solution, and hydrogen bonds are formed between heme and solvent molecules; (b) In glycerol solution, one isopropionate side chain forms HBs with solvent molecules and the other forms HBs with Arg45. (c) In nonpolar atomic solvent the two isopropionate side chains turn inward and form hydrogen bonds with charged or polar side chains of protein residues Arg45, Lys63, Ser92 and Lys96. The heme is depicted using a CPK model.
TABLE 3: Excess Kinetic Energy Relaxation Time Constants in Heme Following Ligand Photolysis Simulated at 10 Ka time constant (ps) solvent
native heme Mb
modified heme Mb
normal water glycerol-h8
7.5 ( 0.3 7.4 ( 0.3
9.3 ( 0.4 10.2 ( 0.4
a The simulation data are well fitted by a single-exponential function in each case.
Figure 6. Excess kinetic energy redistribution in protein for the first 5 ps following ligand photolysis from the heme at 10 K: (a) native myoglobin aqueous solution, (b) modified myoglobin aqueous solution; (c) native myoglobin glycerol solution; (d) modified myoglobin glycerol solution. The results were averaged over ten trajectories. The kinetic energy evolution of the protein is similar in all cases. This observation indicates that the solvent dependent heme cooling rate is not the result of different heme pocket structure and subsequent differences in hemeprotein interaction.
transfer mediated by nonbonded collision between the heme and its surroundings.23,25,26 For native myoglobin in normal water (Figure 6a), at least three phases of protein kinetic energy evolution can be recognized. (i) The nearly instantaneous increase of 13 kcal/ mol followed by rapid decay of ∼4 kcal/mol within 300 fs. During the first 200 fs, the time dependence of the whole protein kinetic energy is determined by the dynamics of residues 5895, where the process is dominated by initial “through projectile” energy transfer to the heme pocket residues His64 and Val68 and rapid localized structural change induced by heme doming (mainly in His93 and its immediate neighbors). For the time period 200-300 fs, there is a significant contribution from residues 36-57 leading to a roughly 2 kcal/mol increase in kinetic energy. The solvent molecules near the heme isopropi-
onate side chains are observed to be heated almost instantaneously due to the “protein quake” following ligand photolysis. Subsequent energy transfer from the heated solvent to residues 36-57 was observed. The other source of heating of residues 36-57 is “through projectile” transfer to residues including Phe43. This channel was found to be less important as rapid relaxation expected following “through projectile” energy transfer was not observed here. The excess kinetic energy of His93 is converted into potential energy and transferred to the nearby solvent molecules (one water molecule was found in the proximal heme pocket for most trajectories), causing a rapid decay. (ii) Following the first process, there is a slow increase of ∼4 kcal/mol between 0.3and ∼1.0 ps, during which the whole protein kinetic energy increase is faster than that of residues 58-95. In addition to the contribution from warming of the C and D helices, the “through space” channel of interactions between the heme side chains and nearby heme pocket residues results in an energy increase in the E and F helices (about 2 kcal/mol) and other residues (about 2 kcal/mol). (iii) For the remainder of the 30 ps trajectory, multiple component relaxation was observed (not shown here). The energy difference between residues 58-95 and the whole protein is nearly constant during this period, indicating that relaxation completely dominated by the energy transfer from E and F helices to the surrounding solvent. Several distinctions were observed for kinetic energy evolution in other solutions. For native myoglobin in glycerol solution, residues 36-57 receive more excess energy (∼1.5 kcal/mol) due to the hydrogen bonds between the heme and Arg45; the energy increase in the whole protein (∼19 kcal/mol) is much higher than that in residues 58-95 (∼11 kcal/mol), which is different from other solution systems. This is due to the “through projectile” energy transfer to Leu29 (of helix B), in addition to His64 and Val68, that occurs independent of solvent. It also can be observed that the protein relaxation is slower in glycerol solutions than in water solutions, which indicates that the kinetic energy transfer between protein and solvent is less efficient between protein and glycerol than that between protein and water. We have assumed that the energy transfer among protein residues is slow. This assumption seems reasonable because, for example, no indication of energy flow from helices E and F to surrounding residues was observed (see Figure 6). This is consistent with previous theoretical studies.26,62 It should be noted that the “through projectile” energy transfer and rapid “through bond” transfer due to heme doming must be considered separately from heme vibrational energy relax-
3248 J. Phys. Chem. B, Vol. 111, No. 12, 2007 ation.25 Although heme pocket structural change is induced by the solvent, the excess vibrational energy transferred from the excited heme into protein residues is the same (∼3 kcal/mol) in all cases and is small relative to the total excess energy relaxation of the heme (>50 kcal/mol). This agrees with the results of previous molecular dynamics studies at 300 K.23,25,26 The above kinetic redistribution analysis was carried out at 10 K, whereas real “biological” process and most experimental measurements occur at room temperature. The direct simulation and comparison of this anisotropic property at high-temperature would be interesting. However, the large thermal fluctuations at room temperature make such observations almost impossible. On the other hand, we note that heme cooling is complete within 10 ps even at this low temperature. In contrast, relaxation of the protein structure (other than heme doming, which was observed in our 10 K simulations) following photolysis happens on a relatively long time scale even at room temperature.11,19,63 Low-temperature simulations provide at least qualitative insight in the relaxation processes occurring at high temperature. And we can conclude that (a) the effect of heme pocket structural change on heme cooling is minimal and that (b) the solvent dependent kinetic energy relaxation rate observed in this work (see Table 1) is the result of the differing heme-solvent interactions in different solvation environments and not due to heme-protein interactions. 3.3. Solvent Dependent Heme Cooling Rate is the Result of the Microscopic Properties of Heme-solvent Interaction. The solvent dependence of vibrational energy relaxation in trans-stilbene was studied previously by monitoring the CdC stretching Raman peak position.64 For the solvents studied, a linear dependence on the thermal diffusivity of the solvent was observed for the vibrational cooling rate. The authors explained this correlation by assuming a two-step energy transfer process. During the first step, the excess energy of the excited mode is rapidly shared with the closest solvent molecules, and in the second step the energy slowly flows from this first solvation shell into the solvent bulk. The observed linear correlation is attributed to this slow process. A similar two-step energy transfer model was also used by Terazima and co-workers to explain the two-phase relaxation process of hydroxybenzophenone (HBP) in hexane solution.65 By assuming a spherical solute in continuous medium, Li and Champion derived a one-boundary model of vibrational energy relaxation,66 in which the thermal temperature of a solute is a function of the surface conductivity between the solute and solvent and the thermal diffusivity of the solvent. This model was used later to explain the solvent dependent behavior of iron(III) meso-tetrakis(4-sulfonatophenyl)porphyrin (FeTSPP),35 for which a single-exponential decay was observed in aqueous solution and a double-exponential decay in methanol/benzene mixtures. The first decay constant in these mixtures was found to be consistent with that in aqueous solution and the second depends on the ratio of the two components. To further explore solvent dependence of the rate of VER, the inverse kinetic energy fluctuation metric ΩKE(0)/ΩKE(t) was computed as a function of time for the solvents studied in present work. The kinetic energy diffusion constants derived are summarized in Table 2. Water demonstrates higher friction than glycerol due to quicker kinetic energy diffusion in bulk water relative to bulk glycerol. This trend agrees qualitatively with observed variation in experimental thermal diffusivity constants. However, our computed rates of kinetic energy relaxation of the heme are faster in glycerol solutions than in water. This
Zhang et al.
Figure 7. Density of states from INM calculations for various solvents. The data have been smoothed for clarity and normalized to have unit area. Observed differences in the density of states between 700 and 1100 cm-1 were found to be important for the V-V energy transfer from the excited heme to solvent.
observation cannot be explained using models based on macroscopic properties of the solvent bulk. To understand the microscopic mechanism of the heme cooling following ligand photolysis, three intermolecular energy transfer mechanisms34,67 were considered here: (i) energy transfer via hydrogen bonds interaction, (ii) direct vibrationvibration (V-V) energy transfer via resonant interaction, and (iii) energy transfer via vibration-translation (V-T), or vibration-rotation (V-R) interaction, or in other words, via thermal collisions. 3.3.1. SolVent Hydrogen Bonding to Heme Side Chains Key to the Relaxation in Excited Heme. The two isopropionate side chains of the native heme are negatively charged in our simulation model. In water and glycerol solutions, hydrogen bonds form with solvent molecules (see Figure 5). However, for the modified heme, the two isopropionate side chains were replaced by aliphatic hydrogen atoms, which make the hydrogen bonds difficult to form and weak. The slower relaxation of kinetic energy for the modified heme cases in these solvents and for native myoglobin in nonpolar solvent can be explained at least partly by this difference. In the native myoglobin glycerol solution, the kinetic energy transfer from the excited heme to protein through the hydrogen bonds between a heme isopropionate side chain and Arg45 was found to be minimaly important (only ∼2 kcal/mol at 10 K). If the energy transferred through these hydrogen bonds were large, it will not accelerate the heme cooling due to a corresponding reduction in hydrogen bonds between the heme and solvent molecules. The situation is quite different in the nonpolar solvent solution. No hydrogen bonds can be formed between the heme and solvent molecules. HBs between the heme and protein residues open an additional heme cooling channel. This additional channel effectively accelerates the heme cooling. 3.3.2. SolVent Deuteration Exposes V-V Energy Transfer Pathway for Heme Cooling. Using our classical models, the equilibrium properties of normal water and normal glycerol solutions and the corresponding deuterated solutions must be similar. The distribution of vibrational frequencies will differ. The instantaneous normal modes were computed for these solvents, and the densities of states (using eq 2) are shown in Figure 7. In water solutions, the native heme myoglobin relaxation time constants show a slowdown of 1.4 ps upon deuteration. This indicates that the V-V pathway is involved in energy transfer
Energy Relaxation of the Myoglobin Heme Moiety
J. Phys. Chem. B, Vol. 111, No. 12, 2007 3249 slowest in nonpolar solvent. For modified heme, the closer heme-solvent contact leads to more effective heme cooling in nonpolar solvent than in normal water or heavy water. The heme has faster relaxation rates in glycerol than nonpolar solvent due to the predominance of the V-V mechanism. 4. Summary and Conclusions
Figure 8. Radial distribution functions showing the structure of solvent heavy atoms around the heme. The results were averaged over all heme atoms.
between the heme and normal water. A significant number of modes observed between 700 and 1100 cm-1 for normal water is absent for deuterated water. We conjecture that these mediumfrequency modes are important to the mechanism of heme cooling. For the modified heme cases, no deuteration dependence was observed. This behavior is similar to that of nonpolar metalloporphyrin NiOEP for which the V-V energy transfer mechanism was excluded.34 For glycerol solutions, a dependence of heme cooling rate on deuteration was observed for both native and modified heme. This suggests that the V-V energy transfer mechanism contributes in both cases. This is supported by the popularity of medium modes in both solvents. The greater number of the medium-frequency modes in deuterated glycerol accelerates the rate of heme cooling relative to normal glycerol. Compared to that in water solutions, the stronger V-V interaction in native myoglobin glycerol solution and additional V-V interaction in modified myoglobin results in faster heme kinetic energy relaxation in glycerol solutions by 0.4 ps for native myoglobin and 1.0 ps for modified myoglobin. 3.3.3. V-T and V-R Energy Transfer Dominates the Heme Cooling Process in the Nonpolar SolVent Solutions. The vibration-translation and vibration-rotation energy transfer, or in other words, thermal collisions, contribute to the energy transfer in all cases considered in this work. The efficiency of the thermal collision mechanism is proportional to the collision frequency, which in turn is proportional to the local density of solvent particles. The radial distribution function was computed for the heme and solvent heavy atoms (see Figure 8). The distributions of nonpolar atoms around the native heme and the modified heme are similar, especially for short distances. The observed differences in the rate of heme cooling in these two systems must be due to the nature of hydrogen bonding between the isopropionate side chains and protein residues in the native heme, as suggested in the previous section. An alternative mechanism for intermolecular energy transfer in nonpolar solvent solutions is V-T and V-R energy transfer as has been demonstrated by Kitagawa34 and Kaiser.67 Because of the smaller molecular radius, the nonpolar solvent has closer contact with the heme than other solvents examined in this study in both native and modified myoglobin. However, for the native heme, this contact factor is secondary to the predominant channels of V-V energy transfer and HB interaction in normal water solution and glycerol solutions. That factor appears to be balanced by the HB interaction in heavy water solution. As a result, native heme cooling was observed to be
This work has investigated the time scales and pathways of vibrational energy relaxation in both wild type and modified heme myoglobin following ligand photolysis. Five different solvent models, including normal water, heavy water, normal glycerol, deuterated glycerol, and a nonpolar solvent, were used in the simulation. A solvent dependent heme cooling rate was observed. A single-exponential decay was found to be characteristic of the excess kinetic energy of the heme in all cases studied. Relaxation of the heme was observed to be complete within 4.3-8.5 ps after photolysis of CO. It was observed to be faster in the native myoglobin than in modified myoglobin in all solvents, consistent with previous experimental and theoretical results.15-17,20,23,25 By analyzing the excess kinetic energy redistribution at 10 K, we found the direct heme-solvent energy transfer to be the most important heme cooling pathway following ligand photolysis. However, this dependence cannot be explained using bulk transport properties of the solvent, including macroscopic thermal diffusion. Three intermolecular energy transfer mechanisms were considered: (i) energy transfer mediated by hydrogen bonds, (ii) direct vibration-vibration energy transfer via resonant interaction, and (iii) energy transfer via vibrationtranslation or vibration-rotation interaction. The dominant microscopic intermolecular interaction mechanism is different upon the modification of the heme side chains and solvent environment. The vibration-vibration interaction and hydrogen bond interaction between the heme and solvent molecules dominates the energy transfer for native myoglobin in aqueous solution and glycerol solutions, resulting in faster heme cooling than in nonpolar solution. For modified myoglobin, the V-V interaction is observed to be effective in glycerol solution. However, the V-V energy transfer mechanism is absent or less effective in aqueous solution. This explains the faster relaxation observed in glycerol solution compared with that in aqueous solution. Thermal collision is the only available energy transfer pathway for modified myoglobin in water solution or for modified myoglobin in a nonpolar environment. The closer contact of nonpolar solvent with the heme results in faster heme kinetic energy relaxation of modified myoglobin in nonpolar solutions than in water. Hydrogen bonds between heme isopropionate side chains and protein residues Arg45, Lys63, Ser92, and Lys96 in native myoglobin in nonpolar solvent serve as an additional kinetic energy dissipation pathway, resulting in faster heme cooling than observed for modified myoglobin in nonpolar solvent. Solvent-induced protein structure alteration was observed in this study. Due to the high viscosity of the glycerol solvent, the initial crystal structure of myoglobin was preserved on the time scale of our simulation. In water solution, the protein was observed to be more flexible; the heme isopropionate side chains were observed to be completely extended into the solvent. The simulated heme pocket structure was observed to be different in these two solvent solutions, consistent with previous experi
3250 J. Phys. Chem. B, Vol. 111, No. 12, 2007 mental results derived from resonance Raman spectroscopy studies.59 Due to the different heme pocket structures, the “through projectile” energy transfer between CO and protein residues following ligand photolysis has different targets in the two solutions. However, the excess vibrational energy dissipated into protein residues “through space” or “through bonds” from the excited heme is the same in the two cases. Acknowledgment. We are grateful for the generous support of this research by the National Science Foundation (Grant No. CHE-0316551) and Boston University’s Center for Computer Science. We thank Prof. D. Tobias, University of California at Irvine, for providing the glycerol parameters. Note Added in Proof. During the proof process, we were informed by Professors T. Kitagawa and Y. Mizutani that their recent experimental results (refs 68 and 69) directly confirm our conclusion that the interactions between the heme isopropionate side chains and the surrounding solvent molecules form the dominant pathway for dissipation of the heme’s excess kinetic energy. References and Notes (1) Asplund, M. C.; Zanni, M. T.; Hochstrasser, R. M. Proc. Natl. Acad. Sci. U.S.A. 2000, 97, 8219-8224. (2) Kholodenko, Y.; Volk, M.; Gooding, E.; Hochstrasser, R. M. Chem. Phys. 2000, 259, 71-87. (3) Muenck, E.; Champion, P. M. Ann. N. Y. Acad. Sci. 1975, 244, 142-162. (4) Sage, J. T.; Paxson, C.; Wyllie, G. R. A.; Sturhahn, W.; Durbin, S. M.; Champion, P. M.; Alp, E. E.; Scheidt, W. R. J. Phys.: Condens. Matter 2001, 13, 7707-7722. (5) Sage, J. T.; Durbin, S. M.; Sturhahn, W.; Wharton, D. C.; Champion, P. M.; Hession, P.; Sutter, J.; Alp, E. E. Phys. ReV. Lett. 2001, 86, 4966-4969. (6) Henry, E. R.; Eaton, W. A.; Hochstrasser, R. M. Proc. Natl. Acad. Sci. U.S.A. 1986, 83, 8982-8986. (7) Elber, R.; Karplus, M. Science 1987, 235, 318-321. (8) Elber, R.; Karplus, M. J. Am. Chem. Soc. 1990, 112, 9161-9175. (9) Miller, R. J. Annu. ReV. Phys. Chem. 1991, 42, 581-614. (10) Straub, J. E.; Karplus, M. Chem. Phys. 1991, 158, 221-248. (11) Li, H.; Elber, R.; Straub, J. E. J. Biol. Chem. 1993, 268, 1790817916. (12) Mizutani, Y.; Kitagawa, T. Science 1997, 278, 443-446. (13) Mizutani, Y.; Kitagawa, T. Chem. Rec. 2001, 1, 258-275. (14) Okazaki, I.; Hara, Y.; Nagaoka, M. Chem. Phys. Lett. 2001, 337, 151-157. (15) Petrich, J. W.; Poyart, C.; Martin, J. L. Biochemistry 1988, 27, 4049-4060. (16) Li, P.; Sage, J. T.; Champion, P. M. J. Chem. Phys. 1992, 97, 32143227. (17) Lim, M.; Jackson, T. A.; Anfinrud, P. A. J. Phys. Chem. 1996, 100, 12043-12051. (18) Kitagawa, T.; Haruta, N.; Mizutani, Y. Biopolymers 2002, 61, 207213. (19) Kruglik, S. G.; Mojzes, P.; Mizutani, Y.; Kitagawa, T.; Turpin, P.-Y. J. Phys. Chem. B 2001, 105, 5018-5031. (20) Lian, T.; Locke, B.; Kholodenko, Y.; Hochstrasser, R. M. J. Phys. Chem. 1994, 98, 11648-11656. (21) Genberg, L.; Heisel, F.; Mclendon, G.; Miller, R. D. J. Phys. Chem. 1987, 91, 5521-5524. (22) Henry, E. R.; Hochstrasser, R. M. Proc. Natl. Acad. Sci. U.S.A. 1987, 84, 6142-6146. (23) Bu, L.; Straub, J. E. J. Phys. Chem. B 2003, 107, 10634-10639. (24) Bu, L.; Straub, J. E. J. Phys. Chem. B 2003, 107, 12339-12345. (25) Sagnella, D. E.; Straub, J. E. J. Phys. Chem. B 2001, 105, 70577063. (26) Sagnella, D. E.; Straub, J. E.; Thirumalai, D. J. Chem. Phys. 2000, 113, 7702-7711. (27) Ma, J.; Huo, S.; Straub, J. E. J. Am. Chem. Soc. 1997, 119, 25412551.
Zhang et al. (28) Cao, W.; Ye, X.; Sjodin, T.; Christian, J. F.; Demidov, A. A.; Berezhna, S.; Wang, W.; Barrick, D.; Sage, J. T.; Champion, P. M. Biochemistry 2004, 43, 11109-11117. (29) Rodriguez, J.; Kirmaier, C.; Holten, D. J. Chem. Phys. 1991, 94, 6020-6029. (30) Rodriguez, J.; Holten, D. J. Chem. Phys. 1989, 91, 3525-3531. (31) Mizutani, Y.; Kitagawa, T. Bull. Chem. Soc. Jpn. 2002, 75, 965971. (32) Mizutani, Y.; Uesugi, Y.; Kitagawa, T. J. Chem. Phys. 1999, 111, 8950-8962. (33) Kruglik, S. G.; Mizutani, Y.; Kitagawa, T. Chem. Phys. Lett. 1997, 266, 283-289. (34) Mizutani, Y.; Kitagawa, T. Bull. Chem. Soc. Jpn. 2002, 75, 623639. (35) Mizutani, Y.; Kitagawa, T. J. Mol. Liq. 2001, 90, 233-242. (36) Kuriyan, J.; Wilz, S.; Karplus, M.; Petsko, G. A. J. Mol. Biol. 1986, 192, 133-154. (37) Brooks, B. R.; Bruccoleri, R. E.; Olafson, B. D.; States, D. J.; Swaminathan, S.; Karplus, M. J. Comput. Chem. 1983, 4, 187-217. (38) MacKerell, A. D. J.; Bashford, D.; Bellott, M.; Dunbrack, R. L.; Evanseck, J. D.; Field, M. J.; Fischer, S.; Gao, J.; Guo, H.; Ha, S.; JosephMcCarthy, D.; Kuchnir, L.; Kuczera, K.; Lau, F. T. K.; Mattos, C.; Michnick, S.; Ngo, T.; Nguyen, D. T.; Prodhom, B.; Reiher, W. E. I.; Roux, B.; Schlenkrich, M.; Smith, J. C.; Stote, R.; Straub, J.; Watanabe, M.; Wiorkiewicz-Kuczera, J.; Yin, D.; Karplus, M. J. Phys. Chem. B 1998, 102, 3586-3616. (39) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. J. Chem. Phys. 1983, 79, 926-935. (40) Tobias, D. Private communication. (41) van Gunsteren, W. F.; Karplus, M. Biochemistry 1982, 21, 22592274. (42) Verlet, L. Phys. ReV. 1967, 1, 98-103. (43) Kachalova, G. S.; Popov, A. N.; Bartunik, H. D. Science 1999, 284, 473-476. (44) Schlichting, L.; Berendzen, J.; Phillips, G. N., Jr.; Sweet, R. M. Nature 1994, 371, 808-812. (45) Meuwly, M.; Becker, O. M.; Stote, R.; Karplus, M. Biophys. Chem. 2002, 98, 183-207. (46) Thirumalai, D.; Mountain, R. D. Phys. ReV. A 1990, 42, 45744587. (47) Straub, J. E.; Rashkin, A. B.; Thirumalai, D. J. Am. Chem. Soc. 1994, 116, 2049-2063. (48) Straub, J. E.; Thirumalai, D. Proc. Natl. Acad. Sci. U.S.A. 1993, 90, 809-813. (49) Sagnella, D. E.; Straub, J. E. Biophys. J. 1999, 77, 70-84. (50) Bu, L.; Straub, J. E. Biophys. J. 2003, 85, 1429-1439. (51) Sagnella, D. E.; Straub, J. E.; Jackson, T. A.; Lim, M.; Anfinrud, P. A. Proc. Natl. Acad. Sci. U.S.A. 1999, 96, 14324-14329. (52) Goodyear, G.; Stratt, R. M. J. Chem. Phys. 1997, 107, 3098-3120. (53) Goodyear, G.; Larsen, R. E.; Stratt, R. M. Phys. ReV. Lett. 1996, 76, 243-246. (54) Ladanyi, B. M.; Stratt, R. M. J. Phys. Chem. A 1998, 102, 10681082. (55) Seeley, G.; Keyes, T. J. Chem. Phys. 1989, 91, 5581-5586. (56) Straub, J. E.; Choi, J.-K. J. Phys. Chem. 1994, 98, 1097810987. (57) Wan, Y.; Stratt, R. M. J. Chem. Phys. 1994, 100, 5123-5138. (58) Kleinert, T.; Doster, W.; Leyser, H.; Petry, W.; Schwarz, V.; Settles, M. Biochemistry 1998, 37, 717-733. (59) Sage, J. T.; Schomacker, K. T.; Champion, P. M. J. Phys. Chem. 1995, 99, 3394-3405. (60) Rodgers, K. R.; Spiro, T. G. Science 1994, 265, 1697-1699. (61) Sato, A.; Mizutani, Y. Biochemistry 2005, 44, 14709-14714. (62) Yu, X.; Leitner, D. M. J. Phys. Chem. B 2003, 107, 16981707. (63) Findsen, E. W.; Scott, T. W.; Chance, M. R.; Friedman, J. M.; Ondrias, M. R. J. Am. Chem. Soc. 1985, 107, 3355-3357. (64) Iwata, K.; Hamaguchi, H.-o. J. Phys. Chem. A 1997, 101, 632637. (65) Okazaki, T.; Hirota, N.; Terazima, M. J. Mol. Liq. 2001, 90, 243249. (66) Li, P.; Champion, P. M. Biophys. J. 1994, 66, 430-436. (67) Scherer, P. O. J.; Seilmeier, A.; Kaiser, W. J. Chem. Phys. 1985, 83, 3948-3957. (68) Gao, Y.; Koyama, M.; El-Mashtoly, S. F.; Hayashi, T.; Harada, K.; Mizutani, Y.; Kitagawa, T. Chem. Phys. Lett. 2006, 429, 239-243. (69) Koyama, M.; Neya, S.; Mizutani, Y. Chem. Phys. Lett. 2006, 430, 404-408.
