Article pubs.acs.org/JPCA

Communication Maps of Vibrational Energy Transport Through Photoactive Yellow Protein Yao Xu and David M. Leitner* Department of Chemistry and Chemical Physics Program, University of Nevada, Reno, Nevada 89557, United States ABSTRACT: We calculate communication maps for Photoactive Yellow Protein (PYP) from the purple phototropic eubacterium Halorhodospira halophile and use them to elucidate energy transfer pathways from the chromophore through the rest of the protein in the ground and excited state. The calculations reveal that in PYP excess energy from the chromophore flows mainly to regions of the surrounding residues that hydrogen bond to the chromophore. In addition, quantum mechanics/molecular mechanics and molecular dynamics (MD) simulations of the dielectric response of the protein and solvent environment due to charge rearrangement on the chromophore following photoexcitation are also presented, with both approaches yielding similar time constants for the response. Results of MD simulations indicate that the residues hydrogen bonding to the chromophore make the largest contribution to the response.

Frequency-resolved communication maps23 identify the network of energy transport pathways in a protein or other macromolecule, one of a number of methods developed in recent years to locate energy transport channels in biomolecules.11,12,22,24−31 Frequency-resolved communication maps, like some of these methods,29,30 provide information about how specific vibrational modes contribute to pathways for energy transport. It is a coarse graining procedure that in the limit of the size of the biomolecule itself goes over to the frequency resolved energy diffusion coefficients from which thermal transport coefficients can be determined. Because most of the methods have only been tested on a few proteins it is useful to compare results of different approaches where possible, which we do in this study by comparing communication maps for PYP with calculations of energy conductivity in the same protein by Ishikura and Yamato.22 We thermally average the frequency resolved maps to compute temperature-dependent communication maps, and sample the communication maps to find pathways for energy transport beginning at any residue of PYP, though we focus mainly on transport channels beginning at the deprotonated p-coumeric acid chromophore. We carry out the analysis in both the ground and excited state of PYP to compare the energy transport channels in each state. We find that in both states the network of hydrogen bonds beginning with the residues hydrogen bonding to the chromophore mediates most of the transport of excess energy from the chromophore across the protein. We find little qualitative difference in the energy transport channels for PYP in the ground and excited states.

1. INTRODUCTION Hydrogen bonds, beyond their role in protein and peptide structure and dynamics, mediate energy flow in peptides and proteins, between biomolecules and water, and through networks of water molecules.1−4 In peptide−water clusters, for example, the kinetics of hydrogen bond rearrangement, i.e., the shuttling of a water molecule from one hydrogen-bonding site on a peptide to another,5 is influenced by energy flow between vibrations of the peptide and hydrogen bonds between the peptide and water.6 Hydrogen bonds mediate thermal transport between organic monolayers and water.7,8 In proteins hydrogen bonds serve to transfer energy efficiently, perhaps best illustrated by the exceptionally high thermal conductivity of spider silk due to the hydrogen bond network supporting the fibers.9 Hydrogen bonds funnel excess energy from protein cofactors, e.g., excess energy following ligand photolysis at heme groups,10−16 energy that can flow along an allosteric pathway between binding sites, e.g., in the dimeric Scapharca hemoglobin.17−20 Energy transport through hydrogen bonds in a PDZ domain and the surrounding water mediate allosteric transitions from the peptide−ligand binding site.21 Calculations of thermal energy conductivity in the ground state of photoactive yellow protein (PYP) reveal the central role that hydrogen bonds play in carrying energy from the deprotonated p-coumeric acid chromophore anisotropically through the protein.22 In this article we present calculations of communication maps in the ground and excited state of the protein to explore the role of hydrogen bond networks in energy transport from the chromophore in PYP. The force field model adopted for the calculation is also used to simulate the dielectric response to photoexcitation of the chromophore and is compared with results of a quantum mechanics/molecular mechanics (QM/MM) calculation of the dielectric response, which yields similar time constants. © 2014 American Chemical Society

Special Issue: Kenneth D. Jordan Festschrift Received: November 16, 2013 Revised: January 16, 2014 Published: February 19, 2014 7280

dx.doi.org/10.1021/jp411281y | J. Phys. Chem. A 2014, 118, 7280−7287

The Journal of Physical Chemistry A

Article

present and discuss our results for energy flow in PYP, and the dielectric response. Concluding remarks are given in section 4.

For the calculation of the communication maps we adopt a force field model that we recently used to calculate the dielectric response to photoexcitation of the chromophore.32 The major focus of this article is the calculation of communication maps of PYP, and as one check on this force field model we compare results of the dielectric response with those we obtain from a semiempirical QM/MM simulation. We find that both models yield similar time constants. Photoactive yellow protein, discovered in the bacterium Halorhodospira halophila, is a water-soluble photoreceptor protein containing 125 amino acids that absorbs blue light, initiating a signaling state.33 The PYP chromophore, a pcoumaric acid anion (pCA), is covalently bonded on one end to a cysteine residue and hydrogen bonds at the other end, the O− phenyl group, to two residues (Figure 1). The chromophore

2. THEORETICAL AND COMPUTATIONAL METHODS A. Frequency-Resolved Communication Maps. To obtain structures for the calculation of communication maps, we begin with the PYP (1OTB) structure (Figure 1) obtained from the Protein Data Bank (PDB). This structure was introduced into a 50 Å3 cubic box of equilibrated TIP3P water model35 and simulated using the AMBER force field36 in the GROMACS program.37 The charges of the chromophore in both ground (S0) and vertical excited state (S1) were obtained from an ab initio CASSCF (10,12)/6-31G** calculation with Gaussian03. 38 Force field parameters adopted for the chromophore in the ground and excited states are provided in ref 32. The particle-mesh Ewald (PME) method39 was employed to account for long-range electrostatic effects and a constant pressure of 1.0 atm was maintained. The SHAKE algorithm40 was employed to constrain all bonds. Periodic boundary conditions were used with standard cutoffs. The excess potential energy due to bad contacts and strain was then reduced using the adopted basis Newton−Raphson (ABNR) minimization method for about 1000 steps so that the norm of the gradient was 5.0 kcal/mol. The molecular dynamics simulation was carried out using the Verlet algorithm. The temperature was increased slowly to 300 K in intervals of 10 K for about 6 ps. Equilibration dynamics followed for about 100 ps by keeping the temperature constant at 300 K with a time step of 1 fs. Of course, PYP undergoes photoisomerization much faster than the equilibration time in the excited state, and in this respect we have more plausibly modeled the dielectric response of PYP analogs where the chromophore does not undergo photoisomerization.41 We calculate communication maps in the excited state for a structure close to the ground state structure to examine if there are significant differences from the ground state results. The system was minimized using the ABNR method for about 50 000 steps to satisfy the norm of the gradient of 1 × 10−6 kcal/mol. Normal mode analysis was then carried out on the systems in both the ground state and the excited state. Using the normal modes for the protein, we computed frequency resolved communication maps. Details concerning the calculation of communication maps have been presented elsewhere and in this subsection we summarize the method. The approach begins with the heat current operator in harmonic approximation4243

Figure 1. Ribbon diagram of PYP showing the chromophore, which is attached to Cys69, and nearby residues Glu46, Tyr42, to which it hydrogen bonds, and Thr50, which also appears close the chromophore anion during MD and QM/MM simulation. Arg52, Thr70, Ser72, Tyr94, and Tyr98 also lie near the chromophore.

undergoes trans → cis photoisomerization within a few picoseconds upon photoexcitation, modifying hydrogen bonding between the chromophore and nearby residues and triggering the PYP photocycle.34 To follow the flow of energy in PYP from the chromophore, Ishikura and Yamato22 derived the interresidue energy conductivity in terms of timecorrelation functions of the interatomic energy flux. They then sampled molecular dynamics (MD) simulation trajectories of the PYP molecule in its aqueous solution environment and evaluated numerically the energy conductivities. As a result of this analysis, they identified several pathways for energy transport in the molecule.22,31 The primary energy transport pathway they found from the chromophore in the ground state extends via a hydrogen bond network out to the N-terminal cap. In the following section we summarize our calculation of frequency-dependent communication maps. We then discuss computation of the dielectric response of the protein and solvent environment to photoexcitation. In section 3 we

S=

∑ Sαβ aα†aβ (1)

α ,β

where α and β are two modes of the protein. The coefficient, Sαβ, corresponding to the protein as a whole can be expressed in terms of the Hessian matrix, H, and eigenmodes, e, of the object. We break the coefficient up into contributions from various regions. For the region, AA′, it is S{αβAA ′} =

iℏ(ωα + ωβ ) 4V ωαωβ

∑

∑

elαHrrll′′(R l − R l ′)elβ′

r , r ′∈ (x , y , z) l , l ′∈ AA ′

(2)

where Rl is the position of atom l and r is a coordinate (x, y, or z). When the region AA′ spans the protein, eq 2 expresses the matrix elements of the heat current operator for the whole system in harmonic approximation given in ref 42. Similarly, we 7281

dx.doi.org/10.1021/jp411281y | J. Phys. Chem. A 2014, 118, 7280−7287

The Journal of Physical Chemistry A

Article

⟨ΔE⟩; and the brackets denote the equilibrium ensemble average. The linear response approximation to the dynamic Stokes shift becomes more questionable on times near and beyond the reaction time of ≈1 ps. We recently examined the subpicosecond response of the protein/solvent environment to charge rearrangement of the chromophore using a molecular mechanics force field model we developed for the PYP chromophore,32 and the purpose of the analysis here is to compare that result with a QM/MM calculation. We adopt the semiempirical PM3 approach, which, though only a rough approximation, is feasible for the length of simulation required to calculate the dielectric response. Our computation of C(t) with eq 5, the dielectric response of the protein and solvent environment to charge redistribution on the chromophore, parallels earlier work by Schulten and co-workers,53 who examined the protein response to photoexcitation of bacteriorhodopsin. For the calculation of the dielectric response, MD simulations were carried out by using the GROMACS package. PYP was minimized for 3000 steps with the steepest descent algorithm using the AMBER03 force field. Then it was solvated in a cubic water box of 60 Å on each side with TIP3P water. The system, which comprised 6499 water molecules, was neutralized by adding 7 Na+ ions and 1 Cl− ion using the genion tool in GROMACS37 and minimized again using the steepest decent algorithm. The system was then equilibrated for 300 ps using the AMBER03 force field with periodic boundary conditions. In the first 100 ps the protein was restrained before it was released in the next 200 ps. The SHAKE algorithm was employed to constrain bonds with hydrogen atoms throughout the simulation. Nonbonded interactions were gradually brought to zero by using a PME switch for electrostatic interactions and a switch function for van der Waals at 10 and 12 Å cutoffs, respectively. A standard leapfrog algorithm with an integration time step of 0.5 fs was used to integrate Newton’s equations of motion. All the MD simulations were done in the microcanonical (NVE) ensemble for 8 ns at energies corresponding to a temperature near 300 K with atomic coordinates stored every 1 fs. As in the communication map calculations described above, the charges of the chromophore (deprotonated cinnamic acid) in anionic form at both ground (S0) and vertical excited state (S1) were obtained from an ab initio CASSCF (10,12)/631G** calculation with Gaussian03.38 Force field parameters adopted for the chromophore in the ground and excited states are listed in ref 32 and the force field is similar to one used by us in earlier MD simulations of PYP.54,55 The semiempirical PM3 QM/MM calculation was carried out following the methods detailed in ref 56. In the QM/MM calculation, the chromophore was carefully built up using the xLEAP module in the AMBER9 simulation package,57 and we followed a similar simulation protocol as in its MD counterpart;32 i.e., we run 8 ns simulations starting with an identical chromophore geometry using the respective charges in the ground and excited state. In both systems (chromophore in solvated PYP and in bulk water), we modeled the chromophore in the QM/MM simulations using the semiempirical PM3 method with the remaining part of the system treated classically.

can sum over all regions to obtain Sαβ for the whole molecule. Considering only the local region AA′ we write the local energy diffusivity in mode α in harmonic approximation as Dα{AA ′} =

πV 2 3ℏ2ωα 2

∑ |S{αβAA ′}|2 δ(ωα − ωβ) β≠α

(3)

When the local region AA′ spans the molecule, eq 3 gives the mode diffusivity, from which the coefficient of thermal conductivity, κ, can be expressed for the whole molecule, κ = ΣαCαDα, where Cα is the heat capacity per unit volume of the molecule for mode α. The coefficient of thermal conductivity has been calculated for proteins44−47 and here we identify the regions of the protein that primarily contribute to thermal transport through the molecule. For a practical calculation on a finite-sized system we substitute a rectangular window of width η for the delta function. The width should be large enough to envelop several vibrational modes, enough so that the mode diffusivity converges. For the results presented here, we calculated the mode diffusivity using η = 15 cm−1, which is large enough to include many modes in the averaging; results for Dα did not change significantly with larger η. Communi′} cation maps can then be constructed by plotting D{AA for all α regions AA′ at frequencies, ωα. Below we use frequency-resolved communication maps to locate energy transport channels that include a particular site on the protein, in this case the chromophore. The procedure can be carried out for any particular vibrational mode of the protein to determine how vibrational modes around a given frequency carry energy through the protein. We carried out a thermal average over the frequency-resolved communication maps for the structures at 300 K, using for the thermal averaging the communication maps constructed at 50 cm−1 and 100 cm−1 and continuing in 50 cm−1 intervals until 400 cm−1, above which the communication maps contribute very little to the ′} thermal average due to the small values of D{AA at higher α frequency. The thermally averaged communication map is obtained by assigning a Boltzmann weight to the frequency resolved communication maps. Each of these frequencyresolved maps is actually an average over maps for four modes closest in frequency to the designated frequency and an average over four protein structures. B. Solvation Dynamics. The response of water to photoexcitation of a solute immersed in it and the corresponding time-dependent Stokes shift (TDSS) of the solute have been explored with the help of MD simulations for some time.48,49 The solvation response function is usually expressed in normalized form as S( t ) =

v (t ) − v (∞ ) v(0) − v(∞)

(4)

where v(t) is the time evolution of the frequency of maximum fluorescence. The Stokes shift for PYP is 2180 cm−1.50 Within the linear response approximation, the normalized solvation response function, S(t), is equal to the equilibrium solvation time correlation function, C(t), which is49,51,52 C(t ) =

⟨δ ΔE(t ) δ ΔE(0)⟩ ⟨(δ ΔE)2 ⟩

3. RESULTS AND DISCUSSION A. Communication Maps. We have computed frequency resolved communication maps for ground and excited state PYP. The thermal average at 300 K for each state is plotted in Figure 2. These maps indicate pairs of residues between which

(5)

where ΔE(t) = Ee − Eg is the difference in the chromophore− bath interaction energy when the chromophore is in the excited state, e, and the ground state, g, at time, t; δΔE(t) = ΔE(t) − 7282

dx.doi.org/10.1021/jp411281y | J. Phys. Chem. A 2014, 118, 7280−7287

The Journal of Physical Chemistry A

Article

Table 1. Top 10 Hot Residues for PYP Based on Ranking the Diffusivity of the Energy out of the Chromophore from Communication Map Calculationa communication map

Yamato

ground

excited

ground

Glu46 Cys69 Thr50 Tyr42 Ala67 Arg52 Phe96 Gln99 Asp97 Phe62

Cys69 Glu46 Thr50 Arg52 Tyr42 Phe62 Tyr98 Ala67 Phe96 Ile49

Thr70 Tyr42 Thr50 Pro68 Glu46 Ala67 Arg52 Tyr94 Tyr98 Ser72

a

Hot residues from the calculation of ref 22 are also listed for comparison.

with those obtained for the ground state by the energy conductivity calculations reported in ref 22. We notice that, although the ordering may be different, most of the residues identified in ref 22 are also identified in our communication maps computed as described in section 2. The force field models used in both studies are the same, apart from the chromophore. Examining Table 1, we find that the main differences between results of this study and those in ref 22 concern residues of the main chain to which the chromophore is attached. For example, Thr70 is identified in ref 22 to lie along the energy transport pathway from the chromophore, whereas in this study Cys69, to which the chromophore is directly bonded, is found to be the gateway to energy transport into the main chain. This difference is likely due to the definition of an extended chromophore as both pCA and Cys69 in ref 22. Apart from this essentially technical difference, the results obtained in this study and ref 22 are quite similar. Both sets of results highlight the importance of the residues that hydrogen bond to the chromophore, which we observe here to be the main gateway for energy transport from the chromophore in both the ground and excited state. To illustrate the location of the residues into which energy mainly flows from the chromophore, we plot in Figure 3 the location of the residues in PYP that are listed in the first two columns of Table 1. We observe in Figure 3 that the neighboring residues hydrogen bonding to the chromophore, which actively transport excess energy away from the chromophore at early times, are Tyr42, Glu46, and Thr50, which initiate vibrational energy transport from the chromophore in both the ground and excited electronic state. In addition to this group of residues there is another energy transport pathway that dissipates the excess energy effectively through Arg52, Tyr98, and Gln99, none of which hydrogen bond to the chromophore. In particular, although Arg52 is not directly in contact with the chomophore (it is about 6.5−7.0 Å away), it can form a counterion to the phenolic oxygen anion of the chromophore, which, though of lower energy diffusivity compared to those residues that can hydrogen bond to the chromophore, nevertheless opens up an energy transport channel from it. Looking beyond the first residues into which energy flows from the chromophore, we examine the next group of residues into which energy flows from the first group. A schematic mapping of the transport pathways is plotted in Figure 4. This

Figure 2. Thermal average communication maps plotted for PYP at 300 K at ground (top) and excited state (bottom). Red, blue, and green symbols correspond to values of D between residues that are greater than 100 Å2 ps−1, between 10 and 100 Å2 ps−1, and between 1 and 10 Å2 ps−1, respectively. In PYP, the chromophore is represented as residue 126.

energy flow is rapid, in particular where red symbols appear. The chromophore is in position 126, in the last column of the map. The communication maps can be used to determine energy transport channels between specific sites on the protein, which are the intersection of the residue pairs between which energy flow is fast. Our focus here is identifying energy transport channels that begin at the chromophore. We examine the communication maps for PYP at 300 K and rank all the residues on the basis of their energy diffusivity value to examine the early steps of energy transfer from the chromophore. In Table 1 we rank the residues by their interaction with the chromophore, i.e., in order of the value of D between that residue and the chromophore at 300 K. We compare the rankings obtained from the communication maps 7283

dx.doi.org/10.1021/jp411281y | J. Phys. Chem. A 2014, 118, 7280−7287

The Journal of Physical Chemistry A

Article

Table 2. Top Hot Residues for PYP Based on Ranking the Diffusivity of the Energy out of the Residues That Hydrogen Bond to the Chromophore (Primary) and Nearby Residues (Secondary) from Communication Map Calculation residue

ground state

excited state Primary

Tyr42 Glu46 Thr50 Asn43 Arg52

Asn43, Gln41, Val57 Ala45, Gly47, CRO126 Ile49, Gly51, CRO126, Glu46 Secondary Tyr42, Ala44, Ala30, Leu23, Glu46, Asp24 Asp53, Gly51, Gln99, Gly47, Tyr98

Gln41, Asn43, Val57 Ala45, Gly47, CRO126 Ile49, Thr50, Gly51, Arg52 Ala44, Tyr42, Ala30, Phe28, Leu23, Asp24 Gly51, Asp53, Gln99, Thr50, Tyr98

chromophore and extends across PYP. We find that after energy flows out of the residues that hydrogen bond to the chromophore, it follows mainly along two pathways. The first one goes from Tyr42 through Asn43 and then to Leu23 and Asp24 at the N-terminal cap; the second pathway is from Thr50 through Arg52 and then to Tyr98 and Gln99. Secondary pathways are also listed in Table 2. Overall, all the residues along the pathways are the same as those found by Ishikura and Yamato,22 who calculated the interresidue energy conductivity in the ground state in terms of time-correlation functions of the interatomic energy flux. Given the consistent results for the energy transport channels from these two methods, and the fact that we observe the same energy transport channels in the ground state and excited state of PYP, we can conclude that Figure 4 describes well the long-range intramolecular signaling of PYP, which evolves from the hydrogen bonding between the chromophore and specific residues that surround it. Recently, we computed mode damping rates of the chromophore due to coupling to the protein and solvent environment by MD simulations.32 We observed that for chromophore modes up to 600 cm−1 the damping rate was controlled almost entirely by coupling to Tyr42, Glu46, and Thr50, which we find from calculation of communication maps to be the gateway residues from the chromophore into energy transport channels to remote parts of the protein. B. Solvation Dynamics. In a recent study we examined the short-time (subpicosecond) dielectric response of the protein and water environment to photoexcitation of the chromophore using force field models developed for the ground and excited states of PYP. The result was decomposed into contributions of the protein and the water, where the latter was found to contribute negligibly, as well as into contributions of individual residues. Those that hydrogen bond to the chromophore were found to respond significantly to changes in charge on the chromophore. Here, we compare the dielectric response obtained using the force field model of our previous study with the dielectric response obtained by QM/MM calculations described in section 2. We consider the response out to 100 ps, even though such times are well beyond the reaction time of PYP. In this respect we are in essence modeling the dielectric response of PYP analogs where the chromophore does not undergo photoisomerization.41 Our purpose is to examine whether the description of the dielectric response obtained with the force field model used in the previous work, and in the calculation of communication maps in the previous section, yields results that are comparable to a QM/MM simulation result, and for that purpose we make the comparison out to relatively long times.

Figure 3. Ten residues with the largest values of D to the chromophore in the ground state (top) and excited state (bottom).

Figure 4. Summary of energy transfer pathways from the chromophore (CRO, yellow) to the N-terminal cap (red). The residues that form the hydrogen bond network with the chromophore and in the helix α3 are shown in green and orange, respectively. Pathways leading into the sheets β4 and β5 are shown in blue.

figure shows the gateway that commences with the residues that hydrogen bond to the chromophore, leading to transport to more distant parts of the molecule and does not include local energy flow through the main chain, which commences at Cys69 and, to a lesser extent, Ala67. We rank the residues near Tyr42, Glu46, and Thr50 in terms of their interresidue energy diffusivity (see Table 2 for details) with these gateway residues to map an energy transport pathway that begins at the 7284

dx.doi.org/10.1021/jp411281y | J. Phys. Chem. A 2014, 118, 7280−7287

The Journal of Physical Chemistry A

Article

obtained by the QM/MM simulations. Though the time constants are similar with the two approaches to the simulations, the contribution of each decay term, aj, j = 1, 2, 3, to the overall decay differs. The MD approach yields a fast decay that weighs more heavily, 74%, to the overall decay than does the QM/MM approach (51%). The latter approach yields a more dominant picosecond time-scale decay (47% vs 21%), and both predict a small contribution for the long-time decay (5% from MD simulations and 2% from QM/MM). We recently examined the role of residues that hydrogen bond to the chromophore in the dielectric response of the protein to charge rearrangement on the chromophore when PYP is excited to the S1 state.32 We found that the response can be accounted for largely by the response of two of the residues, Tyr42 and Glu46, both of which hydrogen bond to the chromophore, and both of which serve as gateway residues for energy transport from the chromophore. Fluorescence decay measurements by Mataga and co-workers reveal coherent oscillations of about 50 and 120 cm−1 in the fluorescence decay,63,64 which they attributed to coupled vibrations of the chromophore and these residues. From all of these perspectives the network of hydrogen bonds between the chromophore and its immediate environment control the flow of excess energy from the chromophore.

The dielectric response to photoexcitation of the chromophore, C(t), computed for solvated PYP at 300 K are plotted in Figure 5 for times out to 100 ps. We find that we can fit the

Figure 5. Dielectric response of the protein and solvent to photoexcitation of the PYP chromophore, calculated using the interaction energy, E, between the chromophore and the rest of the system (jagged curves, green for MD and black for QM/MM). The two results both fit well, other than some oscillations, to a triexponential, eq 6 (smooth curves, blue for MD and red for QM/ MM).

4. CONCLUDING REMARKS Frequency-resolved communication maps provide a coarsegrained view of regions of a protein and water through which vibrational energy flow is facile as a function of the frequency of the vibrations that carry the energy. In this study we have calculated and examined communication maps for PYP from the purple phototropic eubacterium H. halophile to identify the energy transport channels that carry excess energy from the chromophore. Sampling the 300 K communication map for PYP, we find excess energy in both the ground and excited state to be mainly transported to other regions of the protein via neighboring residues that hydrogen bond to the chromophore, specifically Tyr42, Glu46, and Thr50. In addition to these residues, which serve as the gateway for the main energy transport channel from the p-coumeric acid anion chromophore, we observed a secondary energy transport pathway beginning at Arg52, Tyr98, and Gln99, none of which hydrogen bond to the chromophore. Both of the pathways are consistent with those found by Ishikura and Yamato22 for the ground state of PYP on the basis of their calculation on interresidue energy conductivity by using a time-correlation function of the interatomic energy flux. We found similar energy transport channels by the calculation of communication maps for PYP in the ground and excited state. We also carried out QM/MM and MD simulations of the dielectric response of the protein to photoexcitation of the chromophore. The time constants obtained by both simulations span several orders of magnitude, from tens of femtoseconds to hundreds of picoseconds, consistent with the dielectric response of many biomolecules to photoexcitation of a chromophore.51 Of course, PYP undergoes photoisomerization much faster than the longer decay times obtained in the simulation, and in this respect we have more plausibly modeled the dielectric response of PYP analogs where the chromophore does not undergo photoisomerization.41 The aim of this part of the study was to compare time constants obtained using the force field model developed for the study of energy flow, the main focus of the article, with the results of QM/MM

dielectric response, C(t), obtained from the QM/MM and the MD simulations to a triexponential fitting function, Cfit(t ) = a1 exp( −t /τ1) + a 2 exp(−t /τ2) + a3 exp(−t /τ3) (6)

In Table 3 we list the values of the constants in eq 6 that we obtain when we fit this equation to the data plotted in Figure 5. Table 3. Parameters Used To Fit Eq 6 to the Computed Dielectric Response for the Chromophore in Solvated PYP (Figure 5) interactions

a1

τ1/ps

a2

τ2/ps

a3

τ3/ps

chromophore in solvated PYP (MD) chromophore in solvated PYP (QM/ MM)

0.74

0.04

0.21

2.61

0.05

168.56

0.51

0.06

0.47

1.85

0.02

108.91

We find a fast, inertial component on the order of tens of femtoseconds, a slower component on the order of picoseconds, and one or several slower components on the order of tens and hundreds of picoseconds or longer, consistent with time scales for TDSS obtained in many computational and experimental studies on solvated proteins and DNA.51,58−62 We observe in Table 3 that the time constants obtained from the MD simulations using the molecular mechanics force field model for the chromophore in the ground and the excited state compare very well with the time constants we obtain from the QM/MM simulations. The fast component obtained from the MD simulation is 0.04 ps compared to 0.06 ps fit to the results of the QM/MM simulations. The picosecond time scale components also compare very well; we obtain 2.61 ps from the MD simulations and 1.85 ps from the QM/MM simulations. Finally, both approaches yield time constants over 100 ps, specifically 261 ps for the fit to the dielectric relaxation computed with MD simulations and 185 ps for the fit to C(t) 7285

dx.doi.org/10.1021/jp411281y | J. Phys. Chem. A 2014, 118, 7280−7287

The Journal of Physical Chemistry A

Article

In Proteins: Energy, Heat and Signal Flow; Leitner, D. M., Straub, J. E., Eds.; Taylor & Francis Group, CRC Press: Boca Raton, FL, 2009; pp 169−196. (15) Zhang, Y.; Straub, J. E. Directed Energy Funneling in Proteins: From Structure to Function. In Proteins: Energy, Heat and Signal Flow, Leitner, D. M., Straub, J. E., Eds.; Taylor & Francis Group, CRC Press: Boca Raton, FL, 2009; pp 199−228. (16) Zhang, Y.; Fujisaki, H.; Straub, J. E. Molecular Dynamics Study on the Solvent Dependent Heme Cooling Following Ligand Photolysis in Carbonmonoxymyoglobin. J. Phys. Chem. B 2007, 111, 3243−3250. (17) Gnanasekaran, R.; Xu, Y.; Leitner, D. M. Dynamics of Water Clusters Confined in Proteins: A Molecular Dynamics Simulation Study of Interfacial Waters in a Dimeric Hemoglobin. J. Phys. Chem. B 2010, 114, 16989−96. (18) Gnanasekaran, R.; Agbo, J. K.; Leitner, D. M., Communication Maps Computed for Homodimeric Hemoglobin: Computational Study of Water-Mediated Energy Transport in Proteins. J. Chem. Phys. 2011, 135, art. no. 065103. (19) Knapp, J. E.; Pahl, R.; Srajer, V.; Royer, W. E. Allosteric Action in Real Time: Time-Resolved Crystallographic Studies of a Cooperative Dimeric Hemoglobin. Proc. Natl. Acad. Sci. U. S. A. 2006, 103, 7649−7654. (20) Royer, W. E.; Pardanani, A.; Gibson, Q. H.; Peterson, E. S.; Friedman, J. M. Ordered Water Molecules as Key Allosteric Mediators in a Cooperative Dimeric Hemoglobin. Proc. Natl. Acad. Sci 1996, 93, 14526−14531. (21) Buchli, B.; Waldauer, S. A.; Walser, R.; Donten, M. L.; Pfister, R.; Bloechliger, N.; Steiner, S.; Caflisch, A.; Zerbe, O.; Hamm, P. Kinetic Response of a Photoperturbed Allosteric Protein. Proc. Natl. Acad. Sci. U. S. A. 2013, 110, 11725−30. (22) Ishikura, T.; Yamato, T. Energy Transfer Pathways Relevant for Long-Range Intramolecular Signaling of Photosensory Protein Revealed by Microscopic Energy Conductivity Analysis. Chem. Phys. Lett. 2006, 432, 533−537. (23) Leitner, D. M., Frequency-Resolved Communication Maps for Proteins and Other Nanoscale Materials. J. Chem. Phys. 2009, 130, art. no. 195101, 1−9. (24) Ota, N.; Agard, D. A. Intramolecular Signaling Pathways Revealed by Modeling Anisotropic Thermal Diffusion. J. Mol. Biol. 2005, 351, 345−354. (25) Piazza, F.; Sanejouand, Y. H. Long-Range Energy Transfer in Proteins. Phys. Biol. 2009, 6, 046014. (26) Kong, Y.; Karplus, M. Signaling Pathways of PDZ2 Domain: A Molecular Dynamics Interaction Correlation Analysis. Proteins: Struct. Funct. Bioinform. 2009, 74, 145−154. (27) Nguyen, P. H.; Derreumaux, P.; Stock, G. Energy Flow and Long-Range Correlations in Guanine-Binding Riboswitch: A Nonequilibrium Molecular Dynamics Study. J. Phys. Chem. B 2009, 113, 9340−9347. (28) Leitner, D. M.; Straub, J. E. Proteins: Energy, Heat and Signal Flow; Taylor and Francis Press, New York, 2009. (29) Sharp, K.; Skinner, J. J. Pump-Probe Molecular Dynamics as a Tool for Studying Protein Motion and Long Range Coupling. Proteins: Struct. Funcy. Bioinform. 2006, 65, 347−361. (30) Kaledin, M.; Kaledin, A. L.; Brown, A.; Bowman, J. M. Driven Molecular Dynamics for Normal Modes of Biomolecules without the Hessian, and Beyond. In Normal Mode Analysis: Theory and Applications to Biological and Chemical Systems; Cui, Q., Bahar, I., Eds.l CRC Press: Boca Raton, FL, 2005; pp 281−300. (31) Yamato, T. Energy Flow Pathways in Photoreceptor Proteins. Proteins: Energy, Heat and Signal Flow; Leitner, D. M., Straub, J. E., Eds.; Taylor and Francis: New York, 2009; pp 129−147. (32) Gnanasekaran, R.; Leitner, D. M. Dielectric Response and Vibrational Energy Relaxation in Photoactive Yellow Protein: A Molecular Dynamics Simulation Study. Chem. Phys. Lett. 2011, 516, 102−105. (33) Hellingwerf, K. J.; Hendriks, J.; Gensch, T. Photoactive Yellow Protein, a New Type of Photoreceptor Protein: Will This Yellow Lab

simulations, and the two methods produced nearly the same relaxation times. The MD simulations reveal that the dielectric response can be accounted for largely by the response of two residues, Tyr42 and Glu46, both of which hydrogen bond to the chromophore, and both of which, as seen in the communication maps, serve as gateway residues for energy transport when there is excess energy in the chromophore.

■

AUTHOR INFORMATION

Corresponding Author

*D. M. Leitner: e-mail, [email protected]. Notes

The authors declare no competing financial interest.

■ ■

ACKNOWLEDGMENTS Support from NSF CHE-0910669 is gratefully acknowledged. REFERENCES

(1) Jiang, H.; Myshakin, E. M.; Jordan, K. D.; Warzinski, R. P. Molecular Dynamics Simulations of the Thermal Conductivity of Methane Hydrate. J. Phys. Chem. B 2008, 112, 10207−10216. (2) Nguyen, P. H.; Park, S. M.; Stock, G. Nonequilibrium Molecular Dynamics Simulation of the Energy Transfer through a Peptide Helix. J. Chem. Phys. 2010, 132, 025102. (3) Zhang, Y.; Fujisaki, H.; Straub, J. E. Mode Specific Vibrational Energy Relaxation of Amide I and Ii Modes in N-Methylacetamide/ Water Clusters: The Intra- and Inter-Molecular Energy Transfer Mechanisms. J. Phys. Chem. A 2009, 113, 3051−3060. (4) Castonguay, T. C.; Wang, F. Kinetic Monte Carlo Modelling of Chemical Reactions Coupled with Heat Transfer. J. Chem. Phys. 2008, 128, art. no. 124706, 1−8. (5) Clarkson, J. R.; Baquero, E.; Shubert, V. A.; Myshakin, E. M.; Jordan, K. D.; Zwier, T. S. Laser-Initiated Shuttling of a Water Molecule between H-Bonding Sites. Science 2005, 307, 1443−1446. (6) Agbo, J. K.; Leitner, D. M.; Myshakin, E. M.; Jordan, K. D., Quantum Energy Flow and the Kinetics of Water Shuttling between Hydrogen Bonding Sites on Trans-Formanilide (Tfa). J. Chem. Phys. 2007, 127, art. 064315, 1−10. (7) Shenogina, N.; Godawat, R.; Keblinski, P.; Garde, S., How Wetting and Adhesion Affect Thermal Conductance of a Range of Hydrophobic and Hydrophilic Aqueous Interfaces. Phys. Rev. Lett. 2009, 102, art. no. 156101. (8) Acharya, H.; Mozdzierz, N. J.; Keblinski, P.; Garde, S. How Chemistry, Nanoscale Roughness, and the Direction of Heat Flow Affect Thermal Conductance of Solid-Water Interfaces. Ind. Eng. Chem. Res. 2012, 51, 1767−1773. (9) Huang, X.; Liu, G.; Wang, X. New Secrets of Spider Silk: Exceptionally High Thermal Conductivity and Its Abnormal Change under Stretching. Adv. Mater. 2012, 24, 1482−86. (10) Takayanagi, M.; Okumura, H.; Nagaoka, M. Anisotropic Structural Relaxation and Its Correlation with the Excess Energy Diffusion in the Incipient Process of Photodissociated Mbco: HighResolution Analysis Via Ensemble Perturbation Method. J. Phys. Chem. B 2007, 111, 864−869. (11) Sagnella, D. E.; Straub, J. E. Directed Energy “Funneling” Mechanism for Heme Cooling Following Ligand Photolysis or Direct Excitation in Solvated Carbonmonoxy Myoglobin. J. Phys. Chem. B 2001, 105, 7057−7063. (12) Bu, L.; Straub, J. E. Simulating Vibrational Energy Flow in Proteins: Relaxation Rate and Mechanism for Heme Cooling in Cytochrome C. J. Phys. Chem. B 2003, 107, 12339−12345. (13) Fujii, N.; Mizuno, M.; Mizutani, Y. Direct Observation of Vibrational Energy Flow in Cytochrome C. J. Phys. Chem. B 2011, 115, 13057−64. (14) Nagaoka, M.; Yu, I.; Takayanagi, M. Energy Flow Analysis in Proteins Via Ensemble Molecular Dynamics Simulations: TimeResolved Vibrational Analysis and Surficial Kirkwood-Buff Theory. 7286

dx.doi.org/10.1021/jp411281y | J. Phys. Chem. A 2014, 118, 7280−7287

The Journal of Physical Chemistry A

Article

Bring Us Where We Want to Go? J. Phys. Chem. A 2003, 107, 1082− 94. (34) Mizuno, M.; Kamikubo, H.; Kataoka, M.; Mizutani, Y. Changes in the Hydrogen Bond Network around the Chromophore of Photoactive Yellow Protien in the Ground and Excited States. J. Phys. Chem. B 2011, 115, 9306−9310. (35) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. Comparison of Simple Potential Functions for Simulating Liquid Water. J. Chem. Phys. 1983, 79 (2), 926−935. (36) Duan, Y.; Wu, C.; Chowdhury, S.; Lee, M. C.; Xiong, G.; Zhang, W.; Yang, R.; Cieplak, P.; Luo, R.; Lee, T.; Caldwell, J.; Wang, J.; Kollman, P. A Point-Charge Force Field for Molecular Mechanics Simulations of Proteins Based on Condensed-Phase Quantum Mechanical Calculations. J. Comput. Chem. 2003, 24, 1999−2012. (37) Spoel, D. v. d.; Lindahl, E.; Hess, B.; Groenhof, G.; Mark, A. E.; Berendsen, H. J. C. Gromacs: Fast, Flexible and Free. J. Comput. Chem. 2005, 26, 1701−1718. (38) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Montgomery, J. A.; Vreven, T.; Kudin, K. N.; Burant, J. C.; J. M. Millam, S. S. I., Tomasi, J.; Barone, V.; Mennucci, B.; Cossi, M.; Scalmani, G.; Rega, N.; Petersson, G. A.; Nakatsuji, H.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Klene, M.; Li, X.; Knox, J. E.; Hratchian, H. P.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Ayala, P. Y.; Morokuma, K.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Zakrzewski, V. G.; Dapprich, S.; Daniels, A. D.; Strain, M. C.; Farkas, O.; Malick, D. K.; Rabuck, A. D.; Raghavachari, K.; Foresman, J. B.; Ortiz, J. V.; Cui, Q.; Baboul, A. G.; Clifford, S.; Cioslowski, J.; Stefanov, B. B.; Liu, G.; Liashenko, A.; Piskorz, P.; Komaromi, I.; Martin, R. L.; Fox, D. J.; Keith, T.; M. A. Al-Laham, Peng, C. Y.; Nanayakkara, A.; Challacombe, M.; Gill, P. M. W.; Johnson, B.; Chen, W.; Wong, M. W.; Gonzalez, C.; , and Pople, J. A. Gaussian 03; Gaussian, Inc.: Wallingford, CT, 2004. (39) Essmann, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen, L. G.; Smooth, A. Particle Mesh Ewald Potential. J. Chem. Phys. 1995, 103, 8577−8592. (40) Ryckaert, J.-P.; Ciccotti, G.; Berendsen, H. J. C. Numerical Integration of the Cartesian Equations of Motion of a System with Constraints: Molecular Dynamics of N-Alkanes. J. Comput. Phys. 1977, 23, 327. (41) Cordfunke, R.; Kort, R.; Pierik, A.; Gobets, B.; Koomen, G.-J.; Verhoeven, J. W.; Hellingwerf, K. J. Trans/Cis (Z/E) Photoisomerization of the Chromophore of Photoactive Yellow Protein Is Not a Prerequisite for the Initiation of the Photocycle of This Photoreceptor Protein. Proc. Natl. Acad. Sci. U. S. A. 1998, 95, 7396− 401. (42) Allen, P. B.; Feldman, J. L. Thermal Conductivity of Disordered Harmonic Solids. Phys. Rev. B 1993, 48, 12581−12588. (43) Hardy, R. J. Energy-Flux Operator for a Lattice. Phys. Rev. 1963, 132, 168−177. (44) Leitner, D. M. Heat Transport in Molecules and Reaction Kinetics: The Role of Quantum Energy Flow and Localization. Adv. Chem. Phys. 2005, 130B, 205−256. (45) Leitner, D. M. Energy Flow in Proteins. Annu. Rev. Phys. Chem. 2008, 59, 233−259. (46) Yu, X.; Leitner, D. M. Vibrational Energy Transfer and Heat Conduction in a Protein. J. Phys. Chem. B 2003, 107, 1698−1707. (47) Yu, X.; Leitner, D. M. Heat Flow in Proteins: Computation of Thermal Transport Coefficients. J. Chem. Phys. 2005, 122, 054902−1 11. (48) Jimenez, R.; Fleming, G. R.; Kumar, P. V.; Maroncelli, M. Femtosecond Solvation Dynamics of Water. Nature 1994, 369, 471− 472. (49) Bagchi, B. Water Dynamics in the Hydration Layer around Proteins and Micelles. Chem. Rev. 2005, 105, 3197−3219. (50) Changenet-Barret, P.; Plaza, P.; Martin, M. M.; Chosrowjan, H.; Taniguchi, S.; Mataga, N.; Imamoto, Y.; Kataoka, M. Structural Effects

on the Ultrafast Photoisomerization of Photoactive Yellow Protein. Absorption Spectroscopy of Two Point Mutants. J. Phys. Chem. C 2009, 113, 11605−11613. (51) Furse, K. E.; Corcelli, S. A. Molecular Dyamics Simulations of DNA Solvation Dynamics. J. Phys. Chem. Lett. 2010, 1, 1813−1820. (52) Golosov, A. A.; Karplus, M. Probing Polar Solvation Dynamics in Proteins: A Molecular Dynamics Simulation Analysis. J. Phys. Chem. B 2007, 111 (6), 1482−1490. (53) Xu, D.; Martin, C.; Schulten, K. Molecular Dynamics Study of Early Picosecond Events in the Bacteriorhodopsin Photocycle: Dielectric Response, Vibrational Cooling and the J, K Intermediates. Biophys. J. 1996, 70, 453−460. (54) Yu, X.; Leitner, D. M. Chromophore Vibrations During Isomerization of Photoactive Yellow Protein: Analysis of Normal Modes and Energy Transfer. Chem. Phys. Lett. 2004, 391, 181−186. (55) Maisuradze, G. G.; Yu, X.; Leitner, D. M. Normal Mode Analysis and Calculation of the Cooling Rates of the Chromophore Vibrations During Isomerization of Photoactive Yellow Protein. J. Biolog. Phys. Chem 2007, 7, 25−29. (56) Walker, R. C.; Crowley, M. F.; Case, D. A. The Implementation of a Fast and Accurate Qm/Mm Potential Method in Amber. J. Comput. Chem. 2008, 29, 1019−1031. (57) Case, D. A.; Darden, T. A.; T.E. Cheatham, I.; Simmerling, C. L.; Wang, J.; Duke, R. E.; Luo, R.; Merz, K. M.; Pearlman, D. A.; Crowley, M.; Walker, R. C.; Zhang, W.; Wang, B.; Hayik, S.; Roitberg, A.; Seabra, G.; Wong, K. F.; Paesani, F.; Wu, X.; Brozell, S.; Tsui, V.; Gohlke, H.; Yang, L.; Tan, C.; Mongan, J.; Hornak, V.; Cui, G.; Beroza, P.; Mathews, D. H.; Schafmeister, C.; Ross, W. S.; Kollman, P. A. Amber9; University of California: San Francisco, 2006. (58) Halle, B.; Nilsson, L. Does the Dynamic Stokes Shift Report on Slow Protein Hydration Dynamics? J. Phys. Chem. B 2009, 113 (24), 8210−8213. (59) Nilsson, L.; Halle, B. Molecular Origin of Time-Dependent Fluorescence Shifts in Proteins. Proc. Natl. Acad. Sci. U. S. A. 2005, 102 (39), 13867−13872. (60) Li, T.; Hassanali, A. A.; Singer, S. J. Origin of Slow Relaxation Following Photoexcitation of W7 in Myoglobin and the Dynamics of Its Hydration Layer. J. Phys. Chem. B 2008, 112 (50), 16121−16134. (61) Jha, A.; Ishii, K.; Udgaonkar, J. B.; Tahara, T.; Krishnamoorthy, G. Exploration of the Correlation between Solvation Dynamics and Internal Dynamics of a Protein. Biochemistry 2011, 50, 397−408. (62) Pal, S.; Maiti, P. K.; Bagchi, B.; Hynes, J. T. Multiple Time Scales in Solvation Dynamics of DNA in Aqueous Solution: The Role of Water, Counterions, and Cross-Correlations. J. Phys. Chem. B 2006, 110 (51), 26396−26402. (63) Chosrowjan, H.; Taniguchi, S.; Mataga, N.; Unno, M.; Yamauchi, S.; Hamada, N.; Kumauchi, M.; Tounaga, F. LowFrequency Vibrations and Their Role in Ultrafast Photoisomerization Reaction Dynamics of Photoactive Yellow Protein. J. Phys. Chem. B 2004, 108, 2686−98. (64) Mataga, N.; Chosrowjan, H.; Taniguchi, S.; Hamada, N.; Tokunuga, F.; Imamoto, Y.; Kataoka, M. Ultrafast Photoreactions in Protein Nanospaces as Revealed by Fs Fluorescence Dynamics Measurements on Photoactive Yellow Protein and Related Systems. Chem. Phys. Phys. Chem. 2003, 5, 2454−2460.

7287

dx.doi.org/10.1021/jp411281y | J. Phys. Chem. A 2014, 118, 7280−7287