J. Phys. Chem. A 2010, 114, 133–142
133
Reorientations of Aromatic Amino Acids and Their Side Chain Models: Anisotropy Measurements and Molecular Dynamics Simulations Krzysztof Kuczera,†,‡ Jay Unruh,† Carey K. Johnson,† and Gouri S. Jas*,§ Department of Chemistry and Department of Molecular Biosciences, UniVersity of Kansas, Lawrence, Kansas 66045, and Department of Chemistry and Biochemistry, and Institute of Biomedical Studies, Baylor UniVersity, Waco, Texas 76706 ReceiVed: July 31, 2009; ReVised Manuscript ReceiVed: October 5, 2009
We present a study of reorientation dynamics of the three aromatic amino acids and their side chain models in aqueous solution. Experimentally, anisotropy decay measurements with picosecond time resolution were performed for blocked tryptophan, tyrosine, and phenyalanine and model compounds p-cresol and 3-methylindole. Computationally, rotational diffusion was modeled by molecular dynamics simulations for the three aromatic residues and their side chain models: benzene, toluene, phenol, p-cresol, indole, and 3-methylindole in explicit water. Our simulations used the CHARMM protein force field and associated TIP3P water model and tend to underestimate the rotational correlation times. However, the simulations yield several interesting qualitative insights into reorientational motions that complement the experimental measurements. The effects of substituent and temperature on reorientations of the parent compounds are well reproduced computationally. Additionally, simulations indicate strongly anisotropic reorientations for most of the studied compounds and a separation of time scales between conformational dynamics and rotational diffusion. Comparison with continuum hydrodynamic models suggests that we may consider that the blocked amino acids move under stick boundary conditions, while the dynamics for most of the model compounds falls between stick and slip conditions. Our systematic treatment of blocked amino acids, starting from the parent compounds (benzene, phenol, and indole) provides a baseline for understanding the anisotropy decay signals of more complicated peptide systems. Introduction Studies of molecular reorientations in solution provide the most direct access to the rates of nanosecond and subnanosecond molecular motions and interactions with the environment.1-4 To interpret dynamics of larger systems such as peptides and proteins,5,6 it is important to understand motions of the simplest building blocks: individual aromatic amino acids and their sidechain models. The goal of our work is to perform a systematic study of these systems, presenting both basic experimental data to quantify the time scales of the occurring processes and results of molecular dynamics simulations and continuum modeling aimed at providing insights into the microscopic underpinnings of the observations. There are several previously reported measurements of the rotational dynamics of amino acids. Rotational correlation times were measured previously for N-acetylated and C-amidated tyrosine, AcYa2 (59 ps at 5 °C, 37 ps at 20 °C) and N-acetylated and C-amidated tryptophan, acWa7 (47 ps at 20 °C) by Lakowicz and co-workers using frequency-domain methods. We are not aware of any previous measurements of fluorescence anisotropy decays for blocked amino acid residues by time-domain methods. Measurements of the rotational correlation time of the unblocked tryptophan (Trp) amino acid have been reported previously by several laboratories8-10 yielding values of 35-45 ps at temperatures of 20-25 °C. For unblocked tyrosine (Tyr), a rotational correlation time of 36 ps has been reported for Tyr * Corresponding author. E-mail:
[email protected]. † Department of Chemistry, University of Kansas. ‡ Department of Molecular Biosciences, University of Kansas. § Baylor University.
in water at pH 7.11 The aromatic amino acid side chains and their models have also been the subject of several computer simulations. Hu and Fleming, using the CHARMM force field and limited solvation (112-352 TIP3 waters) found reorientation times of 24-28 ps for the La axis and 36-67 ps for the Lb axis of Trp.12 They also studied several substituted indoles, obtaining relaxation times of 6 ps (both La and Lb) for indole, and 14 ps (La) and 20 ps (Lb) for 3-methylindole. Mark and Nilsson performed extensive simulations for zwitterionic Trp in water, employing several water models and different simulation conditions,13 obtaining values of 15-34 ps for the La and 22-45 ps for the Lb axis. They reported that the SPC model, which had the predicted viscosity closest to the water experimental value, also yielded the best predictions of rotational relaxation times. We present a systematic study of aromatic amino acid reorientations in solution. In the experimental part, we performed anisotropy decay measurements for blocked amino acids (Nacetylated and C-amidated phenylalanine, acFa (for the first time), acYa, and acWa) as well as for the model compounds p-cresol and 3-methylindole. In the complementary simulations, we followed molecular dynamics trajectories of the same amino acids and also of their side chain models (benzene, toluene, phenol, p-cresol, indole, and 3-methylindole). See Figure 1. In the simulations we employed the widely used CHARMM protein force field and its associated TIP3P water model. Due to the well-known underestimation of viscosity by the TIP3P water model, our rotational relaxation times tend to be systematically lower than the experimental times. Thus, our simulations have two main goals. The first is providing a baseline of results against which aromatic side chain motions in peptides and
10.1021/jp907382h 2010 American Chemical Society Published on Web 10/29/2009
134
J. Phys. Chem. A, Vol. 114, No. 1, 2010
Kuczera et al. the fluorescence decays were carried out with a nonlinear leastsquares iterative reconvolution software program written inhouse. See Figure 2. The decays of fluorescence polarized parallel and perpendicular to the vertical excitation polarization were collected sequentially. The vertically and horizontally polarized decays were normalized by collecting the same number of counts integrated over a period at time delays long relative to the anisotropy decay times. The decays were globally fit to
1 I|(t)) I(t)[1 + 2r(t)] 3
(1)
1 I⊥(t) ) I(t)[1 - r(t)] 3
(2)
and
Figure 1. Structures and axis systems of the studied molecules. The structure of benzene, toluene, and phenylalanine with acetylated N-terminus and amidated C-terminus are located in the first column (top to bottom). Structures in the second column are phenol, p-cresol, and tyrosine with acetylated N-terminus and amidated C-terminus (top to bottom). In the final column, shown are the structures of indole, 3-methylindole, and tryptophan with acetylated N-terminus and amidated C-terminus (top to bottom). The coordinate system shows the directions of the in-plane axes (x and y) and out-of-plane axis (z).
proteins modeled with CHARMM/TIP3P may be assessed. The second is a qualitative analysis of the influence of temperature, substituents, and conformational dynamics on reorientation times that may be obtained from studies of series of compounds. Finally, we also analyze the shape effects using simple hydrodynamic models. Methods A. Experimental Details. The time-resolved fluorescence studies were carried out with a time-correlated single-photon counting setup described previously11 with a mode-locked, Nd: YAG laser (Coherent Antares) and a synchronously pumped, cavity-dumped Rhodamine 6G dye laser. The dye laser pulses were frequency doubled for excitation at 287 nm (acYa) or 300 nm (acWa). The sample was excited close to the edge of the fluorescence cuvette to avoid significant reabsorption of emission. Emission was collected through a monochromator set at 303 nm (acYa) or 360 nm (acWa). The dye laser pulses were frequency doubled for excitation at 285-287 nm (p-cresol and acYa) or 300 nm (3-methylindole and acWa). Sample concentrations were 30 µM (acWa and 3MI), 100 µM (acYa and p-cresol), and 5 mM (acFa). The sample was excited close to the edge of the fluorescence cuvette to avoid significant reabsorption of emission. Emission was collected through a monochromator set at 303 nm (p-cresol and acYa) or 360 nm (3-methylindole and acWa). The cavity-dumped repetition rate was 7.6 MHz and the count rates were kept below 0.1% of the laser repetition rate to avoid pulse pile-up. acFa was excited with the fourth harmonic of a mode-locked Nd:YAG laser (Coherent Antares) at 266 nm, and the emission was collected at 280 nm. Because of the fast pulse repetition rate (76 MHz), KI (50 mM) was added to the acFa solution to decrease the fluorescence lifetime to ca. 2.5 ns so that fluorescence decay is complete within the 13 ns time window. The instrument response function collected with a dilute colloidal solution of nondairy creamer had a full-width at half-maximum of 100 ps. Nonlinear least-squares fitting of
respectively. The anisotropy r(t) was fit to a single-exponential decay. The fluorescence decay I(t) was fit to single or doubleexponential decays. B. Computational Details. All simulations were preformed with the program CHARMM version 25, using the version 22 protein topology and parameter files.14,15 The blocked amino acids, acFa, acYa, and acWa were constructued with CHARMM in conformations with trans values of the main dihedral angles. The CHARMM residues for side chain models (benzene, toluene, phenol, p-cresol, indole, and 3-methylindole) were created from the protein dictionary, and initial structures were generated from the parameters with CHARMM. The solutes were solvated in a bcc cell on the basis of a cube of side a ) 38 Å. After deletion of the overlapping water, the systems included one solute molecule and ca. 820 TIP3P water molecules. Each system underwent a 50 ps equilibration of solvent with fixed solute, a 50 ps equilibration of the whole system. The final molecular dynamics (MD) trajectories were generated under conditions of constant temperature (either 278 or 298 K) and pressure (1 atm). The trajectories were 12 ns long for the model compounds and 26 ns long for the blocked peptides, so that the simulation lengths were at least 200 times longer than the experimental correlation times. SHAKE constraints were applied to bonds involving hydrogens and a 2 fs integration time step was employed. Coordinates were saved for analysis every 0.2 ps. The van der Waals interactions were truncated at 12 Å, smoothed with a switching function between 10 and 12 Å, while the electrostatic interactions were smoothed by a shift function between at 12 Å. The simulations were performed on an SGI ORIGIN 2400 computer at the Kansas Center for Advanced Scientific Computing. The MD trajectories were used to evaluate the second-order autocorrelation functions C2(t) of the transition dipole axes of the aromatic rings:
C2(t) )
〈
3 cos2 θ(t) - 1 2
〉
(3)
where θ(t) is the axis reorientation angle during time t and the angled brackets denote a trajectory average. For acFa and acYa the transition dipole directions were taken to be in the plane of the benzene rings, in the direction perpendicular to the Cγ-Cζ bond. For acWa the two vectors corresponding to the 1La and 1Lb transitions were employed,
Reorientations of Aromatic Acids
J. Phys. Chem. A, Vol. 114, No. 1, 2010 135
Figure 2. Fluorescence decays for blocked aromatic amino acids measured by time-correlated single-photon counting. Each plot shows the decay of fluorescence polarized parallel to the excitation polarization (with blue fit line), the decay of fluorescence polarized perpendicular to the excitation polarization (with red fit line), and the instrument function. Weighted residuals are shown for fits to the parallel (upper residual plots) and perpendicular (lower residual plots) fluorescence decays. (A) Fluorescence decays of acWa, at 5 °C (left) and 25 °C (right). (B) Fluorescence decays of acYa, at 5 °C (left) and 25 °C (right). (C) Fluorescence decays of acFa at 5 °C.
136
J. Phys. Chem. A, Vol. 114, No. 1, 2010
Kuczera et al.
TABLE 1: Rotational Times of Amino Acid Derivatives sample
r0 (fixed)
φ (ps)
χ2
AcYa, 25 °C AcYa, 5 °C AcWa, 25 °C AcWa, 5 °C p-cresol, 25 °C p-cresol, 5 °C 3-MI, 25 °C + 50 mm acrylamide 3-MI, 5 °C + 50 mM acrylamide AcFa, 5 °C + 50 mM KI
0.32a 0.32a 0.29b 0.29b 0.32a 0.32a 0.29b 0.29b 0.065c
40 ( 5 70 ( 5 48 ( 5 88 ( 7 17 ( 5 22 ( 5 22 ( 6 30 ( 6 98 ( 30
1.11 1.36 1.07 1.08 1.11 1.04 1.07 1.29 3.34
a Lackowicz; et al. Biochemistry 1987, 26, 82.2 b Shen; et al. J. Phys. Chem. 2001, 105, 6260.9 c Duneau; et al. Biophys. Chem. 1998, 73, 109.22
defined as described in Yamamoto and Tanaka.16 The correlation functions were relatively well described by single-exponential decays, although fit quality could be improved by introducing a second exponential. Four approaches were employed to estimate the rotational correlation times from the C2(t) decays (Figure 3). (a) The area under the C2(t) plot, evaluated from t ) 0 to t ) tcut, where the cutoff is defined by the time it takes for the correlation to decay to 1% of its initial value, C2(tcut) ) 0.01. This yields the correlation time for single-exponential decays and the weighted average time for functions described by combinations of several exponentials. (b) The correlation time resulting from fitting C2(t) to a simple one-exponential decay exp(-t/τ) over the first 50 ps. (c) Unrestrained twoexponential fit over the first 50 ps, to a · exp(-t/τ1) + (1 a) · exp(-t/τ2). In this case the value of the second, longer correlation time τ2 was used. (d) Restrained two-exponential fit, with the coefficient a and time τ1 determined from fitting to 0-1 ps the data of, and the time τ2 obtained from fitting the first 50 ps of the decay with a and τ1 kept fixed. Rotational diffusion coefficients Dx, Dy, and Dz around the three molecular axes x, y, and z are simply related to the C1(t) ) 〈cos θ(t)〉 correlation functions, e.g., for the x-axis C1,x(t) ) exp[-(Dy + Dz)t] (Figure 4).17,18 To estimate the diffusion coefficients, we calculated the decay times for the C1 functions of the three axes using the area method. Solution of the three equations (one for each axis) yielded Dx, Dy, and Dz, while error propagation gave the error estimates. The errors are equal for all three coefficients, due to the algebra. A separate axis system was introduced to describe rotations of the whole blocked amino acid residues, rather than the side chains only. In Phe these whole-molecule axes are defined as the CR-Cζ axis (a), the axis connecting the C of the N-terminal and the N of the C-terminal blocking groups (b), and the axis perpendicular to a and b (c, equal to the cross product of a and b). For Tyr the a axis connects CR to the hydroxyl O, and for Trp a is from CR to the Cδ2, with the b and c definitions as for Phe. Results and Discussion Experimental Correlation Times. The experimental data are presented in Table 1. At 5 °C the measured rotational correlation times were 98 ns for AcFa, 70 ps for AcYa, and 88 ns for AcWa. The results for the three blocked amino acids are comparable, roughly coinciding within experimental errors. As expected, reorientations are significantly accelerated at 25 °C, with measured correlation times of 40 ps for AcYa and 48 ns for AcWa. (It was not possible to resolve the anisotropy decay of AcFa at 25 °C because of the longer (ca. 70 ps) pulse widths of the fourth harmonic of the Nd:YAG laser at 266 nm.) The
compounds p-cresol and 3-MI were chosen as models of the Tyr and Trp side chains, respectively. At 5 °C they exhibited correlation times of 22 ns (p-cresol) and 30 ns (3MI), while at 25 °C the times were 17 ns (p-cresol) and 22 ns (3MI). Our measurements for AcYa, AcWa, and the model compounds are in good agreement with previously reported data. The qualitative trends of an approximately 80% increase in rotational rates at 25 °C compared to those at 5 °C approximately agrees with the change in viscosity of water, from 1.52 to 0.89 cP [CRC Handbook of Chemistry and Physics (55th ed., 1974)] in that temperature range. The rotational correlation time of acWa is longer than that of acYa at both 5 and 25 °C, as expected for its increased volume. (Because of the large uncertainty limits for the rotational correlation time of acFa at 5 °C, resulting from the longer pulse widths at the excitation wavelength 266 nm, it is not useful to compare its rotational correlation time with those of acYa and acWa.) Also, faster reorientation of the smaller model compounds compared to that for the corresponding blocked amino acids is in agreement with simple hydrodynamic theory (see below). Simulated Correlation Times. Simulation results for the aromatic side chain models are presented in Tables 2 and 3. In this section, we describe only the reorientations of the transition dipole axes. A more detailed analysis of the reorientations around the different axes is presented in the Anisotropic Rotational Diffusion section. The rotational correlation times of the transition dipole axis (y) are 1.4-1.8 ps for benzene, 2.1-2.3 ps for toluene, and 10.4-11.3 for blocked acFa. Compared to those for benzene, the reorientations of the phenyl ring are about 50% slower in toluene and about 7 times slower in acFa. This suggests that the phenyl side chain is not reorienting by itself but appears to be moving together with the molecule as a whole. The rotational correlation times for the transition dipole axis (y) are 3.8-4.5 ps in phenol, 4.4-5.4 ps in p-cresol, and 14.3-15.1 ps in blocked Tyr. In this case p-methylation slows down the reorientations by about 10-20%, an effect that rises only slightly above the noise level. However, the reorientations of the phenol side chain are significantly slowed down upon incorporation into acYa, by about a factor of 3. The presence of the -OH group, which can hydrogen bond with the solvent water, has a profound influence on the reorientation times. The reorientations are slower by about a factor of 3 in phenol compared to those in benzene, by about a factor of 2 in p-cresol compared to those in toluene and by about 40% in blocked Tyr compared to those in blocked Phe. For indole derivatives (Table 3), we consider reorientations of the two transition dipoles, corresponding to 1La and 1Lb transitions. Interestingly, the rotational correlation times for the two axes are quite similar in the model compounds, 4.3-5.0 ps for indole and 5.4-6.2 ps for 3-methylindole (3MI). 3MI exhibits a systematic tendency for slower reorientations than indole itself; the differences are small, but significant, 20-25%. Two effects are evident in the case of the blocked Trp: slowing down of the reorientations of the axes by a factor of 4-6 compared to the case for indole, and a significantly slower reorientation of the 1Lb axis, by about 50% compared to that for 1La. Again, the slow reorientations in acWa suggest that the indole side chain does not execute independent reorientations, but rather the molecule tumbles as a whole rigid unit. The rotational correlation times of the transition dipole axes of the blocked amino acids are presented in Table 4. At 5 °C the correlation times are 19 ps for acFa, 27 ps for acYa, and 25 and 37 ps for the two axes of acWa. The corresponding values
Reorientations of Aromatic Acids
J. Phys. Chem. A, Vol. 114, No. 1, 2010 137
Figure 3. C2 correlation functions for transition dipole axes from explicit water all atom molecular dynamics simulations: (A) benzene (-), toluene ( · · · ), and blocked phenylalanine, acFa (---); (B) phenol (-), p-cresol ( · · · ), and blocked tyrosine, acYa (---); (C) La axes of indole (-), 3-methylindole ( · · · ), and blocked tryptophan, acWa (---) axes; (D) Lb axes of indole (-), 3-methylindole ( · · · ), and blocked tryptophan, acWa axes (---).
at 25 °C are 13 ps for acFa, 16 ps for acYa, and 20 and 30 ps for acWa. In accord with expectations, the reorientations of each axis speed up upon temperature increase. For the three systems the increase factor is 1.2-1.7. Of the three blocked amino acid residues, reorientations of acFa are the fastest at both temperatures. Rotations of acYa and the 1La axis of AcWa occur at comparable rates, while those of the 1Lb axis of acWa are the slowest. Simulation results are compared with experimental data in Tables 5 and 6. There is a clear systematic effect: all calculated times are lower than measured values. The lowering, corresponding to faster reorientational motion, is by a factor of 3-4 for the model compounds, and a factor of 2-5 for the blocked amino acids. The largest difference occurred for acFa at 5 °C, for which the experimental data have large errors. In the simple hydrodynamic theory of rotational diffusion, which treats the solute as a rigid body and the solvent as a continuum, the rotational correlation times are directly proportional to the viscosity. The predicted viscosity of the TIP3P water model used in our simulations is about 2.5 times lower than the experimental value.13 The deviations of our results from the experimental data are consistent with this property of the water model. The qualitative trends found in the experimental data are correctly, though not perfectly, reproduced by the simulations. Thus, the correlation times tend to increase with system size. At 25 °C the ratio of τ measured for acYa and p-cresol was
2.4, while the calculated one was 3.2. The analogous ratios for acWa and 3MI were: measured, 2.2; calculated, 3.3. and 5.0 (the latter for 1La and 1Lb, respectively). The reorientation times also tended to increase with lowering of the temperature. The measured ratios of τ at 5 and 25 °C were 1.8 for both acYa and acWa. The corresponding simulation ratios were 1.7 for acYa and ca. 1.2 for both axes of acWa. Our results are in generally good agreement with previous calculations for model systems and aromatic acids acids, as described in the Introduction.12,13 Anisotropy of Rotational Diffusion. The calculated rotational diffusion coefficients around the three molecular axes are presented in Table 7. The results are in accord with intuition based on molecular shapes and substituent types. Benzene exhibits the fastest reorientations among all the molecules studied. For this molecule motions around the two in-plane axes, x and y, occur at identical rates, in accord with its high symmetry. Spinning around the out-of-plane axis z is the fastest motion, since it does not require displacement of solvent. In toluene, addition of a methyl group along the x axis leads to slowing down reorientations along y and z, which move the substituent through solvent, while rotation around x is not significantly affected. The fast spinning motion of the benzene ring is abolished in phenol, which contains an -OH group that may hydrogen bond to the solvent water. The presence of the -OH substituent on the x axis in phenol slows down reorientations around y and z compared to the case for benzene.
138
J. Phys. Chem. A, Vol. 114, No. 1, 2010
Kuczera et al.
Figure 4. C1 correlation functions of molecule-fixed axes calculated from MD simulations reveal anisotropic reorientations. Axes defined in Figure 1. (A) Molecule-fixed axes, x (-), y ( · · · ), z (---) of benzene. (B) Molecule-fixed axes, x (-), y ( · · · ), z (---) of toluene. (C) Molecule-fixed axes, x (-), y ( · · · ), z (---) of phenol. (D) Molecule-fixed axes, x (-), y ( · · · ), z (---) of p-cresol. (E) Molecule-fixed axes, x (-), y ( · · · ), z (---) of indole. (F) Molecule-fixed axes, x (-), y ( · · · ), z (---) of 3-methylindole. (G) Molecule-fixed axes, x (-), y ( · · · ), z (---) of acFa. (H) Molecule-fixed axes, x (-), y ( · · · ), z (---) of acYa. (I) Molecule-fixed axes, x (-), y ( · · · ), z (---) of acWa.
Reorientations of Aromatic Acids
J. Phys. Chem. A, Vol. 114, No. 1, 2010 139
TABLE 2: Calculated Rotational Correlation Times for Aromatic Side Chain Modelsa axis x
axis y
axis z
system
a
b
a
b
a
b
benzene toluene AcFa phenol p-cresol AcYa
1.4 2.8 20.4 4.7 7.4 24.7
1.8 3.2 20.2 5.2 8.0 27.0
1.4 2.1 10.4 3.8 4.4 14.3
1.8 2.3 11.3 4.5 5.4 15.1
3.0 3.6 12.4 3.3 4.1 14.3
3.4 3.8 13.2 3.8 4.6 15.1
a Benzene derivatives at 25 °C. Axis x: short in-plane. Axis y: long in-plane. Axis z: ring plane normal. See Methods. Correlation times obtained by single-exponential fit (a) and area method (b). Units: ps. Estimated relative errors: 3-5% for models, 6-8% for amino acids.
TABLE 3: Calculated Rotational Correlation Times for Aromatic Side Chain Modelsa System Indole 3MI AcWa
Axis x
Axis y
Axis z
1
La
1
Lb
a b a b a b a b a b 4.0 4.3 5.0 5.3 6.0 6.2 4.5 5.0 4.3 4.6 5.3 5.9 6.0 7.0 7.7 8.1 5.4 6.0 6.0 6.2 23.3 25.0 19.4 20.9 18.7 19.5 17.2 21.3 28.8 30.1
a Indole derivatives at 25 °C. Axis x: short in-plane. Axis y: long in-plane. Axis z: ring plane normal. 1La ad 1Lb: transition dipoles as in ref 10. See Methods. Correlation times obtained by single-exponential fit (a) and area method (b). Units: ps. 3MI ) 3-methylindole. Statistical errors: 3-5% for models, 6-8% for AcWa.
Surprisingly, the effects of -CH3 and -OH substitution on benzene are qualitatively similar. In p-cresol addition of the -CH3 group along the x axis slows down reorientations around the y and z axes compared to the situation for phenol. In both phenol and p-cresol reorientations are fastest around x, the axis along which the substituents lie, while the spinning motion around z, dragging substituents through solvent, becomes the slowest one. In the case of indole, spinning around the outof-plane z axis is the fastest motion, analogous to the case in benzene. For the indole in-plane axes, reorientations are faster around the long axis x, in agreement with intuition. Methylation of indole at the 3 position slows down reorientations around all three axes. However, both indole and 3MI exhibit similar patterns and rotational diffusion rates, with fastest motion around z and slowest along x. The blocked amino acid residues reorient much more slowly than side chain models themselves. For acFa and acYa, the fastest reorientation is along the x axis, while the reorientations of acWa are almost isotropic. The results clearly show a trend for slowing down of the rotational diffusion upon addition of substituents, as expected. Further, the rotational diffusion of the aromatic rings in the blocked amino acids is markedly slower than in the model compounds. This last observation indicates that the side chains are reorienting together with the whole blocked amino acid residues rather than reorienting independently with respect to the backbone. Dihedral Angle Transitions. Time evolution of the χ1 dihedral angles in the blocked amino acids are shown in Figure 5, while the calculated transition rates are given in Table 8. In all three peptides, only g- and t conformers are found in the MD trajectories. The t form is dominant for acFa and acWa, while g- dominates for acYa. The transition rates are the fastest in acFa, intermediate for acYa, and slowest for acWa. The fastest
transition rate is kg-t ) 2.6 ns-1 in acFa, which corresponds to a relaxation time τr ) 1/kg-t ) 0.4 ns. This is about 10 times slower than the rotational correlation times in Table 6, and about 3-8 times slower than the reorientational rates (Dx, Dy, Dz; Table 7). For acYa the calculated equilibrium constant is close to the measured value. Both the forward and reverse transition rates are overestimated compared to experimental data (by ca. 40% and ca. 80%, respectively). Generally, the potential energy surface for the rotation is well represented by the CHARMM potential: the differences from the experimental data correspond to ca. 0.7 kJ/mol in the free energy differences between the t and g states, and 0.7 and 1.5 kJ/mol for the forward and reverse rates, respectively. Whole Molecule Reorientations in Blocked Amino Acids. For the blocked amino acids, rotational correlation times obtained for axes defined in terms of the chromophores (side chains, see Table 6) are compared with times for axes spanning the whole molecules in Table 9. The results are quite interesting and insightful. In all cases the reorientations of the wholemolecule axes are slower than the side chain-based axes. For almost all axes this difference is quite small, ca. 10-20%, comparable to the statistical errors. The only difference is for acWa, for which the x-axis correlation time is 23 ps, while the a axis time is 33 ps. Thus, the blocked amino acids appear to be moving mostly as whole rigid units, and only along one direction in acWa the reorientation of the side chain is significantly (ca. 30%) faster than reorientation of the whole molecule. Hydrodynamic Modeling: Boundary Conditions. It is interesting to compare the experimental and simulation results to predictions of simple hydrodynamic models, which use Debye-Stokes-Einstein relations to calculate rotational diffusion coefficients for molecules approximated as spheres or ellipsoids. In this approach we model the solvent as a viscous continuum and the solute as a rigid body, either as (a) a sphere with stick boundary conditions or as (b) an ellipsoid with slip boundary conditions. For a sphere with stick boundary conditions the rotational diffusion coefficient D is predicted as
D)
kT 6ηV
(4)
with k the Boltzmann constant, T absolute temperature, η solvent viscosity, and V solute volume. For an ellipsoid with axis ratio a:b:c and slip boundary conditions, we use previously published numerical tables19,20 (with linear interpolation when needed). The results, corresponding to T ) 298 K and water viscosity at that temperature, η ) 0.89 × 10-3 kg m-1 s-1, are presented in Table 10. Typical ratios of ellipsoid long to short axes in our systems are about 2. For this ratio the ellipsoid rotational diffusion with stick conditions is only about 20% slower than for a sphere of the same volume.21 Since taking into account these deviations from spherical shape does not lead to a qualitative change of the model predictions, we limit our discussion to spheres under stick conditions, for the sake of simplicity. Qualitatively, values of the rotational diffusion coefficients D from the stick boundary condition model are the lowest, those of the slip boundary model are the highest, while the MD simulation values fall between the two model predictions. For the model compounds the slip boundary model values are within a factor of 2 above the MD results, while the stick boundary model predictions are within a factor of 5-10 below the MD.
140
J. Phys. Chem. A, Vol. 114, No. 1, 2010
Kuczera et al.
TABLE 4: Calculated Rotational Correlation Times for Transition Dipole Axes of Blocked Amino Acids from MD Simulationsa 5 °C
25 °C
system
a
b
c
d
AcFa AcYa AcWab AcWac
16.9 24.8 24.0 36.5
15.8 22.3 23.5 35.0
21.8 32.4 25.8 37.1
20.2 27.1 26.3 37.9
average 18.9 ( 4.4 26.7 ( 6.9 24.9 ( 2.2 36.6 ( 2.0
a
b
c
d
11.3 15.1 21.3 30.1
10.4 14.3 17.2 28.8
14.7 20.4 21.5 30.0
14.2 14.2 20.4 31.7
average 12.6 ( 3.4 16.0 ( 4.7 20.1 ( 3.2 30.2 ( 1.9
a a ) area method, b ) single-exponential fit, c ) unrestrained two-exponential fit, d ) restrained two-exponential fit (see Methods). Units: ps. Errors are described in Methods. b 1La transition dipole axis. c 1Lb transition dipole axis.
TABLE 8: Rates of χ1 Dihedral Angle Transitions at 25 °C
TABLE 5: Rotational Correlation Times in psa 5 °C system
experiment
p-cresol indole 3MI
22 ( 5
25 °C MD
30 ( 6
MD
experiment
MD
system
17 ( 5 12b 22 ( 6
5(1 5 ( 1/4 ( 1 6 ( 1/6 ( 1
AcFa AcYa AcWa
2.6 0.66 0.54
experiment
ktg-
kg-t
K
1.3 1.8 0.46
1.9 0.36 1.2
kg-t 0.49
ktg1.0
K 0.48
a Experimental: fluorescence anisotropy decay. MD: average rotational correlation times of transition dipole axes from molecular dynamics simulations. b Hu and Fleming, 1991.12
a kg-t: g f t, ktg-: t f g- and equilibrium constant. K ) [t]/[g-]. Units of rates: ns-1 ) 109 s-1. Relative statistical errors are in the 20-30% range for MD data. Expt data: Unruh et al., 2007.23
TABLE 6: Rotational Correlation Times in psa
TABLE 9: Blocked Amino Acidsa
5 °C system experiment AcFa AcYa AcWa
98 (30 70 ( 5 88 ( 7
25 °C
chromophore axes
whole molecule axes
MD
experiment
MD
system
x
y
z
a
b
c
19 ( 5 27 ( 7 25 ( 3/36 ( 2
40 ( 5 48 ( 5
13 ( 4 16 ( 5 20 ( 4/30 ( 2
AcFa AcYa AcWa
20.4 27.7 23.3
16.4 14.3 19.4
12.4 14.3 18.7
22.4 30.1 32.7
17.6 18.5 20.2
13.6 16.7 18.0
a Experimental: fluorescence anisotropy decay. MD: average rotational correlation times of transition dipole axes from molecular dynamics simulations.
a Rotational correlation times from MD simulations, calculated for axes defined for the chromophores (side chains) and for axes defined for the whole molecules. All results from single-exponential fit method, at 25 °C. (See Methods). Units: ps.
TABLE 7: Rotational Diffusion Coefficients for the x, y, and z Axes of Model Compounds and Blocked Amino Acidsa system
Dx
Dy
Dz
error
benzene toluene phenol p-cresol indole 3MI AcFa AcYa AcWa
5.6 6.9 5.0 5.7 2.5 1.4 2.1 1.9 0.8
5.6 2.7 3.7 2.4 3.5 2.7 0.7 0.6 0.7
26. 9.9 3.4 1.6 4.8 3.9 1.1 0.7 0.5
1.2 0.6 0.4 0.3 0.4 0.3 0.2 0.1 0.1
a Calculated from C1 decays; see Methods. Units: 10-2 ps-1 ) 1010 s-1.
The exception here is benzene, for which the slip boundary rotations are predicted to be extremely high, due to the shape of the molecule (see below). In the case of the three blocked amino acids, the slip model values lie as much as 7-15 times above the MD results, while the stick boundary predictions are 2-4 times slower than corresponding MD results. As discussed above for motions of molecular axes, the simulations exhibit reorientation rates about 2-3 times faster than observed experimentally, due to the use of the TIP3P water model. Thus, the stick boundary predictions, which are systematically lower than the MD results, appear to be in better agreement with experimental data. In the case of the slip models, the molecular shapes were approximated by asymmetric ellipsoids. In several cases,
Figure 5. Time evolution of the χ1 dihedral angles of three aromatic residues from molecular dynamics simulations: (A) blocked phenylalanine, acFa; (B) blocked tyrosine, acYa; (C) blocked tryptophan, acWa.
especially for benzene and the three blocked amino acids, the shapes had two axes of similar lengths, making the D predictions quite sensitive to small changes in size. This is because within this model, under slip conditions motions which do not displace solvent proceed without friction, i.e., with extremely high diffusion rates. Overall, comparison of the hydrodynamic models and MD rotation rates leads to the conclusion that the rotation of benzene corresponds to slip boundary conditions, while rotations of the blocked amino acids correspond to a stick boundary. For the
Reorientations of Aromatic Acids
J. Phys. Chem. A, Vol. 114, No. 1, 2010 141
TABLE 10: Comparison of Simulation Results with Hydrodynamic Modelsa system
Vb (A3)
D(stick)c,d
a:b:ce
D(slip)f,d
D(MD)g,d
benzene toluene phenol p-cresol indole 3MI AcFa AcYa AcWa
111 133 120 142 144 166 249 258 283
7.0 5.8 6.4 5.4 5.4 4.6 3.1 3.0 2.7
1.0:0.9:0.5 1.0:0.8:0.5 1.0:0.85:0.45 1.0:0.75:0.4 1.0:0.8:0.4 1.0:0.8:0.4 1.0:0.7:0.6 1.0:0.7:0.6 1.0:0.9:0.6
217 66 89 34 48 42 86 82 106
124 65 40 32 36 27 12 10 6.5
a MD results calculated from C1 decays; see Methods. Molecular volume calculated from CHARMM atomic radii. c For sphere with volume V, using eq 4. d Units of D: ns-1 ) 109 s-1 . e Ratios of molecular size along principal axes of moment of inertia, with largest dimension aligned on x and smallest on z axis. f For ellipsoid with volume V and aspect ratio a:b:c, using table from ref 19. g From our MD simulations (Table 7). b
majority of the models the reorientation rate is intermediate between the stick and slip boundary predictions. These results are in line with trends in the strength of interactions with water, which are expected to be strongest for the relatively polar blocked amino acids and weakest for the solutes of lowest polarity, benzene and toluene. Conclusions We present a study of reorientations of aromatic amino acids and their models in aqueous solution. Experimentally, anisotropy decay measurements with picosecond time resolution are performed for blocked tryptophan, tyrosine, and phenyalanine and for the model compounds p-cresol and 3-methylindole. Computationally, rotational diffusion is modeled by molecular dynamics simulations for the same aromatic peptides and their side chain models (indole, 3-methylindole, phenol, p-cresol, benzene, and toluene) in explicit water. Our systematic treatment of blocked amino acids and models provides a baseline for understanding the anisotropy decay signals of more complicated systems (larger peptides and proteins). Experimentally we find an increase in rotational correlation time in the blocked amino acids compared to those for the sidechain species alone, both for p-cresol compared to acYa and for 3-MI compared to acWa. This result shows that the side chains in the blocked amino acids do not undergo reorientations of large amplitude independent of the backbone to which they are attached, a finding that is confirmed in the MD simulations. The temperature dependence observed for the rotational correlation times accords with the temperature dependence of the viscosity of water. The simulation results include quantities that overlap with the experimental data, as well as properties that are difficult to access experimentally but provide a more complete picture of the nature of the studied reorientational motions. The primary calculation results are the second-order correlation times for the transitional dipole axes, which correspond directly to the measured rotational correlation times. In our calculations we use the CHARMM protein force field and its associated TIP3P water model, which is known to underestimate the viscosity of water by about a factor of 3. The calculated values of rotational correlation times reflect this discrepancy and fall systematically below the measured values. While the simulation results are not in quantitative agreement with experimental data, they provide useful information about several interesting aspects of the reorientational motions of blocked amino acids in solution.
Study of a series of model compounds at two temperatures allows the qualitative description of the influence of substituents on rotational motions. Comparison of the mostly rigid models with the blocked amino acids allows the description of the influence of conformational dynamics on the anisotropy decay signal. The simulation data allow us to follow reorientations around various axes, yielding a description of the anisotropy of the diffusion rates, as well as another view of the effects of conformational flexibility. Predictions of simple hydrodynamic models are compared with the simulation results, providing qualitative information about the nature of the solute-solvent interactions. Overall, our simulations are a baseline against which motions of aromatic side chains in larger peptides and proteins modeled with CHARMM/TIP3P may be assessed. Finally, the systematic deviation from experimental data should provide a stimulus for implementation of improved water models for biomolecular simulations. We see two clear trends in the simulation data that are in accord with the available experimental data and intuition. The first is for a systematic slowing down of reorientations as successively bulkier substituents are added to the parent molecules (benzene and indole). The second is an acceleration of rotations with increasing temperature. These relative trends in rotational dynamics are thus well represented by the models. The simulations indicate that the reorientations of most of the studied molecules are anisotropic, occurring with different rates around different molecular axes. The extreme example here is the smallest studied system, benzene, for which a ratio of 5:1 is predicted for rotations around the ring normal vs in-plane axes. On the other extreme is the bulkiest molecule, acWa, for which the rotations are close to isotropic. These kinds of data are difficult to access experimentally, as typically reorientations for a single transition dipole axis are measured. Perhaps the most interesting question that could be addressed by our studies is the search for the signature of intramolecular conformational dynamics in the anisotropy decay signals of the blocked amino acids. The simulated transition rates of χ1 dihedral angles at 25 °C were in the nanosecond range for all three blocked amino acids, and MD predictions agreed quite well with available measurements for acYa. The conformational transitions were thus much slower than the reorientations, suggesting that the observed bulk anisotropy signal describes an average over contributions from multiple conformers. This finding suggests two interesting conclusions. First, the separation of time scales between conformational transitions and reorientations raises the possibility of observing anisotropy decays of individual conformers in single-molecule experiments. Second, the data indicates that it would be worthwhile to study aromatic reorientations in peptides of moderate size, in which conformational dynamics and overall reorientation time scales would overlap, possibly introducing effects of conformational transitions into polarization anisotropy decay signals. The calculations using simple hydrodynamic models may be used to qualitatively gauge the strength of the solute-solvent interactions. Thus, the motions may be best described by slip boundary conditions for benzene, intermediate conditions for the remaining model compounds, and stick boundary conditions for the blocked amino acids, roughly following the expected strengths of solute-solvent interactions. Overall, our joint experimental-computational study provides a systematic picture of the rotational dynamics of aromatic amino acids and their model compounds in solution. Bringing together measured rotational correlation times with results of molecular dynamics simulations and continuum hydrodynamic
142
J. Phys. Chem. A, Vol. 114, No. 1, 2010
models presents a more complete picture of the reorientational motion than possible by any of the methods separately. The work provides a baseline for analysis of motions in larger peptides and provides suggestions for future studies. Acknowledgment. We thank the Kansas Center for Advanced Scientific Computing for use of computer resources (funded by NSF grant 9627330). G.S.J. thanks Baylor intramural funding support and Shen-En Chen for critical reading. References and Notes (1) Ben-Amotz, D.; Drake, J. M. J. Chem. Phys. 1988, 89, 1010. (2) Lakowicz, J. R.; Laczko, G; Gryczynski, I. Biochemistry 1987, 26, 82. (3) Jas, G. S.; Wang, Y.; Pauls, S. W.; Johnson, C. K.; Kuczera, K. J. Chem. Phys. 1997, 107, 8800. (4) Schaffer, J.; Volkmer, A.; Eggeling, C.; Subramaniam, V.; Striker, G.; Seidel, C. A. M. J. Phys. Chem. A 1999, 103, 331. (5) Rosenberg, S. A.; Quinlan, M. E.; Forkey, J. N.; Goldman, Y. E. Acc. Chem. Res. 2005, 38, 583. (6) Moors, S. L. C.; Jonckheer, A.; Engelborghs, Y.; Ceulemans, A. Current Protein Peptide Sci. 2008, 9, 427. (7) Lakowicz, J. R.; Laczko, G.; Gryczynski, I.; Cherek, H. J. Biol. Chem. 1986, 261, 2240. (8) Chen, L. X. Q.; Engh, R. A.; Fleming, G. R. J. Phys. Chem. 1988, 92, 4811.
Kuczera et al. (9) Shen, X.; Knutson, J. R. J. Phys. Chem. B 2001, 105, 6260. (10) Larsen, O. F. A.; van Stokkum, I. H. M.; Pandit, A.; van Grondelle, R.; van Amerongen, H. J. Phys. Chem. B 2003, 107, 3080. (11) Harms, G. S.; Pauls, S. W.; Hedstrom, J. F.; Johnson, C. K. J. Fluoresc. 1997, 7, 273–282. (12) Hu, Y.; Fleming, G. R. J. Chem. Phys. 1991, 94, 3857. (13) Mark, P.; Nilsson, L. J. Phys. Chem. B 2002, 106, 9440. (14) Brooks, B. R.; Bruccoleri, B. R.; Olafson, B. D.; States, D. J.; Swaminathan, S.; Karplus, M. J. Comput. Chem. 1983, 4, 187. (15) MacKerell, A. D., Jr.; Bashford, D.; Bellott, M.; Dunbrack, R. L., Jr.; Evanseck, J. D.; Field, M. J.; Fischer, S.; Gao, J.; Guo, H.; Ha, S.; Joseph-McCarthy, D.; Kuchnir, L.; Kuczera, K.; Lau, F. T. K.; Mattos, C.; Michnick, S.; Ngo, T.; Nguyen, D. T.; Prodhom, B.; Reiher, W. E.; Roux, B.; Schlenkrich, M.; Smith, J. C.; Stote, R.; Straub, J.; Watanabe, M.; Wiorkiewicz-Kuczera, J.; Karplus, M. J. Phys. Chem. B 1998, 102, 3586. (16) Yamamoto, Y.; Tanaka, J. J. Bull. Chem. Soc. Jpn. 1972, 45, 1362. (17) Ehrenberg, M.; Rigler, R. Chem. Phys. Lett. 1972, 14, 539. (18) Chuang, T. J.; Eisenthal, K. B. J. Chem. Phys. 1972, 57, 5094. (19) Sension, R. J.; Hochstrasser, R. M. J. Chem. Phys. 1992, 98, 2490. (20) Youngren, G. K.; Acrivos, A. J. Chem. Phys. 1975, 63, 3846. (21) Cantor, C. R.; Schimmel, P. R. Biophysical Chemistry: Techniques for the Study of Biological Structure and Function; Freeman: New York, 1980; Vol. II. (22) Duneau, J. P.; Garnier, N.; Cremel, G.; Nulans, G.; Hubert, P.; Genest, D.; Vincent, M.; Gallay, J.; Genest, M. Biophys. Chem. 1998, 73, 109. (23) Unruh, J. R.; Liyanage, M. R.; Johnson, C. K. J. Phys.Chem. B 2007, 111, 5494.
JP907382H