This article was downloaded by: [Rutgers University] On: 12 April 2015, At: 09:22 Publisher: Taylor & Francis Informa Ltd Registered in England and Wales Registered Number: 1072954 Registered office: Mortimer House, 37-41 Mortimer Street, London W1T 3JH, UK
Journal of Biomolecular Structure and Dynamics Publication details, including instructions for authors and subscription information: http://www.tandfonline.com/loi/tbsd20
Using Molecular Dynamics Simulations on Crambin to Evaluate the Suitability of Different Continuum Dielectric and Hydrogen Atom Models for Protein Simulations Rick L. Ornstein
a
a
Pacific Northwest Laboratory Molecular Science Research Center , Richland , Washington , 99352 Published online: 21 May 2012.
To cite this article: Rick L. Ornstein (1990) Using Molecular Dynamics Simulations on Crambin to Evaluate the Suitability of Different Continuum Dielectric and Hydrogen Atom Models for Protein Simulations, Journal of Biomolecular Structure and Dynamics, 7:5, 1019-1041, DOI: 10.1080/07391102.1990.10508543 To link to this article: http://dx.doi.org/10.1080/07391102.1990.10508543
PLEASE SCROLL DOWN FOR ARTICLE Taylor & Francis makes every effort to ensure the accuracy of all the information (the “Content”) contained in the publications on our platform. However, Taylor & Francis, our agents, and our licensors make no representations or warranties whatsoever as to the accuracy, completeness, or suitability for any purpose of the Content. Any opinions and views expressed in this publication are the opinions and views of the authors, and are not the views of or endorsed by Taylor & Francis. The accuracy of the Content should not be relied upon and should be independently verified with primary sources of information. Taylor and Francis shall not be liable for any losses, actions, claims, proceedings, demands, costs, expenses, damages, and other liabilities whatsoever or howsoever caused arising directly or
indirectly in connection with, in relation to or arising out of the use of the Content.
Downloaded by [Rutgers University] at 09:22 12 April 2015
This article may be used for research, teaching, and private study purposes. Any substantial or systematic reproduction, redistribution, reselling, loan, sub-licensing, systematic supply, or distribution in any form to anyone is expressly forbidden. Terms & Conditions of access and use can be found at http://www.tandfonline.com/page/terms-and-conditions
Journal ofBiomolecular Structure & Dynamics, ISSN 0739-1102 Volume 7. Issue Number 5 (1990), ®Adenine Press (1990).
Using Molecular Dynamics Simulations on Crambin to Evaluate the Suitability of Different Continuum Dielectric and Hydrogen Atom Models for Protein Simulations Downloaded by [Rutgers University] at 09:22 12 April 2015
Rick L. Ornstein Pacific Northwest Laboratory* Molecular Science Research Center Richland, Washington 99352 Abstract Molecular dynamics simulations of enzymes with enough explicit waters of solvation to realistically account for solute-solvent interactions can burden the computational resources required to perform the simulation by more than two orders of magnitude. Since enzyme simulations even with an implicit solvation model can be imposing for a supercomputer, it is important to assess the suitability of different continuum dielectric models for protein simulations. A series of 1DO-picosecond molecular dynamics simulations were performed on the X-ray crystal structure of the protein crambin to examine how well computed structures, obtained using seven continuum dielectric and two hydrogen atom models, agreed with the X-ray structure. The best level of agreement between computed and experimental structures was obtained using a constant dielectric of2 and the all-hydrogen model. Continuum dielectric models of 1, 1*r, and 2*r also led to computed structures in reasonably good agreement with the X-ray structure. In all cases, the all-hydrogen model gave better agreementthan the united atom model, although, in one case, the difference was not significant. Dielectric models of 4, 80, and 4*r with either hydrogen model yielded significantly poorer fits. It is especially noteworthy that the observed trends did not semiquantitatively converge until about 50 picoseconds into the simulations, suggesting that validation studies for protein calculations based on energy minimizations or short simulations should be viewed with caution.
Introduction Starting with the pioneering molecular dynamics simulations on pancreatic trypsin inhibitor(l,2), the dramatic progress in understanding the dynamics of proteins has been significantly fueled by computer simulations based on classical mechanics and dynamics (3). Even with today's supercomputers, protein dynamics simulations of any significant duration can totally consume available resources for all but small proteins. Electrostatic interactions are of crucial importance in protein folding and stability as evidenced by their pH dependence (4); the calculation of this component *Pacific Northwest Laboratory is operated for the U.S. Department of Energy by Battelle Memorial Institute under Contract DE-AC06-76RLO 1830.
1019
1020
Ornstein
Downloaded by [Rutgers University] at 09:22 12 April 2015
represents a major scientific and practical problem. Charge-charge interactions are calculated with either a microscopic description of the solvent or a continuum dielectric model. In theory, the microscopic treatment is expected to be superior; but it demands much more computer time (and memory), and current microscopic models do not always reproduce the dielectric properties of even liquid water (5 and references therein). Continuum models are thought to incorporate the essential qualitative electrical properties of the systems being studied and are much more computationally practical. A number of force fields have been proposed and used in protein simulations (6-11 ). Some researchers prefer to use a constant dielectric value in the Coulombic term, while others prefer a distance-dependent dielectric. The theoretical or experimental justification for either choice is subject to debate (6-18). In general, however, force field development and validation are usually applied to small molecular fragments and amino acids and not to proteins. Since proteins possess unique properties such as cooperativity and "buried" residues that their constituent amino acids do not possess, we believe that extensive force field validation studies must also be performed on whole proteins and not just on their constituent moieties. The interior of a protein molecule is generally thought to be a low dielectric medium, while the exterior surface is in a higher dielectric medium (17). Our strategy has been to implicitly account for the solvent by adjusting the effective dielectric value in the Coulombic term of the potential energy and to screen down to zero the net charge on any side chains of Arg, Asp, Glu, and Lys residues, as well as N-and C-termini. This type of approach has also been taken by others (19-22). By screening ionic residues, which tend to be solvated on the protein surface, it is reasonable to anticipate that simple continuum models may overcome known defects (5) involving relative overestimation ofelectrostatic interactions at the protein surface and underestimation of electrostatic interactions at the interior. The secondary structure and hydrogen bonding of crambin in solution have recently been investigated using two-dimensional nuclear magnetic resonance (NMR) (23) and were essentially identical to those of the crystal structure (24,25). Because of its small size and well-known structure, the protein crambin plays an important role as a model system in both experimental and tl: ·etical methods development. Crambin has been used to test potential energy funl-.10ns used in energy minimization (26) and restrained molecular dynamics simulations using simulated NMR distance restraints (27-29). Whitlow and Teeter(26) performed a large series of minimizations on the X-ray crystal structure of crambin with different force field conditions, using the program AMBER (1 0). This noteworthy effort suggested that a continuum dielectric value of 4*r leads to the best agreement between computed and X-ray structures. Because of the existence of a large number of minima in the multidimensional conformational energy hypersurface of proteins (30), energy minimization calculations starting from even an X-ray structure may still locate a minimum significantly less stable energetically than a nearby structure. Such a situation could complicate the interpretation of the type of analysis performed by Whitlow and Teeter (26). Therefore, we reexamined the question of what continuum dielectric and hydrogen-atom
Molecular Dynamics on Crambin
1021
model best accounts for the X-ray structure of crambin based on a set of 100 picosecond molecular dynamics simulations. Results using the DISCOVER program (31) suggest different conclusions and indicate that the earlier findings based on energy minimizations alone are difficult to interpret because of the multiple-minima problem.
Downloaded by [Rutgers University] at 09:22 12 April 2015
Methods The heavy atom coordinates for crambin, with diffraction data to 1.5 A(24), were taken from the October 1988 release version of the Brookhaven Protein Data Bank (32). All hydrogens were added in standard fashion or according to the united-atom approximation. In the latter procedure, nonpolar hydrogen atoms are omitted and the van der Waals radius of the nonhydrogen atoms to which they are bonded is increased accordingly. For proteins, the united-atom approximation reduces the number of atoms by about 30%, reduces the calculation time per iteration by about 50%, and reduces the total number of iterations required to obtain a specific convergence criterion in energy minimization calculations (33). The starting structure for all calculations was the X-ray structure for crambin. Prior to any molecular dynamics simulations, the starting structure was "completely energy minimized" using the following two-step procedure: minimize using the steepest descents method until the maximum derivative was less than 2.0 kcal/Aand then minimize using the conjugate gradients method until the maximum derivative was less than 0.02 kealiA All calculations were performed using the DISCOVER molecular simulation software package version 2.21, running on a 1-by-8 CRAY XMP computer. The default" consistentvalence force field" was used with no cross terms and with a simple harmonic potential for the bond-stretching terms; details of the force field and parameters are available elsewhere (7,33,34, and references therein). All nonbonded interactions were computed, i.e., an infinite cutoff distance was used. The acidic and basic residues as well as the N- and C-termini were assigned neutral partial charges as available within DISCOVER Solvent was not explicitly treated. The root mean square (r.m.s.) deviations reported between computed and X-ray structures were calculated for all nonhydrogen atoms of the backbone between Cys-3 and Cys-40; residues 1, 2, and 41-46 were not included in the r.m.s. calculations to minimize the influence of the more flexible terminal residues. The r.m.s. deviations were computed after least squares fitting to the specified backbone heavy atoms. All simulations used an integration time step of 1 femtosecond. All trajectories were carried out using a constant temperature bath (35). Each 300K dynamics simulation was divided into 20 serial segments 5.0 picoseconds long plus complete minimization. Each 5.0 picosecond segment at 300K consisted of the following sequence of steps: 1) initialize dynamics for0.5 picosecond at 300K with a randomization of the atomic velocities; 2) resume dynamics for an additional4.5 picoseconds at 300K; 3) minimize using steepest descents method until the maximum derivative was less than 2.0 kcal/ A; 4) minimize using conjugate gradients method until the maximum derivative was less than 0.02 kcal/A, and 5) compute r.m.s. deviation between the minimized structure and the X-ray structure. The 100-picosecond simulations using the united-
1022
Ornstein a)
United-Atom/De=1
380
:=0 E
-•
3-5
•..> I
>< 370
2 E
u
.Mi
..•
:
>a 360
'• '•
~
20
b)
40
60
0
E ::::
•u
0
•
Q
0 II! E 100 .;
80
....
United-Atom/De=2
4(
4-
400
.....
>•
~
350
0
c
...v· \ .'.....__,.'''
~ .~
w Downloaded by [Rutgers University] at 09:22 12 April 2015
0
..
;
~:
c::
-.
. -. •> I
>
- 370 a
..• c::
w
2
·;•
•
'
360
;
:'
Q
100
e .;
.,;
350
0
20
c)
40
60
80
Unlted-Atom/De:4
=0 E ::::
•u
•.
370
s,c
360
4
350
.Mi
......
.. •
.·....:, . ..
w
E
0
3 ;0
~~
...
'•
320
.... c
.·~
340
a>- 330 c::
-
-c
6 ..... >
380
2
>• & .,;
310 0
20
40 60 Time (ps)
80
e
1 100 ..:
Figure 1: Dynamics simulations at 300Kforprotein crambin using the united-atom approximation. The energy versus simulation time is plotted with a dashed line, while the r.m.s. deviation of the computed structure from the X-ray structure is plotted with a solid line. Simulations were repeated with different continuum dielectric models: a) 1, b) 2, c) 4, d) 80, e) 1*r, f) 2*r, and g) 4*r.
Molecular Dynamics on Crambin d)
1023
United-Atom/De:SO
360 - - - - - - - - - r - - - . , . . . - - r 6 350 5
340
4
330
..•c,.. w Downloaded by [Rutgers University] at 09:22 12 April 2015
Q
320
~
..>< E -. >-
Ill
0
280
E
:::: Ill
u
2
270
.IC ....,
c .!!
->
>- 260
.. Q
~ w
0
II
•
250
Q
20
0
f)
40
60
80
United-Atom/De=2*r
;(
34o . . . - - - - - - - - - - - - - - r s •
=0
E
>-
..•. ...:·· 1,
...Ill
··..- /\
330
4
>
•
~·
Q
310~~~~--r-~~~-,----+
20
40 Time
...0 c
~
\:• : ~ : \ .... .: ' ....,. 't
0
4
E
'•
>aa
-
United-Atom/De=4*r
60
(ps)
80
1
100
•e .:
Figure 1 continued
atom approximation and the all-hydrogen atom model required about 4.0 and 7.5 CPU hours, respectively.
Results and Discussion I. Effect ofContinuum Dielectric and Hydrogen-Atom Model on r.m.s. Deviation Between Computed Structures and X-Ray Structure
Starting from the X-ray structure of crambin, 100-picosecond molecular dynamics simulations were computed at 300K for seven dielectric models: constant values of 1, 2, 4 and 80 and linear models of 1*r, 2*r, and 4*r. The results obtained using the united-atom approximation are plotted in Figure 1, and the results obtained using the all-hydrogen atom model are shown in Figure 2. The data are shown in the form of double-Y plots with energy and r.m.s. deviations along theY-axis and simulation time along the X-axis. The data points correspond to fully minimized structures (see above for details) at time zero and at 5-picosecond intervals during the simulation. In general, the second half of the simulations (from 50 to 100 picoseconds) tended to have equilibrated. The plots ofheavy atom mainchain r.m.s. deviations (between Xray and computed structures) versus time for the 14 simulations in Figures 1 and 2 are displayed in Figure 3.1t is evident from these data that some potential functions resulted in good agreement with the X-ray structure, while some potential functions did not. The average r.m.s. deviations between computed and X-ray structures over the last 50 picoseconds of each simulation are reported in Table I. For dielectric models of 1, 2,4, 80, 1*r, 2*r, and4*r, the average r.m.s. values were 1.83, 2.84, 4.68, 5.62, 1.75, 3.81, and 4.36, respectively, using the united-atom hydrogen model and 1.44, 1.28, 3.92,
1025
Molecular Dynamics on Crambin a)
450
AII-Hydrogen/De=1
.,...-------------"T 2.0 ~
.
1.8 ~
...... 0 440 E
::::as u
...... ~
...
430
1.2 c 0 ;: 1.0 •
420
'i
0.8 Ci
Downloaded by [Rutgers University] at 09:22 12 April 2015
ILl
410 ..__ 0
__,__,....-.,........,..---r----+ 0.8 20
•u
450
.. .
·.... ~
...E
. I
1.6
\...
)(
·~
0
E
-... -·=• 0
1.4
..
......
.
.;
100
80
> CD
~
>Q
80
40
AII-Hydrogen/De=2
b) 460
...... 0 E ::::
E
...2
1.4
>Q
CD 1:
>(
1.6
c .2
I
1.2
440
CD 1:
1.0 Ci
ILl
430
•
""-----........---.---...-,.-......--.--+ 0.8 E 100 80 40 60 0 20 ..:
c)
480
AII·Hydrogen/De=4
.,...-------------------------~rs ~ >-
1!
...... 0 E :::: 450 u
4 >(
•
...... ~
.
>Q
440
.. •\
CD
•"
c
ILl
430
3
.:\,I...\ .....
!·..
; ·: ·.. ! \ :: ' ......... : ·.
~
.. "l
2
..-~ g
·;:
·;CD
•
Cl
-+-......-..--......---.--...-----..----+o IIi 1 oo 20 40 60 80 0 ..:E Time (ps)
Figure 2: Dynamics simulations at 300K for protein crambin using the all-hydrogen atom model. The energy versus simulation time is plotted with a dashed line, while the r.m.s. deviation of the computed structure from the X-ray structure is plotted with a solid line. Simulations were repeated with different continuum dielectric models: a) I, b) 2, c) 4, d) 80, e) I *r, f) 2*r, and g) 4*r.
Ornstein
1026 d) 440
AII-HydrogeniDe:BO
"T""------------~6
.c -
..• >
-
.....
430 5
0
E
:::: u
•
420
E
4
~ 410
. •
>Q
c
Downloaded by [Rutgers University] at 09:22 12 April 2015
w
3
-
400
2
390
2 c
-• i 0
Q
~o~~~~--r-~~--~~~+1
0
20
40
60
80
100
e)
AII-Hydrogen/De=1 *r
.g
1.8
f!
1.8
>
,...•. .
... .., ... ..
I
>
a
c w
80
.
390
··'
\
.
~:
~
380 +-------.--r---r-------~+ 0 0 20 40 60 80 100
Time (ps) Figure 2 continued
c
.
0 ~
·;: Q
.,;
•
~ ,_
1027
Molecular Dynamics on Crambin
g)
420
-
.
....
·,• •l
0
E :::::: 410
AII-Hydrogen/De=4*r
.......
.•••
•...>-
••
I
•
..•• ••
•u
3
>
- 400 Q
...0
I
•• •• •• ••
...•
...
Downloaded by [Rutgers University] at 09:22 12 April 2015
-
4~
•c w
•• 1 •
390
c
-• 0
>• Q
+-...-......._..--,__,P"""""''r---..--+ 0 ,;. 0 20 40 60 80 100 ~ Time (ps) Figure 2 continued
4.74, 1.71, 1.99, and 3.55, respectively, using the all-hydrogen model. Comparing simulations with the same dielectric model, the .t:\(r.m.s.) values between the unitedatom and all-hydrogen model simulations were 0.39, 1.56, 0.76, 0.88, 0.04, 1.82, and 0.81, respectively, for dielectric models 1, 2, 4, 80, 1*r, 2*r, and 4*r. The following trends are evident: a) for each dielectric model, the r.m.s. deviation between the computed and X-ray structures over the last 50 picoseconds for the allhydrogen model was always less than the r.m.s. deviation for the united-atom hydrogen model; b) increasing the dielectric value beyond 1 for the constant models and beyond 1*r for the linear models resulted in significantly and progressively poorer agreement between computed and X-ray structures, with one exception. The exception was the simulation with a constant dielectric of2 using the all-hydrogen model; this particular simulation is also noteworthy because the r.m.s. deviation was significantly lower than for any other combination of dielectric and hydrogen model. As noted above, the united-atom approximation can result in very significant computational savings in comparison with the all-hydrogen model. It has been suggested, however, that aliphatic hydrogen atoms may play a significant role in protein structure since they may have nonnegligible dipole moments of up to 0.4 debye in a C-H bond (11 ). The results described above suggest that using the united-atom approximation causes an increase in the r.m.s. deviation, ranging from 0.04 to 1.82 Adepending on the dielectric model, in comparison with the corresponding all-hydrogen simulation. Based on current data for crambin and on results for sea snake venom and pancreatic trypsin inhibitor (to be presented elsewhere), caution is advised in using the unitedatom approximation. On the other hand, our data clearly indicates that choosing an inappropriate dielectric value can result in even larger r.m.s. deviations.
1028
--
Downloaded by [Rutgers University] at 09:22 12 April 2015
-c
Ornstein
3.5
.
.,;
.. E
2.5
20
0
60
40
80
100
Time (ps) Figure 3: Comparison of r.m.s. versus time for the 14 simulations considered individually in Figures 1 and 2. Solid line curves are for the all-hydrogen model with different continuum dielectric models, while dashed line curves are for the united-atom model with different continuum dielectric models. The dielectric values used were 1 (curves without symbols), 2 (curves with filled circles), 4 (curves with filled triangles), 80 (curve with crosses), l*r (curves with open squares), 2*r (curves with open circles) and 4*r (curves with open triangles). Table I Average r.m.s. Deviation Between Computed and X-Ray Structures Over the Last 50 Picoseconds of Protein Crambin Simulations Hydrogen Atom Model
Dielectric Model Constant Value 4 2
80
Distance-Dependent l*r 2*r 4*r
United-atom All-hydrogen
1.83 1.44
2.84 1.28
4.68 3.92
5.62 4.74
1.75 1.71
3.81 1.99
4.36 3.55
A(r.m.s.)"
0.39
1.56
0.76
0.88
0.04
1.82
0.81
•Each A(r.m.s.) value is obtained by subtracting from the united-atom r.m.s. the all-hydrogen r.m.s. for a particular dielectric model.
1029
Molecular Dynamics on Crambin
Table II r.m.s. Deviation Between Lowest Energy-Minimized Structure and X-Ray Structure at Different Times During 100 Picosecond Molecular Dynamics Simulations at 300K for Crambin with Various Continuum Dielectric Models and Hydrogen-Atom Models
Downloaded by [Rutgers University] at 09:22 12 April 2015
r.m.s. Deviation (A) Duration of Simulation (picoseconds) 25 50 75
Dielectric Value
0
1 2 4 80 l*r 2*r 4*r
0.85 1.02 1.16 1.13 0.83 1.11 1.15
United-Atom Approximation 1.25 1.25 1.22 1.72 2.19 1.26 3.98 3.92 2.14 3.94 4.19 4.04 1.14 1.22 2.17 1.45 3.61 1.45 3.28 1.31 3.55
1.25 2.19 4.91 5.75 1.67 3.98 4.08
1.25 3.23 5.10 5.65 1.67 4.06 4.08
1 2 4 80 l*r 2*r 4*r
0.78 0.83 0.95 1.10 0.77 0.88 0.96
All-Hydrogen Model 1.33 1.23 1.19 1.19 2.09 1.63 2.84 1.10 1.08 1.21 1.58 1.46 1.62 1.62
1.65 1.19 2.73 3.76 1.71 1.36 3.29
1.65 1.42 2.73 4.58 1.74 2.43 3.70
10
1.65 1.19 2.73 3.76 1.52 1.36 3.04
200
We now consider the simulation data in a second way. For each simulation, the r.m.s. deviations between the computed and X-ray structures at 0, 10, 25, 50, 75, and 100 picoseconds are reported in Table II and Figure 4. The reported computed structure at each time period is the most stable energy-minimized structure to that point in the simulation resulting from complete minimization at time zero and every 5 picoseconds thereafter. For simulations using the united-atom approximation (Figure 4a), only simulations using dielectric values of 1 and 1*r resulted in r.m.s. differences ofless than 2.0 A throughout the 100 picosecond interval. The next smallest r.m.s. deviation was 3.2 for the dielectric constant of 2. For simulations using the all-hydrogen model (Figure 4b), only simulations using dielectric values of 1, 2, and 1*r resulted in r.m.s. deviations ofless than 2.0 Athroughout the 100 picosecond interval. The other constant value models (4 and 80) resulted in significantly and progressively larger r.m.s. deviations. The linear models 2*r and 4*r also resulted in large r.m.s. deviations. The poorest agreement between computed and the X-ray structures occurred with a dielectric value of 80. We believe that the general progressive worsening of the r.m.s. deviations as the dielectric value increased indicates the internal consistency, if not the physical plausibility, of the computational results.
2. Effect ofContinuum Dielectric on Dihedral Angle Difference Between Computed Structures and X-Ray Structure Compared with the united-atom model, the all-hydrogen model always leads to a smaller r.m.s. deviation between the heavy atoms ofthe backbone of the computed structures and the X-ray structure for each continuum dielectric model. The results
1030
Ornstein a) United-Atom
~
..•
Approximation
...
,)1·------------
5.5
,• ,,,'
> I
)(
,,
,,,"' ,,... ........ ---•' #>........
4.5
Jl-------··-·•=: .............-
E
2
Downloaded by [Rutgers University] at 09:22 12 April 2015
0
.;
1.5
E
.:
·-------------
-------- ••• ....*-~-:8--------··
--
!
3.5
c
2.5
~-"
l ..,as••"• ,l ~-------····;,, : l l ,' ,,,, •: l l ,l , ,,' , . ,,'' : l ,l ,,' , ..•' : £ / ,, •• 4::-----------·/ ,'' l' _,~.: ... --:::-·· --........___ _ :
t
I
= • -; a•
--.
: ,, ,' ,.........,, --·-
-"'11-------------
~.,, --~~---o •'' P"-:!i!!====••·a•···---······------··-------····---------
;.s 0.5 +-----r--.....- - r - -.....-T"""---r----+ 0
20
40
60
80
, 00
Time (ps)
b)
All-Hydrogen Model
~ >
• a:
3.75
I
)(
-.. E 0
2.75
c .!
-• ~
.;
1.75
..~ 0.75 -9-----r--.....- - r - -.....- - r - -.....-T"""---+ 0 20 40 80 80 100
Time (pe) Figure 4: r.m.s. deviation between computed lowest energy-minimized structure and X-ray structure at different times during 100 picosecond molecular dynamics simulations for the protein crambin. Solid line curves are for the all-hydrogen model with different continuum dielectric models, while dashed line curves are for the united-atom model with different continuum dielectric models. The dielectric values used were 1 (curves without symbols), 2 (curves with filled circles), 4 (curves with filled triangles), 80 (curves with crosses), 1*r (curves with open squares), 2*r (curves with open circles) and 4*r (curves with open triangles).
1031
Molecular Dynamics on Crambin a) Backbone Dihedral Angles
'iii a. + .c S;
Downloaded by [Rutgers University] at 09:22 12 April 2015
-
10
s:.
0
""
5
Time (ps)
b) Total Backbone Hydrogen Bonds 30
25 liD "C
c 0
.a c
CD
01
...0 "C
20
15
>-
-"" :::t
10
0
5
Time (ps) Figure 6: Average change in the number of hydrogen bonds during molecular dynamics simulations for crambin with the all-hydrogen model and continuum dielectric values of l (solid line without symbols), 2 (solid line and filed circles), 4 (solid line and filled triangles), 80 (solid line and crosses), l *r (dashed line without symbols), 2*r (dashed line and open circles), and 4*r (dashed line and open triangles).
1035
Molecular Dynamics on Crambin c) Specific Sidechain Hydrogen Bonds
fit
10
'1:1
c: 0
.a
c: G) Downloaded by [Rutgers University] at 09:22 12 April 2015
Q
...0 '1:1
5
:r:>-
-
0
0
-5~--~--,---~--,---~--,---~--,---~--~
0
20
60
40
80
100
Time (ps)
d) Total Sidechain Hydrogen Bonds 20~------------------------------------,-
/ ,'
fit
'1:1
....
15
c:
,..., l\ '', l \. ,-.., 'v' \,''
0
.a
c: G) Q
... '1:1 0
10
:r:>-
-
0
5
0~--~--~--~--,---~--,---~--~--~--+
0
20
40
60
Time (ps) Figure 6 continued
80
100
1036
Ornstein
Downloaded by [Rutgers University] at 09:22 12 April 2015
During the last 50 picoseconds, simulations with continuum dielectric models of 1, 2, 1*r, and 2*r on average maintained 67% to 82% of the specific hydrogen bonds noted in the energy-minimized X-ray structure. The corresponding percent for continuum models 4, 80, and 4*r ranged from 12% to l'l'Al. The total numberofhydrogen bonds maintained during the last 50 picoseconds, during simulations with the former group of continuum models, exceeded 92% of the number in the energy-minimized X-ray structure. The total number was less than 46% for the latter group of continuum model simulations, suggesting that the backbone hydrogen bonds are relatively conserved only for simulations with the continuum dielectric models 1, 2, 1*r, and 2*r. The specific sidechain hydrogen bonds were not well maintained for any continuum model simulation. The two best cases were for continuum models of2 and 2*r in which 39% and 40% of the sidechain hydrogen bonds were maintained throughout the last 50 picoseconds. On the other hand, the total number of sidechain hydrogen bonds during the last 50 picoseconds was 1 to 4 greater than the number observed in the energy-minimized X-ray structure. This suggests that sidechain hydrogen bonds are exchanging easily with one another and that the implicit solvent calculation results in more intramolecular hydrogen bonds involving sidechains than would be expected in the crystal structure (and presumably in aqueous solution). On the basis of the above hydrogen bond analysis, continuum dielectric models of 1, 2, 1*r, and 2*r seem about equally able to maintain the hydrogen bonding pattern observed in the X-ray structure.
4. Convergence of the Simulations and Implications The time zero r.m.s. comparison noted above (and in Table II) corresponds to the type of analysis performed by Whitlow and Teeter (26). In the Whitlow and Teeter study, the X-ray structure was essentially completely energy minimized and the computed structures and the X-ray structure were compared to suggest which dielectric value led to the closest agreement with the X-ray structure. The data presented above indicates, however, that the r.m.s. deviations obtained at various times during the molecular dynamics simulations using the program DISCOVER do not converge or permit semiquantitative comparison of the effect of using different dielectric values' until at least 50 picoseconds into the simulations. This finding implies that drawing conclusions based on energy minimization alone may not yield trends that have converged with those based on physically more plausible molecular dynamics simulations. Results from adequately lengthy molecular dynamics simulations are expected to be more realistic than simple energy minimization calculations because of the multiple-minima problem in protein calculations (30). Therefore, we suggest caution in using energy minimization calculations to empirically establish the most appropriate force field parameters. A comparison of the AMBER results after adequate simulation time, however, would be interesting to determine if the dielectric trends observed using DISCOVER are force field dependent.
5. Significance of the Simulation with a Dielectric Constant of 2 The above simulations suggest that a continuum dielectric model of 2 is the best
Downloaded by [Rutgers University] at 09:22 12 April 2015
Molecular Dynamics on Crambin
1037
choice for crambin simulations when using the prescription described. In order to better assess the significance of the r.m.s. difference afforded by the simulation with the all-hydrogen model and a dielectric constant of2 compared with the simulation with the all-hydrogen model and a dielectric constant of 1 (the next best continuum model), we compared the r.m.s. deviations for the backbone atoms in the secondary structure domains of crambin (Figure 7). The three secondary structure domains considered were the alpha-helix domain spanning residues 7 to 19 (Figure 7a), the alpha-helix domain between residues 23 to 30 (Figure 7b ), and the beta-sheet domain involving residues 1 to 4 and 32 to 35 (Figure 7c). Over the last 50 picoseconds of each simulation, the average r.m.s. deviation for the simulation with a dielectric constant of 2 was smaller for secondary structure domains 1, 2, and 3 by 0.26, 0.51, and 0.12 A, respectively. A dielectric value of2 also resulted in the smallest average deviation of the dihedral angles between the computed and X-ray structures. In addition, a continuum dielectric model of 2 did as good a job in preserving the backbone and sidechain hydrogen bonding pattern in the X-ray structure as any other continuum model. Figure 8 shows stereo views of the 21 energy-minimized structures at 5 picosecond intervals along the 100 picosecond trajectory using a continuum dielectric model of 2 and the all hydrogen model force-field. The structure in Figure Sa is the energyminimized X-ray structure (time 0) and indicates the location of the disulfide bridges. The alpha-carbon backbone of the structures between 5 and 100 picoseconds were least-squares fitted to the alpha-carbon positions of the energy-minimized X-ray structure (Figure 8b ). It is (visually) evident that the helical and beta secondary structure domains show considerably less conformational variability then the various turn and coil-like regions. The r.m.s. deviation for mainchain heavy atoms, over lengths of three residues, and averaged over the last 50 picoseconds of the simulation was compared to the secondary structure observed in the X-ray structure. This comparison is presented in Figure 9. The range of r.m.s. deviations for the six alpha-helical three-residue units is 0.14 to 0.36 with an average value of0.21. The corresponding average values for the two beta and six turn regions are 0.47 and 0.93, respectively. TheN-terminus and C-terminus three-residue units have r.m.s. deviations of0.44 and 1.35, respectively. It is noteworthy that theN-terminus region is involved in beta sheet formation, while the C-terminus region adopts a slightly different conformation in the solution NMR study of crambin (23). Thus it appears that the implicit continuum dielectric model2 simulation, with all-hydrogens included, can reasonable account for the structural properties of crambin. The observation that a dielectric model of2 works best may at first glance seem surprising since the DISCOVER force field was originally developed for use with a dielectric constant of 1 (7 and references therein). However given that the force field parameterization was performed on molecular fragments and amino acids and not on whole proteins and since we have chosen to use a neutral charge set for the ionic residues with an implicit water model, it is reasonable to anticipate that a somewhat different continuum model may lead to better agreement between computed and X-ray structures. In addition, the present parameterization studies were relatively lengthy molecular dynamics simulations coupled with "complete" minimizations that are expected to result in excellent confidence concerning convergence of the data.
1038
Ornstein a) Helix 1 : residues 7-19 1.4.....-----------------.-
--
..,.z:
I ll
E
..:
1.2 1.0
I
r ' ',.
I
0.8
/" '-""
,--.._J
'~
I
I I I I
0.6 0.4 0.2
Downloaded by [Rutgers University] at 09:22 12 April 2015
0
20
40
60
80
100
Time (picoseconds)
b) Helix 2: residues 23-30 1.4....---------------r t'
--
I
1.2
I
..,.z: 1.0
a;
E
..:
I
0.8
~-----,
1\ 1\ ,,/
I I
'I
/\:
\ .. J
',/
I
0.6 0.4 0.2 0
20
40
60
80
100
Time (picoseconds)
c) Sheet 1: residues 1-4 + 32-35
-
1.0
~
0.9
ui
0.8
E
0.7
..:
0.6
20
40
60
80
100
Time (picoseconds) Figure 7: r.m.s. difference in the three secondary structure domains of crambin determined by molecular dynamics simulations using all-hydrogen model and a dielectric constant of I (dashed curves) or 2 (solid line curves).
1039
Downloaded by [Rutgers University] at 09:22 12 April 2015
Molecular Dynamics on Crambin
Figure 8: Stereo views of the alpha carbon trace of crambin. (a) X-ray structure, (b) 21 energy-minimized structures at 5 picosecond intervals along the 100 picosecond trajectory (least-squares superimposed) using a continuum dielectric model of 2 and the all-hydrogen model force-field .
.44
.39
.18 .18
.14 .36
.76
.19 .18
.97 .49 1.17
.95
1.35
rms
_I_I __ I_I __ I__ IXXXI_I_L_IXXXI __ L_I_IXXXI___ I__ I 3
li
6
9
H
12
H
H
15 18 H
19
22
I
25
H
28 H
29
helix average
0.21
bela average
0.47
turn average
0.93
N-lerminus
0.44
C-lerminus
1.35
32 35 38 40 I fi I
43
46
residue
Figure 9: r.m.s. deviation in three-residue segments and comparison with secondary structure for crambin. The r.m.s. values are averaged over the last 50 picoseconds of the 100 picosecond simulation using a continuum dielectric value of 2 and the all-hydrogen model. The secondary structure is indicated as li (beta), t (tum or coil) and H (helical).
1040
Ornstein
Based on the present analysis, empirically choosing a continuum dielectric model seems a reasonable and practical alternative to an explicit solvent simulation that may be more than two orders of magnitude more computationally demanding. By doing this empirically, one can expect to reduce the protein specific error introduced by complex phenomena that cannot easily be computed, such as secondary structure cooperativity, the effects ofburied residues and the effects ofburied solvent. Similar studies that are in progress for more hydrophilic proteins indicate that, it is presently best to consider the choice of continuum dielectric model as a variable to be empirically obtained for each protein system.
Downloaded by [Rutgers University] at 09:22 12 April 2015
Conclusions The results discussed in this article suggest that 100 picosecond molecular dynamics simulations using the DISCOVER program with a continuum dielectric constant of 2 can represent the X-ray structure of the protein crambin better than simulations with several other continuum models. All-hydrogen model simulations resulted in lower r.m.s. deviations then the united-atom simulations for all continuum models tested. It should be noted, however, that cram bin is a rather small and hydrophobic protein; therefore, our results may have limited applicability for other proteins. With this in mind, we are currently performing similar investigations for hydrophilic proteins.
Acknowledgments The author thanks Dr. Mark Paulsen for critically reading the manuscript and Gail Sanders for excellent technical assistance. Thanks are due to Gail Sanders and Dr. Greg Schenter for helping to write the least-squares superimposing program. References and Footnotes
1. J.A McCammon, B.R. Gel in and M. Karp! us, Nature 267, 585 (1977). 2. M. Karp! us and J.A McCammon, Nature 277, 578 (1979). 3. J.A McCammon and S.C. Harvey, in Dynamics ofProteins and Nucleic Acids, Cambridge University Press (1987). 4. M.F. Perutz, Science 201, 1187 (1978). 5. MK Gilson and B.H. Honig, Proteins 3, 32 ( 1988). 6. M. Levitt and S. Lifson,J Mol. Bioi. 46, 269 (1969). 7. AT. Hagler, in The Peptides 7, 213 (1985). 8. J. Hermans, H.J.C. Berendsen, W.F. van Gunsteren and J.P.M. Postma, Biopolymers 23, 1513 (1984). 9. B.R. Brooks et al., J Comp. Chern. 4, 187 (1983). 10. S.P. Weiner eta/., JAm. Chern. Soc. 106,765 (1984). 11. P. Ahlstrom, 0. Teleman, B. Jonsson and S. Forsen,J Am. Chern. Soc. 109, 1541 (1987). 12. S. Lifson, NATO Adv. Study Inst./FEBS Adv. Course No. 78, (1981 ). 13. J.B. Matthew, Ann. Rev. Biophys. Biophys. Chern. 14, 387 (1985). 14. D. van Belle, I. Couplet, M. Prevost and S.J. Wodak,J Mol. Bioi. 198,721 (1987). 15. A Warshel, Nature 330, 15 (1987). 16. M.J.E. Sternberg et al., Nature 330, 86 (1987). 17. MK Gilson and B.H. Honig, Biopolymers 25, 2097 (1986). 18. MK Gilson and B.H. Honig, Nature 330, 84 ( 1987). 19. J. Aqvist, W.F. van Gunsteren, M. Leijonmarck and 0. Tapia,! Mol. Bioi. 183,461 (1985).
Downloaded by [Rutgers University] at 09:22 12 April 2015
Molecular Dynamics on Crambin
1041
20. 21. 22. 23. 24. 25. 26. 27. 28.
J. Aqvist eta/., J Mol. Bioi. 192, 593 (1987). P. Kruger eta/., Eur. Biophys. J 14,449 (1987). J. Aqvist. M. Leijonmarck and 0. Tapia, Eur. Biophys. J 16, 327 ( 1989). R.M.J.N. Lamerichs eta/., Eur. J Biochem. 171, 307 ( 1988). W.A Hendrickson and M.M. Teeter, Nature 290, 107 (1981 ). M.M. Teeter, Proc. Nat/. Acad. Sci. USA 81,6014 (1984). M. Whitlow and M.M. Teeter.J Am. Chern. Soc. 108,7163 (1986). AT. Brunger et al.. Science 235, 1049 (1987). AT. Brunger, G.M. Clore, AM. Gronenbom and M. Karplus, Proc. Nat/. Acad. Sci. USA 83,3801
29. 30. 31. 32. 33. 34. 35.
(1986). G.M. Clore, AT. Brunger, M. Karplus and AM. Gronenbom.J Mol. Bioi. 191,523 (1986). G. Nemethy and H.A Scheraga, Q. Rev. Biophys. 10,239 (1977). DISCOVER is a product of Biosym Technologies, Inc., San Diego, CA. T.F. Bernstein eta/., J Mol. Bioi. 112,535 (1977). P. Dauber-Osguthorpe et al., Proteins 4, 31(1988). J.R. Maple, U. Dinur and AT. Hagler, Proc. Nat/. Acad. Sci. USA 85, 5350 (1988). H.J.C. Berendsen et al.. J Chern. Phys. 81, 3684 (1984).
Date Received: July 12, 1989
Communicated by the Editor Wayne Hendrickson
Journal of Biomolecular Structure and Dynamics Publication details, including instructions for authors and subscription information: http://www.tandfonline.com/loi/tbsd20
Using Molecular Dynamics Simulations on Crambin to Evaluate the Suitability of Different Continuum Dielectric and Hydrogen Atom Models for Protein Simulations Rick L. Ornstein
a
a
Pacific Northwest Laboratory Molecular Science Research Center , Richland , Washington , 99352 Published online: 21 May 2012.
To cite this article: Rick L. Ornstein (1990) Using Molecular Dynamics Simulations on Crambin to Evaluate the Suitability of Different Continuum Dielectric and Hydrogen Atom Models for Protein Simulations, Journal of Biomolecular Structure and Dynamics, 7:5, 1019-1041, DOI: 10.1080/07391102.1990.10508543 To link to this article: http://dx.doi.org/10.1080/07391102.1990.10508543
PLEASE SCROLL DOWN FOR ARTICLE Taylor & Francis makes every effort to ensure the accuracy of all the information (the “Content”) contained in the publications on our platform. However, Taylor & Francis, our agents, and our licensors make no representations or warranties whatsoever as to the accuracy, completeness, or suitability for any purpose of the Content. Any opinions and views expressed in this publication are the opinions and views of the authors, and are not the views of or endorsed by Taylor & Francis. The accuracy of the Content should not be relied upon and should be independently verified with primary sources of information. Taylor and Francis shall not be liable for any losses, actions, claims, proceedings, demands, costs, expenses, damages, and other liabilities whatsoever or howsoever caused arising directly or
indirectly in connection with, in relation to or arising out of the use of the Content.
Downloaded by [Rutgers University] at 09:22 12 April 2015
This article may be used for research, teaching, and private study purposes. Any substantial or systematic reproduction, redistribution, reselling, loan, sub-licensing, systematic supply, or distribution in any form to anyone is expressly forbidden. Terms & Conditions of access and use can be found at http://www.tandfonline.com/page/terms-and-conditions
Journal ofBiomolecular Structure & Dynamics, ISSN 0739-1102 Volume 7. Issue Number 5 (1990), ®Adenine Press (1990).
Using Molecular Dynamics Simulations on Crambin to Evaluate the Suitability of Different Continuum Dielectric and Hydrogen Atom Models for Protein Simulations Downloaded by [Rutgers University] at 09:22 12 April 2015
Rick L. Ornstein Pacific Northwest Laboratory* Molecular Science Research Center Richland, Washington 99352 Abstract Molecular dynamics simulations of enzymes with enough explicit waters of solvation to realistically account for solute-solvent interactions can burden the computational resources required to perform the simulation by more than two orders of magnitude. Since enzyme simulations even with an implicit solvation model can be imposing for a supercomputer, it is important to assess the suitability of different continuum dielectric models for protein simulations. A series of 1DO-picosecond molecular dynamics simulations were performed on the X-ray crystal structure of the protein crambin to examine how well computed structures, obtained using seven continuum dielectric and two hydrogen atom models, agreed with the X-ray structure. The best level of agreement between computed and experimental structures was obtained using a constant dielectric of2 and the all-hydrogen model. Continuum dielectric models of 1, 1*r, and 2*r also led to computed structures in reasonably good agreement with the X-ray structure. In all cases, the all-hydrogen model gave better agreementthan the united atom model, although, in one case, the difference was not significant. Dielectric models of 4, 80, and 4*r with either hydrogen model yielded significantly poorer fits. It is especially noteworthy that the observed trends did not semiquantitatively converge until about 50 picoseconds into the simulations, suggesting that validation studies for protein calculations based on energy minimizations or short simulations should be viewed with caution.
Introduction Starting with the pioneering molecular dynamics simulations on pancreatic trypsin inhibitor(l,2), the dramatic progress in understanding the dynamics of proteins has been significantly fueled by computer simulations based on classical mechanics and dynamics (3). Even with today's supercomputers, protein dynamics simulations of any significant duration can totally consume available resources for all but small proteins. Electrostatic interactions are of crucial importance in protein folding and stability as evidenced by their pH dependence (4); the calculation of this component *Pacific Northwest Laboratory is operated for the U.S. Department of Energy by Battelle Memorial Institute under Contract DE-AC06-76RLO 1830.
1019
1020
Ornstein
Downloaded by [Rutgers University] at 09:22 12 April 2015
represents a major scientific and practical problem. Charge-charge interactions are calculated with either a microscopic description of the solvent or a continuum dielectric model. In theory, the microscopic treatment is expected to be superior; but it demands much more computer time (and memory), and current microscopic models do not always reproduce the dielectric properties of even liquid water (5 and references therein). Continuum models are thought to incorporate the essential qualitative electrical properties of the systems being studied and are much more computationally practical. A number of force fields have been proposed and used in protein simulations (6-11 ). Some researchers prefer to use a constant dielectric value in the Coulombic term, while others prefer a distance-dependent dielectric. The theoretical or experimental justification for either choice is subject to debate (6-18). In general, however, force field development and validation are usually applied to small molecular fragments and amino acids and not to proteins. Since proteins possess unique properties such as cooperativity and "buried" residues that their constituent amino acids do not possess, we believe that extensive force field validation studies must also be performed on whole proteins and not just on their constituent moieties. The interior of a protein molecule is generally thought to be a low dielectric medium, while the exterior surface is in a higher dielectric medium (17). Our strategy has been to implicitly account for the solvent by adjusting the effective dielectric value in the Coulombic term of the potential energy and to screen down to zero the net charge on any side chains of Arg, Asp, Glu, and Lys residues, as well as N-and C-termini. This type of approach has also been taken by others (19-22). By screening ionic residues, which tend to be solvated on the protein surface, it is reasonable to anticipate that simple continuum models may overcome known defects (5) involving relative overestimation ofelectrostatic interactions at the protein surface and underestimation of electrostatic interactions at the interior. The secondary structure and hydrogen bonding of crambin in solution have recently been investigated using two-dimensional nuclear magnetic resonance (NMR) (23) and were essentially identical to those of the crystal structure (24,25). Because of its small size and well-known structure, the protein crambin plays an important role as a model system in both experimental and tl: ·etical methods development. Crambin has been used to test potential energy funl-.10ns used in energy minimization (26) and restrained molecular dynamics simulations using simulated NMR distance restraints (27-29). Whitlow and Teeter(26) performed a large series of minimizations on the X-ray crystal structure of crambin with different force field conditions, using the program AMBER (1 0). This noteworthy effort suggested that a continuum dielectric value of 4*r leads to the best agreement between computed and X-ray structures. Because of the existence of a large number of minima in the multidimensional conformational energy hypersurface of proteins (30), energy minimization calculations starting from even an X-ray structure may still locate a minimum significantly less stable energetically than a nearby structure. Such a situation could complicate the interpretation of the type of analysis performed by Whitlow and Teeter (26). Therefore, we reexamined the question of what continuum dielectric and hydrogen-atom
Molecular Dynamics on Crambin
1021
model best accounts for the X-ray structure of crambin based on a set of 100 picosecond molecular dynamics simulations. Results using the DISCOVER program (31) suggest different conclusions and indicate that the earlier findings based on energy minimizations alone are difficult to interpret because of the multiple-minima problem.
Downloaded by [Rutgers University] at 09:22 12 April 2015
Methods The heavy atom coordinates for crambin, with diffraction data to 1.5 A(24), were taken from the October 1988 release version of the Brookhaven Protein Data Bank (32). All hydrogens were added in standard fashion or according to the united-atom approximation. In the latter procedure, nonpolar hydrogen atoms are omitted and the van der Waals radius of the nonhydrogen atoms to which they are bonded is increased accordingly. For proteins, the united-atom approximation reduces the number of atoms by about 30%, reduces the calculation time per iteration by about 50%, and reduces the total number of iterations required to obtain a specific convergence criterion in energy minimization calculations (33). The starting structure for all calculations was the X-ray structure for crambin. Prior to any molecular dynamics simulations, the starting structure was "completely energy minimized" using the following two-step procedure: minimize using the steepest descents method until the maximum derivative was less than 2.0 kcal/Aand then minimize using the conjugate gradients method until the maximum derivative was less than 0.02 kealiA All calculations were performed using the DISCOVER molecular simulation software package version 2.21, running on a 1-by-8 CRAY XMP computer. The default" consistentvalence force field" was used with no cross terms and with a simple harmonic potential for the bond-stretching terms; details of the force field and parameters are available elsewhere (7,33,34, and references therein). All nonbonded interactions were computed, i.e., an infinite cutoff distance was used. The acidic and basic residues as well as the N- and C-termini were assigned neutral partial charges as available within DISCOVER Solvent was not explicitly treated. The root mean square (r.m.s.) deviations reported between computed and X-ray structures were calculated for all nonhydrogen atoms of the backbone between Cys-3 and Cys-40; residues 1, 2, and 41-46 were not included in the r.m.s. calculations to minimize the influence of the more flexible terminal residues. The r.m.s. deviations were computed after least squares fitting to the specified backbone heavy atoms. All simulations used an integration time step of 1 femtosecond. All trajectories were carried out using a constant temperature bath (35). Each 300K dynamics simulation was divided into 20 serial segments 5.0 picoseconds long plus complete minimization. Each 5.0 picosecond segment at 300K consisted of the following sequence of steps: 1) initialize dynamics for0.5 picosecond at 300K with a randomization of the atomic velocities; 2) resume dynamics for an additional4.5 picoseconds at 300K; 3) minimize using steepest descents method until the maximum derivative was less than 2.0 kcal/ A; 4) minimize using conjugate gradients method until the maximum derivative was less than 0.02 kcal/A, and 5) compute r.m.s. deviation between the minimized structure and the X-ray structure. The 100-picosecond simulations using the united-
1022
Ornstein a)
United-Atom/De=1
380
:=0 E
-•
3-5
•..> I
>< 370
2 E
u
.Mi
..•
:
>a 360
'• '•
~
20
b)
40
60
0
E ::::
•u
0
•
Q
0 II! E 100 .;
80
....
United-Atom/De=2
4(
4-
400
.....
>•
~
350
0
c
...v· \ .'.....__,.'''
~ .~
w Downloaded by [Rutgers University] at 09:22 12 April 2015
0
..
;
~:
c::
-.
. -. •> I
>
- 370 a
..• c::
w
2
·;•
•
'
360
;
:'
Q
100
e .;
.,;
350
0
20
c)
40
60
80
Unlted-Atom/De:4
=0 E ::::
•u
•.
370
s,c
360
4
350
.Mi
......
.. •
.·....:, . ..
w
E
0
3 ;0
~~
...
'•
320
.... c
.·~
340
a>- 330 c::
-
-c
6 ..... >
380
2
>• & .,;
310 0
20
40 60 Time (ps)
80
e
1 100 ..:
Figure 1: Dynamics simulations at 300Kforprotein crambin using the united-atom approximation. The energy versus simulation time is plotted with a dashed line, while the r.m.s. deviation of the computed structure from the X-ray structure is plotted with a solid line. Simulations were repeated with different continuum dielectric models: a) 1, b) 2, c) 4, d) 80, e) 1*r, f) 2*r, and g) 4*r.
Molecular Dynamics on Crambin d)
1023
United-Atom/De:SO
360 - - - - - - - - - r - - - . , . . . - - r 6 350 5
340
4
330
..•c,.. w Downloaded by [Rutgers University] at 09:22 12 April 2015
Q
320
~
..>< E -. >-
Ill
0
280
E
:::: Ill
u
2
270
.IC ....,
c .!!
->
>- 260
.. Q
~ w
0
II
•
250
Q
20
0
f)
40
60
80
United-Atom/De=2*r
;(
34o . . . - - - - - - - - - - - - - - r s •
=0
E
>-
..•. ...:·· 1,
...Ill
··..- /\
330
4
>
•
~·
Q
310~~~~--r-~~~-,----+
20
40 Time
...0 c
~
\:• : ~ : \ .... .: ' ....,. 't
0
4
E
'•
>aa
-
United-Atom/De=4*r
60
(ps)
80
1
100
•e .:
Figure 1 continued
atom approximation and the all-hydrogen atom model required about 4.0 and 7.5 CPU hours, respectively.
Results and Discussion I. Effect ofContinuum Dielectric and Hydrogen-Atom Model on r.m.s. Deviation Between Computed Structures and X-Ray Structure
Starting from the X-ray structure of crambin, 100-picosecond molecular dynamics simulations were computed at 300K for seven dielectric models: constant values of 1, 2, 4 and 80 and linear models of 1*r, 2*r, and 4*r. The results obtained using the united-atom approximation are plotted in Figure 1, and the results obtained using the all-hydrogen atom model are shown in Figure 2. The data are shown in the form of double-Y plots with energy and r.m.s. deviations along theY-axis and simulation time along the X-axis. The data points correspond to fully minimized structures (see above for details) at time zero and at 5-picosecond intervals during the simulation. In general, the second half of the simulations (from 50 to 100 picoseconds) tended to have equilibrated. The plots ofheavy atom mainchain r.m.s. deviations (between Xray and computed structures) versus time for the 14 simulations in Figures 1 and 2 are displayed in Figure 3.1t is evident from these data that some potential functions resulted in good agreement with the X-ray structure, while some potential functions did not. The average r.m.s. deviations between computed and X-ray structures over the last 50 picoseconds of each simulation are reported in Table I. For dielectric models of 1, 2,4, 80, 1*r, 2*r, and4*r, the average r.m.s. values were 1.83, 2.84, 4.68, 5.62, 1.75, 3.81, and 4.36, respectively, using the united-atom hydrogen model and 1.44, 1.28, 3.92,
1025
Molecular Dynamics on Crambin a)
450
AII-Hydrogen/De=1
.,...-------------"T 2.0 ~
.
1.8 ~
...... 0 440 E
::::as u
...... ~
...
430
1.2 c 0 ;: 1.0 •
420
'i
0.8 Ci
Downloaded by [Rutgers University] at 09:22 12 April 2015
ILl
410 ..__ 0
__,__,....-.,........,..---r----+ 0.8 20
•u
450
.. .
·.... ~
...E
. I
1.6
\...
)(
·~
0
E
-... -·=• 0
1.4
..
......
.
.;
100
80
> CD
~
>Q
80
40
AII-Hydrogen/De=2
b) 460
...... 0 E ::::
E
...2
1.4
>Q
CD 1:
>(
1.6
c .2
I
1.2
440
CD 1:
1.0 Ci
ILl
430
•
""-----........---.---...-,.-......--.--+ 0.8 E 100 80 40 60 0 20 ..:
c)
480
AII·Hydrogen/De=4
.,...-------------------------~rs ~ >-
1!
...... 0 E :::: 450 u
4 >(
•
...... ~
.
>Q
440
.. •\
CD
•"
c
ILl
430
3
.:\,I...\ .....
!·..
; ·: ·.. ! \ :: ' ......... : ·.
~
.. "l
2
..-~ g
·;:
·;CD
•
Cl
-+-......-..--......---.--...-----..----+o IIi 1 oo 20 40 60 80 0 ..:E Time (ps)
Figure 2: Dynamics simulations at 300K for protein crambin using the all-hydrogen atom model. The energy versus simulation time is plotted with a dashed line, while the r.m.s. deviation of the computed structure from the X-ray structure is plotted with a solid line. Simulations were repeated with different continuum dielectric models: a) I, b) 2, c) 4, d) 80, e) I *r, f) 2*r, and g) 4*r.
Ornstein
1026 d) 440
AII-HydrogeniDe:BO
"T""------------~6
.c -
..• >
-
.....
430 5
0
E
:::: u
•
420
E
4
~ 410
. •
>Q
c
Downloaded by [Rutgers University] at 09:22 12 April 2015
w
3
-
400
2
390
2 c
-• i 0
Q
~o~~~~--r-~~--~~~+1
0
20
40
60
80
100
e)
AII-Hydrogen/De=1 *r
.g
1.8
f!
1.8
>
,...•. .
... .., ... ..
I
>
a
c w
80
.
390
··'
\
.
~:
~
380 +-------.--r---r-------~+ 0 0 20 40 60 80 100
Time (ps) Figure 2 continued
c
.
0 ~
·;: Q
.,;
•
~ ,_
1027
Molecular Dynamics on Crambin
g)
420
-
.
....
·,• •l
0
E :::::: 410
AII-Hydrogen/De=4*r
.......
.•••
•...>-
••
I
•
..•• ••
•u
3
>
- 400 Q
...0
I
•• •• •• ••
...•
...
Downloaded by [Rutgers University] at 09:22 12 April 2015
-
4~
•c w
•• 1 •
390
c
-• 0
>• Q
+-...-......._..--,__,P"""""''r---..--+ 0 ,;. 0 20 40 60 80 100 ~ Time (ps) Figure 2 continued
4.74, 1.71, 1.99, and 3.55, respectively, using the all-hydrogen model. Comparing simulations with the same dielectric model, the .t:\(r.m.s.) values between the unitedatom and all-hydrogen model simulations were 0.39, 1.56, 0.76, 0.88, 0.04, 1.82, and 0.81, respectively, for dielectric models 1, 2, 4, 80, 1*r, 2*r, and 4*r. The following trends are evident: a) for each dielectric model, the r.m.s. deviation between the computed and X-ray structures over the last 50 picoseconds for the allhydrogen model was always less than the r.m.s. deviation for the united-atom hydrogen model; b) increasing the dielectric value beyond 1 for the constant models and beyond 1*r for the linear models resulted in significantly and progressively poorer agreement between computed and X-ray structures, with one exception. The exception was the simulation with a constant dielectric of2 using the all-hydrogen model; this particular simulation is also noteworthy because the r.m.s. deviation was significantly lower than for any other combination of dielectric and hydrogen model. As noted above, the united-atom approximation can result in very significant computational savings in comparison with the all-hydrogen model. It has been suggested, however, that aliphatic hydrogen atoms may play a significant role in protein structure since they may have nonnegligible dipole moments of up to 0.4 debye in a C-H bond (11 ). The results described above suggest that using the united-atom approximation causes an increase in the r.m.s. deviation, ranging from 0.04 to 1.82 Adepending on the dielectric model, in comparison with the corresponding all-hydrogen simulation. Based on current data for crambin and on results for sea snake venom and pancreatic trypsin inhibitor (to be presented elsewhere), caution is advised in using the unitedatom approximation. On the other hand, our data clearly indicates that choosing an inappropriate dielectric value can result in even larger r.m.s. deviations.
1028
--
Downloaded by [Rutgers University] at 09:22 12 April 2015
-c
Ornstein
3.5
.
.,;
.. E
2.5
20
0
60
40
80
100
Time (ps) Figure 3: Comparison of r.m.s. versus time for the 14 simulations considered individually in Figures 1 and 2. Solid line curves are for the all-hydrogen model with different continuum dielectric models, while dashed line curves are for the united-atom model with different continuum dielectric models. The dielectric values used were 1 (curves without symbols), 2 (curves with filled circles), 4 (curves with filled triangles), 80 (curve with crosses), l*r (curves with open squares), 2*r (curves with open circles) and 4*r (curves with open triangles). Table I Average r.m.s. Deviation Between Computed and X-Ray Structures Over the Last 50 Picoseconds of Protein Crambin Simulations Hydrogen Atom Model
Dielectric Model Constant Value 4 2
80
Distance-Dependent l*r 2*r 4*r
United-atom All-hydrogen
1.83 1.44
2.84 1.28
4.68 3.92
5.62 4.74
1.75 1.71
3.81 1.99
4.36 3.55
A(r.m.s.)"
0.39
1.56
0.76
0.88
0.04
1.82
0.81
•Each A(r.m.s.) value is obtained by subtracting from the united-atom r.m.s. the all-hydrogen r.m.s. for a particular dielectric model.
1029
Molecular Dynamics on Crambin
Table II r.m.s. Deviation Between Lowest Energy-Minimized Structure and X-Ray Structure at Different Times During 100 Picosecond Molecular Dynamics Simulations at 300K for Crambin with Various Continuum Dielectric Models and Hydrogen-Atom Models
Downloaded by [Rutgers University] at 09:22 12 April 2015
r.m.s. Deviation (A) Duration of Simulation (picoseconds) 25 50 75
Dielectric Value
0
1 2 4 80 l*r 2*r 4*r
0.85 1.02 1.16 1.13 0.83 1.11 1.15
United-Atom Approximation 1.25 1.25 1.22 1.72 2.19 1.26 3.98 3.92 2.14 3.94 4.19 4.04 1.14 1.22 2.17 1.45 3.61 1.45 3.28 1.31 3.55
1.25 2.19 4.91 5.75 1.67 3.98 4.08
1.25 3.23 5.10 5.65 1.67 4.06 4.08
1 2 4 80 l*r 2*r 4*r
0.78 0.83 0.95 1.10 0.77 0.88 0.96
All-Hydrogen Model 1.33 1.23 1.19 1.19 2.09 1.63 2.84 1.10 1.08 1.21 1.58 1.46 1.62 1.62
1.65 1.19 2.73 3.76 1.71 1.36 3.29
1.65 1.42 2.73 4.58 1.74 2.43 3.70
10
1.65 1.19 2.73 3.76 1.52 1.36 3.04
200
We now consider the simulation data in a second way. For each simulation, the r.m.s. deviations between the computed and X-ray structures at 0, 10, 25, 50, 75, and 100 picoseconds are reported in Table II and Figure 4. The reported computed structure at each time period is the most stable energy-minimized structure to that point in the simulation resulting from complete minimization at time zero and every 5 picoseconds thereafter. For simulations using the united-atom approximation (Figure 4a), only simulations using dielectric values of 1 and 1*r resulted in r.m.s. differences ofless than 2.0 A throughout the 100 picosecond interval. The next smallest r.m.s. deviation was 3.2 for the dielectric constant of 2. For simulations using the all-hydrogen model (Figure 4b), only simulations using dielectric values of 1, 2, and 1*r resulted in r.m.s. deviations ofless than 2.0 Athroughout the 100 picosecond interval. The other constant value models (4 and 80) resulted in significantly and progressively larger r.m.s. deviations. The linear models 2*r and 4*r also resulted in large r.m.s. deviations. The poorest agreement between computed and the X-ray structures occurred with a dielectric value of 80. We believe that the general progressive worsening of the r.m.s. deviations as the dielectric value increased indicates the internal consistency, if not the physical plausibility, of the computational results.
2. Effect ofContinuum Dielectric on Dihedral Angle Difference Between Computed Structures and X-Ray Structure Compared with the united-atom model, the all-hydrogen model always leads to a smaller r.m.s. deviation between the heavy atoms ofthe backbone of the computed structures and the X-ray structure for each continuum dielectric model. The results
1030
Ornstein a) United-Atom
~
..•
Approximation
...
,)1·------------
5.5
,• ,,,'
> I
)(
,,
,,,"' ,,... ........ ---•' #>........
4.5
Jl-------··-·•=: .............-
E
2
Downloaded by [Rutgers University] at 09:22 12 April 2015
0
.;
1.5
E
.:
·-------------
-------- ••• ....*-~-:8--------··
--
!
3.5
c
2.5
~-"
l ..,as••"• ,l ~-------····;,, : l l ,' ,,,, •: l l ,l , ,,' , . ,,'' : l ,l ,,' , ..•' : £ / ,, •• 4::-----------·/ ,'' l' _,~.: ... --:::-·· --........___ _ :
t
I
= • -; a•
--.
: ,, ,' ,.........,, --·-
-"'11-------------
~.,, --~~---o •'' P"-:!i!!====••·a•···---······------··-------····---------
;.s 0.5 +-----r--.....- - r - -.....-T"""---r----+ 0
20
40
60
80
, 00
Time (ps)
b)
All-Hydrogen Model
~ >
• a:
3.75
I
)(
-.. E 0
2.75
c .!
-• ~
.;
1.75
..~ 0.75 -9-----r--.....- - r - -.....- - r - -.....-T"""---+ 0 20 40 80 80 100
Time (pe) Figure 4: r.m.s. deviation between computed lowest energy-minimized structure and X-ray structure at different times during 100 picosecond molecular dynamics simulations for the protein crambin. Solid line curves are for the all-hydrogen model with different continuum dielectric models, while dashed line curves are for the united-atom model with different continuum dielectric models. The dielectric values used were 1 (curves without symbols), 2 (curves with filled circles), 4 (curves with filled triangles), 80 (curves with crosses), 1*r (curves with open squares), 2*r (curves with open circles) and 4*r (curves with open triangles).
1031
Molecular Dynamics on Crambin a) Backbone Dihedral Angles
'iii a. + .c S;
Downloaded by [Rutgers University] at 09:22 12 April 2015
-
10
s:.
0
""
5
Time (ps)
b) Total Backbone Hydrogen Bonds 30
25 liD "C
c 0
.a c
CD
01
...0 "C
20
15
>-
-"" :::t
10
0
5
Time (ps) Figure 6: Average change in the number of hydrogen bonds during molecular dynamics simulations for crambin with the all-hydrogen model and continuum dielectric values of l (solid line without symbols), 2 (solid line and filed circles), 4 (solid line and filled triangles), 80 (solid line and crosses), l *r (dashed line without symbols), 2*r (dashed line and open circles), and 4*r (dashed line and open triangles).
1035
Molecular Dynamics on Crambin c) Specific Sidechain Hydrogen Bonds
fit
10
'1:1
c: 0
.a
c: G) Downloaded by [Rutgers University] at 09:22 12 April 2015
Q
...0 '1:1
5
:r:>-
-
0
0
-5~--~--,---~--,---~--,---~--,---~--~
0
20
60
40
80
100
Time (ps)
d) Total Sidechain Hydrogen Bonds 20~------------------------------------,-
/ ,'
fit
'1:1
....
15
c:
,..., l\ '', l \. ,-.., 'v' \,''
0
.a
c: G) Q
... '1:1 0
10
:r:>-
-
0
5
0~--~--~--~--,---~--,---~--~--~--+
0
20
40
60
Time (ps) Figure 6 continued
80
100
1036
Ornstein
Downloaded by [Rutgers University] at 09:22 12 April 2015
During the last 50 picoseconds, simulations with continuum dielectric models of 1, 2, 1*r, and 2*r on average maintained 67% to 82% of the specific hydrogen bonds noted in the energy-minimized X-ray structure. The corresponding percent for continuum models 4, 80, and 4*r ranged from 12% to l'l'Al. The total numberofhydrogen bonds maintained during the last 50 picoseconds, during simulations with the former group of continuum models, exceeded 92% of the number in the energy-minimized X-ray structure. The total number was less than 46% for the latter group of continuum model simulations, suggesting that the backbone hydrogen bonds are relatively conserved only for simulations with the continuum dielectric models 1, 2, 1*r, and 2*r. The specific sidechain hydrogen bonds were not well maintained for any continuum model simulation. The two best cases were for continuum models of2 and 2*r in which 39% and 40% of the sidechain hydrogen bonds were maintained throughout the last 50 picoseconds. On the other hand, the total number of sidechain hydrogen bonds during the last 50 picoseconds was 1 to 4 greater than the number observed in the energy-minimized X-ray structure. This suggests that sidechain hydrogen bonds are exchanging easily with one another and that the implicit solvent calculation results in more intramolecular hydrogen bonds involving sidechains than would be expected in the crystal structure (and presumably in aqueous solution). On the basis of the above hydrogen bond analysis, continuum dielectric models of 1, 2, 1*r, and 2*r seem about equally able to maintain the hydrogen bonding pattern observed in the X-ray structure.
4. Convergence of the Simulations and Implications The time zero r.m.s. comparison noted above (and in Table II) corresponds to the type of analysis performed by Whitlow and Teeter (26). In the Whitlow and Teeter study, the X-ray structure was essentially completely energy minimized and the computed structures and the X-ray structure were compared to suggest which dielectric value led to the closest agreement with the X-ray structure. The data presented above indicates, however, that the r.m.s. deviations obtained at various times during the molecular dynamics simulations using the program DISCOVER do not converge or permit semiquantitative comparison of the effect of using different dielectric values' until at least 50 picoseconds into the simulations. This finding implies that drawing conclusions based on energy minimization alone may not yield trends that have converged with those based on physically more plausible molecular dynamics simulations. Results from adequately lengthy molecular dynamics simulations are expected to be more realistic than simple energy minimization calculations because of the multiple-minima problem in protein calculations (30). Therefore, we suggest caution in using energy minimization calculations to empirically establish the most appropriate force field parameters. A comparison of the AMBER results after adequate simulation time, however, would be interesting to determine if the dielectric trends observed using DISCOVER are force field dependent.
5. Significance of the Simulation with a Dielectric Constant of 2 The above simulations suggest that a continuum dielectric model of 2 is the best
Downloaded by [Rutgers University] at 09:22 12 April 2015
Molecular Dynamics on Crambin
1037
choice for crambin simulations when using the prescription described. In order to better assess the significance of the r.m.s. difference afforded by the simulation with the all-hydrogen model and a dielectric constant of2 compared with the simulation with the all-hydrogen model and a dielectric constant of 1 (the next best continuum model), we compared the r.m.s. deviations for the backbone atoms in the secondary structure domains of crambin (Figure 7). The three secondary structure domains considered were the alpha-helix domain spanning residues 7 to 19 (Figure 7a), the alpha-helix domain between residues 23 to 30 (Figure 7b ), and the beta-sheet domain involving residues 1 to 4 and 32 to 35 (Figure 7c). Over the last 50 picoseconds of each simulation, the average r.m.s. deviation for the simulation with a dielectric constant of 2 was smaller for secondary structure domains 1, 2, and 3 by 0.26, 0.51, and 0.12 A, respectively. A dielectric value of2 also resulted in the smallest average deviation of the dihedral angles between the computed and X-ray structures. In addition, a continuum dielectric model of 2 did as good a job in preserving the backbone and sidechain hydrogen bonding pattern in the X-ray structure as any other continuum model. Figure 8 shows stereo views of the 21 energy-minimized structures at 5 picosecond intervals along the 100 picosecond trajectory using a continuum dielectric model of 2 and the all hydrogen model force-field. The structure in Figure Sa is the energyminimized X-ray structure (time 0) and indicates the location of the disulfide bridges. The alpha-carbon backbone of the structures between 5 and 100 picoseconds were least-squares fitted to the alpha-carbon positions of the energy-minimized X-ray structure (Figure 8b ). It is (visually) evident that the helical and beta secondary structure domains show considerably less conformational variability then the various turn and coil-like regions. The r.m.s. deviation for mainchain heavy atoms, over lengths of three residues, and averaged over the last 50 picoseconds of the simulation was compared to the secondary structure observed in the X-ray structure. This comparison is presented in Figure 9. The range of r.m.s. deviations for the six alpha-helical three-residue units is 0.14 to 0.36 with an average value of0.21. The corresponding average values for the two beta and six turn regions are 0.47 and 0.93, respectively. TheN-terminus and C-terminus three-residue units have r.m.s. deviations of0.44 and 1.35, respectively. It is noteworthy that theN-terminus region is involved in beta sheet formation, while the C-terminus region adopts a slightly different conformation in the solution NMR study of crambin (23). Thus it appears that the implicit continuum dielectric model2 simulation, with all-hydrogens included, can reasonable account for the structural properties of crambin. The observation that a dielectric model of2 works best may at first glance seem surprising since the DISCOVER force field was originally developed for use with a dielectric constant of 1 (7 and references therein). However given that the force field parameterization was performed on molecular fragments and amino acids and not on whole proteins and since we have chosen to use a neutral charge set for the ionic residues with an implicit water model, it is reasonable to anticipate that a somewhat different continuum model may lead to better agreement between computed and X-ray structures. In addition, the present parameterization studies were relatively lengthy molecular dynamics simulations coupled with "complete" minimizations that are expected to result in excellent confidence concerning convergence of the data.
1038
Ornstein a) Helix 1 : residues 7-19 1.4.....-----------------.-
--
..,.z:
I ll
E
..:
1.2 1.0
I
r ' ',.
I
0.8
/" '-""
,--.._J
'~
I
I I I I
0.6 0.4 0.2
Downloaded by [Rutgers University] at 09:22 12 April 2015
0
20
40
60
80
100
Time (picoseconds)
b) Helix 2: residues 23-30 1.4....---------------r t'
--
I
1.2
I
..,.z: 1.0
a;
E
..:
I
0.8
~-----,
1\ 1\ ,,/
I I
'I
/\:
\ .. J
',/
I
0.6 0.4 0.2 0
20
40
60
80
100
Time (picoseconds)
c) Sheet 1: residues 1-4 + 32-35
-
1.0
~
0.9
ui
0.8
E
0.7
..:
0.6
20
40
60
80
100
Time (picoseconds) Figure 7: r.m.s. difference in the three secondary structure domains of crambin determined by molecular dynamics simulations using all-hydrogen model and a dielectric constant of I (dashed curves) or 2 (solid line curves).
1039
Downloaded by [Rutgers University] at 09:22 12 April 2015
Molecular Dynamics on Crambin
Figure 8: Stereo views of the alpha carbon trace of crambin. (a) X-ray structure, (b) 21 energy-minimized structures at 5 picosecond intervals along the 100 picosecond trajectory (least-squares superimposed) using a continuum dielectric model of 2 and the all-hydrogen model force-field .
.44
.39
.18 .18
.14 .36
.76
.19 .18
.97 .49 1.17
.95
1.35
rms
_I_I __ I_I __ I__ IXXXI_I_L_IXXXI __ L_I_IXXXI___ I__ I 3
li
6
9
H
12
H
H
15 18 H
19
22
I
25
H
28 H
29
helix average
0.21
bela average
0.47
turn average
0.93
N-lerminus
0.44
C-lerminus
1.35
32 35 38 40 I fi I
43
46
residue
Figure 9: r.m.s. deviation in three-residue segments and comparison with secondary structure for crambin. The r.m.s. values are averaged over the last 50 picoseconds of the 100 picosecond simulation using a continuum dielectric value of 2 and the all-hydrogen model. The secondary structure is indicated as li (beta), t (tum or coil) and H (helical).
1040
Ornstein
Based on the present analysis, empirically choosing a continuum dielectric model seems a reasonable and practical alternative to an explicit solvent simulation that may be more than two orders of magnitude more computationally demanding. By doing this empirically, one can expect to reduce the protein specific error introduced by complex phenomena that cannot easily be computed, such as secondary structure cooperativity, the effects ofburied residues and the effects ofburied solvent. Similar studies that are in progress for more hydrophilic proteins indicate that, it is presently best to consider the choice of continuum dielectric model as a variable to be empirically obtained for each protein system.
Downloaded by [Rutgers University] at 09:22 12 April 2015
Conclusions The results discussed in this article suggest that 100 picosecond molecular dynamics simulations using the DISCOVER program with a continuum dielectric constant of 2 can represent the X-ray structure of the protein crambin better than simulations with several other continuum models. All-hydrogen model simulations resulted in lower r.m.s. deviations then the united-atom simulations for all continuum models tested. It should be noted, however, that cram bin is a rather small and hydrophobic protein; therefore, our results may have limited applicability for other proteins. With this in mind, we are currently performing similar investigations for hydrophilic proteins.
Acknowledgments The author thanks Dr. Mark Paulsen for critically reading the manuscript and Gail Sanders for excellent technical assistance. Thanks are due to Gail Sanders and Dr. Greg Schenter for helping to write the least-squares superimposing program. References and Footnotes
1. J.A McCammon, B.R. Gel in and M. Karp! us, Nature 267, 585 (1977). 2. M. Karp! us and J.A McCammon, Nature 277, 578 (1979). 3. J.A McCammon and S.C. Harvey, in Dynamics ofProteins and Nucleic Acids, Cambridge University Press (1987). 4. M.F. Perutz, Science 201, 1187 (1978). 5. MK Gilson and B.H. Honig, Proteins 3, 32 ( 1988). 6. M. Levitt and S. Lifson,J Mol. Bioi. 46, 269 (1969). 7. AT. Hagler, in The Peptides 7, 213 (1985). 8. J. Hermans, H.J.C. Berendsen, W.F. van Gunsteren and J.P.M. Postma, Biopolymers 23, 1513 (1984). 9. B.R. Brooks et al., J Comp. Chern. 4, 187 (1983). 10. S.P. Weiner eta/., JAm. Chern. Soc. 106,765 (1984). 11. P. Ahlstrom, 0. Teleman, B. Jonsson and S. Forsen,J Am. Chern. Soc. 109, 1541 (1987). 12. S. Lifson, NATO Adv. Study Inst./FEBS Adv. Course No. 78, (1981 ). 13. J.B. Matthew, Ann. Rev. Biophys. Biophys. Chern. 14, 387 (1985). 14. D. van Belle, I. Couplet, M. Prevost and S.J. Wodak,J Mol. Bioi. 198,721 (1987). 15. A Warshel, Nature 330, 15 (1987). 16. M.J.E. Sternberg et al., Nature 330, 86 (1987). 17. MK Gilson and B.H. Honig, Biopolymers 25, 2097 (1986). 18. MK Gilson and B.H. Honig, Nature 330, 84 ( 1987). 19. J. Aqvist, W.F. van Gunsteren, M. Leijonmarck and 0. Tapia,! Mol. Bioi. 183,461 (1985).
Downloaded by [Rutgers University] at 09:22 12 April 2015
Molecular Dynamics on Crambin
1041
20. 21. 22. 23. 24. 25. 26. 27. 28.
J. Aqvist eta/., J Mol. Bioi. 192, 593 (1987). P. Kruger eta/., Eur. Biophys. J 14,449 (1987). J. Aqvist. M. Leijonmarck and 0. Tapia, Eur. Biophys. J 16, 327 ( 1989). R.M.J.N. Lamerichs eta/., Eur. J Biochem. 171, 307 ( 1988). W.A Hendrickson and M.M. Teeter, Nature 290, 107 (1981 ). M.M. Teeter, Proc. Nat/. Acad. Sci. USA 81,6014 (1984). M. Whitlow and M.M. Teeter.J Am. Chern. Soc. 108,7163 (1986). AT. Brunger et al.. Science 235, 1049 (1987). AT. Brunger, G.M. Clore, AM. Gronenbom and M. Karplus, Proc. Nat/. Acad. Sci. USA 83,3801
29. 30. 31. 32. 33. 34. 35.
(1986). G.M. Clore, AT. Brunger, M. Karplus and AM. Gronenbom.J Mol. Bioi. 191,523 (1986). G. Nemethy and H.A Scheraga, Q. Rev. Biophys. 10,239 (1977). DISCOVER is a product of Biosym Technologies, Inc., San Diego, CA. T.F. Bernstein eta/., J Mol. Bioi. 112,535 (1977). P. Dauber-Osguthorpe et al., Proteins 4, 31(1988). J.R. Maple, U. Dinur and AT. Hagler, Proc. Nat/. Acad. Sci. USA 85, 5350 (1988). H.J.C. Berendsen et al.. J Chern. Phys. 81, 3684 (1984).
Date Received: July 12, 1989
Communicated by the Editor Wayne Hendrickson
