J. Phys. Chem. 1992, 96, 1915-1921
(as observed), it cannot account for the experimentally observed SEDOR decays. Partially segregated samples would result for all compositions in bimodal 31P-113CdSEDOR decays, the faster decaying component of which should reveal 31P-113Cddipoledipole couplings that are stronger than the calculated average. This kind of behavior is not a t all observed in the experiments. An alternative possibility, which is qualitatively consistent with all of the findings here, is shown in Figure 14. This model invokes disorder on both the anion and the cation sublattices. It is then assumed that the phosphorus atoms are, on average, coordinated to more than two germanium atoms and less than two cadmium atoms, whereas the As atoms, conversely, prefer cadmium over germanium coordination. While this tendency could be quantified in principle, we will refrain from doing so here, in view of the simplifying inherent assumptions made in the analysis. Specifically, the contributions from indirect 31P-113Cdspin-spin interactions have been neglected in the SEDOR analysis. Previous 203*205Tl NMR work has shown that, in covalently bonded systems involving isotopes heavier than 13Cdnuclei, the strength of the indirect interactions can exceed that of the direct dipolar couling.^^ Justification for our 113Cd-31PSEDOR approach comes from the experimental observation that this method produces the correct 113Cd-31Pdipolar second moment in crystalline CdP2.24 Furthermore, the isotropic IJ(31P-1I3Cd) value as extracted from the 31PMAS-NMR spectra of the corresponding isotopomers is only 330 Hz. Magnitude and orientation of the anisotropic component of the indirect 31P-113Cdinteraction are currently not known. Such knowledge requires the measurement of singlecrystal rotation patterns, which are planned for the future. (32) Bloembergen, N.; Rowland, T. J. Phys. Reu. 1955, 97, 1679.
1915
Conclusions The atomic distribution in the 11-IV-V2 semiconductor alloys Z I I G ~ A ~ ~ and - ~ PCdGeAs2-,Px , has been studied by X-ray powder diffraction and l13Cd and 31PMAS and spin-echo and 31P-113Cd spin-echo double-resonance NMR. In contrast to the results obtained for 11-VI semiconductors, the MAS-NMR spectra of the present systems show a monotonic chemical shift trend and do not permit resolution of individual peaks for specific nearest neighbor environments. Site-selective 13Cd-31Pspin-echo double-resonance experiments confirm that the chemical shifts in these samples are dominated by nonlocal contributions. The 31Pspin echoes decay somewhat more rapidly than calculated from the crystal structure with randomly occupied sites in the anionic sublattice. This result cannot solely be attributed to un-refocusable 75As31P interactions because the I13Cdspin-echo decays, in conjunction with accompanying simulations, suggest that the 75Asspins fluctuate so rapidly on the time scale of the 31Pspin-echo experiments, that their influence on the 31Pspin echoes is too weak. Conversely, the 31P-113CdSEDOR results in the CdGeAs,P, alloys suggest that the 31P-113Cdinteractions are substantially weaker than calculated from the above model. A possible explanation for these results is a structure involving both cation and anion disorder with a stronger tendency of P-Ge and Cd-As bond formation as compared to P-Cd and As-Ge bond formation.
Acknowledgment. Thanks are due to David A. Lathrop and Francis A. Dempsey for assistance with computer program development. Financial support of this research by N S F Grant DMR-89-13738 is gratefully acknowledged. Acknowledgment is also made to the donors of the Petroleum Research Fund, administered by the American Chemical Society.
Molecular Dynamics Simulations of the Conformational Dynamics of Tryptophan H. L. Gordon,*'+ H. C. Jarrell, A. G. Szabo, K. J. Willis, and R. L. Somorjai Institute for Biological Sciences, National Research Council of Canada, Ottawa, Ontario, Canada Kl A OR6 (Received: July 3, 1991; In Final Form: October 22, 1991)
Molecular dynamics simulations of a single tryptophan molecule were performed using CHARMm-parametrized potential functions, where both explicit molecular and dielectric continuum models were used for the solvent. When all hydrogens were modeled explicitly, rotations about the x2 dihedral were more frequent than about the x1dihedral. The presence of a molecular solvent damped librational motion about both dihedral angles. Conversely, when a united atom representation for hydrogen was used, rotations about the xIdihedral were more frequent than about the x2 dihedral. We show that conflicting rotamer models for the biexponential fluorescence of tryptophan zwitterion can be supported, depending on the model for hydrogen representation in the simulation. We also show that the predicted relative populations of tryptophan rotamers are not consistent with the available experimental data. Simulators are thus warned not to impute universal applicability to empirical potential energy functions found in biomacromolecular molecular mechanics packages.
Introduction A great deal of effort has gone into the development of computational methodologies for molecules of biological interest, such as proteins. Well-known packages, such as CHARMM,] ECEPP,2-4 AMBER,S,6and GROMOS,' have been developed specifically for the computational investigation of biomacromolecules. All of these programs rely on the use of empirical intermolecular potential energy functions. These have been parametrized to reproduce experimental results of selective structural and thermodynamic properties of small molecules and subsequently adjusted for the description of macromolecules. These potential functions and parameters have been designed with the intent that they be widely transferable, so that the energy of any arbitrary *To whom correspondence should be addressed. Canadian Government Laboratory Visiting Fellow, 1990-1991.
protein, for example, can be described given the atom types, their connectivities, and their coordinates. However, it would be unreasonable to expect that these empirically determined potentials (1) Brooks, B. R.; Bruccoleri, R. E.; Olafson, B. D.; States, D. J.; Swaminathan, S.; Karplus, M. J . Compur. Chem. 1983, 4 , 187-217. (2) Momany, F. A.; McGuire, R. F.; Burgess, A. W.; Scheraga, H. A. J . Phys. Chem. 1975, 79, 2361-2381. (3) Dunfield, L. G.; Burgess, A. W.; Scheraga, H. A. J . Phys. Chem. 1978,
82, 2609-26 16. (4) VBsquez, M.; NCmethy, G.; Scheraga, H. A. Macromolecules 1983, 16, 1043-1049. ( 5 ) Weiner, S. J.; Kollman, P. A.; Case, D. A,; Singh, U. C.; Ghio, C.; Alagona, G.; Profeta, S., Jr.; Weiner, P. J . Am. Chem. Soc. 1984, 106, 765-784. (6) Weiner, S. J.; Kollman, P. A.; Nguyen, D. T.; Case, D. A. J . Comput. Chem. 1986, 7 , 23C-252. (7) van Gunsteren, W. F.; Berendsen, H. J. C. Groningen Molecular Sim-
ulation (GROMOS), University of Groningen, Groningen, The Netherlands.
0022-365419212096-1915%03.00/0 0 1992 American Chemical Society
1916 The Journal of Physical Chemistry, Vol. 96, No. 4, 1992
Gordon et al. TABLE I: Six Minimum Energy Rotamers for Tryptophan” rotamer Xlb X,b llp EC
0
\
\
PerPX2 g+ perp-x2 t Perp-x2 ganti-x2 g+ anti-x2 t anti-x2 g-
\
Figure 1. x , , x2, and # dihedral angles in tryptophan: xI N-C,-C& dihedral angle; x2 CuCB-C7-C62dihedral angle; # O-C-C,-C, dihedral angle.
can reproduce or predict all types of experimental results for any given molecule. A concrete example of limitations in the applicability of empirically determined potential functions is given by results of computer simulations of distributed point charge models of water. Although these models have been parametrized to reproduce certain structural and thermodynamic properties of liquid water,8 some models completely fail to “predict” the large experimental value of the bulk dielectric constant of water at rcom temperature? Such deficiencies do not entirely rule out the use of these models, but caution is necessary if they are used in studies of, say, electrolyte solutions, where the dielectric screening properties of water are important. We issue a caveat that simulations employing different representations of hydrogen atoms may not agree, even qualitatively. We use, as a specific example, the contradictory results when two different representations of hydrogen atoms are used for the description of tryptophan. Under most circumstances, both the explicit and the united atom representations of hydrogen give very similar results when used in studies of proteins and peptide~.~J However, results presented here agree with those of ref 10, in that the united atom representation is found to provide a very different pattern of conformational dynamics of trytophan than does the explicit hydrogen representation.
Rotamer Models of Tryptophan Time-resolved fluorescence spectroscopy of proteins and polypeptides can often given useful structural and dynamic information.”i12 The naturally occurring fluorescent amino acid tryptophan is widely used as the fluorescent probe in these studies. Multiexponential fluorescence decay kinetics are frequently observed for tryptophan, its derivatives, and for individual tryptophyl residues in proteins. This complexity is thought to originate from the ground-state heterogeneity arising due to the presence of rotamers of the side ~ h a i n . I ~However, -~~ recent data for tryp(8) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. J. J . Chem. Phys. 1983, 79, 926-935. (9) Neumann, M. J. Chem. Phys. 1985,82, 5663-5672. (10) Axelsen, P. H.; Gratton, E.; Prendergast, F. G. Biochemistry 1991, 30, 1173-1 179. (11) Beechem, J. M.; Brand, L. Annu. Reo. Biochem. 1985, 54, 43-71. (12) Eftink, M. In Methods of Biochemical Analysis, Volume 35: Protein Structure Determination; Suelter, C. H., Ed.; John Wiley: New York, 1991; pp 127-205. (13) Donzel, B.; Gauduchon, P.; Wahl, Ph. J. Am. Chem. Soc. 1974, 96, 80 1-808. (14) Gauduchon, P.; Wahl, Ph. Biophys. Chem. 1978,8, 87-104. (15) Fleming, G. R.; Morris, J. M.; Robbins, R. J.; Woolfe, G. J.; Thistlethwaite, P. J.; Robinson, G. W. Proc. Narl. Acad. Sei. U.S.A. 1978, 75, 4652-4656. (16) Szabo, A. G.; Rayner, D. M. J . Am. Chem. SOC.1980,102,554-563. (17) Robbins, R. J.; Fleming, G. R.; Beddard, G. S.; Robinson, G. W.; Thistlethwaite, P. J.; Woolfe, G. J. J. Am. Chem. Soc. 1980, 102,6271-6279.
60.1 -178.6 -60.4 60.8 -117.5 -60.6
91.9 85.3 91.5 -90.3 -93.4 -85.4
-68.8 65.4 65.8 -62.9 69.9 65.8
Conditions described under Calculations. Potential energy (in kcal/mol).
1.228 1.528 1.682 1.347 1.325 1.724
See Figure 1.
tophan zwitterion suggest that excited-state reaction may also be important.27 Individual rotamers are expected to decay monoexponentially, with lifetimes dependent on the position of the side chain relative to the indole ring. In order to develop these rotamer models, it is clearly important to identify the preferred conformations and determine their relative populations and interconversion frequencies. Conformational energy calculations suggest six low-energy (x,, x2) (see Figure 1) rotamers for t r y p t ~ p h a n . ~ N~ M , ~R~stud,~~ ies3w33and surveys of the conformation of tryptophyl residues in proteins (from X-ray structural data)34,35identify similar preferred rotamers. Six unassigned rotamers were also found for the neutral form of tryptophan in cold supersonic molecular beams.3b38The interconversion frequencies between these various rotamers are uncertain. Insights into the conformational dynamics of tryptophan, both as a single molecule and as a constituent amino acid of proteins, have been provided by molecular dynamics (MD) simulations. 10*22,39 In order to rationalize some recent experimental data on tryptophan fluorescence, we attempted to reproduce and extend some previous MD results.22 We found that the frequency of dihedral transitions of tryptophan zwitterion depended upon what representation was used for the hydrogen atoms in the simulation. As a consequence, results from MD simulations on tryptophan (18) Petrich, J. W.; Chang, M. C.; McDonald, D. B.; Fleming, G. R. J. Am. Chem. SOC.1983, 105, 3824-3832. (19) Gudgin, E.; Lopez-Delgado, R.; Ware, W. R. J . Phys. Chem. 1983, 87, 1559-1565. (20) Gudgin-Templeton, E. F.; Ware, W. R. J. Phys. Chem. 1984, 88, 4626-4631. (21) James, D. R.; Ware, W. R. Chem. Phys. Lett. 1985,120,450-454. (22) Engh, R. A.; Chen, L. X.-Q.; Fleming, G. R. Chem. Phys. Lett 1986, 126, 365-372. (23) Boens, N.; Janssens, L. D.; De Schryver, F. C. In Ultrafast Phenom-
ena V& Yajima, T., Yoshihara, K., Harris, C. B., Shionoya, S., a s . ; Springer-Verlag: Berlin, 1988; pp 486-488. (24) Tilstra, L.; Sattler, M. C.; Cherry, W. R.; Barkley, M. D. J . Am. Chem. SOC.1990, 112, 9176-9182. (25) Colucci, W. J.; Tilstra, L.; Sattler, M. C.; Fronczek, F. R.; Barkley, M. C. J. Am. Chem.Soc. 1990, 112, 9182-9190. (26) Chen, R. F.; Knutson, J. R.; Ziffer, H.; Porter, D. Biochemistry 1991, 30, 5184-5195. (27) Willis, K. J.; Szabo,A. G.; Krajcarski, D. T. Chem. Phys. Lett. 1991, 182, 614-616. (28) Zimmerman, S. S.; Pottle, M. S.; NBmethy, G.; Scheraga, H.A. Macromolecules 1977, 10, 1-9. (29) Anderson, J. S.; Bowitch, G. S.; Brewster, R. L. Biopolymers 1983, 22,2459-2476. (30) Cavanaugh, F. R. J . Am. Chem. SOC.1970, 92, 1488-1493. (31) Baici, A.; Rizzo, V.; Skrabal, P.; Luisi, P.L. J. Am. Chem. Soc. 1979, 101, 5170-5178. (32) Kobayashi, J.; Higashijima, T.; Sekido, S.; Miyazawa, T. Inr. J. Pept. Protein Res. 1981, 17, 486-494. (33) Dezube, B.; Dobson, C. M.; Teague, C. E. J . Chem. SOC.,Perkin Trans. 2 1981, 73C-735. (34) Janin, J.; Wodak, S.; Levitt, M.; Maigret, B. J . Mol. Biol. 1978, 125, 357-386. (35) Bhat, T. N.; Sasisekharan, V.; Vijayan, M. Int. J. Pept. Protein Res. 1979. 13. 17C-184. (36) k, T. R.; Park, Y. D.; Peteanu, L. A.; Levy, D. M. J. Chem. Phys. 1986,84, 2534-2541. (37) Rizzo, T. R.; Park, Y. D.; Levy, D. H. J. Chem. Phys. 1986, 85, 6945-6951. (38) Philips, L. A.; Webb, S. P.; Martinez, S. J.; Fleming, G. R.; Levy, D. H. J. Am. Chem.Soc. 1988, 110, 1352-1355. (39) Hu, Y.; Fleming, G. R. J . Chem. Phys. 1991, 94, 3857-3866.
The Journal of Physical Chemistry, Vol. 96, No. 4, 1992 1917
Conformational Dynamics of Tryptophan t
9-
0 CY
I
0
o?
TABLE 11: Details of Molecular Dynamics Simulations" starting solvent heating equilib collectn run conforme9 model period: ps period: ps Deriod. DS 3 10 1O W 1 perp-x,g+ continuumC 3 100 l00og 2 perp-x2 t continuumC 3 500 1 OOOh 3 perp-x2 g- continuumC TIP3pd 3 50 500' 4 perp-x2 t "Verlet algorithm; MD time step, 1 fs for all runs; SHAKE constrained covalent bonds to hydrogen; no potential cutoff. bSee Table I. c t = c j i j ; e,, = 78. dReference 8. '5 K of kinetic energy (taken from a Gaussian distribution of velocities) assigned every 50 MD time steps to bring system from 0 to 300 K. /Atom velocities scaled every 20 MD time steps if I(average temperature - 300 K)I > 5 K. # N o scaling of velocities. Atom velocities scaled every 100 MD time steps if I(average temperature - 300 K)( > 10 K. 'Atom velocities scaled every 50 MD time steps if ((average temperature - 300 K)( > 10 K.
0
180 -120 -60 0 x1 Figure 2. Potential energy surface for (x,, x2) rotamers of tryptophan. Contours are plotted in levels of 1 kcal/mol. Regions of potential minima indicated by dashed contours.
60
120
can be used to support contradictory rotamer models for the multiexponential fluorescence decay.
Calculations All our potential energy minimizations and MD simulations were done using the CHARMm package.'-@ The full hydrogen representation was chosen for modeling tryptophan; that is, all hydrogen atoms were included explicitly (27 atoms total). The SHAKE a l g ~ r i t h mwas ~ ' ~used ~ ~ to constrain only covalent bonds to hydrogen. Where the solvent was modeled as a dielectric continuum, a dielectric function E = toruwas selected, where rij is the atom pair separation. The chosen value of E, was 78, the bulk value of water at 300 K."3 (We could have chosen a smaller value of E, given the linear riidependence of the dielectric function. However, as we were interested in qualitative rather than quantitative results, and given that a linear dielectric function is only a crude approximate representation of the electrostatic screening afforded by the solvent, the exact value chosen for E, was not considered to be critical.) The TIP3P water model* was used when a molecular solvent model was needed. (Note that the value of the bulk dielectric constant E, for TIP3P water is unknown, but it is likely to be similar to those of the analogous TIP4P and SPC water models, which have E, = 53 f 2 a t 293 K4"46 and 65 f 9 a t 300 K,45-47respectively). All nonbonded interactions were included in the computation of the potential energy of tryptophan. In addition, when molecular solvent was modeled explicitly, all nonbonded interactions were computed between tryptophan and TIP3P and between TIP3P and TIP3P molecules. The minimum energy structures for the six ( x l , xz) rotamers of tryptophan (with the solvent modeled as a dielectric continuum) were found by using an adopted basis Newton-Raphson algorithm' with an energy difference convergence criterion of 0.0001 kcal/mol. Three local minima were found for each of the six rotamers, corresponding to different orientations of the COOgroup. Table I contains the ( x l , x2) values and potential energies (40) CHARMm version 21.2, Polygen Corp., Waltham, MA, 1990. (41) Ryckaert, J. P.; Ciccotti, G.; Berendsen, H. J. C. J . Compuf. Phys. 1977, 23, 327-341. (42) van Gunsteren, W. F.; Berendsen, H. J. C. Mol. Phys. 1977, 34, 13 11-1 327. (43) Uematsu, M.; Franck, E. U. J . Phys. Chem. ReJ Data 1980, 9, 1291-1306. (44) Neumann, M. J . Chem. Phys. 1986, 85, 1567-1580. (45) Gordon, H.; Goidman, S. Mol. Simul. 1989, 2, 177-187. (46) Alper, H. E.; Levy, R. M. J . Chem. Phys. 1989, 91, 1242-1251. (47) Watanabe, K.; Ferrario, M.; Klein, M. L. J . Phys. Chem. 1988, 92, 8 1 9-82 1.
TABLE III: Results from Molecular Dynamics Simulations" x2 transitionsd proportion of starting x1 perp anti collection timee
- -
run 1 2 3 4
conformerb transitionsC anti perP-x2g+ 0 3 0 3 perp-x2 t 0 2 perp-x, g0 1 perp-x, t
- - --
perp
2 2 1 0
perp-x2 anti-x2 0.90 0.10 0.28 0.72 0.81 0.19 0.20 0.80
"See Table 11. *See Table I. CNumberof transitions during colt, t g-, g' g-. "Number of transitions lection period: g' during collection period: perp-x2 anti-x,. e Approximate proportion of simulation time spent in conformation during collection run.
of the lowest energy COO- conformer for each of the six rotamers. Figure 2 is a contour plot of the minimum potential energy surface as a function of x1 and x2conformations of tryptophan in a continuum solvent. Both x1and xz angles were varied in 30° increments and were constrained during the subsequent energy minimization with respect to all other internal coordinates. (By constrained, we mean that the xI and xzdihedral angles were held close to their initial values by application of a sufficiently large restraining force.) Minima were located initially using an adopted basis Newton-Raphson algorithm with an energy difference convergence criterion of 0.0001 kcal/mol. One to three local minima were located for each fixed set of (xl, x2), corresponding to different orientations of the COO- group. This was ascertained by local minimization of each (xl, x2) conformer, starting from 12 different orientations of rC, that were separated by 30° intervals ( x l and x2 were constrained to their initial values during the minimizations; $was not). Figure 2 was then constructed using the lowest energy COO- conformer found for each (xl. x2). Starting from each putative local minimum, full Newton-Raphson minimization was performed to help ensure that the critical points found for the potential were true minima rather than saddle points. Table I1 lists the details of three MD simulations for tryptophan, using the continuum solvent representation. The g+, t , and g(see Figure 2) conformers in the perpendicular x2 (perp-x2) orientation, as described in Table I, were used as the starting conformations for the three simulations. The equations of motion were solved using the Verlet algorithm. The MD time step was 1 fs. In each case, kinetic energy was added to these minimized structures over a period of 3 ps to heat the system to 300 K. After an initial equilibration period during which the total energy was stabilized by scaling the atom velocities, data were collected over a period of 1 ns. Velocity scaling was applied only during the data collection period of the simulation of the g- rotamer (run 3). In this case, the total system energy was not conserved without scaling of the velocities even if the time step was reduced to 0.1 fs. Table 11also contains the details of one MD simulation starting from the perpx, t conformation immersed in 109 TIP3P water molecules contained in a cubic box of length 15.55 A. Periodic boundary conditions were used in order to mimic an infinite system. The potential cutoff for nonbonded interactions was
Gordon et al.
1918 The Journal of Physical Chemistry, Vol. 96, No. 4, 1992
t - 0
X S
r
i
TABLE I V Relative Populations of Tryptophan Rotamers Pg+ p, p*total
This Work' perp-xz
anti-X2 total
0.54 perp-X2
anti-x, total
anti-X2
pep-X2
0
50
100
150
perp-X2
0
50
100
150
t 200
Time (ps) Time (ps) Figure 3. Molecular dynamics trajectories of x1and x2dihedral angles for t conformer of tryptophan, given continuum and molecular models for solvent. Values of x, and x2 plotted every 0.1 ps over initial 200 ps of collection period for runs 2 and 4 (Table 11). (a) x1versus simulation time where continuum solvent model used; (b) xI versus simulation time where molecular solvent model used, (c) x2Venus simulation time where continuum solvent model used; (d) x2 versus simulation time where molecular solvent model used. Values of the dihedral angles in the range -180' to 0" are plotted in the equivalent range 180' to 360".
Figure 4. Normalized Boltzmann population for (xl, x2)conformers of tryptophan. Computed from (xl, x2) potential energy surface and eq 1 at 300 K.
chosen to be the same as the box length. The MD time step was 1 fs. The nonbonded lists were updated every 50 MD time steps. Kinetic energy was added to the system over 3 ps in order to heat the system to 300 K. Atomic velocities were scaled throughout the 50-ps equilibration and 500-ps data collection periods.
Discussion of Computational Results Our results are presented in Figures 2-4 and in Tables I11 and IV. In Figure 2, there are six clearly defined potential wells corresponding to regions close to the expected22*28-35 six.(xl, x2) rotamers of tryptophan. The potential energy surface indicates that all the barrier heights between perp-x2 and antiperpendicular x2 (anti-Xz) forms of g+, 1, and g- rotamers are lower (3-5 kcal/mol) than those between different x, rotamers having the same x2 orientations (5-7 kcal/mol). These relative barrier heights would imply that, for this system, x2 interconversions would occur more frequently than x I interconversions. Figure 2 is substantively different from Figure 2 in ref 22. Engh et aLZ2used the GROMOS package and a united atom repre-
0.22 0.19 0.41
0.14 0.22 0.36
0.13
0.49
0.10 0.23
0.51 1 .oo
Reference 29b 0.30 0.16
Reference 33, Table 4c 0.05 f 0.05 0.05 i 0.05 0.70 f 0.05 0.80 f 0.15 0.15 f 0.05 0.05 f 0.05 0.00 f 0.05 0.20 f 0.15 0.20 f 0.10 0.10 f 0.10 0.70f 0.10 1.00
Reference 33, Table 2 from 0.32 f 0.10 0.15 f 0.10 0.53 f 0.10 from NOE 0.15 f 0.15 0.10 f 0.10 0.75 f 0.10
c!k
"Computed using eqs 1 and 2: ' 0 Ixl(g+)5 120°, 120" 5 xl(t) 5 240°, 240' I xl(g-) I360'; 0" I x2(perp) 5 180°, 180° 5 x2(anti) I360'. bReference 29 used an ECEPP potential, an explicit representation of hydrogen, and an effective dielectric constant of 2. 'Determined from lanthanide shift, relaxation data, and NOE data. sentation of hydrogen, whereby only polar hydrogens are explicitly included in the model for tryptophan (19 atoms total). Their (x,, x2) potential energy contours imply a more intricate potential surface, especially in regions separating the six potential wells. Since only 144 points were used to construct our Figure 2, the appearance of the potential energy surface is likely to be smoother than if a finer grid search had been performed. Additional data would reveal more local variability in the equipotential energy contours. (The smoothness of contour lines depends both on the fitting routines used to produce them and on the level of detail of the original data. The lack of true periodicity a t the borders of Figure 2 would also be rectified by using a more finely detailed grid search.) In addition, differences between this work and that in ref 22 are probably attributable to the different parametrizations of CHARMm and GROMOS and/or the use of a full hydrogen rather than a united atom approach for treating hydrogens. In particular, recent simulations of a single tryptophan residue in phospholipase A2losuggest that the use of a united atom representation for hydrogen, as was used in ref 22, gives rise to simulation results that are not in agreement with fluorescence anisotropy experiments. The CHARMm MD simulations in ref 10, which compared the united and full atom representation of hydrogen atoms, found that use of the former resulted in too rapid rotation of the tryptophan-3 side chain. Axelsen et a1.I0 attributed this to "an inherently bland representation of aromatic ring electrostatics" in the united atom model, so that interactions between the indole ring and explicitly modeled water (solvent) molecules are not successfully described. However, we have found that the conformational dynamics of tryptophan depend on the representation used for hydrogen, even when a continuum solvent is used. We carried out a few sample computations on tryptophan in a continuum solvent using the united atom representation of hydrogen within CHARMm. These indicated that the potential barriers separating the six rotamers were very different from those shown in Figure 2. For example, the potential barrier height between p e r p x 2 and anti-X2 conformations of g+ was -6 kcal/mol for the united atom representation and -3 kcal/mol for the full hydrogen representation of tryptophan. The barriers between perp-xzg+ and t conformations were, respectively, -5.5 and -4 kcal/mol for the full and united atom representations of hydrogen. The barriers between perpxz t and g- rotamers were, respectively, - 5 and -1 kcal/mol for full and united atom representations of hydrogen. A separation of the components of the total intramolecular potential revealed that, as a function of the xz dihedral, the sum of dihedral energies for the united atom representation had much larger potential barriers than that of the full hydrogen representation. This was not the case as x1 was varied for a fixed value of x2;instead, the van der Waals, bond,
The Journal of Physical Chemistry, Vol. 96, No. 4, 1992 1919
Conformational Dynamics of Tryptophan and angle energies of the united atom representation had lower potential barriers than those of the full hydrogen representation, especially in the region of the g conformer. It appears that, in the case of a continuum solvent model, the differences between potential energies for the full and united atom representations of tryptophan are due mainly to nonCoulombic interactions and may be partitioned between various nonbonded and bonded components, not necessarily the dihedral potential. Consequently, the CHARMm potential energy surfaces for the full and united atom representations of tryptophan will be quite different. This can account, at least in part, for the discrepancies between Figure 2 in ref 22 and our Figure 2. As a consequence of the different descriptions of the potential surface, trajectories of MD simulations of tryptophan using the full versus united atom representation of hydrogen will be different. The above comparison of barrier heights implies that where the united atom representation for hydrogen is used to describe tryptophan, xI interconversions will occur more frequently than x2 interconversions. This prediction was confirmed by a short 100-psMD simulation. Starting from the perp-Xz t rotamer, many interconversions between the t and g- conformers (the barrier height is only -1 kcal/mol) and no x2 interconversions (the barrier heights between perp-Xz and anti-xz forms o f t and g are 4-6 kcal/mol) were found. These results are consistent with those reported in ref 22, in which many x1 and few or no x2 transitions occurred during a typical 100-ps MD simulation. However, it is notable that Engh et a1.22reported many more t g+ than t g interconversions, while our short simulation had many t g and no t g+ interconversions. Thus, it is likely that the relative barrier heights separating xI rotamers are described differently by the united atom potentials of GROMOS and CHARMm. Discrepancies between our results and ref 22 also are due to different parametrizations of the GROMOS and CHARMm potentials, not only to the different representations of nonpolar hydrogens. To summarize, the full versus united atom representations of hydrogen provide very different descriptions of the conformational dynamics of tryptophan. Thus, simulation results should always be checked to ascertain consistency with experiment. Subsequent work by Hu and Fleming39indicated that simulations using the full hydrogen representation of CHARMm gave better agreement with rotational reorientation results from experiment than did simulations using GROMOS and a united hydrogen atom model, at least when a molecular solvent was employed. Table I11 contains results from the three MD simulations employing a continuum dielectric model for the solvent. All interconversions between rotamers involved changes between perp-Xz and anti-;yzconformations; no x1transitions were ever observed. The result that x2 interconversions are more frequent than xI interconversions is consistent with the relative barrier heights separating the six rotamers shown in Figure 2. Interconversions between x2 conformations did follow the lowest energy pathways as indicated by Figure 2; that is, the x2 transitions varied between -*90' through a 180° conformation and with little change in
--
-
-
XI.
Our simulations were not long enough to establish the relative populations of each of the six rotamers, given that no x1transitions were observed and that too few x2transitions were seen for any reliable statistical analysis. However, Table I11 indicates that, for a given xI conformer, the molecule spent more time in the region of the lower energy xz conformation. Table I11 also contains the results from the simulation using a molecular solvent. Figure 3 compares the values of ( x l , xz) dihedral angles over a 200-ps segment of two simulations (runs 2 and 4 in Table 11) started from the same perp-x2 t rotamer but where the two different representations of the solvent were used. Librational motion about both xI and x2torsional angles is damped by the presence of the molecular solvent. This conforms with expectations that the concomitant movement of the surrounding solvent molecules will hinder internal motion of the solute. No x, transitions were observed during this simulation. It is not possible to say whether the use of a molecular solvent reduces the
frequency of xz transitions since only one transition was observed during this relatively short MD simulation. The time taken for a x2transition was longer (- 10 ps) in the molecular solvent than when the dielectric continuum model was used (1-2 ps). The potential well in the region of the anti-Xz t rotamer, shown in Figure 2, is broad and elongated in the x2 direction. Figure 3 shows that, in both the continuum and molecular solvent simulations, the anti-x2 t conformer makes frequent excursions to a conformationhaving xzvalues in the range -120O Ixz I-180'. A shallow local minimum was located (within the continuum solvent representation),having (xl,xz,$) = (178.4°,-145.10,63.80), this being the lowest enegy (1.605 kcal/mol) conformer of three COO- rotamers. Similar results were found in the simulation starting from the perp-Xz g- rotamer, where the libration about x2was not symmetric about -90°, but showed frequent excursions to a conformation with 144O Ix2 I -144'. A shallow local minimum was located, with (x1,x2,+) = (-57.0°,144.80,66.10) being the lowest energy (2.827 kcal/mol) conformer of three COO- rotamers. The existence of these shallow local minima is likely to depend strongly on the particular parameters used to construct the potential energy of the system. Relative populations of the six rotamers would best be found from simulation results, since all (except for the constrained covalent bonds to hydrogen) degrees of freedom in the intramolecular potential are exercised during MD simulations. However, as an approximation, a Boltzmann population analysis of the two-dimensional (xl, x2) potential surface shown in Figure 2 can be used. Figure 4 is a three-dimensional plot of the (scaled) normalized Boltzmann populations calculated at 300 K of the (xl, x2) potential energy surface in Figure 2. It was computed by P, =
expl-E(Xli,Xzj)/Rfl
EE exp(-E(Xli,Xzj) /RT) i j
(1)
where Pijis the Boltzmann population, given potential energy E at conformation (xli, xZj) of the molecule, R is the ideal gas constant, and T i s the temperature. Integration over each peak gives the population fraction of each rotamer. Dezube et used NMR coupling constants and relative nuclear Overhauser effects (NOE) to compute rotamer populations. We have used their method of analysis to compute putative coupling constants and NOES, given the rotamer populations shown in Figure 4. If rotation about x2 is much more rapid than rotation about x1with respect to the NMR time scale, then NMR measurements can be analyzed in terms of Boltzmann weighted averages over the three xI (g+, t, and g)rotamers. The fractional populations for these three rotamers are then calculated as
where x; min and x'; m x are the minimum and maximum values of x, designated as describing rotamer r. Following ref 33, the observed coupling constants and between the C, proton Ha,and the two protons, H,, and Hb2,on C, are given by
(3) In eq 3 Ja,ddki) = KI COS2 4 k i
+ K2 COS dki + K3
(4)
where dki is the ith value of the N-C,-C,-HBk dihedral angle. As in ref 33, we used (K,,Kz,K3)= (9.4,-1.4,1.6) in eq 4. Note that JaBk is a function of C#Jk and hence of x1 but not of x2. To calculate the ratio of NOE effects on H,, and H,, as a consequence of saturation of the Haresonance from the potential energy surface in Figure 2, we use
1920 The Journal of Physical Chemistry, Vol. 96, No. 4 , 1992
-=
. .
where ruSkis the distance from H, to HBk,and (rapk), is the Boltzmann weighted average of r4ekfor the rth rotamer. Equation 5 assumes that librational motion within the potential wells is rapid with respect to the overall tumbling of the molecule and that the interchange between the g+, t , and g- rotamers is rapid with respect to relaxation time Ti. Table IV contains the populations of the six tryptophan rotamers as computed from the potential energy surface in Figure 2 and from lanthanide shift data.33 Table IV also contains the relative populations of the three g+, t , and g- rotamers, assuming that x2 interconversion is rapid with respect to the N M R time scale; determined from experimental values of .I$ and , and zapl/ z,82.33Although there are large error estimates attached to the experimentally determined population fractions, our results are not even in qualitative agreement with experiment. Both the experimental coupling constants and NOE effects predict that the g- rotamer is the most populated and that the t rotamer is the least populated. The potential surface in Figure 2 predicts that the g+ rotamer is the most populated and the g- is the least populated. in ref The experimentally determined values of E:l and 33 were 8.2 f 0.2 and 4.8 f 0.2 Hz, resrtively. As computed from eqs 3-5, our values of and JDns2 were 7.0and 5.4Hz, not within the experimental error. However, given the very different population fractions, our results are not totally unreasonable: the coupling constants are not very sensitive functions of P,. Note that eq 4 is an even function of f#Jkj and so cannot discriminate easily between g+ and g rotamers, whose magnitudes of x1are very similar. Using eq 6,our result of 1.1 1 for Z , , ~ / Z , , ~ is in poor agreement with the experimentally determined value of 0.48 f 0.10. The NOE effects are very sensitive functions of the relative populations, due to the 1/# dependence in eq 5. Although we have mentioned that the M D simulations were not of sufficient length to give rotamer populations with any statistical reliability, we can briefly contrast the MD results to those from the Boltzmann population analysis. The last two columns in Table I11 give the approximate proportion of the simulation run spent in the perp-x2 and anti-x2 orientations for a given xI rotamer. The molecule appears to have spent more time in the energetically more stable xZconformer than is predicted by the Boltzmann population analysis (Table IV): the ratios of perpx2/anti-x2 for g+, t , and g- are 9,0.4,and 4 and 1.6,0.64, and 1.3 for simulation and Boltzmann population results, respectively. This may mean that the other degrees of freedom in the intramolecular potential are important in establishing the stability (and hence populations) of various rotamers. Longer MD simulations would be needed to ascertain whether this Boltzmann analysis of the ( x l , x2) energy surface adequately reflects the actual rotameric populations and what degree of approximation is made by reducing the dimensionality of the intramolecular potential. Nevertheless, the crude results from the MD simulations also do not agree, even qualitatively, with the results interpreted from lanthanide shift N M R experiments in ref 33. Thus, the details of this CHARMm intramolecular potential, while in agreement with a xi rotamer model (i.e., assuming rapid x2 rotation) of tryptophan conformational dynamics used by Dezube et al.33 to interpret their N M R data, do not lead to agreement with their experimentally determined relative rotamer populations. Anderson et performed calculations on tryptophan, using modified ECEPP potentials, and computed Boltzmann populations for x1rotamers. These results, shown in Table IV, are in qualitative agreement with ours, despite the use of a different molecular mechanics potential, a much lower dielectric constant (to = 2), and many fewer conformations contributing to Pg, P,,and Pr (eight in total). Both refs 29 and 33 note that the structure of free tryptophan in crystals is predominantly in the perp-x2 g+ conformation, as found by computations here and by Anderson
c$
Gordon et al. et al.29 However, in crystals of peptides and proteins, the gconformer of tryptophan is favored,34 as was found from the solution N M R results.33 It may be that the parametrizations of empirical potential energy functions such as CHARMm and ECEPP are not sufficiently flexible to span the range of environments in which tryptophan is found. These results have important implications for rotamer models of tryptophan photophysics. In the rotamer models, multiexponential decay kinetics are proposed to result from the presence of a number of ground-state rotamers, some of which do not interconvert on the fluorescence timescale (typically 3-5 ns). Individual rotamers are assumed to exhibit monoexponential decay kinetics. The correlation of the fluorescence decay components with individual rotamers is largely ~ p e c u l a t i v e . l ~Tryptophan -~~ zwitterion, for example, shows biexponential fluorescence decay kineticsi6and probably has six rotamers in the ground state (vide supra). As pointed out by Tilstra et al.,24“for six rotamers to give rise to only two lifetimes, some rotamers must have similar lifetimes or must interconvert on the fluorescence time scale”. The rotamer model originally applied to tryptophan zwitterionI6 stressed the importance of x1 (C,-C,) rotamers in determining the decay kinetics and assumed that these rotamers did not interconvert during the excited-state lifetime. An alternative rotamer model was suggested by Engh et a1.,z2in which x2 (C,-C,) rotamers were correlated with the fluorescence decay and xi rotamers were assumed to interconvert rapidly such that they were equilibrated on the fluorescence time scale. The latter model was formulated on the basis of MD simulations, employing a united atom representation of hydrogen, which implied rapid x1interconversion. 22 Our MD simulation results, in which an explicit hydrogen atom representation was used, did not indicate rapid xi interconversion and were therefore consistent with the xl-based rotamer model. However, we have shown that MD simulations on tryptophan can yield results that support either rotamer model depending on whether the explicit or united atom representation of hydrogen is chosen. Hence, it is not possible to use our simulation results to unambiguously prove or disprove the lifetime decay models. Two recent s t ~ d i e s ~have ~ , ) ~found that the use of the explicit hydrogen, rather than the united atom, representation for simulations of tryptophan conformational dynamics gave better agreement with certain experimental parameters. The x2-based rotamer model can be rejected on other grounds. If the xi rotamers are equilibrated on the fluorescence time scale, due to rapid rotation around the C,-C, bond, the chirality at C, would be effectively destroyed. Both the perp-xz and anti-x2 conformations of the indole ring (which has planar symmetry) would therefore experience the same “averaged” side chain, and monoexponential decay kinetics would be observed. MD simulations of tryptophan conformational dynamics reported to date have used a model for tryptophan designed to mimic its ground state. It is not clear what effect the changes in, for example, geometry and electron density (see ref 48) that accompany excitation will have on the rotational dynamics of the side chain. MD simulations of molecules in the excited state have and work is in progress to explore the recently been proposal that dipole moment changes upon excitation induce rotamer interconver~ion.~~*~’
Conclusions The results of this work suggest that the ( x l , x2) potential energy surface of tryptophan is modeled quite differently by the explicit and united atom representations of hydrogen within CHARMm. There is disagreement as to the description of the relative barrier heights separating the six local minima representing different (xl, x2) rotamers. As a consequence, the results of MD simulations suggest that xz conformer interconversions occur more rapidly than xi interconversions for an explicit hydrogen representation of a single, ground-state tryptophan molecule. This was ~
~~~~
(48) Callis, P. R. J . Chem. Phys., in press. (49) Levy,R. M.; Kitchen, D. B.; Blair, J. T.; Krogh-Je-spersen, K. J . Phys. Chem. 1990, 94, 4470-4476.
J. Phys. Chem. 1992,96,1921-1932
the case whether continuum or molecular descriptions of the solvent were used. However, MD simulations of tryptophan using a united atom representation of hydrogen suggest the exact opposite, that x, interconversions are more frequent than x2 interconversions. The united atom representation of hydrogen was introduced into packages such as CHARMm and GROMOS because it significantly reduces the computational requirements of most problems pertaining to biomacromolecules. This model has been shown to present a satisfactory representation of certain properties such as structural minima. But empirical potentials are parametrized to reproduce properties near local/global minima; there is no guarantee that they will reflect the true nature of the potential surface far from local minima. Thus, there can be relatively good agreement between the explicit and the united atom descriptions of hydrogen, viz. the (x,, x2) geometries of the six rotamers of tryptophan, yet disagreement in regions separating these minima. The nature of the potential surface separating minima strongly influences trajectories of molecular dynamics simulations. Yet it is only with the advent of greater computer resources that some discrepancies, such as that illustrated here, may be exposed by sufficiently long simulations. Simulation times for proteins typically are of the order of tens of picoseconds, much shorter than the 1-ns simulations on the single tryptophan molecule described here. Even so, our simulations were not long enough to establish,
1921
with any statistical reliability, transition rates between rotamers on our potential energy surface. It therefore behooves users of M D simulations to carefully examine their results to establish consistency with experimental data whenever available. Otherwise, completely wrong interpretation and microscopic modeling of the experiments could ensue. For example, in this work, computational results using the CHARMm intramolecular potential (and full hydrogen representation) for tryptophan agree with the xi rotamer model used to interpret some N M R data.33 However, as found p r e v i o u ~ l y ,there ~ ~ is not even qualitative agreement between computational and experimental results concerning the predicted relative rotamer populations. Potential energy surfaces and M D simulations can offer experimentalists a rationale upon which to extend interpretations of their data. For example, it is conventional to interpret N M R coupling constants as arising from rotamers having one fixed geometry, as was done in ref 29. Figure 4 emphasizes that contributions to such experimentally measurable quantities actually come from a distribution of conformer geometries. Moreover, these distributions may not be narrow or symmetric (as is the case for perp-x2 g and anti-x2 t rotamers). Therefore, theoretical calculations can provide a framework upon which data analyses are developed that do incorporate a treatment of the librational motion that contributes to NMR relaxation times and fluorescence anisotropies.
Simulations on the Counterion and Solvent Dlstribution Functions around Two Simple Models of a Polyelectrolyte Heather L. Gordont and Saul Goldman* Department of Chemistry and Biochemistry, and the Guelph- Waterloo Center for Graduate Work in Chemistry, University of Guelph, Guelph, Ontario, Canada, Nl G 2 WI (Received: August 12, 1991)
We present our results from a series of long canonical ensemble Monte Carlo simulations on the counterion and solvent distribution functions in regions close to a high charge density negatively charged polyelectrolyte. The simulations are all for the "salt free" solution; i.e., a neutralizing number of counterions are present, but there is no additional supporting electrolyte. We considered two models for the solvent-SPC water and a dielectric continuum whose static dielectric constant equals that of SPC water-and two models for the polyion-a uniform line of charge along the central axis of a soft rigid cylinder and a helical necklace of charges spiraling around the surface of a soft cylinder. These models were used to carry out three sets of simulations: continuum solvent with axially charged polyion, discrete solvent with axially charged polyion, and discrete solvent with polyion containing the helical necklace of negative charges. Our results with the continuum solvent were consistent with those from a number of theoretical and simulation studies done in the past with similar models. This model produced a high narrow peak for the counterion concentration close to the surface of the polyion, which decayed rapidly with distance from the polyion. On the other hand, our simulations using (molecular) SPC water with a uniformly charged polyion p r o d u d a surprise. This time the counterions avoided regions near the polyion surface. This result was attributed to an unphysically high degree of solvent orientational polarization in this system. The highly polarized solvent particles successfully displaced all the counterions from regions near the polyion. This unphysical degree of solvent polarization turns out to be due to the use of a uniform line of charge to represent the charge distribution of the polyion. This was proven by our final set of simulations, where the charge distribution of the polyion was taken to be the spiraling helical necklace of discrete charges. In these runs the counterions were found to cluster near the polyion surface, and the extent of orientational polarization of the solvent was very much reduced from what was found using the axially charged polyion.
1. Introduction A polyelectrolyte is a macromolecule with regularly spaced ionizable groups. The purpose of this study is to learn about how a rigid, linear polyelectrolyte, in a water-like solvent, influences the distribution of solvent molecules and counterions near its surface. As with previous work on polyelectrolyte^,'-^^ the incentive to do this stems partly from the problem itself, which is 'Present address: Bioinformatics Section, Institute For Biological Sciences, Building M54,Montreal Road Campus, National Research Council of Canada, Ottawa, Ontario, Canada, K1A OR6.
0022-3654/92/2096- 1921$03.00/0
intrinsically interesting, and partly from the DNA connection. Double-stranded DNA in water is a high charge-density polyelectrolyte, whose conformational, thermodynamic, and binding properties are influenced by the presence and distribution of dissolved electrolyte. This is a large field to which a number of groups have made contributions. To put our own interest in context we provide a brief outline of the history of previous work on this problem. The emphasis and direction taken by previous workers is a reflection both of their scientific taste, and the period during which the work was done. The earliest theoretical studies of which we 0 1992 American Chemical Society
