Conformational Sampling Using High-Temperature Molecular Dynamics ROBERT E. BRUCCOLERI * and MARTIN KARPLUS Departments of Biochemistry & Molecular Biology and Chemistry, Harvard University, Cambridge, Massachusetts 021 38

SYNOPSIS

High-temperature molecular dynamics as a method for conformational search was explored on the antigen combining site of McPC 603, a phosphorylcholine binding immunoglobulin. Simulations a t temperatures of 500, 800, and 1500 K were run for 111.5, 101.7, and 76.3 ps, respectively. The effectiveness of the search was assessed using a variety of methods. For the shorter hypervariable loops, molecular dynamics explored an appreciable fraction of the conformational space as evidenced by a comparison to a simple theoretical model of the size of the conformational space. However, for the longer loops and the antigen combining site as a whole, the simulation times were too short for a complete search. The simulations at 500 and 800 K both generated conformations that minimized to energies 200 kcal/mole lower than the crystal structure. However, the 1500 K simulation produced higher energy structures, even after minimization; in addition, this highest temperature run had many cis-trans peptide isomerizations. This suggests that 1500 K is too high a temperature for unconstrained conformational sampling. Comparison of the results of high temperature molecular dynamics with a direct conformational search method, [ R. E. Bruccoleri & M. Karplus (1987) Biopolymers 26,137-1681. showed that the two methods did not overlap much in conformational space. Simple geometric measures of the conformational space indicated that the direct method covered more space than molecular dynamics a t the lower temperature, but not at 1500 K. The results suggest that high-temperature molecular dynamics can aid in conformational searches.

INTRODUCTION One approach to the prediction of the native structure of a protein or a portion of a protein is based on the “generate and test” paradigm. One generates a set of plausible conformations for a polypeptide chain segment and then evaluates them using a function that approximates the free energy. Such an approach is necessary because the empirical energy functions currently in use, have many local minima, and standard minimization methods are unlikely to find the global minimum when starting with an arbitrary structure.

’-’

0 1990 John Wiley & Sons, Inc. CCC 0006-3525/90/141847-16 $04.00 Biopolymers, Vol. 29, 1847-1862 (1990) * Current address: Department of Macromolecular Modeling, The Squibb Institute for Medical Research, P.O. Box 4000, Princeton, N J 08543.

Given this strategy, it is necessary to have methods for generating plausible conformations. Monte Carlo techniques and simulated annealing13 provide possible generation methods. The distance geometry technique of Purisima and Scheraga14 is another method for generating a set of starting conformations. In a previous paper,15 we described a conformational search procedure ( CONGEN ) that uniformly samples the conformational space of short peptide segments within a protein. Energetic criteria, such as Ramachandran maps of the backbone16and van der Waals interactions with the rest of the protein, were used to limit the conformational space while still allowing the global minimum to be found. Because the sampling is combinatorial in character, the number of conformations generated increased exponentially with the number of amino acids, and a complete search is not feasible for more than about eight or nine residues. 1847

1848

BRUCCOLERI AND KARPLUS

Here, we consider high-temperature molecular is the framework whose sequence and structure are dynamics as an approach for generating conforwell preserved in different immunoglobulins. Thus, mations. In accordance with the philosophy from we can restrict the search to six short polypeptide which simulated annealing was derived, l3 molecular segments whose end points and local protein envidynamics provides a natural method for searching ronment are approximately known. conformational space. The use of high temperatures The hypervariable loops of McPC 603 were simprovides the energy necessary to surmount energetic ulated by molecular dynamics a t temperatures of barriers so that many adjacent local minima can be 500, 800, and 1500 K for periods 111.5, 101.7, and explored. Further, the time spent in particular con76.3 ps, respectively. Although such high temperaformations provides an indication of their entropies. tures are not physically meaningful, they are of inFinally, the computer time necessary for a given terest for exploring conformational space. In the simulation time increases at most quadratically with analysis of the three simulations, we used a number the number of residues, so relatively large systems of methods to assess how thoroughly conformational can be explored. However, this does not necessarily space is searched. We also compare the results to direct searches of conformational space generated solve the search problem because the simulation time necessary to adequately sample the conforby the conformational search procedure, CONGEN. mational space may grow at an exponential rate.17 The next section in this paper describes the moThe central part of a molecular dynamics lecular dynamics simulation. The energetics of the generated conformations and the effectiveness of the simulation l8 is the potential energy which is an evsearch simulation are presented in the succeeding erywhere continuous and differentiable function of section. Finally, we discuss the results from the the coordinates of all the atoms in the system. The viewpoint of the utility of high temperature dynamgradient of the potential yields the forces on the ics as a search method. atoms, which are then integrated to generate the atomic positions as a function of time, i.e., a trajectory. Because the required characteristics of a potential function are very general, and because the DESCRIPTION OF THE SIMULATION searching capability of the simulation are dependent only on the temperatures required to surmount barThe coordinates for the antigen binding (Fab) fragriers, this search method can be applied to functions ment of McPC 603 were obtained from David Davies that are not based on a physical model for the po(personal communication) a t an intermediate stage tential energy. In particular, if we have a differenof refinement. In particular, the coordinates were tiable function that provides an estimate of the free refined with respect to the electron density but not energy, molecular dynamics can be used to explore stereochemistry, so the energy of the coordinates the potential of mean force surface. Alternatively, was high. The coordinates were derived from 2.7 experimental constraints ( Nuclear Overhauser data and had a crystallographic R factor of 24%. enhancement^'^ or x-ray intensities20,21)can be inFor the test simulations, we considered only the troduced into the potential to aid in convergence. antigen combining site, which is located in the variIn this paper, we test high-temperature molecular able domain. We took the variable domain to be the dynamics as a method of conformational search for heavy chain residues 1-122 and light chain residues a portion of a protein of known structure. Because 1-113, where the numbering is strictly sequential of its intrinsic interest, we use the antigen binding site of the immunoglobulin, McPC 603, which binds from the N-terminus. Table I shows that sequence phosphorylcholine.22The immunoglobulin system is of this domain. Because the initial coordinates for a dynamics a very useful system for the exploration of protein simulation should not be subjected to large forces folding algorithms because the major structural variation between different immunoglobulins is due to an improper geometry, the crystallographic confined to the antigen combining site.23-25The coordinates were minimized by 100 steps of Adopted variable domain of the immunoglobulins is conBasis Newton-Raphson minimization ( ABNR) .l Harmonic constraints of 5 kcal/ (mole atom A2) structed from two chains, each of which folds into two @-barrelsthat are noncovalently a ~ s o c i a t e d . ~ ~ were used during the minimization to reduce the shift from the crystal coordinates.26The initial enSix turns between sheet segments form the antigen ergy of the variable domain was 7.85 X lo5 kcal/ binding site, and these turns, the hypervariable mole, and the final energy was -142 kcallmole; the ~ O O P S ,can ~ ~ vary greatly from one immunoglobulin to the next. The remainder of the variable domain root mean square (rms) shift between the two was

a

SAMPLING BY MOLECULAR DYNAMICS

Table I

1849

Seauence of McPC 603 and Hypervariable Loopsa Heavy Chain

Glu Leu Ser ASP Pro Ser Glu Ile Leu Glu Asn Val Ser

LYS Gln Ala TY LYS Asn Ser Ser Leu Thr TYr GlY

Leu Pro Thr Met A% LYS Ala Arg Gln Ala Gly Ala

Val GlY Ser Glu Leu Gly Ser ASP Met Ile Ser G~Y

Glu G~Y GlY TrP Glu Asn Val Thr Asn TYr Thr Thr

Ser Ser Phe Val Trp LYS LYS Ser Ala TYr T ~ P Thr

GlY Leu Thr A% Ile TYr G~Y Gln Leu CYS TYr Val

GlY Arg Phe Gln Ala Thr Arg Ser Arg Ala Phe Thr

GlY Leu Ser Pro Ala Thr Phe Ile Ala Arg ASP Val

10 20 30 40 50 60 70 80 90 100 110 120 122

Ser Glu Gln Asn GlY Ala Arg Phe Glu ASP GlY

Pro ‘4% Ser Phe Gln Ser Phe Thr ASP His Thr

Ser Val Leu Leu Pro Thr Thr Leu Leu Ser LYS

Ser Thr Leu Ala Pro Arg GlY Thr Ala Tyr Leu

10 20 30 40 50 60 70 80 90 100 110 113

Light Chain ASP Leu Met Asn Trp LYS Glu Ser Ile Val Pro Glu a

Ile

Val Val CYS G~Y Gln Leu GlY Ser Ser TYr Thr Lvs

Met Ser LYS Asn Gln Ile Val GlY Val CYS Phe

Thr Ala Ser Gln LYS TYr Pro Thr Gln Gln GlY

Gln GlY Ser LYS Pro Gly ASP ASP Ala Asn Ala

The hypervariable loops used in the dynamical simulation are in boldface.

0.472 A. An C, plot of the variable domain is shown as the background in Figure 1 (a-c ) . The molecular dynamics simulation was performed using a developmental copy of version 16 of the program, CHARMM.l Only the hypervariable loops (see Table I ) were permitted to move, and the coordinates for the fixed atoms were set at the minimized positions determined above. The potential energy computed by CHARMM was used to determine the forces on the atoms. The polar hydrogen representation was used. In this representation, only the hydrogens attached to nitrogen, oxygen, and sulfur atoms that can be involved in hydrogen bonding are included explicitly; nonpolar hydrogens are treated as part of extended atoms. The parameter file was similar to that described in Brooks et al.,’ with the differences listed in the appendix. The shifted dielectric potential was used to compute the electrostatic interaction, and the nonbonded selection cutoff, “off” switching, and “on” switching distances were set to 7.0, 6.0, and 5.0 A, respectively. This relatively short distance was used

to speed up the search procedure. The hydrogenbond interaction was computed with an attractive exponent i ,equal to 12; a repulsive exponent j ,equal to 10; an exponent of 4 on the cos (8A.H.D) factor and no cos (8AA.A.H) factor. The selection, off switching, and on switching distances for the hydrogen-bond energy were 6.0, 5.0, and 3.5 A, respectively. The selection, off switching, and on switching angles were 99.9’, 75.0°, and 45.0°, respectively. Because of the presence of hydrogen atoms in the simulation, the SHAKE27 algorithm was used to make possible a suitable time step. It was applied to the bonds involving hydrogens with a tolerance of With SHAKE, a step size of 0.978 X s was used. In all the simulations, the ratio of fluctuations in kinetic energy to fluctuations in total energy was generally greater than 50 for short time ( 10-step) intervals. Three simulations were performed at different temperatures. The first, a t 500 K, was run for 111.5 ps. The run was started using the minimized coordinates with random velocities for all the atoms, and

1860

BRUCCOLERI AND KARPLUS

Figure 1 This figure contains 3 stereo C, plots comparing the conformations of the antigen combining site at the end of each simulation with the starting coordinates. (a-c) Final coordinates for the trajectories at temperatures of 500,800, and 1500 K. The starting coordinates are drawn with solid lines; light lines are used for the framework residues and heavy lines are used for the antigen combining site. The final simulation coordinates are drawn with heavy dashed lines. The residues in the middle of the sequence for each hypervariable loop are labeled.

SAMPLING BY MOLECULAR DYNAMICS

1851

Figure 1. ( Continued from the previous page. )

thereafter, no energy was added or removed from the system. The second, at 800 K, was run for 101.7 ps and was started the same way as the 500 K run. Again, no energy was added or removed during the simulation. The third run at 1500 K started with the coordinates and velocities from the end of the 800 K run. The velocities were scaled in 10 K increments every 0.0978 ps until a temperature of 1500 K was achieved. Thereafter, the temperature was adjusted to 1500 K if the average temperature for a 0.0978-ps interval deviated by more than 20 K in either direction. This run continued for 76.3 ps. The total energy did not drift much at 500 K, but it did drift downward a t both 800 K and even more so at 1500 K (ignoring the jumps in energy caused by temperature adjustments). These “violations” of the conservation of energy were indicative of a step size that is too large for good numerical integration of the equations of motion. Since we are interested only in searching the conformational space, this was not important, but it would have been if we were studying dynamical properties. For all the runs, the coordinates were saved every 10 steps, yielding 29,600 coordinate sets that occupied approximately 2.5 X 10’ bytes of storage. The simulations required 56 CPU days on a VAX 11/780 without a Floating Point Accelerator.

RESULTS The results presented in this section are used to analyze the utility of high-temperature molecular dy-

namics as a method for conformational search. Ideally, a search method should yield the low-energy conformations that cover the conformational space. Energetics of the Conformations

From each trajectory, we selected three time points, each separated by at least 2 ps, which had the lowest potential energies (see Table 11). As can be seen, with increasing temperature, the potential energy of the system increased corresponding to the existence of more distorted structures. Table I1 also includes the final coordinates in each simulation and the crystallographic coordinates. Each of these coordinates was minimized using ABNR minimization until convergence (i.e., the ABNR step size dropped below 1.0 X Table I1 shows the results of the minimizations. Up to a temperature of 800 K, the dynamics yielded conformations that could be minimized to energies significantly lower than the crystal structure. The structures generated at 1500 K did not minimize to the energies obtained at lower temperatures, although the energies were close to the result from the x-ray structure. The rms shift necessary to converge increased with temperature. However, the shift was small relative to the difference between the dynamics and x-ray structures. A comparison of the torsion angles between samples from the 1500 K run and the crystal structures showed a number of cis-peptide bonds. For example, after 75.8 ps into

1852

BRUCCOLERI AND KARPLUS

Table I1 Minimization of Dynamics Coordinates"

T (K)

Timeb (PS)

Initial Potential Energy

Final Potential Energy

Steps'

0 500

0.0 62.1 78.2 103.7 111.5 55.3 75.6 95.3 101.7 8.4 27.3 75.8 76.3

-8.65 114.89 113.23 111.05 235.22 718.67 699.87 708.57 833.43 2184.61 2172.92 2163.77 2407.74

-627.35 -864.86 -848.95 -858.91 -841.40 -810.46 -835.49 -803.17 -823.38 -648.35 -623.91 -582.78 -610.61

1463 375 316 482 511 940 780 987 652 1292 1361 1455 1396

800

1500

Gradientd

Initial rms Shift'

Final rms Shiftf

Minimum rms Shiftg

0.0383 0.0648 0.0826 0.0183 0.0738 0.0566 0.0152 0.0185 0.0250 0.0391 0.0614 0.0286 0.0903

0.000 4.267 4.328 4.453 4.587 5.932 6.747 6.661 6.407 6.778 6.804 10.327 10.245

1.438 4.300 4.312 4.434 4.607 5.991 6.813 6.685 6.568 6.503 6.565 10.039 9.865

1.438 0.598 0.559 0.599 0.506 1.218 0.922 1.216 0.920 1.696 1.753 1.732 1.888

a The coordinates used for ABNR minimization were selected by either being the crystal coordinates, the last coordinates in the simulation, or the three lowest potential energy coordinates in the simulation separated by more than 2 ps. The minimizations all terminated when the ABNR step size decreased below 1.0 X The time point in the simulation for the coordinates selected. ' The number of ABNR steps needed for convergence. The final value for the rms energy gradient evaluated over all the atomic coordinates. The shift in Angstroms from the minimized x-ray coordinates to the coordinates listed here. The shift from the minimized x-ray coordinates to the coordinates after minimization. The shift that occurs during minimization of the dynamics coordinates.

the 1500 K run; Arg 52, Lys 54, Tyr 62, Ala 64, Ser 65, Thr 106, Tyr 108 in the heavy chain, and Ser 28, Leu 30, Asp 97, End Pro 101 in the light chain, had cis peptides. Clearly, for this temperature, the energy was so high that transitions between cis and trans peptides were possible. Therefore, a temperature of 1500 K is too hot, unless special constraining potentials are added to forbid cis-trans isomerization. Such potentials have been used in simulated annealing employed for x-ray Table I11 shows the energy distribution for a minimized structure from the 800 and the 1500 K

dynamics. Most of the energy differences are in the nonbonded terms. Figure 1 shows the final conformation for each trajectory. As the temperature was increased, the hypervariable loops became more disrupted, and many of the large loops adopted a random coil configuration away from the framework. At a temperature of 1500 K, kT is about 3 kcall mole, which is larger than individual van der Waals, electrostatic, and hydrogen-bond interactions. Thus, structural features of folded proteins, such as secondary structures and the compact arrangement of atoms, are destroyed in a few picoseconds at this

Table I11 Comparison of Energy Contributions*

800 1500

75.6 8.4

11.4 13.1

296.6 304.4

136.7 120.8

23.5 31.0

-686.6 -597.4

-347.3 -286.8

-269.7 -224.9

-835.5 -639.8

a A display of the various contributions to the energy for minimized coordinates from two different temperature runs. E b is the bond energy, EBis the bond angle energy, E, is the torsion angle energy, Ei is the improper torsion angle energy, Evdwis the van der Waals energy, Eeleeis electrostatic energy, E h b is the hydrogen-bond energy, and E is the total energy.

SAMPLING BY MOLECULAR DYNAMICS

2600. 2500. 2400. 2300. 2200.

F 8

2 d

p

1853

II E

d

v

s

5

A

0

1000. 900. ~m (, 800. w2 E 700. + 500. Y r: 5 400. s 0 5: o 300. 0

.A

a

II E

200. 100.

I

I

I

I

I

I

I

I

I

0.

10.

20.

30.

40.

50.

60.

70.

80.

I

1

90. 100. 110.

Time (ps)

Figure 2 The potential energy is plotted vs time at the same scale for all three simulations.

At the start of the 500 and 800 K runs, the average potential energy dropped (see Figure 2 ) . In effect, molecular dynamics served to minimize the structures; this was not true of the 1500 K run. From the lower energies after the long minimizations (see Table 11),we can also see that the structures were being improved energetically. Figure 3 shows the potential energy of coordinates

temperature. The minimizations listed in Table I1 had no significant effect on restoring these prominent features. Quenched molecular dynamics ( 1.5 ps to cool to 0 K ) applied to the structure at 8.4 ps did not affect the secondary structure or the cistrans state of any peptides, although it shifted the coordinates 1.760 8,with respect to the 8.4 ps structure.

I

I

I

I

I

I

I

I

I

I

I

I

-200. a

A

8

-300.

d 0

x

v

A

-400.

F

4

-500.

d

3

$ -600.

4

a .u -700. N

4'2

s

-800. -900.

' 0.

I

I

10.

20.

30. 40.

I

50.

I

60.

I

70.

I

80.

I

90.

I

I

100. 110.

Time (ps)

Figure 3 Coordinates from each of the simulations were selected every 1.956ps and then ABNR minimized for 100 steps. The energies of these coordinate sets were then plotted against time all on the same axes.

1854

BRUCCOLERI AND KARPLUS

from the trajectory that had been minimized for 100 ABNR cycles. This short minimization was insufficient for convergence, but removed much of the strain present in the coordinates from the simulation. Here again, we see the minimization effect of dynamics in the 500 and 800 K simulations. More importantly, all of the structures found at 500 K had energies much lower than the starting coordinates, indicating that the simulation found many low-energy conformations.

I

1

1

I

1

I

I

7.

6.

5. 4.

3. h

4

a 2.

+ 2

T = 1500°K

1.

x

Effectiveness of Search a

9

5.

3. : * 2. 1.

4. 3. 2. 1.

* r’

T

I

I

I

I

I

=

800°K

1

I

I

I

I

I

I

1

I

I

20.

40.

60.

80.

100.

T

At = 1500°K

I

I

6.

4.

3d

T

I

8.

v c

To characterize the effectiveness of molecular dynamics as a search method, we used a number of criteria. We also compared the present results with a uniform sampling of the conformational space

I

=

500°K

(P4

Figure 5 The minimum, average, and maximum values of the rms shift time correlation function f ( t ,t + A t ) are plotted vs At for only the backbone atoms (N , H, C,, C and 0 ) . The other details of the plot are given in the caption to Figure 4.

based on CONGEN.15 A simple analysis of the conformational space surveyed by the simulations involved computing the rms differences between all pairs of coordinates obtained by periodically sampling the trajectory. This defines the displacement autocorrelation function, f ( t , t A t ) , of the form,

+

T 20.

40.

60.

80.

=

500’K 100.

At (P4

Figure 4 The minimum, average, and maximum values of the rms shift time correlation function, f ( t , t At), are plotted vs At for the three different simulation temperatures. The average values are plotted with heavy lines, the minimum values is plotted with light lines under the average, and the maximum values are plotted in light lines over the average.

+

where i varies over the x , y , and z coordinates of all the atoms of the antigen combining site, and n is the number of atoms in the site. We computed the average, minimum, and maximum values o f f ( t , t A t ) as a function of At for the entire combining site (Figure 4) and for just the backbone atoms (Figure 5 1. Examination of the minimum values for f ( t , t A t ) shows that the minimum rms differences be-

+

+

SAMPLING BY MOLECULAR DYNAMICS

tween any pair of conformations separated by a short time are always comparable, and that no conformations are revisited during any of simulations. As the temperature of the simulation was increased, the roughness of the minimum of f ( t , t A t ) increased, consistent with greater motion in the simulation. The average values have similar features and highlight the consistent differences between coordinates in the trajectory that were separated by more than a few picoseconds. The curves for the maximum values of the autocorrelation function for the 500 K simulation reach a plateau of 4.5 A. The appearance of this plateau suggests that all the conformational space accessible at 500 K was within a 4.5 A conformational hypersphere. Of course, this limit may well be exceeded during a longer simulation. A similar, although rougher, plateau was seen at 1500 K with a maximum around 9 A. At 800 K, it appears that a plateau is reached after about 80 ps, although the region over which the maximum distance is constant is too short to be certain. The general behavior observed in these high-temperature simulations is similar to that found in a detailed analysis of a room-temperature simulation of myoglobin.28 The magnitude of the autocorrelation functions can be compared to the deviations between bovine pancreatic trypsin inhibitor and early attempts to predict its f ~ l d i n gIn . ~this ~ ~system, ~~ which has a similar number of residues (58) to the present work ( 7 0 ) , globular but incorrect models of the trypsin

+

p)

2

10.0 8.0 6.0

2 b 4

11

4.0

E

2.0

g

0

3

;.

& vl

inhibitor had rms deviations on the order of 6 A. Thus, the high rms displacements seen at 1500 K reflect the disruption of the compact structure of the hypervariable loops. A comparison of the average and minimum values for f ( t , t A t ) for all the atoms in the combining site (Figure 4) and for just the backbone atoms (Figure 5) shows very similar curves with slightly smaller displacements for the backbone only plot. This demonstrates that the motion in the trajectories involved backbone as well as side-chain motion. The displacement between points on the trajectories and their dynamical average (Figure 6 ) shows a decrease during the first part of the 500 and 800 K simulations. This is consistent with the notion that the simulation performed a minimization in the early part of the run. No such decrease is seen at 1500 K. Another approach to analyzing the conformational space available to the peptide backbone is to examine the backbone torsion angles, 4 and rc/. This should be a relatively good approximation and permits a large reduction in the information needed to describe a conformation, i.e., from 3 coordinates per atom for 696 atoms (2088 values) to 140 torsions. To analyze the conformations, a grid that divided each torsion into a small number of intervals was introduced; i.e., the torsion angle space was divided into many small hypercubes. For each time point of the trajectory, each torsion was assigned to one of the hypercubes. We then counted the number of

+

10.0

-

6.0 II

4.0

E

2.0

10.0 8.0 5: 6.0 II 4.0 E 2.0

g

0.

1855

10. 20. 30. 40. 50. 60. 70. 80. 90. 100. 110.

Time (ps)

Figure 6 The rms deviation of the 696 atoms in the antigen combining site from the average coordinates for each simulation are plotted vs time. All scales are the same.

1856

BRUCCOLERI AND KARPLUS

times each torsion hypercube was found in the trajectory. This procedure was applied to the entire antibody combining site (HV) as well as to the individual hypervariable loops. The results are summarized in Table IV. The first observation is that almost all the conformations were different; i.e., for the entire site ( H V ) , the trajectory sampled the same conformational hypercube only twice, even with a resolution of only 120". Clearly, the dynamical simulation was not stuck in a small local minimum. The second observation is that for shorter loops, many conformations were sampled repeatedly in the simulation; e.g., for the H1 loop at 500 K, some conformational hypercubes were sampled more than 1% of the time. For the H1 and L2 loops, the number of 120" hypercubes visited during the simulations is close to following theoretical model: Since these loops are fixed at each end point, only 2n - 6 torsions angles are independent for each loop, where n is the number of residues per Since the condition of chain closure can have up to about 8 solutions (if planarity of the peptide o angle is assumed) ,31,32 the total number of possible 120" hypercubes for a loop of n residues is 8 3 2 ( n p 3 For ) . the H1 loop, this limit is 648, and for the L2 loop, this limit is 52488.

-

Table IV

Summary of Torsion Angle Histograms" T

Loopb

Finally, we compared the molecular dynamics simulation to another conformational search, CONGEN." CONGEN generates conformations for a segment of polypeptide chain using a uniform sampling of the torsion angle space. Since this process has space and time requirements that grow exponentially with the number of residues, we are practically limited to the H1 and L2 loops. We used a backbone and side-chain grid of 30" for the torsions, a peptide backbone energy cutoff of 2 kcal/mole, and a maximum van der Waals contact of 100 kcall mole. Since the comparisons were with the backbone torsions only, we used the "First" method of sidechain generation, '' which finds the first side-chain conformation that fits in the local environment rather than a complete search of the side-chain torsion angle space for the lowest energy side-chain conformation. CONGEN generated 44 conformations for H1 and 835 conformations for L2. The number of conformations is much smaller than that found with high-temperature molecular dynamics, presumably due to the energy restrictions used in CONGEN. However, what is of greater interest is the extent of the conformational space covered by the two methods.

Grid (ded

Number of Conf'

=

500 K

T

Max Hitsd

Number of Conf'

=

800 K

T

Max Hitsd

=

Number of Conf"

1500 K

Max Hitsd _______

HV (70) H I (5) H2 (21) H3 (11) L l (17) L2 (7) L3 (9)

120 120 120 120 120 120 120

11399 252 10288 4708 9971 878 1878

2 864 7 43 6 689 258

10399 502 9910 7169 9617 1462 4338

2 361 6 15 5 399 84

7800 829 7572 6249 7459 4301 5043

1 207 6 10 5 20 22

HV (70) H1 ( 5 ) H2 (21) H3 (11) L1 (17) L2 (7) L3 (9)

30 30 30 30 30 30 30

11400 3003 11388 10866 11352 9604 10271

1 260 2 4 2 7 7

10400 6001 10397 10207 10393 9694 10029

1 39 2 3 2 5 5

7800 6913 7800 7771 7800 7579 7725

1 8 1 2 1 3 3

~~

* The torsion angles for every saved coordinate set in the three trajectories are classified into torsion angle hypercubes whose size is given by the grid column. The number of residues is given in parentheses. 'T h e number of different hypercubes found for the trajectory. The maximum number of hypercubes is equal to the number of coordinates saved from each trajectory, which are 11400,10400, and 7800, for the 500,800, and 1500 K trajectories, respectively. T h e number of conformations in the most frequently occupied hypercube.

SAMPLING BY MOLECULAR DYNAMICS

Table V

1857

Closest Approach of CONGEN Conformations to the Dynamicsa Number of CONGEN Conformations with Minimum rms Torsion Difference Between Total

Run

Loop

0"-20"

20"-40"

4Oo-6O0

6Oo-8O0

T = 500 K T = 800 K T = 1500 K

H1

3 2 3

2 3 24

12 16 16

26 20 1

1 3 0

44

T = 500 K T = 800 K T = 1500 K

L2

0 0 0

23 4 48

185 190 573

517 599 214

110 42 0

835

80"-100"

a The distribution of closest distances of each CONGEN conformation to any dynamics conformation is listed. The distances are computed by taking the rms difference in torsion angles. For example, for the L2 loop, there were 4 CONGEN conformations that were between 20' and 40" rms distance to a t least one of the dynamical conformations a t 800 K.

We compared the two sets of conformations on the basis of the backbone torsion angles, C#J and $. The comparison took one set of conformations, A, and determined the minimum torsion angle rms distance to any conformation of other set, B. If a large percentage of the conformations in set A were close to conformations in set B, this would indicate that a large percentage of the space covered by A was in set B. When we calculated how close each CONGEN conformation was to any dynamics-generated conformation (Table V) , we saw very little matching. The most likely explanation for the limited matching was that molecular dynamics permitted variations in bond lengths and especially bond angles that CONGEN could not permit, and this greatly increased the number of possible conformation^.^^'^^ The rms fluctuations in bond angles computed from trajectory samples every 0.978 ps were 4.6" at 500 Table VI

K, 5.8" at 800 K, and 7.8" at 1500 K. One peculiar aspect of the distribution of matches was the improvement seen when the temperature reaches 1500 K. The higher temperatures appear to cover conformational space more thoroughly and thereby approached closer to the CONGEN generated conformations. Table VI shows the converse comparison, namely how close each dynamics conformation was to some CONGEN-generated conformation. In the case of the H1 loop at 500 K, most dynamics coordinates were close to some CONGEN-generated conformations, and thus the space spanned by the dynamics conformations for H1 at 500 K was a subspace of the CONGEN conformations. However, for the other loop and for higher temperatures, the overlap was much smaller, and the results are similar to those in Table V. It is not clear why the dynamics

Closest Approach of Dynamics Conformations to CONGEN" Number of Dynamics Conformations with Minimum rms Torsion Difference Between

Run

Loop

T

= 500 K T = 800 K T = 1500 K

H1

T = 500 K

L2

T = 800 K T = 1500 K

0'-20"

80'-100

"

Total

20"-40"

40"-60"

60"-80"

1831 358 9

8521 302 1577

1048 2962 3149

0 6778 3033

0 0 32

11400 10400 7800

7 6 0

500 24 269

3332 2346 3062

7561 8024 4316

0 0 153

11400 10400 7800

a The distribution of closest distances of each dynamics conformation to some CONGEN conformation is listed. The distances are computed by taking the rms difference in torsion angles. For example, for the H1 loop, there were 3149 dynamics conformations at 1500 K that were between 40' and 60' rms distance of some CONGEN conformation. The table is analogous to Table V except that dynamics and CONGEN have been reversed.

1858

BRUCCOLERI AND KARPLUS

conformations of the H1 loop a t 500 K span a subspace of the CONGEN-generated conformations, but that the L 2 loop did not. The cis-trans peptide isomerization contributed some error into the interpretation of these two tables. Our comparison of just the 4 and II/ torsion angles involved the assumption that if these torsion angles were similar for two conformations, the conformations were similar. However, if a peptide was cis in one and trans in the other, then this assumption could be violated. Thus, at 1500 K, the comparisons to the CONGEN conformations may not be valid. Given the view of conformational space as a multidimensional torsion angle space, we have developed several measures of the size of the conformational space explored during the simulation. In all of these approaches, we assumed we could treat the concept of distance in a Euclidean fashion even though torsion angle space is closed in all dimensions as angles have a period of 360". This simplifying assumption about the computation of distance permits analysis using straightforward methods. Also, for the 1500 K runs, we did not account for the cis peptides, which resulted in an underestimate of the distance. Given these caveats, our definition of distance is

where ti and tj are various times, and k varies over the torsions in the set of interest. For example, if we are concerned with the H1 loop in McPC 603, then k ranges from $H-31 to II/H-35. The first and second measures of the conformational space are based on the matrix of second moments of the torsion angle distribution. The elements of this symmetric matrix are given by the following expression:

where i and j are indices for particular backbone torsions, and t k specifies a particular time in the simulation. The averages over time were computed by summing the sine and cosine of each torsion angle over the trajectory, and then taking the arc tangent of the sums to recover the average. This procedure avoided the difficulties associated with averaging a periodic variable. When computing the differences, the angles were adjusted modulo 360" so as to limit the magnitude of the differences to less than 180". The eigenvectors of this matrix specify a rotation

that, when applied to the standard basis vectors of the torsion space, will result in zero cross variances (off-diagonal matrix elements) when computed by Eq. ( 2 ) . The division of the eigenvalues of this matrix by the number of time points and then computation of the square root yielded the principal values of the distribution. These principal values defined a hyperspace ellipsoid that measures the conformational space in a way analogous to the calculation of the moment of inertia ellipsoid as a measure of the mass distribution of a physical object. The first measure of the conformational space is the volume of this ellipsoid. The maximum possible value of this volume is determined by the maximum possible values for the principal values that would occur with a random torsion angle distribution. Comput,er trials suggested that the largest principal values would be around 100" each. The second measure is the square root of the sum of the square of the principal values, which is the diagonal length of a hyperspace rectangle whose sides are given by the principal values. Since the computation of the principal values (which represent angular differences) yielded a few values greater than 180", these measures were slight overestimates. Its maximum possible value also depends on the principal values as above. Given maximum principal values of loo", then the maximum of this measure would be 100 X ( 2 N )'I2 where N is the number of residues. The third measure of the size of the conformational space is the maximum distance between two conformations. Since the number of conformations is so large, it was impractical to compute this by considering every pair. Thus, we estimated the largest distance by a simple iterative algorithm. The first step in the algorithm was the selection of an "anchor" conformation. Then we searched for the most distant conformation from the anchor conformation. Next, we exchanged the anchor conformation with the most distant one, and we repeated the search until there was no further increase in the maximum distance. Since the choice of the initial anchor could affect the maximum distance, we tried 11different starting points selected roughly uniformly over the time of the runs, and the largest distance was used. We verified this heuristic approach by a brute force computation of the maximum distance for the two CONGEN runs as well as the H1 and L2 loops at 1500 K. In all four cases, the heuristic procedure found the correct distance. Table VII shows these three measures applied to the entire antibody combining site as well as the

SAMPLING BY MOLECULAR DYNAMICS

Table VII

1859

Sizes of Conformational Space"

Loop

Number of Residues

HV

70

H1

5

H2

21

H3

11

Ll

17

L2

7

L3

9

Run

Volumeb

Diagonal Length"

Maximum Distanced

T = 500 K T = 800 K T = 1500 K CONGEN" T = 500 K T = 800 K T = 1500 K T = 500 K T = 800 K T = 1500 K T = 500 K T = 800 K T = 1500 K T = 500 K T = 800 K T = 1500 K CONGEN' T = 500 K T = 800 K T = 1500 K T = 500 K T = 800 K T = 1500 K

6.4 x 101°4 2.0 x 1.3 x 1 0 1 ~ 5 2.4 x loi3 1.2 x 1012 1.1 x 1014 5.5 x 1016 3.1 x 1.3 x 1.3 X 1.4 x 1 0 2 ~ 2.3 x 9.5 x 6.0 x 1o4l 6.4 x 1 0 4 ~ 3.2 x 1.5 X loz2 1.4 X 10" 5.1 X 10" 2.8 x 1oZ5 3.9 x 1020 6.2 x 1 0 2 ~ 2.5 X 10''

431.4 612.5 824.3 161.9 77.3 118.6 184.3 205.1 339.6 472.5 148.0 209.8 335.4 294.7 365.7 415.4 240.8 121.4 110.1 277.4 120.5 236.4 245.5

994.3 1109.9 1238.8 400.9 264.4 331.0 440.0 553.0 683.9 791.8 452.4 562.5 616.1 687.2 693.1 725.0 481.5 404.0 395.5 540.1 413.4 501.7 531.5

Three different measures of the size of the conformational space for the individual loops in the antigen combining site as well as the entire site are listed (see text and below). T h e volume is a multidimensional ellipsoid whose semiaxes are given by the principal values of the torsion space distribution. The units of this measure are degrees raised to the power of the number of residues times 2. ' The diagonal length is the square root of the sum of the squares of the principal values and has units of degrees. T h e largest torsion distance found between two conformations in the space which also has units of degrees. Conformations generated by the CONGEN program, see text for more details.

individual hypervariable loops. In addition, we list the measures as applied to the CONGEN calculations. In general, the size measurements increased with increasing temperature and the number of residues in the loop. However, there were some exceptions to this. The diagonal length and maximum length measurement for the L1 loop a t 500 and 800 K were bigger than the H2 loop, which is longer. Presumably, this was due to the flexibility of the L1 loop, as evidenced by its very low electron density in the crystallographic Second, the conformational space of the L2 loop at 500 and 1500 K had nearly the same diagonal and maximum length measurement as the L3 loop, which is two residues longer. A likely cause for this disparity is that the L3 loop is more buried than the L2 loop. The average accessible surface36for the L2 loop is 5.99 Az/atom

whereas for the L3 loop it is 2.89 A2/atom. Also, another contributing factor was the presence of a proline in the L3 loop, which restricted the freedom of the 4 torsion in that residue. Third, the measures of the conformational space size and the rms deviations from the crystal for the L2 loop were approximately the same for both 500 and 800 K. Tables V and VI show that the comparison to CONGEN-generated conformations were roughly the same for the two temperatures. Apparently, this loop was in a potential well that required temperatures over 800 K to surmount on the time scale of these simulations. CONGEN generated larger conformational spaces for the loops to which it was applied than dynamics at 500 K, but not at higher temperatures. If the CONGEN runs had been made with the

1860

BRUCCOLERI AND KARPLUS

neighboring loops removed so as to account for the simultaneous motion of the other loops, its coverage of the conformational space of these loops would have been larger still.

CONCLUDING DISCUSSION The failure to see any repeats in the conformational space of the hypervariable loops either by torsion angle (Table IV) or rms measurements (Figures 4 and 5) suggests that the simulation was not run long enough for a thorough search of the space. The method of Dill,17 for estimating the number of low free energy states of a globular protein with n residues, 1.4",yields 1.7 X l0'"or the 70 hypervariable residues. Although it is an underestimate for short loops, it is consistent with the idea that our simulations were too short for searching the entire hypervariable loop system. However, for the individual loops, repeats in the conformational space were seen, with more repeats occurring in the smaller loops, and the number of sampled torsion angle hypercubes was close to theoretical expectations for complete coverage. Thus, this method would be suitable for conformational search for short polypeptide segments, on the order of a dozen residues. Although increasing the temperature of the simulation increased the size of the conformational space searched, the simulation temperature was limited by the requirement of maintaining secondary structure and other protein features. In particular, 1500 K is too hot, as the binding site loses its compactness and many cis peptides form. Constraints would be necessary to maintain these features at such temperatures. The molecular dynamics simulations at 500 and 800 K generated a large number of conformations that could be minimized to energies below that of the starting coordinates. Thus, this method can surmount the barriers surrounding local minima. Even a t room temperature, molecular dynamics can be used to improve the energetics of a macromolecular The comparison of the sizes of the conformational spaces sampled by CONGEN to the molecular dynamics a t lower temperatures illustrate that the more direct conformational search methods explore more of the space, whereas at 1500 K, they do not. However, the use of CONGEN is not feasible for complete searches of systems larger than eight or nine residues, although partial searches of larger systems can be p e r f ~ r m e d . ~ ~

It is important to note that the difficulties in using molecular dynamics to thoroughly sample conformational space are due largely to the shape of the potential surface. There are many local minima, all of which have similar energies. The number of such local minima is much greater than can be sampled in a reasonable time. However, other potential surfaces can be effectively searched using this method e.g., molecular dynamics is a very useful tool in refinement of NMR1' or x-ray crystallographic structures, 20,21 where additional constraints are present. In summary, high-temperature molecular dynamics is a reasonable conformational search method when it is restricted to peptide segments on the order of a dozen residues and perhaps more-for a cyclic system. Given present day serial computer technology, it would take too long for an adequate search of a system the size of the antigen combining site or of a small protein. However, massively parallel machines could be very useful for such conformational searches.

APPENDIX: ENERGY PARAMETERS The following tables present the parameters used in the dynamics simulations which differ from the set given in Ref. 1. (The parameter set name is PARAM1.) Bond Parameters Eb= kb(r- ro)2 kB

Bond C-CH1E C-CH2E C-CH3E C-N C-NH1 C-NH1E C-NH2 C-NH2E C-NP C-NR

c-0 c-oc

CHlE-N CHlE-NH1 CHlE-NH1E CHlE-NH2 CHlE-NH2E CHlE-NH3 CHlE-NH3E CH2E-N CH2E-NH1 CH2E-NH1E

(kcal/mole-A2) 405.0 405.0 405.0 471.0 471.0 471.0 471.0 471.0 471.0 471.0 580.0 580.0 422.0 422.0 422.0 422.0 422.0 422.0 422.0 422.0 422.0 422.0

ro

(A)

1.53 1.53 1.53 1.305 1.305 1.305 1.305 1.305 1.305 1.305 1.215 1.22 1.49 1.49 1.49 1.49 1.49 1.49 1.49 1.49 1.49 1.49

SAMPLING BY MOLECULAR DYNAMICS

Bond Parameters Eb= kb(r- ro)2(Continued)

1861

Bond Angle Parameters EB= ka(O- O0)' (Continued)

kB

Bond

(kcal/mole-A*)

ro

(A) Bond Angle

CH2E-NH2 CH2E-NH2E CH2E-NH3 CH2E-NH3E CT-N CT-NC2 CT-NH1 CT-NH2 CT-NH3

422.0 422.0 422.0 422.0 422.0 422.0 422.0 422.0 422.0

Bond Angle Parameters Ea

=

1.49 1.49 1.49 1.49 1.49 1.49 1.49 1.49 1.49

CH3E-C-NH1 CH3E-C-NH1E CH3E-C-0 CH3E-NHl-H CRlE-C-NH1 CRlE-C-NH1E NH1-C-0 NHlE-C-0 NH2-C-0 NH2E-C-0

20.0 20.0 85.0 35.0 65.0 65.0 65.0 65.0 65.0 65.0

116.0 116.0 120.0 110.0 116.5 116.5 121.5 121.5 122.0 122.0

ka(OTorsion Angle Parameters E@= 1 k, I

klJ

Bond Angle

(kcal/mole-radian')

Oo (degrees) Torsion Angle

C-CHlE-N C-CHlE-NH1 C-CHlE-NH1E C-CHlE-NH2 C-CHlE-NH2E C-CHlE-NH3 C-CHlE-NH3E C-CH2E-NH1 C-CH2E-NH1E C-CH2E-NH2 C-CH2E-NH2E C-CH2E-NH3 C-CHBE-NH3E C-N-CH1E C-N-CH2E C-NH1-CH1E C-NH1-CH2E C-NH1-CH3E C-NH1-CR1E C-NH1-H C-NHlE-CH1E C-NHlE-CH2E C-NHlE-CH3E CHlE-C-N CHlE-C-NH1 CHlE-C-NH1E CHlE-C-0 CHlE-C-OC CHlE-NH1-H CH2E-C-N CH2E-C-NH1 CH2E-C-NH1E CHPE-C-NH2 CHPE-C-NH2E CH2E-C-0 CH2E-NH1-H CH3E-C-N

60 (degrees)

k8 (kcal/mole-radian')

45.0 45.0 45.0 45.0 45.0 45.0 45.0 70.0 70.0 70.0 70.0 70.0 70.0 80.0 80.0 77.5 77.5 77.5 60.0 30.0 77.5 77.5 77.5 20.0 20.0 20.0 85.0 85.0 35.0 20.0 20.0 20.0 20.0 20.0 85.0 35.0 20.0

116.0 111.5 111.5 109.5 109.5 107.5 107.5 116.5 116.5 109.5 109.5 109.5 109.5 121.0 127.5 115.0 116.0 116.0 108.0 110.0 115.0 116.0 116.0 115.5 115.5 115.5 120.5 118.5 110.0 116.0 116.5 116.5 119.0 119.0 119.5 110.0 120.0

CRlE-C-C-CRlE CRlE-C-C-C CRlE-C-C-NH1 CRlE-C-C-NH1E X-CHlE-N-X X-CHlE-NH1-X X-CHlE-NHlE-X X-CHlE-NH2-X X-CHlE-NH3-X X-CHlE-OH1-X X-CH2E-N-X X-CH2E-NH1-X X-CH2E-NHlE-X X-CH2E-NH2-X X-CH2E-NH3-X X-CH2E-OH1-X

-

k, (kcal/mole) 35.0 10.0 10.0 10.0 1.8 1.8 1.8 1.8 1.8 1.8 1.8 1.8 1.8 1.8 -1.8 1.8

k,cos(n4)

n

2.0 2.0 2.0 2.0 3.0 3.0 3.0 2.0 3.0 3.0 3.0 3.0 3.0 2.0 3.0 3.0

We thank Dr. Edgar Haber at Massachusetts General Hospital (MGH) for his generous support; John Newel1 and Dave Murphy for the use of the Cardiac Computer Center a t MGH; Dr. David States and Robert Beach for the use of their graphics programs; and Dr. Jiri Novotny for suggestions, discussions, and the contour program. We also thank Dr. Bill Hoffman for suggesting the comparison of the dynamics and CONGEN results. We are grateful to Dr. David Davies for providing the coordinates of McPC 603. This work was funded in part by grants from the National Science Foundation and the National Institutes of Health.

REFERENCES 1. Brooks, B. R., Bruccoleri, R. E., Olafson, B. D., States, D. J., Swaminathan, S. & Karplus, M. ( 1983) J. Comput. Chem. 4, 187-217.

1862

BRUCCOLERI AND KARPLUS

2. Warshel, A. & Levitt, M. (1974) QCPE 11, 247. 3. Browman, M. J., Carruthers, L. M., Kashuba, K. L., Momany, F. A., Pottle, M. S., Rosen, S. P., Rumsey, S. M., Endres, G. F. & Scheraga, H. A. (1975) QCPE 11,286. 4. Huler, E., Sharon, R. & Warshel, A. (1977) QCPE 11,325. 5. Potenzone, R. Jr., Cavicchi, E., Weintraub, H. J. R. & Hopfinger, A. J. (1977) Comput. Chem. 1,187-194. 6. Browman, M. J., Burgess, A. W., Dunfield, L. G., Rumsey, S. M., Endres, G. F. & Scheraga, H. A. (1978) QCPE 11, 361. 7. Williams, D. E. (1979) QCPE 11,373. 8. Allinger, N. L. & Yuh, Y. A. (1981) QCPE 13, 395. 9. Weiner, P. K. & Kollman, P. A. ( 1981) J . Comput. Chem. 2, 287-303. 10. Metropolis, N., Rosenbluth, A., Rosenbluth, M., Teller, A. & Teller, E. (1953) J . Chem. Phys. 21, 1087-1092. 11. Pangali, C., Rao, M. & Berne, B. J. (1978) Chem. Phys. Lett. 55,413-417. 12. Fine, R. M., Wang, H., Shenkin, P. S., Yarmush, D. L. & Levinthal, C. (1986) Proteins 1, 342-362. 13. Kirkpatrick, S., Gelatt, C. D. & Vecchi, M. P. (1983) Science 220, 671-680. 14. Purisima, E. 0. & Scheraga, H. A. (1986) Proc. Natl. Acad. Sci. USA 83, 2782-2786. 15. Bruccoleri, R. E. & Karplus, M. (1987) Biopolymers 26,137-168. 16. Ramachandran, G. N., Ramakrishnan, C. & Sasisekharan, V. (1963) J . Mol. Biol. 7,95-99. 17. Dill, K. A. ( 1985) Biochemistry 24, 1501-1509. 18. McCammon, J. A., Gelin, B. R. & Karplus, M. (1977) Nature 267, 585-590. 19. Clore, G. M., Gronenborn, A. M., Briinger, A. T. & Karplus, M. (1985) J. Mol. Biol. 186,435-455. 20. Briinger, A. T., Kuriyan, J. & Karplus, M. ( 1987) Science 235,458-460. 21. Briinger, A. T., Karplus, M. & Petsko, G. A. (1989) Acta. Cryst. A 45, 50-61. 22. Segal, D. M., Padlan, E. A., Cohen, G. H., Rudikoff,

S., Potter, M. & Davies, D. R. (1974) Proc. Natl. Acad. Sci. USA 71,4298-4302. 23. Novotny, J., Bruccoleri, R. E., Newell, J., Murphy, D., Haber, E. & Karplus, M. (1983) J. Biol. Chem. 258,14433-14437. 24. Kabat, E. A., Wu, T. T., Reid-Miller, M., Perry, H. M. & Gottesman, K. S. (1987) Sequences of Proteins of Immunological Interest, U.S. Department of Health and Human Services, Public Health Service, National Institutes of Health. 25. Chothia, C., Lesk, A. M., Levitt, M., Amit, A. G., Mariuzza, R. A., Phillips, S. E. V. & Poljak, R. J. (1986) Science 233, 755-758. 26. Bruccoleri, R. E. & Karplus, M. (1985) J . Comput. Chem. 7, 165-175. 27. van Gunsteren, W. F. &'Berendsen, H. J. C. (1977) Mol. Phys. 34, 1311-1327. 28. Elber, R. & Karplus, M. (1987) Science 235, 318321. 29. Levitt, M. & Warshel, A. (1975) Nature 253, 694698. 30. Hagler, A. T. & Honig, B. (1978) Proc. Natl. Acad. Sci. USA 75, 554-558. 31. Go, N. & Scheraga, H. A. (1970) Macromolecules 3, 178-187. 32. Bruccoleri, R. E. & Karplus, M. (1985) Macromolecules 18, 2767-2773. 33. van Gunsteren, W. F. & Karplus, M. (1981) Nature 293, 677-678. 34. Kushick, J. N. & Karplus, M. (1981) Macromolecules 14, 325-332. 35. Satow, Y., Cohen, G. H., Padlan, E. A. & Davies, D. R. (1986) J . Mol. B i d . 190,593-604. 36. Lee, B. & Richards, F. M. (1971) J . Mol. Biol. 55, 379-400. 37. Reid, R. H., Hooper, C. A. & Brooks, B. R. (1989) Biopolymers 28,525-530. 38. Bruccoleri, R. E., Haber, E. & Novotny, J. ( 1988) Nature 335, 564-568. Received June 23, 1989 Accepted November 10, 1989