Mass-Weighted Molecular Dynamics Simulation of Cyclic Polypeptides BORYEU MAO,” C. M. MACCIORA, and K. C. CHOU
U p j o h n Research Laboratories, Kalamazoo, Michigan 49001
SYNOPSIS
A modified molecular dynamics ( M D ) method in which atomic masses are weighted was developed previously for studying the conformational flexibility of neuroregulating tetrapeptide Phe-Met-Arg-Phe-amide ( FMRF-amide) . The method has now been applied to longer and constrained molecules, namely a disulfide-linked cyclic hexapeptide, c [ CYFQNC] , and its linear and “pseudo-cyclic’’ analogues. The sampling of dihedral conformational space of the linear hexapeptide in mass-weighted MD simulations was found to be improved significantly over conventional MD simulations, as in the case of the shorter FMRF-amide molecule studied previously. In the cyclic hexapeptide, the internal constraint of the molecule due to the intramolecular disulfide bond (hence the absence of free terminals in the molecule) does not adversely affect the significant improvement of conformational sampling in mass-weighted MD simulations over normal MD simulations. The pseudocyclic polypeptide is identical to the linear CYFQNC molecule in amino acid sequence (i.e., side chains of the cysteine residues are reduced), but the positions of its two terminal heavy atoms were held fixed in space such that the molecule has a nearly cyclic conformation, For this molecule, the mass-weighted MD simulation generated a wide range of polypeptide backbone conformations covering the internal dihedral degrees of freedom; moreover, the physical space of the pseudo-cyclic structure was also sampled in a complete revolution of the entire molecular fragment about the two fixed termini during the simulation. These characteristics suggest that mass-weighted MD can also be an extremely useful method for conformational analyses of constrained molecules and, in particular, for modeling loops on protein surfaces.
INTRODUCTION For the neuroregulating tetrapeptide Phe-Met-ArgPhe-amide (FMRF-amide) found in mollusks, the sampling of backbone dihedral conformational space of the molecule was found to become much more effective when the atomic masses in molecular dynamics ( MD ) simulations are “weighted” algebraically by a numerical factor.’ The trajectory from the mass-weighted molecular dynamics ( MWMD ) simulation has been analyzed for investigating conformational flexibility of the molecular backbone of FMRF-amide.’ The essence of the MWMD simulation as implemented is summarized as follows. Biopolymers, Vol. 31, 1077-1086 (1991) 0 1991 John Wiley & Sons, Inc.
CCC 0006-3525/91/091077-l0$04.00
* Address correspondenceto this author.
Atomic masses in the molecular simulation system were multiplied uniformly by a numerical weighting factor w ( w > 1.0) ,and rigid internal degrees of freedom of the system (i.e., covalent bonds, bond angles, amide planes, and tetrahedral carbon centers) were maintained by constraints. Potential barriers in the “soft” degrees of freedom (e.g., dihedral rotations) are thereby overcome more readily due to large effective momenta of the particles in the otherwise standard MD simulation. The mass-weighting scheme has an advantage over the traditional hightemperature scheme in that atomic velocities can remain in the regime in which the MD integration maintains numerical stability. A detailed comparison of the mass-weighting approach with normaland high-temperature MD simulations will be reported elsewhere; the physics of the mass-weighted 1077
1078
MAO, MAGGIORA, AND CHOU
system as well as the effects of constraints imposed on rigid internal degrees of freedom and variable mass-weighting factors will be d e ~ c r i b e d . ~ A natural extension of the application of this method would be the conformational study of larger and structurally more complex molecules. This extension is especially relevant since many naturally occurring polypeptides are longer than the tetrapeptide FMRF-amide, and they may be cyclic due to disulfide and other intramolecular linkages. For example, the pituitary hormones vasopressin and oxytocin, which regulate lactation and water balance, respectively, are ring-containing nonapeptides. These and other cyclic polypeptide molecules have previously been studied by computer simulation^.^-^ Herein we report the successful application of the MWMD method to a prototypical constrained molecular system, the cyclic hexapeptide Cys-Tyr-PheGln-Asn-Cys cyclized by an intramolecular disulfide linkage. This disulfide-linked cyclic hexapeptide, denoted as c [ CYFQNC] ,is the ring moiety of bovine vasopressin and is also known as pressinoic acid; bovine vasopressin contains a tripeptide tail, ProArg-Gly-amide, at the carboxy-terminus of pressinoic acid. Results from normal-mass molecular dynamics (NMMD) and MWMD simulations of pressinoic acid c [ CYFQNC ] are compared. From the distribution of numerical values of an individual dihedral angle within its full 360" range in an MD simulation, two numerical quantities will be computed as indicators of the extent of sampling the dihedral conformation space. The comparison of these two indicators obtained from MWMD and NMMD simulations shows that the conformational sampling of backbone dihedral angles (4's and $'s) and disulfide dihedral angles ( X ' S ) of c [ CYFQNC] is increased significantly in MWMD simulations. For comparison with the earlier studies of the tetrapeptide FMRF-amide, conformational sampling of the linear hexapeptide CYFQNC by MWMD was also investigated. The results show that the method is scalable at least to polypeptide molecules of six amino acid residues in length; a possible limit on the molecular size for the application of this method is currently under investigation. As a further indication of the capability of the MWMD method for sampling polypeptide conformations, simulations of another constrained molecular system, a "pseudo-cyclic'' analogue of CYFQNC, are also reported. In this pseudo-cyclic molecule, denoted as s [ CYFQNC] ,the two terminal atoms of the polypeptide backbone were held fixed in spatial positions such that the molecular fragment has a nearly cyclic conformation. Results from the
MWMD simulation of this pseudo-cyclic analogue show that, in addition to the improved sampling of internal dihedral degrees of freedom compared to NMMD, the additional rotational degree of freedom of the molecular fragment around the fixed points was also sampled. The advantage of atomic massweighting demonstrated here for the pseudo-cyclic molecule suggests that such a procedure can generate a wide range of plausible conformations of protein surface loops that could be used in crystallographic structural refinement and loop modeling.'
METHODS Structure generation, energy minimization, and MD calcu1.ations are performed within the macromolecular program system CHARMM'; in particular, constraints and atomic mass-weighting were implemented with standard procedures and options of the program package. Except where noted explicitly, parameters in the mass-weighting scheme of MWMD simulations are those described earlier',' : the mass-weighting factor w was set at 10.0 and the equivalent system temperature was set at 600 K (i.e., the atomic velocity distribution corresponds to that in a normal-mass system at 600 K ) . Constraining potentials placed on rigid internal degrees of freedom have force constants of 10000 kcal/mole-A' and 4000 kcal/mole-rad' for covalent bonds and valence angles, respectively; the force constants for improper dihedral angles, which maintain planar and tetrahedral geometries for amides and tetrahedral carbon centers, respectively, are 400 kcal/ mole-rad' for c [ CYFQNC] , and 1500.0 kcallmole-rad' for the linear hexapeptide CYFQNC and the pseudo-cyclic molecule s [ CYFQNC 1. For c [ CYFQNC] , shown schematically in Figure 1, an initial conformation was first generated from a set of standard amino acid residue conformations averaged from crystal structures of known proteins ( Mao, unpublished data). The distance between the S, atoms of two cysteine residues is 15.75 A in this initial structure; it is shortened to 2.02 A after adjustments of the $ and $ dihedral angles of the polypeptide backbone and subsequent energy minimization of the full cyclic structure. MD simulations of c [ CYFQNC] were started from this energy-minimized structure. The disulfide bond in the energyminimized c [ CYFQNC] structure was then removed and positions of cysteine side-chain atoms were readjusted by additional local energy minimization (cf. Figure 1); the resulting structure served
MWMD SIMULATION OF CYCLIC POLYPEPTIDES
I
C
I
N
Figure 1. Schematic drawing of the cyclic hexapeptide c [ CYFQNC 1. Each arrow indicates a dihedral angle that is allowed to freely rotate during simulations. Dihedral angles 1-10 are the 6 and $ angles of the polypeptide backbone, and dihedral angles 11-15 are side-chain x's of the disulfide linkage between the two cysteine residues.
a s the starting conformation in MD simulations of CYFQNC and s[CYFQNC] molecules. For c [ CYFQNC] and CYFQNC, all atoms are allowed to move during the simulations, whereas for s [ CYFQNC] , the two terminal atoms of the hexapeptide backbone, the nitrogen atom of the aminoterminal cysteine and the carbonyl carbon atom of carhoxv-terminal cysteine (cf. Figure 1 ) , are held fixed in space a t their respective starting positions. For MD simulations, the Verlet algorithm was used for integration, with a time step of 0.001 ps. The dielectric constant is set a t 80.0; lists of atom pairs in nonbonded and hydrogen-bond interactions were generated a t 20-step intervals, with a cutoff distance of 8.0 A for nonbonded interactions and 4.5 A f'or hydrogen-bonded interactions. T h e temperature was raised t o the final value in 20 increments, followed by 20 cycles of reequilibriation a t the final temperature; for the heating and equilibriation. a Gaussian distribution a t the given temperature was assigned to atomic velocities, at intervals of 0.5 ps. For both NMMD and MWMD, data from 100-ps simulations are analyzed and reported. In order to characterize the extent of sampling the individual dihedral conformational space in a quantitative manner, numerical values of each dihedral angle a t all time steps during either a normal-
1079
mass or a mass-weighted simulation were obtained from the trajectory, and from the distribution of these values in the full 360" range two suitable indicators are defined and computed as follows. The full 360" range is divided into 360 bins, each of 1" width, distributed from 1" to 360". T h e number of times during the simulation that the dihedral angle has a value t h a t falls within each bin is counted; the counts are then plotted for all 360 bins. In this distribution, nonempty bins indicate regions in the dihedral angle space that have been reached during the simulation, and the number of such nonempty bins is then a n indicator of the range of conformational sampling ( R ) during the simulation. Another characteristic of this distribution to be recognized as the second indicator of the conforma,which is defined tional sampling is the coverage ( C ) as the fraction of the rectangular area bounded by 1and 360 in the x direction and 0 and the maximum count in the y direction that is located under the distribution curve. According to these definitions, R has a minimum value of 1 and a maximum value of 360; the larger the value of R for a given distribution, the larger the region in the 360" space has been reached during a simulation. On the other hand, C has a minimum value of 0.278% ( 1/360 ) for a n extremely narrow distribution of only one filled bin and a maximum value of 100.0% for a n entirely uniform distribution. These two parameters are intrinsically independent indicators of the conformational sampling in that a distribution of large R may cover most of the 360" space only sparsely, resulting in a relatively small value for C. Thus, in general, high values of both R and C indicate that not only is a large region of the 360" space reachable but that potential wells within the range are also sampled more truthfully. Whereas other details of the distribution (e.g., symmetry in peak locations and heights) may be characterized by additional parameters, R and C as defined here provide straightforward descriptions of the extent of sampling the conformational space of a n individual dihedral angle. For analyses of the conformation of s [ CYFQNC] , a coordinate system is defined for which the origin is located a t the nitrogen atom of the amino-terminal cysteine in the energy-minimized structure; the x axis points in the direction of the carbonyl carbon atom of the carboxy-terminal cysteine and t h e y axis points in a direction such that the pseudo-cyclic ring lies approximately in the xy plane. In this coordinate system, fluctuations of backbone atoms with respect to the plane of the initial ring conformation during MD simulations are characterized by the amplitudes of their motion in the z direction ( z amplitudes).
1080
MAO, MAGGIORA, AND CHOU 600
400 200 0 600
400 200 0
600 400
200 0 0'
90
270'
100'
360
Figure 2. Distribution of the numerical values in the 360' range for three of the dihedral angles in c [ CYFQNC] during MD simulations: ( a ) dihedral angle 8, ( b ) 9, and ( c ) 14. For each dihedral angle, the 360' space is divided into 1' bins and number of counts of its numerical values falling in each bin during MD simulations is plotted against the 360 bins. Curves in thin lines are obtained from NMMD simulation and curves in thick lines are obtained from MWMD simulation ( o= 10.0, T, = 600 K ) .
RESULTS AND DISCUSSION c[CYFQNC J
Figure 2 shows the distributions of numerical values of several representative dihedral angles calculated from MD simulations of c [ CYFQNC] . Narrow and broad peaks are observed in NMMD (dihedral angles 14 and 9 respectively, cf. Figure 1), and multiple peaks (dihedral angle 8 ) are also observed. In MWMD, not only is the entire 360" range of each dihedral angle sampled, but potential wells within the 360' range are covered more thoroughly. For example, the sampling range of dihedral angle 14 is extended from a narrow peak at 114" to the full 360" range (Figure 2c); on the other hand, regions of the two potential wells for dihedral angle 8 are covered more thoroughly (Figure 2a). These two characteristics of conformational sampling are quantified numerically by the range ( R ) and the coverage ( C ) indicators described in Methods. In Table I, values of R and C for each of the 15 dihedral angles in c [ CYFQNC] are listed for NMMD and MWMD. The improvement in the conformational sampling in MWMD over conventional NMMD simulations is measured quantitatively by the two ratios also shown in Table I, the R ratio (column 5 ) and the C ratio (column 6 ) . For dihedral angle 14 described above, the R ratio is 3.16, indicating a
threefold increase in the range of dihedral space that is reached during the simulation (cf. Figure 2c); for dihedral angle 8, the C ratio is 2.23, indicating a twofold increase in the overall coverage of potential Table I The R and the C Indicators for the Conformational Sampling of the 15 Dihedral Angles in c[CYFQNC] Shown in Figure 1 in Normal-Mass and Mass-Weighted MD Simulations NMMD
R
C
MWMD
R
Ratios
C
R
C
~
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
360 171 360 325 277 173 359 279 239 179 209 263 360 114 257
24.32% 15.89% 23.01% 11.15% 14.38% 17.77% 32.31% 15.07% 17.15% 18.26% 10.43% 16.28% 22.06% 9.05% 16.74%
360 360 360 360 360 360 360 360 360 360 360 360 360 360 360
38.27% 30.69% 49.19% 30.34% 53.05% 36.02% 61.46% 33.67% 42.46% 26.58% 44.76% 39.77% 46.66% 47.42% 42.80%
1.00 2.11 1.00 1.11 1.30 2.08 1.00 1.29 1.51 2.01 1.72 1.37 1.00 3.16 1.40
1.57 1.93 2.14 2.72 3.69 2.03 1.90 2.23 2.48 1.46 4.29 2.44 2.11 5.24 2.56
MWMD SIMULATION OF CYCLIC POLYPEPTIDES
1081
Figure 3. Molecular structures of c [ CYFQNC] before (solid lines) and during MWMD simulation (dashed lines). The dihedral angles 1-15 (cf. Figure 1) are -175.2, -68.8,139.8, 60.7, -154.9, -75.7, 113.0, 63.7, 170.0, -86.3, -55.8, -61.8, 96.9, 91.0, -82.3 in the initial structure, and are 141.0, -157.5, 115.2, 114.4, -89.0, -76.4, 69.7, 14.5, 71.9, 54.5, -67.8, -143.0, 165.8, 54.0, 13.7 in the structure extracted from MWMD trajectory.
wells for the dihedral angle (cf. Figure 2a). As shown in Table I, the MWMDimprovesonly the coverage for dihedral angles 1, 3, 7, and 13 since the entire 360" range for each angle was already accessible in the NMMD. The MWMD method improves both R and C for all the remaining dihedral angles. T h e symmetric distribution of the two peaks for dihedral angle 14 about the disulfide bond suggests, moreover, that the sampling of dihedral conformational space, for this dihedral angle a t least, is near convergence. Figure 3 shows the energy-minimized structure of c [ CYFQNC] and a conformation randomly extracted from the MWMD trajectory. The comparison of these two structures indicates that the structural integrity of the molecule is maintained in MWMD and that conformations significantly different from the starting structure are sampled during the simulation.
s [ CYFQNC]
As expected, the two termini of the backbone of the linear CYFQNC molecule have considerable positional fluctuations; in NMMD the RMS positional fluctuations of the nitrogen atom of the amino-terminal cysteine and the carbonyl carbon atom of the carboxy-terminal cysteine are 11.1 and 9.2 A, respectively, and increased slightly to 12.9 and 10.3 A, respectively, in the MWMD simulation. In contrast, in c[CYFQNC] they are 5.9 and 6.0 A in NMMD, and 7.8 and 7.2 A in MWMD. The N-toC distance is 6.6 A in the starting structure of CYFQNC, and varies from 3.2 to 16.9 A in NMMD Table I1 The R and the C Indicators for the Conformational Sampling of Backbone Dihedral Angles in CYFQNC in Normal-Mass and Mass-Weighted MD Simulations NMMD
CYFQNC
For comparison with the earlier MWMD study of the tetrapeptide FMRF-amide, the conformational sampling of the linear molecule CYFQNC in NMMD and MWMD simulations was investigated, and the results are shown in Table 11. R and C from the MWMD simulation of the linear hexapeptide are both higher in comparison with the NMMD results, indicating that the dihedral conformational sampling is improved in MWMD, as was shown in the case of the shorter FMRF-amide studied earlier.],'
1 2 3 4 5 6 7 8 9 10
MWMD
R
C
R
C
360 202 357 306 359 188 259 335 360 222
21.75% 16.19% 29.73% 26.61% 25.43% 19.81% 19.49% 24.28% 40.69% 21.70%
360 359 360 360 360 360 360 360 360 341
53.26% 29.86% 50.51% 51.70% 47.33% 34.60% 38.83% 39.83% 52.73% 33.38%
Ratios
R 1.00 1.78 1.01 1.18 1.00 1.91 1.39 1.07 1.00 1.54
C 2.45 1.84 1.70 1.94 1.86 1.75 1.99 1.64 1.30 1.54
1082
MAO, MAGGIORA, AND CHOU KcaVmole
220.
215.
210.
-
205.
200. 20.0
I
I
I
40.0
60.0
80.0
100.0 ps
Figure 4. Fluctuation of total energy during conventional MD simulation (lower curve) and MWMD simulation (upper curve) of the pseudo-cyclic hexapeptide s [ CYFQNC ] ; the first 20.0 ps of the simulations was devoted to heating and equilibration, and is omitted from the figure. For the MWMD simulation, the total energy is reduced by a factor of 10.0 for plotting on the same scale.
and from 3.4 t o 18.3 A in MWMD; for c [ CYFQNC] , the distance varies from 3.2 t o 8.3 A in NMMD and from 3.1 to 9.1 A in MWMD. Thus, while the extensive conformational sampling of dihedral angle space in linear and cyclic polypeptides during MWMD simulations encourages a similar procedure for generating conformations of surface loops of proteins, the fluctuation of the distance between terminal atoms observed in these molecules is too large for these molecular fragments t o represent loops that are anchored on a protein molecule and have more or less constant N-to-C distances. It is possible, however, to hold the two terminal atoms of the hexapeptide fixed in space in a pseudocyclic conformation representative of a loop; this is achieved readily with options available in CHARMM. Figure 4 shows the fluctuations of the total energy during NMMD and MWMD of the pseudo-cyclic hexapeptide s [ CYFQNC] . The total energy of the system decreases during the 80 ps of NMMD simulation; this instability in the MD integration is due t o the fixed terminal atoms. In contrast, the total energy in MWMD remains essentially constant in the same time period; the larger kinetic energy due to the greater atomic masses increases the magnitude of fluctuations in the total energy, but the balance of the kinetic and potential
energies improves the long term numerical stability of molecular dynamics integration. Table I11 shows the improvements in the R and C indicators for the sampling of dihedral angle conformational space of s [ CYFQNC] in MWMD over NMMD. In Table IV, results from MWMD simulations with a higher mass-weighting factor and/or Table I11 The R and the C Indicators for the Conformational Sampling of Backbone Dihedral Angles in s[CYFQNC] in Normal-Mass and MassWeighted MD Simulations
NMMD R 1 2 3 4 5 6 7 8 9 10
360 175 343 309 284 232 324 352 360 194
MWMD
C
R
C
22.40% 17.47% 18.41% 21.90% 18.87% 23.63% 21.28% 22.73% 38.60% 22.46%
360 353 360 360 348 345 360 360 360 359
39.12% 24.73% 39.18% 40.19% 34.83% 28.26% 35.87% 36.52% 56.10% 31.65%
Ratios
R 1.00 2.02 1.05 1.17 1.23 1.49 1.11 1.02 1.00 1.85
C 1.75 1.42 2.13 1.84 1.85 1.20 1.69 1.61 1.45 1.41
MWMD SIMULATION OF CYCLIC POLYPEPTIDES
1083
Table IV Effects of Varying Mass Weight Factor w and/or Equivalent System Temperature T, on the R and the C Indicators of the Conformational Sampling of Backbone Dihedral Angles of s[CYFQNC] in MWMD Simulations (Both Sets of Ratios are Calculated with Respect to the Simulation at w = 10.0 and T,= 600 K) 10.0/600 K
1 2 3 4 5 6 7 8 9 10
30.0/600 K
30.0/1000 K
Ratios
Ratios
R
C
R
C
R
C
R
C
K
C
360 353 360 360 348 345 360 360 360 359
39.12% 24.73% 39.18% 40.19% 34.83% 28.26% 35.87% 36.52% 56.10% 31.65%
360 360 360 360 360 360 360 360 360 360
44.09% 42.66% 44.69% 50.41% 52.00% 45.30% 47.85% 45.93% 58.98% 37.34%
1.00 1.02 1.oo 1.00 1.03 1.04 1.00 1.oo 1.00 1.00
1.13 1.72 1.14 1.25 1.49 1.60 1.33 1.26 1.05 1.18
360 360 36G 360 360 360 360 360 360 360
61.89% 35.68% 59.65% 43.22% 47.08% 43.51% 50.32% 51.90% 52.94% 52.11%
1.00 1.02 1.00
1.58 1.44 1.52 1.08 1.35 1.54 1.40 1.42 0.94 1.65
a higher equivalent system temperature are compared. The improved sampling of the dihedral (internal) degrees of freedom of s [ CYFQNC] in MWMD is also accompanied by an improved sampling of the physical (external) space of the molecule around the fixed termini. Figure 5 shows the RMS positional fluctuations of the 18 backbone atoms, and Figure 6 shows their maximum amplitudes above and below the initial plane of the molecule. Effects of different mass-weighting factors and/or different equivalent
1.00
1.03 1.04 1.00 1.00 1.00 1.00
system temperatures are also shown in these figures. The improved sampling of the external space by various MWMD simulations is indicated by the clearly larger positional fluctuations and z amplitudes. The largest z amplitudes were observed for the simulation at a mass-weighting factor ( w ) of 30.0 and at an equivalent system temperature (T,) of
A 7. 5.
3.
A
1.
-1.
-3.
1
0.
Y
,
I
I
I
1
5
8
11
14
Y 18
Figure 5. Position fluctuations of atoms in the backbone of s[CYFQNC]. Curve 1 is calculated from the NMMD simulation, and curves 2-4 are calculated respectively from MWMD simulations at 10.0/600 K ( w / T e ) ,30.0/600 K, and 30.0/1000 K.
5
8
11
14
18
Figure 6. Fluctuations of the z coordinate of backbone atoms in s [ CYFQNC] . The ring in the starting conformation for MD simulations is tilted in the xy plane (curve 0 ) . The z amplitudes of each atom during the four simulations described in Figure 5 are plotted; for example, curves l a and l b are respectively the positive z amplitudes and the negative z amplitudes from the normal-mass simulation. Only in the simulation at w = 30.0 and T, = 1000 K are the z amplitudes (4a and 4b) symmetric with respect to the xy plane.
1084
MAO, MAGGIORA, AND CHOU
1000 K; the magnitude of ca. 7 A is comparable to the size of the pseudo-cyclic ring, which is approximately 7.5 A, measured from the fixed termini to the amide linkage between Phe-3 and Gln-4. These large z amplitudes, and their approximate symmetry with respect to the xy plane ( z = O.O), are in fact due to a complete revolution of the pseudo-cyclic ring about its fixed termini. The trajectory of the centroid of the ring, shown in two projections in Figure 7, indicates that the additional external degree of freedom of the molecule is sampled in a complete revolution about the axis of fixed termini dur-
A (z-axis)
A (z-axis) 5.
ing MWMD (Figure 7a); in comparison, the centroid remains stationary in the NMMD simulation (Figure 7b). The sampling of the physical space of s [ CYFQNC] shown in the MWMD simulation is an important and desired capability of any method for generating plausible conformations of a protein surface loop. Two representative conformations in the MWMD simulation that have large z amplitudes are compared with the initial structure in Figure 8; it should be noted that the MWMD structure in Figure 8a is twisted from the nearly planar initial conformation.
9% 5'
3.
3.
1.
1.
-1.
-1.
-3.
-3.-
-5.
I
-5.
I
5.
3.
1.
A (y-axis)
-3.
I
I
I
-1.
1.
3.
5.
:I
5.
(c)
-1.
-1.
-3.
-3.
-1.
I
-5.
A (y-axis)
A (y-axis) 5.
-5.
::
-5. -LA
-1.
-3.
-
I
1.
I
1
I
3.
5.
7.
-5.
A (x-axis)
-1.
I
I
I
,
1.
3.
5.
7.
Figure 7. Trajectories of the ring centroid of s [ CYFQNC] during MD simulations, projected onto the yz plane [ ( a ) for MWMD and ( b ) for NMMD] and onto the xy plane [ ( c ) for MWMD and ( d ) for NMMD]. During the MWMD simulation at w = 30.0 and T, = 1000 K, the ring makes over one revolution about the n axis ( a ) , whereas it remains in the first quadrant during the normal MD simulation ( b ).
A (y-axis)
MWMD SIMULATION OF CYCLIC POLYPEPTIDES
1085
Figure 8. Backbone conformation of s [ CYFQNC] in the starting structure (solid lines) and in structures extracted from the MWMD trajectory at w = 30.0 and T, = 1000 K (dashed lines). The origin of the Cartesian coordinate system is located at the nitrogen atom of the amino-terminal cysteine (labeled N in the figure) and the x axis points to the carbonyl carbon atom of the carboxy-terminal cysteine (labeled C ). One structure from the MWMD trajectory is shown to have a large positive z amplitude ( a ) and another a large negative z amplitude ( b ) .
SUMMARY We demonstrate here that the MWMD method can significantly improve the conformational sampling of dihedral angles in different polypeptide structures. As implemented, the mass-weighting scheme considerably improves the capability of conformational sampling by an MD simulation, and thus reduces the dependence on the starting structure selected for a given simulation. From the distribution of numerical values of a given dihedral angle during a simulation, two parameters ( range and coverage ) are defined as indicators of the extent of confor-
mational sampling of the internal coordinate space. They are measures of the accessibility of the full 360" range of the dihedral angle and the nonsparsity of the distribution. The effects of variable temperatures and variable mass-weighting factors are also described (for the pseudo-cyclic analogue s [ CYFQNC] ) and a systematic study of these variables of the MWMD method will be reported elsewhere. Judging from results on the linear hexapeptide reported here, the method has the potential of scaling up to longer polypeptide chains. For cyclic molecules and other constrained systems, procedures
1086
MAO, MAGGIORA, AND CHOU
for conformation searching that are based on standard MD simulation have the advantage over other methods in that the structural constraints (e.g., “ring closure” ) is satisfied without additional work; the capability for conformational sampling of these procedures, however, is somewhat limited due to the presence of structural constraints (Figure 2 ) T h e results on constrained molecular systems reported here show that this limitation can be greatly removed by the mass-weighting scheme; the absence of free termini in internally constrained cyclic molecules does not adversely affect conformational sampling in MWMD simulations. Whereas the extensive sampling of individual dihedral conformational space of pressinoic acid in MWMD has been quantified by increases in the R and C indicators, the sampling of molecular conformations requires further analyses by procedures similar to those described in the study of a linear polypeptide molecule2;such analyses will be carried out and results compared with experimental data.” T h e procedures should also enhance methods of determining molecular structures that must satisfy experimentally determined intramolecular constraints.13 Finally, the sampling of external physical space of the pseudo-cyclic s [ CYFQNC ] in a complete revolution about a fixed axis observed in the MWMD simulation suggests that the method can have important applications for generating loop conformations in proteins. Due to the presence of neighboring groups in a protein, a polypeptide loop anchored a t protein surface will generally have a smaller space for libration about their anchor points compared to the complete rotational freedom for the isolated pseudo-cyclic peptide studied in this report. Thus, given t h a t the MWMD simulation is capable of sampling the external space available for s [ CYFQNC] in a complete revolution of the molecular fragment, it is expected that the method will efficiently sample the relatively smaller space around a protein surface loop. T h e potential application of
the MWMD method in generating plausible strqctures for loop modeling and for crystallographic refinement of protein surface loops will be investigated. The authors would like to acknowledge helpful discussion with D. J. Duchamp, and thank M. Karplus for making the Harvard macromolecular computation program CHARMM available.
REFERENCES 1. Mao, B. & Friedman, A. R. (1990) Biophys. J. 58, 803-805. 2. Mao, B. (1990) Biophys. J . 60,611-622. 3. Mao, B., in preparation. 4. Hagler, A. T., Osguthorpe, D. J., Dauber-Osguthorpe, P. & Hempel, J. C. (1985) Science 2 2 7 , 1309-1315. 5. Somoza, J. R. & Brady, J. W. (1988) Biopolymers 27, 939-956. 6. Al-Obeidi, F., Hadley, M. E., Pettitt, B. M. & Hruby, V. J. (1989) J. Am. Chem. SOC.111, 3413-3416. 7. Briinger, A., Kuriyan, J. & Karplus, M. (1987) Science 235,458-460. 8. Chothia, C., Lesk, A. M., Tramontano, A., Levitt, M., Smith-Gill, S. J., Air, G., Sheriff, S., Padlan, E. Q., Davies, D., Tulip, W. R., Colman, P. M., Spinelli, S., Alzari, P. M. & Poljak, R. J. (1989) Nature 342,877883. 9. Brooks, B. R., Bruccoleri, R. E., Olafson, B. D., States, D. J., Swaminathan, S. & Karplus, M. (1983) J. Comp. Chem. 4, 187-217. 10. Bruccoleri, R. E. & Karplus, M. ( 1987) Biopolymers 26,137-168. 11. Peishoff, C. E., Dixon, J. S. & Kopple, K. D. (1990) Biopolymers 30, 45-46. 12. Langs, D. A., Smith, G. D., Stezowski, J. J. & Hughes, R. E. (1986) Science 232, 1240-1242. 13. Constantine, K. L., De Marco, A., Madrid, M., Brooks, C. L. & Llinas, M. (1990) Biopolymers 30, 239-256. Received December 11, 1990 Accepted May 6, 2991
U p j o h n Research Laboratories, Kalamazoo, Michigan 49001
SYNOPSIS
A modified molecular dynamics ( M D ) method in which atomic masses are weighted was developed previously for studying the conformational flexibility of neuroregulating tetrapeptide Phe-Met-Arg-Phe-amide ( FMRF-amide) . The method has now been applied to longer and constrained molecules, namely a disulfide-linked cyclic hexapeptide, c [ CYFQNC] , and its linear and “pseudo-cyclic’’ analogues. The sampling of dihedral conformational space of the linear hexapeptide in mass-weighted MD simulations was found to be improved significantly over conventional MD simulations, as in the case of the shorter FMRF-amide molecule studied previously. In the cyclic hexapeptide, the internal constraint of the molecule due to the intramolecular disulfide bond (hence the absence of free terminals in the molecule) does not adversely affect the significant improvement of conformational sampling in mass-weighted MD simulations over normal MD simulations. The pseudocyclic polypeptide is identical to the linear CYFQNC molecule in amino acid sequence (i.e., side chains of the cysteine residues are reduced), but the positions of its two terminal heavy atoms were held fixed in space such that the molecule has a nearly cyclic conformation, For this molecule, the mass-weighted MD simulation generated a wide range of polypeptide backbone conformations covering the internal dihedral degrees of freedom; moreover, the physical space of the pseudo-cyclic structure was also sampled in a complete revolution of the entire molecular fragment about the two fixed termini during the simulation. These characteristics suggest that mass-weighted MD can also be an extremely useful method for conformational analyses of constrained molecules and, in particular, for modeling loops on protein surfaces.
INTRODUCTION For the neuroregulating tetrapeptide Phe-Met-ArgPhe-amide (FMRF-amide) found in mollusks, the sampling of backbone dihedral conformational space of the molecule was found to become much more effective when the atomic masses in molecular dynamics ( MD ) simulations are “weighted” algebraically by a numerical factor.’ The trajectory from the mass-weighted molecular dynamics ( MWMD ) simulation has been analyzed for investigating conformational flexibility of the molecular backbone of FMRF-amide.’ The essence of the MWMD simulation as implemented is summarized as follows. Biopolymers, Vol. 31, 1077-1086 (1991) 0 1991 John Wiley & Sons, Inc.
CCC 0006-3525/91/091077-l0$04.00
* Address correspondenceto this author.
Atomic masses in the molecular simulation system were multiplied uniformly by a numerical weighting factor w ( w > 1.0) ,and rigid internal degrees of freedom of the system (i.e., covalent bonds, bond angles, amide planes, and tetrahedral carbon centers) were maintained by constraints. Potential barriers in the “soft” degrees of freedom (e.g., dihedral rotations) are thereby overcome more readily due to large effective momenta of the particles in the otherwise standard MD simulation. The mass-weighting scheme has an advantage over the traditional hightemperature scheme in that atomic velocities can remain in the regime in which the MD integration maintains numerical stability. A detailed comparison of the mass-weighting approach with normaland high-temperature MD simulations will be reported elsewhere; the physics of the mass-weighted 1077
1078
MAO, MAGGIORA, AND CHOU
system as well as the effects of constraints imposed on rigid internal degrees of freedom and variable mass-weighting factors will be d e ~ c r i b e d . ~ A natural extension of the application of this method would be the conformational study of larger and structurally more complex molecules. This extension is especially relevant since many naturally occurring polypeptides are longer than the tetrapeptide FMRF-amide, and they may be cyclic due to disulfide and other intramolecular linkages. For example, the pituitary hormones vasopressin and oxytocin, which regulate lactation and water balance, respectively, are ring-containing nonapeptides. These and other cyclic polypeptide molecules have previously been studied by computer simulation^.^-^ Herein we report the successful application of the MWMD method to a prototypical constrained molecular system, the cyclic hexapeptide Cys-Tyr-PheGln-Asn-Cys cyclized by an intramolecular disulfide linkage. This disulfide-linked cyclic hexapeptide, denoted as c [ CYFQNC] ,is the ring moiety of bovine vasopressin and is also known as pressinoic acid; bovine vasopressin contains a tripeptide tail, ProArg-Gly-amide, at the carboxy-terminus of pressinoic acid. Results from normal-mass molecular dynamics (NMMD) and MWMD simulations of pressinoic acid c [ CYFQNC ] are compared. From the distribution of numerical values of an individual dihedral angle within its full 360" range in an MD simulation, two numerical quantities will be computed as indicators of the extent of sampling the dihedral conformation space. The comparison of these two indicators obtained from MWMD and NMMD simulations shows that the conformational sampling of backbone dihedral angles (4's and $'s) and disulfide dihedral angles ( X ' S ) of c [ CYFQNC] is increased significantly in MWMD simulations. For comparison with the earlier studies of the tetrapeptide FMRF-amide, conformational sampling of the linear hexapeptide CYFQNC by MWMD was also investigated. The results show that the method is scalable at least to polypeptide molecules of six amino acid residues in length; a possible limit on the molecular size for the application of this method is currently under investigation. As a further indication of the capability of the MWMD method for sampling polypeptide conformations, simulations of another constrained molecular system, a "pseudo-cyclic'' analogue of CYFQNC, are also reported. In this pseudo-cyclic molecule, denoted as s [ CYFQNC] ,the two terminal atoms of the polypeptide backbone were held fixed in spatial positions such that the molecular fragment has a nearly cyclic conformation. Results from the
MWMD simulation of this pseudo-cyclic analogue show that, in addition to the improved sampling of internal dihedral degrees of freedom compared to NMMD, the additional rotational degree of freedom of the molecular fragment around the fixed points was also sampled. The advantage of atomic massweighting demonstrated here for the pseudo-cyclic molecule suggests that such a procedure can generate a wide range of plausible conformations of protein surface loops that could be used in crystallographic structural refinement and loop modeling.'
METHODS Structure generation, energy minimization, and MD calcu1.ations are performed within the macromolecular program system CHARMM'; in particular, constraints and atomic mass-weighting were implemented with standard procedures and options of the program package. Except where noted explicitly, parameters in the mass-weighting scheme of MWMD simulations are those described earlier',' : the mass-weighting factor w was set at 10.0 and the equivalent system temperature was set at 600 K (i.e., the atomic velocity distribution corresponds to that in a normal-mass system at 600 K ) . Constraining potentials placed on rigid internal degrees of freedom have force constants of 10000 kcal/mole-A' and 4000 kcal/mole-rad' for covalent bonds and valence angles, respectively; the force constants for improper dihedral angles, which maintain planar and tetrahedral geometries for amides and tetrahedral carbon centers, respectively, are 400 kcal/ mole-rad' for c [ CYFQNC] , and 1500.0 kcallmole-rad' for the linear hexapeptide CYFQNC and the pseudo-cyclic molecule s [ CYFQNC 1. For c [ CYFQNC] , shown schematically in Figure 1, an initial conformation was first generated from a set of standard amino acid residue conformations averaged from crystal structures of known proteins ( Mao, unpublished data). The distance between the S, atoms of two cysteine residues is 15.75 A in this initial structure; it is shortened to 2.02 A after adjustments of the $ and $ dihedral angles of the polypeptide backbone and subsequent energy minimization of the full cyclic structure. MD simulations of c [ CYFQNC] were started from this energy-minimized structure. The disulfide bond in the energyminimized c [ CYFQNC] structure was then removed and positions of cysteine side-chain atoms were readjusted by additional local energy minimization (cf. Figure 1); the resulting structure served
MWMD SIMULATION OF CYCLIC POLYPEPTIDES
I
C
I
N
Figure 1. Schematic drawing of the cyclic hexapeptide c [ CYFQNC 1. Each arrow indicates a dihedral angle that is allowed to freely rotate during simulations. Dihedral angles 1-10 are the 6 and $ angles of the polypeptide backbone, and dihedral angles 11-15 are side-chain x's of the disulfide linkage between the two cysteine residues.
a s the starting conformation in MD simulations of CYFQNC and s[CYFQNC] molecules. For c [ CYFQNC] and CYFQNC, all atoms are allowed to move during the simulations, whereas for s [ CYFQNC] , the two terminal atoms of the hexapeptide backbone, the nitrogen atom of the aminoterminal cysteine and the carbonyl carbon atom of carhoxv-terminal cysteine (cf. Figure 1 ) , are held fixed in space a t their respective starting positions. For MD simulations, the Verlet algorithm was used for integration, with a time step of 0.001 ps. The dielectric constant is set a t 80.0; lists of atom pairs in nonbonded and hydrogen-bond interactions were generated a t 20-step intervals, with a cutoff distance of 8.0 A for nonbonded interactions and 4.5 A f'or hydrogen-bonded interactions. T h e temperature was raised t o the final value in 20 increments, followed by 20 cycles of reequilibriation a t the final temperature; for the heating and equilibriation. a Gaussian distribution a t the given temperature was assigned to atomic velocities, at intervals of 0.5 ps. For both NMMD and MWMD, data from 100-ps simulations are analyzed and reported. In order to characterize the extent of sampling the individual dihedral conformational space in a quantitative manner, numerical values of each dihedral angle a t all time steps during either a normal-
1079
mass or a mass-weighted simulation were obtained from the trajectory, and from the distribution of these values in the full 360" range two suitable indicators are defined and computed as follows. The full 360" range is divided into 360 bins, each of 1" width, distributed from 1" to 360". T h e number of times during the simulation that the dihedral angle has a value t h a t falls within each bin is counted; the counts are then plotted for all 360 bins. In this distribution, nonempty bins indicate regions in the dihedral angle space that have been reached during the simulation, and the number of such nonempty bins is then a n indicator of the range of conformational sampling ( R ) during the simulation. Another characteristic of this distribution to be recognized as the second indicator of the conforma,which is defined tional sampling is the coverage ( C ) as the fraction of the rectangular area bounded by 1and 360 in the x direction and 0 and the maximum count in the y direction that is located under the distribution curve. According to these definitions, R has a minimum value of 1 and a maximum value of 360; the larger the value of R for a given distribution, the larger the region in the 360" space has been reached during a simulation. On the other hand, C has a minimum value of 0.278% ( 1/360 ) for a n extremely narrow distribution of only one filled bin and a maximum value of 100.0% for a n entirely uniform distribution. These two parameters are intrinsically independent indicators of the conformational sampling in that a distribution of large R may cover most of the 360" space only sparsely, resulting in a relatively small value for C. Thus, in general, high values of both R and C indicate that not only is a large region of the 360" space reachable but that potential wells within the range are also sampled more truthfully. Whereas other details of the distribution (e.g., symmetry in peak locations and heights) may be characterized by additional parameters, R and C as defined here provide straightforward descriptions of the extent of sampling the conformational space of a n individual dihedral angle. For analyses of the conformation of s [ CYFQNC] , a coordinate system is defined for which the origin is located a t the nitrogen atom of the amino-terminal cysteine in the energy-minimized structure; the x axis points in the direction of the carbonyl carbon atom of the carboxy-terminal cysteine and t h e y axis points in a direction such that the pseudo-cyclic ring lies approximately in the xy plane. In this coordinate system, fluctuations of backbone atoms with respect to the plane of the initial ring conformation during MD simulations are characterized by the amplitudes of their motion in the z direction ( z amplitudes).
1080
MAO, MAGGIORA, AND CHOU 600
400 200 0 600
400 200 0
600 400
200 0 0'
90
270'
100'
360
Figure 2. Distribution of the numerical values in the 360' range for three of the dihedral angles in c [ CYFQNC] during MD simulations: ( a ) dihedral angle 8, ( b ) 9, and ( c ) 14. For each dihedral angle, the 360' space is divided into 1' bins and number of counts of its numerical values falling in each bin during MD simulations is plotted against the 360 bins. Curves in thin lines are obtained from NMMD simulation and curves in thick lines are obtained from MWMD simulation ( o= 10.0, T, = 600 K ) .
RESULTS AND DISCUSSION c[CYFQNC J
Figure 2 shows the distributions of numerical values of several representative dihedral angles calculated from MD simulations of c [ CYFQNC] . Narrow and broad peaks are observed in NMMD (dihedral angles 14 and 9 respectively, cf. Figure 1), and multiple peaks (dihedral angle 8 ) are also observed. In MWMD, not only is the entire 360" range of each dihedral angle sampled, but potential wells within the 360' range are covered more thoroughly. For example, the sampling range of dihedral angle 14 is extended from a narrow peak at 114" to the full 360" range (Figure 2c); on the other hand, regions of the two potential wells for dihedral angle 8 are covered more thoroughly (Figure 2a). These two characteristics of conformational sampling are quantified numerically by the range ( R ) and the coverage ( C ) indicators described in Methods. In Table I, values of R and C for each of the 15 dihedral angles in c [ CYFQNC] are listed for NMMD and MWMD. The improvement in the conformational sampling in MWMD over conventional NMMD simulations is measured quantitatively by the two ratios also shown in Table I, the R ratio (column 5 ) and the C ratio (column 6 ) . For dihedral angle 14 described above, the R ratio is 3.16, indicating a
threefold increase in the range of dihedral space that is reached during the simulation (cf. Figure 2c); for dihedral angle 8, the C ratio is 2.23, indicating a twofold increase in the overall coverage of potential Table I The R and the C Indicators for the Conformational Sampling of the 15 Dihedral Angles in c[CYFQNC] Shown in Figure 1 in Normal-Mass and Mass-Weighted MD Simulations NMMD
R
C
MWMD
R
Ratios
C
R
C
~
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
360 171 360 325 277 173 359 279 239 179 209 263 360 114 257
24.32% 15.89% 23.01% 11.15% 14.38% 17.77% 32.31% 15.07% 17.15% 18.26% 10.43% 16.28% 22.06% 9.05% 16.74%
360 360 360 360 360 360 360 360 360 360 360 360 360 360 360
38.27% 30.69% 49.19% 30.34% 53.05% 36.02% 61.46% 33.67% 42.46% 26.58% 44.76% 39.77% 46.66% 47.42% 42.80%
1.00 2.11 1.00 1.11 1.30 2.08 1.00 1.29 1.51 2.01 1.72 1.37 1.00 3.16 1.40
1.57 1.93 2.14 2.72 3.69 2.03 1.90 2.23 2.48 1.46 4.29 2.44 2.11 5.24 2.56
MWMD SIMULATION OF CYCLIC POLYPEPTIDES
1081
Figure 3. Molecular structures of c [ CYFQNC] before (solid lines) and during MWMD simulation (dashed lines). The dihedral angles 1-15 (cf. Figure 1) are -175.2, -68.8,139.8, 60.7, -154.9, -75.7, 113.0, 63.7, 170.0, -86.3, -55.8, -61.8, 96.9, 91.0, -82.3 in the initial structure, and are 141.0, -157.5, 115.2, 114.4, -89.0, -76.4, 69.7, 14.5, 71.9, 54.5, -67.8, -143.0, 165.8, 54.0, 13.7 in the structure extracted from MWMD trajectory.
wells for the dihedral angle (cf. Figure 2a). As shown in Table I, the MWMDimprovesonly the coverage for dihedral angles 1, 3, 7, and 13 since the entire 360" range for each angle was already accessible in the NMMD. The MWMD method improves both R and C for all the remaining dihedral angles. T h e symmetric distribution of the two peaks for dihedral angle 14 about the disulfide bond suggests, moreover, that the sampling of dihedral conformational space, for this dihedral angle a t least, is near convergence. Figure 3 shows the energy-minimized structure of c [ CYFQNC] and a conformation randomly extracted from the MWMD trajectory. The comparison of these two structures indicates that the structural integrity of the molecule is maintained in MWMD and that conformations significantly different from the starting structure are sampled during the simulation.
s [ CYFQNC]
As expected, the two termini of the backbone of the linear CYFQNC molecule have considerable positional fluctuations; in NMMD the RMS positional fluctuations of the nitrogen atom of the amino-terminal cysteine and the carbonyl carbon atom of the carboxy-terminal cysteine are 11.1 and 9.2 A, respectively, and increased slightly to 12.9 and 10.3 A, respectively, in the MWMD simulation. In contrast, in c[CYFQNC] they are 5.9 and 6.0 A in NMMD, and 7.8 and 7.2 A in MWMD. The N-toC distance is 6.6 A in the starting structure of CYFQNC, and varies from 3.2 to 16.9 A in NMMD Table I1 The R and the C Indicators for the Conformational Sampling of Backbone Dihedral Angles in CYFQNC in Normal-Mass and Mass-Weighted MD Simulations NMMD
CYFQNC
For comparison with the earlier MWMD study of the tetrapeptide FMRF-amide, the conformational sampling of the linear molecule CYFQNC in NMMD and MWMD simulations was investigated, and the results are shown in Table 11. R and C from the MWMD simulation of the linear hexapeptide are both higher in comparison with the NMMD results, indicating that the dihedral conformational sampling is improved in MWMD, as was shown in the case of the shorter FMRF-amide studied earlier.],'
1 2 3 4 5 6 7 8 9 10
MWMD
R
C
R
C
360 202 357 306 359 188 259 335 360 222
21.75% 16.19% 29.73% 26.61% 25.43% 19.81% 19.49% 24.28% 40.69% 21.70%
360 359 360 360 360 360 360 360 360 341
53.26% 29.86% 50.51% 51.70% 47.33% 34.60% 38.83% 39.83% 52.73% 33.38%
Ratios
R 1.00 1.78 1.01 1.18 1.00 1.91 1.39 1.07 1.00 1.54
C 2.45 1.84 1.70 1.94 1.86 1.75 1.99 1.64 1.30 1.54
1082
MAO, MAGGIORA, AND CHOU KcaVmole
220.
215.
210.
-
205.
200. 20.0
I
I
I
40.0
60.0
80.0
100.0 ps
Figure 4. Fluctuation of total energy during conventional MD simulation (lower curve) and MWMD simulation (upper curve) of the pseudo-cyclic hexapeptide s [ CYFQNC ] ; the first 20.0 ps of the simulations was devoted to heating and equilibration, and is omitted from the figure. For the MWMD simulation, the total energy is reduced by a factor of 10.0 for plotting on the same scale.
and from 3.4 t o 18.3 A in MWMD; for c [ CYFQNC] , the distance varies from 3.2 t o 8.3 A in NMMD and from 3.1 to 9.1 A in MWMD. Thus, while the extensive conformational sampling of dihedral angle space in linear and cyclic polypeptides during MWMD simulations encourages a similar procedure for generating conformations of surface loops of proteins, the fluctuation of the distance between terminal atoms observed in these molecules is too large for these molecular fragments t o represent loops that are anchored on a protein molecule and have more or less constant N-to-C distances. It is possible, however, to hold the two terminal atoms of the hexapeptide fixed in space in a pseudocyclic conformation representative of a loop; this is achieved readily with options available in CHARMM. Figure 4 shows the fluctuations of the total energy during NMMD and MWMD of the pseudo-cyclic hexapeptide s [ CYFQNC] . The total energy of the system decreases during the 80 ps of NMMD simulation; this instability in the MD integration is due t o the fixed terminal atoms. In contrast, the total energy in MWMD remains essentially constant in the same time period; the larger kinetic energy due to the greater atomic masses increases the magnitude of fluctuations in the total energy, but the balance of the kinetic and potential
energies improves the long term numerical stability of molecular dynamics integration. Table I11 shows the improvements in the R and C indicators for the sampling of dihedral angle conformational space of s [ CYFQNC] in MWMD over NMMD. In Table IV, results from MWMD simulations with a higher mass-weighting factor and/or Table I11 The R and the C Indicators for the Conformational Sampling of Backbone Dihedral Angles in s[CYFQNC] in Normal-Mass and MassWeighted MD Simulations
NMMD R 1 2 3 4 5 6 7 8 9 10
360 175 343 309 284 232 324 352 360 194
MWMD
C
R
C
22.40% 17.47% 18.41% 21.90% 18.87% 23.63% 21.28% 22.73% 38.60% 22.46%
360 353 360 360 348 345 360 360 360 359
39.12% 24.73% 39.18% 40.19% 34.83% 28.26% 35.87% 36.52% 56.10% 31.65%
Ratios
R 1.00 2.02 1.05 1.17 1.23 1.49 1.11 1.02 1.00 1.85
C 1.75 1.42 2.13 1.84 1.85 1.20 1.69 1.61 1.45 1.41
MWMD SIMULATION OF CYCLIC POLYPEPTIDES
1083
Table IV Effects of Varying Mass Weight Factor w and/or Equivalent System Temperature T, on the R and the C Indicators of the Conformational Sampling of Backbone Dihedral Angles of s[CYFQNC] in MWMD Simulations (Both Sets of Ratios are Calculated with Respect to the Simulation at w = 10.0 and T,= 600 K) 10.0/600 K
1 2 3 4 5 6 7 8 9 10
30.0/600 K
30.0/1000 K
Ratios
Ratios
R
C
R
C
R
C
R
C
K
C
360 353 360 360 348 345 360 360 360 359
39.12% 24.73% 39.18% 40.19% 34.83% 28.26% 35.87% 36.52% 56.10% 31.65%
360 360 360 360 360 360 360 360 360 360
44.09% 42.66% 44.69% 50.41% 52.00% 45.30% 47.85% 45.93% 58.98% 37.34%
1.00 1.02 1.oo 1.00 1.03 1.04 1.00 1.oo 1.00 1.00
1.13 1.72 1.14 1.25 1.49 1.60 1.33 1.26 1.05 1.18
360 360 36G 360 360 360 360 360 360 360
61.89% 35.68% 59.65% 43.22% 47.08% 43.51% 50.32% 51.90% 52.94% 52.11%
1.00 1.02 1.00
1.58 1.44 1.52 1.08 1.35 1.54 1.40 1.42 0.94 1.65
a higher equivalent system temperature are compared. The improved sampling of the dihedral (internal) degrees of freedom of s [ CYFQNC] in MWMD is also accompanied by an improved sampling of the physical (external) space of the molecule around the fixed termini. Figure 5 shows the RMS positional fluctuations of the 18 backbone atoms, and Figure 6 shows their maximum amplitudes above and below the initial plane of the molecule. Effects of different mass-weighting factors and/or different equivalent
1.00
1.03 1.04 1.00 1.00 1.00 1.00
system temperatures are also shown in these figures. The improved sampling of the external space by various MWMD simulations is indicated by the clearly larger positional fluctuations and z amplitudes. The largest z amplitudes were observed for the simulation at a mass-weighting factor ( w ) of 30.0 and at an equivalent system temperature (T,) of
A 7. 5.
3.
A
1.
-1.
-3.
1
0.
Y
,
I
I
I
1
5
8
11
14
Y 18
Figure 5. Position fluctuations of atoms in the backbone of s[CYFQNC]. Curve 1 is calculated from the NMMD simulation, and curves 2-4 are calculated respectively from MWMD simulations at 10.0/600 K ( w / T e ) ,30.0/600 K, and 30.0/1000 K.
5
8
11
14
18
Figure 6. Fluctuations of the z coordinate of backbone atoms in s [ CYFQNC] . The ring in the starting conformation for MD simulations is tilted in the xy plane (curve 0 ) . The z amplitudes of each atom during the four simulations described in Figure 5 are plotted; for example, curves l a and l b are respectively the positive z amplitudes and the negative z amplitudes from the normal-mass simulation. Only in the simulation at w = 30.0 and T, = 1000 K are the z amplitudes (4a and 4b) symmetric with respect to the xy plane.
1084
MAO, MAGGIORA, AND CHOU
1000 K; the magnitude of ca. 7 A is comparable to the size of the pseudo-cyclic ring, which is approximately 7.5 A, measured from the fixed termini to the amide linkage between Phe-3 and Gln-4. These large z amplitudes, and their approximate symmetry with respect to the xy plane ( z = O.O), are in fact due to a complete revolution of the pseudo-cyclic ring about its fixed termini. The trajectory of the centroid of the ring, shown in two projections in Figure 7, indicates that the additional external degree of freedom of the molecule is sampled in a complete revolution about the axis of fixed termini dur-
A (z-axis)
A (z-axis) 5.
ing MWMD (Figure 7a); in comparison, the centroid remains stationary in the NMMD simulation (Figure 7b). The sampling of the physical space of s [ CYFQNC] shown in the MWMD simulation is an important and desired capability of any method for generating plausible conformations of a protein surface loop. Two representative conformations in the MWMD simulation that have large z amplitudes are compared with the initial structure in Figure 8; it should be noted that the MWMD structure in Figure 8a is twisted from the nearly planar initial conformation.
9% 5'
3.
3.
1.
1.
-1.
-1.
-3.
-3.-
-5.
I
-5.
I
5.
3.
1.
A (y-axis)
-3.
I
I
I
-1.
1.
3.
5.
:I
5.
(c)
-1.
-1.
-3.
-3.
-1.
I
-5.
A (y-axis)
A (y-axis) 5.
-5.
::
-5. -LA
-1.
-3.
-
I
1.
I
1
I
3.
5.
7.
-5.
A (x-axis)
-1.
I
I
I
,
1.
3.
5.
7.
Figure 7. Trajectories of the ring centroid of s [ CYFQNC] during MD simulations, projected onto the yz plane [ ( a ) for MWMD and ( b ) for NMMD] and onto the xy plane [ ( c ) for MWMD and ( d ) for NMMD]. During the MWMD simulation at w = 30.0 and T, = 1000 K, the ring makes over one revolution about the n axis ( a ) , whereas it remains in the first quadrant during the normal MD simulation ( b ).
A (y-axis)
MWMD SIMULATION OF CYCLIC POLYPEPTIDES
1085
Figure 8. Backbone conformation of s [ CYFQNC] in the starting structure (solid lines) and in structures extracted from the MWMD trajectory at w = 30.0 and T, = 1000 K (dashed lines). The origin of the Cartesian coordinate system is located at the nitrogen atom of the amino-terminal cysteine (labeled N in the figure) and the x axis points to the carbonyl carbon atom of the carboxy-terminal cysteine (labeled C ). One structure from the MWMD trajectory is shown to have a large positive z amplitude ( a ) and another a large negative z amplitude ( b ) .
SUMMARY We demonstrate here that the MWMD method can significantly improve the conformational sampling of dihedral angles in different polypeptide structures. As implemented, the mass-weighting scheme considerably improves the capability of conformational sampling by an MD simulation, and thus reduces the dependence on the starting structure selected for a given simulation. From the distribution of numerical values of a given dihedral angle during a simulation, two parameters ( range and coverage ) are defined as indicators of the extent of confor-
mational sampling of the internal coordinate space. They are measures of the accessibility of the full 360" range of the dihedral angle and the nonsparsity of the distribution. The effects of variable temperatures and variable mass-weighting factors are also described (for the pseudo-cyclic analogue s [ CYFQNC] ) and a systematic study of these variables of the MWMD method will be reported elsewhere. Judging from results on the linear hexapeptide reported here, the method has the potential of scaling up to longer polypeptide chains. For cyclic molecules and other constrained systems, procedures
1086
MAO, MAGGIORA, AND CHOU
for conformation searching that are based on standard MD simulation have the advantage over other methods in that the structural constraints (e.g., “ring closure” ) is satisfied without additional work; the capability for conformational sampling of these procedures, however, is somewhat limited due to the presence of structural constraints (Figure 2 ) T h e results on constrained molecular systems reported here show that this limitation can be greatly removed by the mass-weighting scheme; the absence of free termini in internally constrained cyclic molecules does not adversely affect conformational sampling in MWMD simulations. Whereas the extensive sampling of individual dihedral conformational space of pressinoic acid in MWMD has been quantified by increases in the R and C indicators, the sampling of molecular conformations requires further analyses by procedures similar to those described in the study of a linear polypeptide molecule2;such analyses will be carried out and results compared with experimental data.” T h e procedures should also enhance methods of determining molecular structures that must satisfy experimentally determined intramolecular constraints.13 Finally, the sampling of external physical space of the pseudo-cyclic s [ CYFQNC ] in a complete revolution about a fixed axis observed in the MWMD simulation suggests that the method can have important applications for generating loop conformations in proteins. Due to the presence of neighboring groups in a protein, a polypeptide loop anchored a t protein surface will generally have a smaller space for libration about their anchor points compared to the complete rotational freedom for the isolated pseudo-cyclic peptide studied in this report. Thus, given t h a t the MWMD simulation is capable of sampling the external space available for s [ CYFQNC] in a complete revolution of the molecular fragment, it is expected that the method will efficiently sample the relatively smaller space around a protein surface loop. T h e potential application of
the MWMD method in generating plausible strqctures for loop modeling and for crystallographic refinement of protein surface loops will be investigated. The authors would like to acknowledge helpful discussion with D. J. Duchamp, and thank M. Karplus for making the Harvard macromolecular computation program CHARMM available.
REFERENCES 1. Mao, B. & Friedman, A. R. (1990) Biophys. J. 58, 803-805. 2. Mao, B. (1990) Biophys. J . 60,611-622. 3. Mao, B., in preparation. 4. Hagler, A. T., Osguthorpe, D. J., Dauber-Osguthorpe, P. & Hempel, J. C. (1985) Science 2 2 7 , 1309-1315. 5. Somoza, J. R. & Brady, J. W. (1988) Biopolymers 27, 939-956. 6. Al-Obeidi, F., Hadley, M. E., Pettitt, B. M. & Hruby, V. J. (1989) J. Am. Chem. SOC.111, 3413-3416. 7. Briinger, A., Kuriyan, J. & Karplus, M. (1987) Science 235,458-460. 8. Chothia, C., Lesk, A. M., Tramontano, A., Levitt, M., Smith-Gill, S. J., Air, G., Sheriff, S., Padlan, E. Q., Davies, D., Tulip, W. R., Colman, P. M., Spinelli, S., Alzari, P. M. & Poljak, R. J. (1989) Nature 342,877883. 9. Brooks, B. R., Bruccoleri, R. E., Olafson, B. D., States, D. J., Swaminathan, S. & Karplus, M. (1983) J. Comp. Chem. 4, 187-217. 10. Bruccoleri, R. E. & Karplus, M. ( 1987) Biopolymers 26,137-168. 11. Peishoff, C. E., Dixon, J. S. & Kopple, K. D. (1990) Biopolymers 30, 45-46. 12. Langs, D. A., Smith, G. D., Stezowski, J. J. & Hughes, R. E. (1986) Science 232, 1240-1242. 13. Constantine, K. L., De Marco, A., Madrid, M., Brooks, C. L. & Llinas, M. (1990) Biopolymers 30, 239-256. Received December 11, 1990 Accepted May 6, 2991
