J. Phys. Chem. 1996, 100, 2637-2645
2637
Motion and Conformation of Side Chains in Peptides. A Comparison of 2D Umbrella-Sampling Molecular Dynamics and NMR Results Thomas C. Beutler, Tobias Bremi, Richard R. Ernst, and Wilfred F. van Gunsteren* Laboratorium fu¨ r Physikalische Chemie, Eidgeno¨ ssische Technische Hochschule, 8092 Zu¨ rich, Switzerland ReceiVed: June 20, 1995; In Final Form: August 28, 1995X
An extension of the technique of umbrella sampling to two dimensions is presented and applied to the study of the kinetic and thermodynamic properties of side chains in peptides. Four two-dimensional potential-ofmean-force surfaces for the two degrees of freedom χ1 and χ2 of the four phenylalanine residues in antamanide at 300 K are calculated using the stochastic dynamics simulation method with the molecular force field GROMOS. From these surfaces, time constants for the transitions between the favored conformations along χ1 and χ2 are calculated using transition state theory. The results are compared with experimental values determined by NMR. In addition, motionally averaged 3J-coupling constants are calculated from the rotamer probability distributions derived from the χ1, χ2 potential-of-mean-force surfaces using the generalized Karplus equations. The average values are compared to the experimental NMR results. It is shown that the twodimensional umbrella-sampling approach is capable of providing detailed information on the motional behavior of side chains and can be used to interpret and complement the experimental findings. Furthermore, by comparing the simulation results to results obtained with energy minimization, the importance of sampling a representative set of conformations is demonstrated. Although the GROMOS force field and the stochastic treatment of solvent effects produce transition rates of the correct order of magnitude, the experimental differences between the rate constants of the four residues are not reproduced. The root-mean-square difference between simulated and experimental averaged 3J-coupling constants ranges between 1.92 Hz (Phe10) and 3.22 Hz (Phe5).
1. Introduction It is well-known that proteins and peptides are not rigid.1 The intramolecular motion of proteins and peptides has been the subject of many theoretical and experimental studies within the past few years,2-4 recognizing that motion is of importance for their biological function, such as ligand binding or enzymatic catalysis.5-7 It is likely that not only the large-scale backbone rearrangements are relevant but also the librational motions and conformational transitions of the side chains, on which this article focuses. Side chains buried inside a protein mostly show conformational rigidity or very slow dynamics unless located in regions such as active-site pockets with less steric hindrance. The conformation of these interior side chains is relevant for the tertiary structure. Their motion is restricted to small-angle libration within the favored potential energy well. Rotamer libraries have been developed using data from X-ray crystallographic studies to predict the most probable side-chain conformations.8,9 Side chains located on the surface of a protein, on the other hand, often exhibit dynamics between several conformations. The resulting multimodal distributions of the dihedral angles may influence the shape and the distribution of polar and apolar regions and thus affect the interactions with other molecules, e.g. with the solvent. With the advent of NMR techniques that allow one to measure homo- and heteronuclear multiple-bond J-coupling constants, which are directly related to dihedral-angle degrees of freedom, and with the possibility to measure relaxation parameters of side chains, much can be learned about the probability distribuX
Abstract published in AdVance ACS Abstracts, October 15, 1995.
0022-3654/96/20100-2637$12.00/0
tions of the dihedral angles and their dynamical behavior. Since the measurement only delivers time-averaged values of Jcoupling and relaxation parameters, a motional model is required to interpret the experimental results at the atomic level. Such a model can be provided by molecular dynamics (MD) or stochastic dynamics (SD) simulations based on a particular atomic force field. The comparison of simulated with experimental averaged values can be used to (in)validate the motional model or force field. Molecular motions occur on various time scales. Smallamplitude librations within one conformational energy minimum are in the range of femto- to picoseconds. The transitions between different conformations usually occur in the range of pico- to microseconds. Many of the slower interwell motions are still beyond the scope of atomic level molecular dynamics simulations. To counter this problem, time-saving techniques, such as e.g. umbrella sampling,10 have been developed.11 In this paper, an improved umbrella-sampling simulation method is presented to investigate the conformational and dynamic properties of side chains with two degrees of freedom. By combining one-dimensional10 and two-dimensional umbrellasampling techniques, a free energy (potential-of-mean-force12) map for the side chains, depending on two degrees of freedom, was calculated and used for comparison to experiment. The molecular dynamics simulations were carried out applying the stochastic dynamics (SD) method,13 where the influence of the nonpolar solvent molecules, in this case chloroform, was taken into account via the frictional and the stochastic force terms of the Langevin equation of motion.14 This technique has been shown to be a good approximation for nonpolar solvent effects.15 Being far less time-consuming than an explicit treatment of solvent degrees of freedom, it is attractive for the extensive calculations required. © 1996 American Chemical Society
2638 J. Phys. Chem., Vol. 100, No. 7, 1996
Beutler et al.
As a test system, the cyclic decapeptide antamanide (Figure 1) has been chosen, as its side-chain conformations and dynamics have been thoroughly investigated.16-24 It contains four inequivalent phenylalanine residues with different motional properties regarding the rotations about the CRCβ and CβCγ bonds with the rotation angles χ1 and χ2. On the basis of the simulated (χ1, χ2) potentials of mean force for the four phenylalanine side chains and using the generalized Karplus equations, vicinal J-coupling constants are computed and compared to experimentally determined values. The rotamer populations extracted from the simulations are also compared to those derived from measurements using a different motional model.16,17 Additionally, the simulation results were used to determine transition probabilities for the χ1 and χ2 conformational transitions by transition-state theory and were compared to rate constants derived from experimental 13C relaxation data.24 For comparison with the stochastic dynamics simulations, energy-minimization techniques were applied to determine the surfaces of minimal potential energy for selected χ1 and χ2 values (adiabatic potential energy surfaces25) of the phenylalanine side chains. These surfaces were also used to estimate rotamer populations and energy barriers. The quality of this computationally less demanding approach is discussed by comparison to the stochastic molecular dynamics results. The essential physical quantities and model assumptions are briefly described in section 2. In section 3, the computational methodology is outlined, while the details of the twodimensional umbrella-sampling technique are given in the Appendix. The results are discussed in section 4. Section 5 summarizes the results. 2. System and Theory The side-chain conformations of the phenylalanine residues depend on the two dihedral angles χ1 and χ2, defined in Figure 1b and c. The three-bond J-coupling constants, measured by NMR, depend on either χ1 or χ2, and the observed average values 〈3Jkl〉 are determined by the one-dimensional probability densities F1(χ1) or F2(χ2), respectively. The one-dimensional probability densities are projections of the two-dimensional probability density of the side-chain conformations F(χ1,χ2):
F1(χ1) dχ1 ) (∫0 F(χ1,χ2) dχ2) dχ1 2π
(1)
and similarly for F2(χ2). In stochastic dynamics (SD) simulations, the integration of the Langevin equations of motion yields a trajectory from which F(χ1,χ2) can be directly estimated, provided that the conformational states are not separated by high-energy barriers. In the presence of high barriers (greater than a few kBT), the conformational transitions become too infrequent to be statistically meaningful within the computationally feasible time frame of dynamics simulations. However, with umbrella-sampling techniques, it is possible to determine the potential of mean force also for degrees of freedom plagued by high energy barriers.26-28 Thus, umbrella sampling allows one to calculate the populations of different rotamer states and the heights of the free energy barriers separating them. Below, we use a generalization to two dimensions of the standard one-dimensional umbrellasampling technique. Since the generalization mainly involves technical issues, it is discussed in the Appendix. For the discussion of the side-chain dynamics reflected in the NMR relaxation measurements, it is desirable to have available the potential of mean force, i.e. the free energy hypersurface w(χ1,χ2) and the corresponding one-dimensional profiles w1(χ1) and w2(χ2). They can be obtained from F(χ1,χ2)
a
b
c
Figure 1. (a) Primary structure of antamanide. The phenylalanine residues are marked in bold face. (b) Geometry of the phenylalanine residue with stereospecific assignments and with the designations of the angles used in the text. (c) Newman projections of the three rotamers I, II, and III of phenylalanine.
and Fi(χi) by inversion of the following two relations:
F(χ1,χ2) dχ1 dχ2 ) exp[-w(χ1,χ2)/(kBT)] dχ1 dχ2
∫0 ∫02π exp[-w(χ′1,χ′2)/(kBT)] dχ′1 dχ′2 2π
Fi(χi) dχi )
exp[-wi(χi)/(kBT)] dχi
∫0
2π
(2)
(3)
exp[-wi(χ′i)/(kBT)] dχ′i
Here T is the absolute temperature and kB the Boltzmann constant. A zero Kelvin approximation of the potential-of-mean-force surface w(χ1,χ2) is given by the surface-of-minimum-energy configurations at fixed χ1 and χ2, umin(χ1,χ2). This adiabatic potential energy surface can be obtained at low computational costs directly by minimizing the potential energy with respect to all molecular degrees of freedom for fixed χ1 and χ2 values. It is of interest to compare molecular parameters deduced from umin(χ1,χ2) with those computed from w(χ1,χ2) to decide whether the efforts of calculating w(χ1,χ2) are worthwhile. The corresponding probability density Fu(χ1,χ2) is obtained from an expression equivalent to eq 2, replacing w(χ1,χ2) by umin(χ1,χ2). From the one-dimensional probability densities Fu,i(χi), obtained by projection in analogy to eq 1, it is then possible to compute the one-dimensional profiles of minimum energy umin,i(χi) by inversion of an equation analogous to eq 3 or, combining the
Motion and Conformation of Side Chains in Peptides
J. Phys. Chem., Vol. 100, No. 7, 1996 2639
three steps, by
umin,1(χ1) )
[2π1 ∫
-kBT ln
2π
0
]
exp[-umin(χ1,χ2)/(kBT)] dχ2 (4)
and analogously for umin,2(χ2) by integration over χ1. Following the original work by Karplus,29,30 the dependence of a homonuclear or heteronuclear coupling constant 3Jkl(χi) between the vicinal nuclei k and l is given by the relation29-31 3
Jkl(χi) ) A cos2(χi - χkl) + B cos(χi - χkl) + C
(5)
A, B, and C are empirical coefficients, specific for the considered molecular fragment and the considered type of nuclei k and l. The angle χkl denotes a phase shift such that the argument of the trigonometric functions is the dihedral angle between the vicinal nuclei k and l. Since the inverse function of eq 5 is multiple-valued, it is necessary to measure several 3Jkl-coupling constants between different nuclei k and l to uniquely determine χi even in the case that only one single rotamer χi is populated. Moreover, the approximate nature of eq 5 may further complicate the derivation of conformational probabilities from measured 3Jkl values. For rapid conformational dynamics, the average value of the coupling constant,
〈3Jkl〉 ) ∫0 Fi(χi) 3Jkl(χi) dχi 2π
(6)
can be calculated with the help of the one-dimensional probability densities Fi(χi) extracted from simulations. They can be compared to J-couplings derived from simplified motional models16,17,32 and to the experimental ones. In order to obtain dynamical information from simulations in the presence of high energy barriers, transition state theory (and higher order corrections) can be applied.33-35 The rate constant ki,nfm of a transition from a conformational state n over a barrier, located at χqi , to a neighboring state m can be written as
ki,nfm )
F(χiq) 1 ) κi,nfm〈|χ˙ i|〉 τi,nfm 2 ∫nF(χi) dχi 1
(7)
where τi,nfm denotes the time constant of the process. The index n indicates that the integration extends over the region of χi belonging to the state n. 〈|χ˘ i|〉 stands for the average absolute value of the velocity of crossing the barrier evaluated at the transition state coordinate value χqi . κi,nfm is the transmission coefficient which takes the recrossing rate into account and is a higher order correction of transition state theory. If κi,nfm is set equal to unity, eq 7 gives the standard transition state theory rate constant, which will be used in this article for the analysis of transitions along the side-chain degrees of freedom χ1 and χ2. The motions along χ1 and χ2 are of different complexity. The potential energy landscape is periodic in χ2 with a period of 180°, and there is usually one single well per period. Accordingly, the determination of the kinetic properties is straightforward. By contrast, along χ1, there can be three or even more different potential energy wells and the system is periodic with a period of 360°. As there are five independent rate constants in the case of three stable rotamers about the CRCβ bond (I, II, and III; see Figure 1c), the motional behavior is more complex. fast , is If the time constant of an intramolecular process, τint fast within the window τc/10 e τint e 10τc, centered at the time constant τc of the overall tumbling of the molecule in the solvent,
TABLE 1: Average Absolute Velocities 〈|χ3 i|〉 for the Two Degrees of Freedom χ1 and χ2 at the Various Transition States As Obtained from 50 ps Simulations with Harmonic Biasing Potential Energy Functions (eq A-1, k1 ) k2 ) 150 kJ mol-1 rad-2), Centered at the Transition States Located at the Saddle Points of the Surfaces of Minimum Energy umin(χ1,χ2) (See Figure 2a) I T II Phe5 Phe6 Phe9 Phe10
323 365 303 401
〈|χ˘ 1|〉 (deg/ps) I T III II T III 378 354 376 357
367 421 332 294
〈|χ˘ 2|〉 (deg/ps) ITI 526 456 516 531
the relaxation rate constants observed by NMR are sensitive to fast τfast int , which allows for determination of τint . Sometimes also the type of the motion can be determined by comparing the measured relaxation rate constants to simulated ones based on various motional models. The side-chain dynamics of the phenylalanine residues of antamanide have been shown to involve librational motions within this observation window for the actual experimental conditions.24 In addition to these fast librational motions, some residues exhibit interwell dynamics on a slower time scale not affecting NMR T1 relaxation. These dynamics with time constants in the range of micro- to milliseconds are reflected in motionally averaged 3Jkl-coupling constants which do not allow the determination of the time slow but yield information about the rotational probconstant τint ability distribution Fi(χi). 3. Computational Methods and Simulations The stochastic dynamics simulations were performed using the GROMOS program library and the GROMOS force field.36 The friction coefficients γi for the solute atoms were set to be
γi ) ωiγ
(8)
with ωi being a weight factor, chosen to approximate the relative surface area of each solute atom that is accessible to the solvent; γ denotes the solvent friction coefficient, which was set to 19 ps-1. The function representing the accessible-surface-area weight factor ωi was chosen as follows.15 It is zero in the presence of six neighboring atoms within a distance of 0.3 nm of atom i and increases linearly for a decreasing number of neighbors until it reaches the value 1 in the absence of any neighbor within that distance. The bond lengths were kept fixed within a relative tolerance of 10-4 using the SHAKE method.37 The simulations were started from the configuration obtained after a 500 ps simulation at 300 K with identical parameters during which a stable, low-energy backbone conformation was reached.22 The efficiency and quality of the umbrella simulations depend on (i) the choice of the umbrella potential energy functions and (ii) the start configuration used for the simulations. These are discussed in the Appendix together with the construction of the minimum energy surface umin(χ1, χ2). The average absolute velocities on the saddle points, 〈|χ˘ i|〉, depend on the effective masses of the corresponding degrees of freedom. They are determined by the kinetic energy terms of the Hamiltonian but are independent of the potential energy terms. For this reason, they can be obtained from simulations with a biasing potential energy term that focuses sampling on the saddle point and thus allows one to obtain converged results. In the simulations, a biasing potential energy function of parabolic form (eq A-1) was applied with force constants k1 ) k2 ) 150 kJ mol-1 rad-2 and (χi,χj) corresponding to the saddle points of the surface umin(χ1, χ2). The resulting average
Figure 2. (a, left) Surfaces of minimum potential energy, umin(χ1,χ2) (adiabatic potential energy surfaces, after energy minimization for fixed χ1 and χ2), of the four phenylalanine residues in antamanide. Contour levels are drawn at 2, 4, 6, 8, 10, 15, 20, 30, 40, 50, 75, and 100 kJ/mol. (b, right) Potential-of-mean-force surfaces w(χ1,χ2) for the phenylalanine residues in antamanide at 300 K. Contour levels are drawn at 2, 4, 6, 8, 10, 15, 20, 30, 40, and 50 kJ/mol. Regions where the sampling was unreliable are shaded.
2640 J. Phys. Chem., Vol. 100, No. 7, 1996 Beutler et al.
Motion and Conformation of Side Chains in Peptides
Figure 3. (a, top) Potential energy profiles along χ1, umin,1(χ1), as obtained by weighted averaging of the surfaces of minimum energy, umin(χ1,χ2), of the four phenylalanine residues in antamanide at 300 K. The offsets of the potentials were chosen such that all minima near χ1 ) -60° are at zero energy. (b, bottom) Potentials of mean force along χ1, w1(χ1), of the four phenylalanine residues in antamanide at 300 K. The offsets of the potentials were chosen such that the minima near χ1 ) -60° are at zero potential.
velocities are all similar (see Table 1), and thus, for the calculation of the rate constants of the four phenylalanine residues, identical average values were chosen for all transitions along χ1, 〈|χ˘ i|〉 ) 350 deg/ps, and for the transitions along χ2, 〈|χ˘ 1|〉 ) 500 deg/ps. 4. Results and Discussion Potential-of-Mean-Force and Minimum Energy Surfaces. In Figure 2a the surfaces of minimum potential energy, umin(χ1,χ2), are shown. Because of the mentioned symmetry of the phenyl ring, the surfaces are only given in the range from 0° to 180° along χ2. In Figure 2b the potential-of-mean-force surfaces w(χ1,χ2) are shown. The regions where reliable sampling is inhibited by the shape of the underlying energy landscape are shaded. Appropriate sampling was obtained in all the wells of the potential of mean force and at the low saddle points. Significant differences are apparent between the potentialof-mean-force surfaces w(χ1,χ2) and the surfaces of minimum potential energy umin(χ1,χ2). In particular, it can be seen that some fine-structure features in umin(χ1,χ2) are not reproduced in w(χ1,χ2). Also, the shapes and the relative depths of the energy wells differ, obvious particularly for Phe9, demonstrating the importance of sampling a representative set of configurations in computing w(χ1,χ2). In Figures 3b and 4b, the one-dimensional potentials of mean force for the phenylalanine degrees of freedom χ1 and χ2 as determined from umbrella sampling are shown. The transition states at χ1 ≈ 120° lie significantly higher than the other χ1 transition states (difference > kBT), as is also the case for the surfaces of minimum energy shown in Figures 3a and 4a. The isomerization dynamics will be dominated by the motion over
J. Phys. Chem., Vol. 100, No. 7, 1996 2641
Figure 4. (a, top) Potential energy profiles along χ2, umin,2(χ2), as obtained from the surfaces of minimum energy, umin(χ1,χ2), of the four phenylalanine residues in antamanide at 300 K. The offsets of the potentials were chosen such that the minima near χ1 ) 90° are at zero energy. (b, bottom) Potentials of mean force along χ2, w2(χ2), of the phenylalanine residues in antamanide at 300 K. The offsets of the potentials were chosen such that the minima near χ2 ) 90° are at zero potential.
the lower barriers. Accordingly, the motion over the highest barrier χ1 ≈ 120° can be neglected, as will be done in the following. The importance of sampling can again be seen by comparing the one-dimensional potentials of mean force w1(χ1) and w2(χ2) (Figures 3b and 4b) with the corresponding onedimensional minimum energy profiles umin,1(χ1) and umin,2(χ2) (Figures 3a and 4a). In particular, the minimum energy profiles umin,2(χ2) seem to contain minima irrelevant for higher temperatures at χ2 ≈ 0° and for Phe9 also at χ2 ≈ 80°. The observed deviations between umin and w are on the order of kBT and larger and thus have a strong impact on quantities derived from simulations that can be compared to experimental results. The deviations show that the quality of molecular force fields for the simulation of biomolecules can only be judged if experimental results are compared to simulation results that take dynamic effects properly into account. Nevertheless, the surfaces of minimum energy approximate qualitatively the potential-of-mean-force surfaces and may be useful for crude qualitative arguments. As expected, the χ1 potentials of mean force of the individual phenylalanine residues vary significantly among each other. Their general qualitative behavior is determined by the dihedralangle potential energy due to steric repulsions. The highest barrier is found between χ1 ) 60° and χ1 ) 180°, and the lowest, between χ1 ) -60° and χ1 ) 180°. Nevertheless, the cyclic peptide gives each phenylalanine a unique environment that influences its properties considerably. Rotamer Distribution along χ1. The rotamer distribution along χ1, obtained from molecular dynamics simulation, is given by F1(χ1), the one obtained by energy minimization by Fu,1(χ1). The two rotamer distribution functions are used to compute the dynamically averaged coupling constants 〈3Jkl〉 according to eq 6 for the five J-couplings 3JR,β2, 3JR,β3, 3JC′,β2, 3JC′,β3, and 3JR,Cγ,
2642 J. Phys. Chem., Vol. 100, No. 7, 1996
Figure 5. Comparison of the motionally averaged 3J values of antamanide at 300 K: (1) calculated from the potential of mean force w(χ1); (2) calculated from the potential energy profile umin(χ1); (3) experimental values according to refs 17 and 38-40; (4) fit of experimental values by three-state model of ref 16; (5) fit of experimental values by three-state model of ref 17. The parameters A, B, and C for the generalized Karplus equation were taken from eq 16 in ref 31 for 3J(HRHβ), from ref 42 for 3J(C′Hβ), and from ref 43 for 3J(H C ). R γ
which all depend on the dihedral angle χ1. The two sets of results are represented in graphical form in Figure 5 together with the experimental 3Jkl values and two back-calculated sets from a simple three-state analysis given in refs 16 and 17. In these studies, sets of 3Jkl-coupling constants depending on the same χ1 angle were used to determine the populations pi of the three discrete rotamers I, II, and III,32 assuming a jump model with rapid interconversions between the rotamers leading to
〈3Jkl〉 ) pI3Jkl(χ1)-60°) + pII3Jkl(χ1)180°) + pIII3Jkl(χ1)60°) (9) Root-mean-square deviations (RMSD’s) between experiment and MD calculation for the five J couplings in each of the four residues are given in Table 2. Considering the fact that the experimental coupling constants have errors of less than (1 Hz, the agreement between the experimental and the calculated values is less perfect than expected on the basis of the computational accuracy of the potentials of mean force which should be better than kBT. While the best rationalization in terms of an optimally adapted threestate model by Schmidt17 gives RMSD’s of less than 1.06 Hz, the RMSD’s for the MD simulations based on w1(χ1) reach values up to 3.22 Hz and those based on umin,1(χ1), even up to 6.22 Hz. Especially for Phe9 the agreement is rather poor for 3JR,β2, 3JC′,β3, and 3JR,Cγ, while for Phe6 and Phe10 the agreement is reasonable. Indeed, the reproduction of experimental 〈3Jkl〉 values is a sensitive test for the force field and the representation of the solvent interactions used in the simulation. To obtain better insight into the deviations between the threestate model data evaluation and the simulated potentials of mean force used here, the populations of the three minimum energy conformations were also calculated from umin,1(χ1) and w1(χ1) and are given in Table 2 together with the populations of the three-state model analyses. The agreement is good for Phe10
Beutler et al. for the potential of mean force as well as for the potential energy profile. Rotamer I is preferentially populated. Rather good agreement is also found for Phe6, where NMR and MD simulation both predict a dynamic equilibrium with all three conformations contributing to the observed 〈3Jkl〉 values. For Phe5 and Phe9, the populations obtained from the simulations agree less with the populations obtained by fitting the threestate model to the experimental data. While for Phe5 both the potential of mean force and the potential energy profile still predict a preferred rotamer I and a small population of II, for Phe9 a weighting of rotamer III with respect to the rotamers I and II is found that disagrees with the three-state interpretation of the experimental data,16,17 while the relative weighting of the latter two is correctly predicted. Dynamics of the χ1 Rotation. Using transition-state theory, the various time constants (τ1,nfm ) 1/k1,nfm) for the interconversions along χ1 between the conformations m, n ) I, II, and III (see Figure 1c) through the relevant transition states were calculated and are given in Table 3. When these values are compared with the correlation times τ1 determined by NMR relaxation measurements (τ1 ) 180, 150, 200, and 1550 ps for Phe5, Phe6, Phe9, and Phe10, respectively),24 it must be taken into account that the NMR values have been interpreted in terms of a dominant contribution from intrawell dynamics within a segment characterized by a restriction angle 2χmax 1 . The computational results allow conclusions to be drawn on the possible contribution of interwell dynamics to the rate constants observed by NMR. For Phe5 and Phe10, the simulations indicate an almost exclusive population of rotamer I. The dynamical processes between state I and the scarcely populated state II are very rapid (τ1,IfII e 125 ps) while the transfer to state III is much slower (τ1,IfIII e 5 µs). However, due to the low population of states II and III, interwell dynamics hardly contribute to the NMR relaxation, and the measured NMR correlation times can be safely attributed to intrawell librational motion. The small ) 42° for Phe10 is in agreement with restriction angle 2χmax 1 this conclusion while for the explanation of the larger value ) 64° of Phe5, a rather wide single-well potential is 2χmax 1 required which is not clearly reflected in the profiles of Figure 3. In the case of Phe6, all three rotamers are appreciably populated according to Table 2. While the time constants involving state III are in microseconds on the basis of the umin,1(χ1) potential, they are in nanoseconds on the basis of the more reliable w1(χ1) potential. The fastest processes occur among states I and II. All these interwell processes can appreciably contribute to the observed NMR relaxation. It thus agrees well with the interpretation in terms of exclusive intrawell ) motion that leads to a rather large restriction angle 2χmax 1 80°. For Phe9, the potential of mean force, w1(χ1), and the potential energy profile, umin,1(χ1), differ significantly. From the potential of mean force, a behavior similar to the one of Phe5 and Phe10 is expected with little contribution of interwell dynamics to relaxation, but here the preferentially populated state is rotamer III. According to the potential energy profile, Phe9 should behave more like Phe6 with all three states appreciably populated. Only the transitions I T II are rapid and can conceivably contribute to relaxation. However the processes are so rapid (τITII ) 10 ps and τIITI ) 1 ps) that they would mainly lead to an averaging of the intrawell behavior in the two wells. The intermediate restriction angle of 64°, determined by NMR, does not allow an unambiguous conclusion in favor of either model.
Motion and Conformation of Side Chains in Peptides
J. Phys. Chem., Vol. 100, No. 7, 1996 2643
TABLE 2: Populations of the Three Minimum Energy χ1 Rotamers I, II, and III of the Phenylalanine Residues in Antamanide (See Figure 1) and the Root-Mean-Square Deviations (RMSD’s) from the Experimental 3Jkl Constants (3Jr,β2, 3Jr,β3, 3JC′,β2, 3JC′,β3, and 3Jr,Cγ) Determined by Different Methods Phe5 pI (%) umin,1(χ1)a w1(χ1)b Kessler et al.16 c Schmidt17 c
99 99 33 56
Phe6
pII (%)
pIII (%)
RMSD (Hz)
1 1 33 44
0 0 33 0
3.22 3.22 2.04 1.06
pI (%) 72 35 33 30
pII (%)
pIII (%)
RMSD (Hz)
8 7 33 22
20 58 33 48
1.53 2.32 1.36 0.77
pII (%)
pIII (%)
RMSD (Hz)
1 1 13 16
0 0 9 0
1.77 1.92 1.09 1.01
Phe9 pI (%) umin,1(χ1)a w1(χ1)b
55 2 91 88
Kessler et al.16 c Schmidt17 c a
Phe10
pII (%)
pIII (%)
RMSD (Hz)
9 1 14 12
36 97 -5 0
6.22 2.88 0.61 0.55
pI (%) 99 99 78 84
umin,1(χ1): energy minimization. b w1(χ1): umbrella sampling. c Kessler et al.16 and Schmidt17: three-state analysis of NMR data.
TABLE 3: Time Constants τnfm ) (knfm)-1 for the Transitions between Rotamers n, m ) I, II, and III (See Figure 1) in Antamanide at 300 K Determined with Transition State Theory from the Potential Energy Profiles umin,1(χ1) or the Potential of Mean Force w1(χ1) τITII Phe5 50 ps Phe6 nca Phe9 10 ps Phe10 125 ps Phe5 Phe6 Phe9 Phe10
100.3 ps 83.2 ps 37.6 ps 56.4 ps
τIITI
τITIII
τIIITI
τnfm from umin(χ1) 0.3 ps 5 µs 500 ps nc 33 µs 11 µs 1 ps 16 µs 11 µs 1 ps 1 µs 125 ps tnfm from w(χ1) 0.5 ps 4.3 ns 3.3 ps 16.6 ps 4.9 ns 8.2 ns 12.3 ps 20.8 ps 1.0 ns 0.4 ps 13 ns 2.2 ps
τIITIIIb
τIIITIIb
400 µs 150 µs 3 ms 76 ms
9 µs 350 µs 12 ms 2 ms
261.5 ns 2.8 ns 798 ps 128.6 ns
43.3 ns 23.9 ns 118 ns 2.9 ns
a nc (not calculated) indicates that a transition state cannot be defined, as there is no maximum between the two rotamers. b Because of the poor statistics for the potential of mean force at the saddle points between rotamers II and III, these values are only rough estimates.
Dynamics of χ2 Rotation. Transition state theory, based on the profiles umin,2(χ2) and w2(χ2), leads to the time constants of the phenyl ring 180° flips listed in Table 4. The NMR relaxation data in ref 24 have been interpreted in terms of rapid interwell dynamics, represented by a two-site jump model. The computed and the NMR data of Table 4 are therefore directly comparable. The time constants are calculated according to the convention used in NMR: τ2 ) 1/(2k2). The time scale of this motion is reproduced by the values derived from the potentials of mean force; they lie in the range τ2 ) 15-227 ps with an average value of τ2 ) 82.5 ps, compared to τ2 ) 37-194 ps with an average value of τ2 ) 102.5 ps given by NMR. The time constants derived from the potential energy profiles indicate (except for Phe9) a much slower motion (τ2 ) 1314 ps) than measured by NMR, which again shows the importance of sampling a representative set of configurations. While the proper time range of the motion is obtained, the differences between the motional behavior of the four residues are not reproduced by the SD simulations using the GROMOS force field and a stochastic model for solvent effects.
TABLE 4: Time Constants τ2 ) (2k2)-1 for the Phenylalanine Ring Flips in Antamanide at 300 K Using Transition State Theory Combined with the Potential Energy Profiles, umin,2(χ2) or the Potentials of Mean Force, w2(χ2) τ2 ) 1/(2k2) (ps) Phe5 Phe6 Phe9 Phe10 a
from umin,2(χ2)
from w2(χ2)
from NMRa
660 900 125 3570
15 20 227 68
113 37 66 194
The NMR values have been taken from ref 24.
K prove to be a demanding test system. It turns out that the qualitative rotational behavior about χ1 and χ2 is simulated in agreement with the NMR data. The slow interwell motional behavior about χ1, giving rise to the averaged 3J-coupling constants, cannot always be predicted with sufficient accuracy, since it is highly sensitive to even small differences in the relative depths of the potential energy wells, which seem to be beyond the accuracy of the force field and the stochastic treatment of solvent effects. Nevertheless, the study shows that some residues (Phe5 and Phe10) are well modeled. Even on the slow time scale their simulated behavior is consistent with the experimental results. Furthermore, it is demonstrated that the quality of the molecular force field and solvent effects can only be assessed if in the simulations a Boltzmann-weighted set of configurations is sampled. The energy minimization approach is too crude, giving significantly different results. In particular, different populations of the side-chain rotamers are obtained and the fine structure in the potential hypersurface, which can be observed in the energy minimization results, disappears when averaging over sampled molecular configurations. Thus, entropic effects must not be neglected in the refinement of parameters of a molecular force field. The present study concentrated on simulations of a peptide in an apolar solvent. In polar solvents, where entropic effects are even more pronounced, representative conformational sampling and correspondingly long simulation periods are therefore a must.
5. Conclusions The combined calculation of one- and two-dimensional potentials of mean force for side chains by means of umbrellasampling techniques can be utilized to gain insight into the conformational preferences and the motional behavior along two degrees of freedom in amino-acid side chains in proteins and peptides. The four phenylalanine residues in antamanide at 300
Acknowledgment. This research has been supported by the Swiss National Science Foundation. We would like to thank Dr. Matthias Ernst for many stimulating discussions and Dr. Roger M. Brunne for making trajectories available to us. We also acknowledge the generous financial support of Zeneca Ltd., U.K.
2644 J. Phys. Chem., Vol. 100, No. 7, 1996
Beutler et al.
6. Appendix Two-Dimensional Umbrella-Sampling Procedure. The force fields used in molecular simulations define a complex energy hypersurface. In general, the hypersurface contains many local energy minima and energy wells which are separated from each other by energy barriers. For equilibrium properties, energy wells contribute if their energy does not lie higher than a few kBT above the global energy minimum. In a molecular simulation it is possible that the system is trapped in an energy well due to high barriers surrounding it; i.e., it is possible that the system cannot relax to other probable configurations within the available simulation period. In such simulations, not all of the configurational space is sampled that is representative for the system at the temperature of interest. By applying biasing potential energy functions during the simulations, the system can be forced to overcome barriers and thus it is possible to explore a wider configurational space within the available simulation period. Accordingly, biasing potential energy functions (so-called umbrella potential energy functions) in combination with energy minimization techniques and molecular dynamics can be used to obtain configurations within all relevant energy wells. The choice of the (umbrella) biasing potential energy functions is crucial for achieving representative sampling efficiently. A two-step procedure was applied. In the first step (i) successively two different sets of biasing potential energy functions and energy minimization techniques were used to determine the surfaces of minimum energy at fixed χ1 and χ2, umin(χ1,χ2), for each of the phenylalanine side chains. Subsequently (ii) these surfaces were used to decide upon the umbrella simulation setup. Furthermore, the surfaces were used to assess the quality of the static (energy minimization) approximation by comparison to the potential-of-mean-force surfaces w(χ1,χ2), which are obtained from the computationally more demanding umbrella simulations. (i) Surface of Minimum Energy. To determine the surfaces of minimum energy (adiabatic energy surfaces), umin(χ1,χ2), the following protocol was applied. To begin with, a grid of 30° spacing was defined for the χ1, χ2 space. Harmonic biasing potential energy functions
1 1 Ur,s(χ1,χ2) ) k1(χ1 - χr)2 + k2(χ2 - χs)2 2 2
(A-1)
at first with force constants k1 ) k2 ) 75 kJ mol-1 rad-2 were defined. In a series of 20 ps stochastic dynamics simulations, χr and χs were fixed successively at the values of the grid points. The first simulation of the series was started from the end configuration of a previous 500 ps trajectory.22 Subsequently, the end configurations of the previous simulations in the series were used as the start configurations for the next simulations with χr and χs set to the values of neighboring grid points. Next, a denser grid of 5° spacing and more restricting harmonic biasing potential energy functions Ur,s(χ1,χ2) with force constants k1 ) k2 ) 1000 kJ mol-1 rad-2 were defined. In a series of simulations, χr and χs were set successively to the values of the grid points. Every simulation comprised 5 ps simulated annealing followed by 400 steps steepest descent energy minimization. The final χ1 and χ2 angles were determined and the resulting potential energies were corrected for the biasing potential energy contribution. By repeating the simulations starting from different configurations stemming from the simulations with k1 ) k2 ) 75 kJ mol-1 rad-2 (30° spaced grid) and taking the lowest potential energy value as an estimate for umin(χ1,χ2), trapping in local energy minima was avoided.
The large biasing potential energy force constant ensured that angles close to the grid points were obtained. The dense spacing (5°) allowed us to obtain by interpolation the corresponding surfaces given in Figure 2a. (ii) Free Energy Calculations. To obtain the potentials of mean force, a one-dimensional grid of 12° spacing was defined along the χ1 axis. A harmonic biasing potential
1 Ur(χ1) ) k1(χ1 - χr)2 2
(A-2)
was chosen as the umbrella potential energy function. The force constant k1 was set to 150 kJ mol-1 rad-2. In a series of 100 ps simulations, χr was successively set to the values of the grid points. As start configurations, the end configurations from the simulations with the two-dimensional biasing potential energy functions on the 30° spaced grid (eq A-1) with χs ) 90° were taken. The umbrella potential energy functions ensured sampling in a narrow area around the grid points χr. In the area where the system has been sampled, the potential of mean force w1,r(χ1), obtained from the simulated probability distribution F1,r(χ1) by
w1,r(χ1) ) -kBT ln(F1,r(χ1)) + Ur(χ1) + cr
(A-3)
is identical to the potential of mean force w1(χ1) of the system. The cr’s denote (still undetermined) parameters which shift the segments w1,r(χ1) to obtain a continuous function w1(χ1). Since in every simulation of the series sampling is focused on a small region along χ1, a second-order regression of w1,r(χ1) above and below the probability maximum accurately gives the local behavior of w1(χ1). After fitting the w1,r(χ1) segment by segment, the parameters cr were adjusted such that
w1,r(ξ) ) w1,r+1(ξ)
(A-4)
with ξ being the value of χ1 halfway between two subsequent probability maxima. This procedure of joining the segments ensures that the joints are placed in regions of reliable sampling. An alternative would be to use overlapping parts of the sampled distributions, which is, however, less preferable since in these overlapping parts the sampling is restricted by the umbrella potential. The resulting curves are shown in Figure 3b. If in an umbrella simulation only a narrow distribution along χ1 is sampled, the corresponding distribution along χ2, F2,r(χ2), probes directly the potential-of-mean-force surface at fixed χ1, w(χ1)χr,χ2).41 Thus one can write
wr(χ1)χr,χ2) ≈ -kBT ln(F2,r(χ2)) + cr′
(A-5)
The symbol cr′ denotes again parameters that shift the individual functions wr(χ1)χr,χ2) to obtain w(χ1)χr,χ2). The parameters cr′ can be determined with the help of w1(χ1), since one can write
(
w1(χ1) ) -kBT ln
[
] )
w (χ ,χ ) - c ′ 1 2π ∫ exp - r 1 kB2T r dχ2 2π 0
(A-6)
By determining all parameters cr′, the one-dimensional potential of mean force, w1(χ1), can be extended into the second dimension to yield the potential of mean force surface w(χ1,χ2) in the region of the χ1,χ2-space where the system was sampled, that is, in the area of the free energy wells corresponding to the rotameric states and also along the isomerization reaction pathway for transitions between the rotamers along χ1, where the umbrella potential energy function enforces sampling. However, no information is obtained on the regions of the saddle
Motion and Conformation of Side Chains in Peptides points relevant for transitions along χ2. From the surfaces of minimum energy, umin(χ1,χ2) (Figure 2a), one sees that transitions along χ2 can proceed over different saddle points (χ1 ) (60°, χ2 ) 0°/(180°) that are separated by barriers. Thus, in a simulation with restrained motion along χ2, ergodic sampling cannot be expected. However, in two independent series of simulations with harmonic umbrella potential energy functions acting along χ2, it is possible to determine the potentials of mean force for the subspaces defined by 0° < χ1 < 120° and 120° < χ1 < 360°. These potentials of mean force will, in the following, be denoted by w+(χ2) and w-(χ2), respectively. After extension of w+(χ2) and w-(χ2) into their two-dimensional analogs, the two curves can then be merged to yield w2(χ2) by calculating
w2(χ2) )
(
[
]
w+(χ1,χ2) - c+ dχ1 + kBT w-(χ1,χ2) - c2π (4π/3)-1∫2π/3exp dχ1 (A-7) kBT
-kBT ln (2π/3)-1∫0
2π/3
exp -
[
] )
where c+ and c- denote the parameters that calibrate the independently obtained w+(χ1,χ2) and w-(χ1,χ2) to w(χ1,χ2). These parameters can be determined, since, in the region of the wells that correspond to the rotameric states, w(χ1,χ2) is already known from the simulations with the umbrella potential energy function acting along χ1. Figure 4b shows the resulting w2(χ2) for the four phenylalanine residues in antamanide. Figure 2b displays the two-dimensional potentials of mean force w(χ1,χ2) obtained after the merging. References and Notes (1) Frauenfelder, H.; Sligar, S.G.; Wolynes, P.G. Science 1991, 254, 1598. (2) Brooks, C. L., III; Karplus, M.; Pettitt, B. M. Proteins: A theoretical perspectiVe of dynamics, structures, and thermodynamics; John Wiley and Sons: New York, 1987. (3) van Gunsteren, W. F.; Berendsen, H. J. C. Angew. Chem., Int. Ed. Engl. 1990, 29, 992. (4) Tycko, R., Ed. Nuclear magnetic resonance probes of molecular dynamics; Kluwer Academic Publishers: Amsterdam, 1994. (5) Debrunner, P. G.; Frauenfelder, H. Annu. ReV. Phys. Chem. 1982, 33, 283. (6) Janin, J.; Wodak, S. J. Prog. Biophys. Mol. Biol. 1983, 42, 21. (7) Ringe, D.; Petsko, G. Prog. Biophys. Mol. Biol. 1985, 45, 197. (8) Janin, J.; Wodak, S. J.; Levitt, M.; Maigret, B. J. Mol. Biol. 1987, 125, 357. (9) Dunbrack, R. L.; Karplus, M. J. Mol. Biol. 1993, 230, 543. (10) Torrie, G. M.; Valleau, J. P. J. Comput. Phys. 1977, 23, 187.
J. Phys. Chem., Vol. 100, No. 7, 1996 2645 (11) van Gunsteren, W. F. Computer simulation of biomolecular systems: Overview of time-saving techniques. AdVances in Biomolecular Simulations, American Institute of Physics (A.I.P.) Conference Proceedings Obernai (France); Lavery, R., Rivail, J.-L., Smith, J., Eds.; New York, 1991; Vol. 239, p 131. (12) Hill, T. L. Statistical Mechanics; McGraw-Hill: New York, 1956. (13) van Gunsteren, W. F.; Berendsen, H. J. C. Mol. Sim. 1988, 1, 173. (14) McQuarry, D. A. Statistical Mechanics; Harper and Row: New York, 1976. (15) Shi, Y. Y.; Wang, L.; van Gunsteren, W. F. Mol. Sim. 1988, 1, 369. (16) Kessler, K.; Griesinger, C.; Wagner, K. J. Am. Chem. Soc. 1987, 109, 6927. (17) Schmidt, J. M. To appear in Supramol. Struct. Funct. (18) Ma´di, Z. L.; Griesinger, C.; Ernst, R. R. J. Am. Chem. Soc. 1990, 112, 2908. (19) Bru¨schweiler, R.; Roux, B.; Blackledge, M.; Griesinger, C.; Karplus, M.; Ernst, R. R. J. Am. Chem. Soc. 1992, 114, 2285. (20) Bru¨schweiler, R.; Blackledge, M.; Ernst, R. R. J. Biomol. NMR 1991, 1, 3. (21) Blackledge, M. J.; Bru¨schweiler, R.; Griesinger, C.; Schmidt, J. M.; Xu, P.; Ernst, R. R. Biochemistry 1993, 32, 10960. (22) Brunne, R. M.; van Gunsteren, W. F.; Bru¨schweiler, R.; Ernst, R. R. J. Am. Chem. Soc. 1993, 115, 4764. (23) Schmidt, J. M.; Bru¨schweiler, R.; Ernst, R. R.; Dunbrack, R. L., Jr.; Joseph, D.; Karplus, M. J. Am. Chem. Soc. 1993, 115, 8747. (24) Bremi, T.; Ernst, M.; Ernst, R. R. J. Phys. Chem. 1994, 98, 9322. (25) Gelin, B. R.; Karplus, M. Biochemistry 1979, 18, 1256. (26) Hermans, J.; Yun, R. H.; Anderson, A. G. J. Comput. Chem. 1992, 13, 429. (27) Straatsma, T. P.; McCammon, J. A. J. Chem. Phys. 1994, 101, 5032. (28) Hooft, R. W. W.; van Eijck, B. P.; Kroon, J. J. Chem. Phys. 1992, 97, 6690. (29) Karplus, M. J. Chem. Phys. 1959, 30, 11. (30) Karplus, M. J. Am. Chem. Soc. 1963, 85, 2870. (31) Bystrov, V. F. Prog. Nucl. Magn. Reson. Spectrosc. 1976, 10, 41. (32) Pachler, K. G. R. Spectrochim. Acta 1963, 19, 2085. (33) Chandler, D. J. Chem. Phys. 1978, 68, 2959. (34) McCammon, J. A.; Karplus, M. Proc. Natl. Acad. Sci. U.S.A. 1979, 76, 3585. (35) Northrup, S. H.; Pear, M. R.; Lee, C.-Y.; McCammon, J. A.; Karplus, M. Proc. Natl. Acad. Sci. U.S.A. 1982, 79, 4035. (36) van Gunsteren, W. F.; Berendsen, H. J. C. Groningen Molecular Simulation Package and Manual (GROMOS); Biomos: University of Groningen, Groningen, 1987. (37) Ryckaert, J.-P.; Ciccotti, G.; Berendsen, H. J. C. J. Comput. Phys. 1977, 23, 327. (38) Kessler, H.; Mu¨ller, A.; Pook, K.-H. Liebigs Ann. Chem. 1989, 903. (39) Kessler, H.; Mu¨ller, A.; Pook, K.-H. Liebigs Ann. Chem. 1989, 913. (40) Griesinger, C.; Sørensen, O. W.; Ernst, R. R. J. Magn. Reson. 1987, 75, 474. (41) Beutler, T. C.; van Gunsteren, W. F. Chem. Phys. Lett. 1995, 237, 308. (42) Fischman, A. J.; Live, D. H.; Wyssbrod, H. R.; Agosta, W. C.; Cowburn, D. J. Am. Chem. Soc. 1979, 102, 2533. (43) Wasylishen, R.; Schaefer, T. Can. J. Chem. 1972, 50, 2710.
JP951713K
