J. Am. Chern. Soc. 1990, 112, 2773-2781
`
`2773
`
`Molecular Dynamics of Proteins with the OPLS Potential
`Functions. Simulation of the Third Domain of Silver Pheasant
`Ovomucoid in Water
`
`Julian Tirado-Rives and William L. Jorgensen*
`
`Contribution from the Department of Chemistry, Purdue University,
`West Lafayette, Indiana 47907. Received May 15, 1989
`
`Abstract: A molecular dynamics simulation using the OPLS non bonded potential functions has been carried out for the third
`domain of silver pheasant ovomucoid in aqueous solution. Insights have been obtained on the quality of the force field, the
`convergence of such calculations, differences in the protein's structure in the crystal and in aqueous solution, protein hydration,
`and the dynamics of water molecules near a protein. The simulation covered 100 ps at 25 °C, which allowed complete equilibration
`prior to averaging and analysis of the results. Continuous monitoring of the potential energy, root-mean-square deviations
`from the crystal structure, and other properties indicated that convergence to a stable structure was achieved after 30-40 ps.
`The RMS deviation of the instantaneous structure from the crystal structure after 100 psis 1.43 A for the backbone atoms
`of residues 8-56 and 1.61 A for all residues. There is substantial reorganization of hydrogen bonds that do not involve secondary
`structure in comparing the crystal and solution structures, though in the simulation Ala-44 is displaced from the a-helix and
`Lys-29, Thr-30, Tyr-31, and Gly-32 are moved out of hydrogen-bonding distance in the triple-stranded anti parallel !3-sheet.
`Analyses of the protein-water hydrogen bonding were also carried out and are compared with results from previous simulations
`and NMR experiments.
`
`The rational design or modification of biomolecules, including
`the development of selective inhibitors for enzymes, requires
`detailed knowledge of the structure, dynamics, and corresponding
`energetics. Importantly, the continuous improvement of crys(cid:173)
`tallographic techniques has made possible the precise determi(cid:173)
`nation of the structures of many proteins, as reflected in the more
`than 300 entries now deposited in the Brookhaven Protein Data
`Bank. 1 However, the structures obtained in this fashion represent
`an ordered crystalline state, while biological processes normally
`occur in solution. Furthermore, the data obtained from crystal
`structures are static in nature, although some dynamic information
`can be obtained from the temperature (B) factors. 2 Recent
`advances in nuclear magnetic resonance spectroscopy (NMR),
`particularly the use of 2-D nuclear Overhauser effects, have been
`very valuable, since sets of distance constraints are obtained that
`can be transformed into three-dimensional structures of proteins
`in solution.3 Although this methodology allows the direct ob(cid:173)
`servation of proteins in their native solution state, the structures
`obtained reflect conformational averaging and are not unique
`solutions for the data sets. Nevertheless, the combination of the
`experimental structural results with molecular dynamics (MD)
`calculations is proving to be a powerful approach to the detailed
`characterization of the structure, dynamics, and energetics of
`proteins.4
`Since the pioneering MD simulation of bovine pancreatic trypsin
`inhibitor (BPTI) in vacuo,5 there have been numerous molecular
`dynamics calculations of proteins.4 However, the aqueous medium
`has rarely been represented in molecular detail. Some exceptions
`are for BPT1,6- 8 avian pancreatic polypeptide hormone (APP),9
`
`(I) Bernstein, F. C.; Koetzle, T. F.; Williams, G. J. B.; Meyer, E. F., Jr.;
`Brice, M. D.; Rodgers, J. R.; Kennard, 0.; Shimanouchi, T.; Tasumi, M. J.
`Mol. Bioi. 1977, I I 2, 535.
`(2) Petsko, G. A.; Ringe, D. Annu. Rev. Blophys. Bioeng. 1984,13,331.
`Fraunfelder, H.; Hartmann, H.; Karplus, M.; Kuntz, I. D., Jr.; Kuriyan, J.;
`Parak, F.; Petsko, G. A.; Ringe, D.; Tilton, R. A., Jr.; Connolly, M. L.; Max,
`N. Biochemistry 1987, 26, 254.
`(3) Wuthrich, K. NMR of Proteins and Nucleic Acids; John Wiley &
`Sons: New York, 1986.
`(4) For recent reviews, see: (a) Karplus, M.; McCammon, J. A. CRC Crit.
`Rev. Biochem. 1986, 9, 293. (b) McCammon, J. A.; Harvey, S.C. Dynamics
`of Proteins and Nucleic Acids; Cambridge University Press: Cambridge,
`England, 198 7.
`(5) McCammon, J. A.; Gelin, B. R.; Karplus, M. Nature 1977, 267, 585.
`(6) van Gunsteren, W. F.; Berendsen, H. J. C. J. Mol. Bioi. 1984, 176, 559.
`(7) Ghosh, 1.; McCammon, J. A. J. Phys. Chern. 1987, 91, 4878.
`(8) Levitt, M.; Sharon, R. Proc. Nat!. Acad. Sci. U.S.A. 1988, 85, 7557.
`
`and parvalbumin 10 in aqueous solution and BPTI in its full
`crystalline environment. 11 Other recent calculations have been
`more focused toward the modeling of enzyme-inhibitor complexes
`in water, including trypsin-benzamidine12 and thermolysin(cid:173)
`phosphonamidate.13
`With the exceptions of the parvalbumin 10 and one of the BPTI
`calculations,8 all of these simulations were run for very short times.
`Total times including the equilibrium periods were 15-30 ps for
`some of the BPTI and the APP simulations and 45 ps for the
`trypsin-benzamidine complex, as compared to 106 and 210 ps
`for the parvalbumin and the longest of the BPTI calculations,
`respectively. It is not clear that the shorter simulation times allow
`the systems to achieve equilibrium and to remove the biases from
`the starting conformation, typically obtained from a crystal
`structure. In this setting, the present study was undertaken to
`follow a molecular dynamics simulation for a protein in water for
`a long enough time to assess the convergence issue, to further test
`the performance of the OPLS force field, 14 and to obtain insights
`on protein hydration and possible differences between the solution
`and crystal structures. The protein selected for the study was the
`third domain of silver pheasant ovomucoid (OMSVP3). 15
`Ovomucoids make up about 10% of the protein in avian egg
`whites, in which they are the dominant inhibitors of serine pro(cid:173)
`teases. They are members of the Kazal family of pancreatic
`secretory inhibitors, which are generally important in controlling
`the premature activation of pancreatic zymogens. The complete
`ovomucoid consists of three homologous, tandem domains, each
`
`(9) Kruger, P.; Stra/3burger, W.; Wollmer, A.; van Gunsteren, W. F. Eur.
`Biophys. J. 1985, I 3, 77.
`(10) (a) Ahlstrom, P.; Teleman, 0.; Jonsson, B.; Forsen, S. J. Am. Chern.
`Soc. 1987, 109, 1541. (b) Ahlstrom, P.; Teleman, 0.; Jonsson, B. J. Am.
`Chern. Soc. 1988, 110, 4198.
`(II) (a) van Gunsteren, W. F.; Karplus, M. Biochemistry 1982, 21, 2259.
`(b) Swaminathan, S.; Ichiye, T.; van Gunsteren, W. F.; Karplus, M. Bio(cid:173)
`chemistry 1982, 21, 5230. (c) van Gunsteren, W. F.; Berendsen, H. J. C.;
`Hermans, J.; Hoi, W. G. J.; Postma, J.P. M. Proc. Nat!. Acad. Sci. U.S.A.
`1983, 80, 4315.
`(12) Wong, C. F.; McCammon, J. A.lsr. J. Chern. 1986, 27, 211; J. Am.
`Chern. Soc. 1986, 108, 3830.
`(13) Bash, P. A.; Singh, U. C.; Brown, F. K.; Langridge, R.; Kollman, P.
`A. Science 1987, 235, 574.
`(14) Jorgensen, W. L.; Tirado-Rives, J. J. Am. Chern. Soc. 1988, 1/0,
`1657.
`(15) A preliminary report on this work was provided in the Proceedings
`of the 1988 Nobel Symposium: Jorgensen, W. L.; Tirado-Rives, J. Chern. Scr.
`1989, 29A, 191.
`
`0002-7863/90/1512-2773$02.50/0
`
`© 1990 American Chemical Society
`
`

`
`2774 J. Am. Chern. Soc., Vol. 112, No. 7, 1990
`
`Tirado- Rives and Jorgensen
`
`of which may inhibit a protease. The third domain contains 56
`amino acid residues that can be detached by controlled proteolysis.
`These fragments are typically at least as active as the full ovo(cid:173)
`mucoid.
`Sequences have now been determined for ovomucoid third
`domains from over I 00 avian species by Laskowski and co(cid:173)
`workers.16 They have also determined association constants for
`many of these inhibitors with a-chymotrypsin (AC), elastase
`(HLE), subtilisin, and Streptomyces griseus proteases A and B
`(SGPA and SGPB). 17•18 Many of the sequences differ by only
`one or a few residue changes, so an unusually complete struc(cid:173)
`ture/activity data base is being constructed.
`Besides the association and hydrolytic constants, some structural
`data are also available. Crystal structures have been determined
`for complexes of turkey ovomucoid third domain with SGPB, 19
`HLE,20 and AC,21 for two isolated third domains, silver pheasant22
`and Japanese quail,23 and for their corresponding hydrolyzed
`forms. 24 However, structural data in solution are limited to two
`studies by 2-D NMR, specifically, for turkey ovomucoid third
`domain in both its native25 and hydrolyzed forms. 26
`The typical structure of an ovomucoid third domain contains
`three disulfide bridges, a I 0-11 residue long a-helix, and a tri(cid:173)
`ple-stranded antiparallel ~-sheet. The combination of small size,
`a large body of experimental data, and interesting structure makes
`the ovomucoids unusually attractive for a series of molecular
`dynamics investigations. Such simulations could provide detailed
`structural and thermodynamic information, which would be of
`great assistance in the interpretation of the biophysical data and
`the development of selective inhibitors. The third domain of silver
`pheasant ovomucoid was chosen for this initial study owing to its
`greater sequence homology with other ovomucoids and the better
`resolution of its crystal structure· than that of Japanese quail.
`
`Computational Procedure
`The entire simulation was conducted on Sun-4 computers using the
`AMBER 3.0 program,27 with minor local modifications to improve its
`use of the UNIX environment. The OPLS non bonded parameters14 were
`used for the protein atoms in conjuction with the TIP3P model for
`water.28 As specified in the OPLS model, the dielectric constant was
`kept fixed at 1.0, and the scaling factors for the I ,4-nonbonded interac(cid:173)
`tions were 8.0 for the Leonard-Jones and 2.0 for the electrostatic inter(cid:173)
`actions.14 The energetics for angle bending and torsional motion were
`described with the AMBER united-atom force field. 27 During the sim(cid:173)
`ulation, all bond lengths and the H-H distances in water were kept
`constant by using the SHAKE algorithm29 with a tolerance of 0.0004 A,
`
`(16) Laskowski, M., Jr.; Kato, I.; Ardelt, W.; Cook, J.; Denton, A.; Empie,
`M. W.; Kohr, W. J.; Park, S. J.; Parks, K.; Schatz1ey, B. L.; Schoenberger,
`0. L.; Tashiro, M.; Vichot, G.; Whatley, H. E.; Wieczorek, A.; Wieczorek,
`M. Biochemistry 1987, 26, 202. Kato, I.; Kohr, W. J.; Laskowski, M., Jr.
`Biochemistry 1987, 26, 193.
`(17) Empie, M. W.; Laskowski, M., Jr. Biochemistry 1987, 21, 2274.
`(18) Laskowski, M., Jr.; Tashiro, M.; Empie, M. W.; Park, S. J.; Kato, 1.;
`Ardelt, W.; Wieczorek, M. In Proteinase Inhibitors: Medical and Biological
`Aspects; Katunuma, N., Umezawa, H., Hozer, H., Eds.; Japan Scientific
`Societies Press: Tokyo, 1983; p 55.
`(19) Fujinaga, M.; Read, R.; Sielecki, A. R.; Ardelt, W.; Laskowski, M.,
`Jr.; James, M. N. G. Proc. Nat/. Acad. Sci. U.S.A. 1982, 79, 4868. Read,
`R.; Fujinaga, M.; Sie1ecki, A. R.; James, M. N. G. Biochemistry 1983, 22,
`4420.
`(20) Bode, W.; Wei, A.-Z.; Huber, R.; Meyer, E.; Travis, J.; Neumann,
`S. EMBO J. 1986, /0, 2453.
`(21) Read, R.; Fujinaga, M.; Sielecki, A. R.; Arde1t, W.; Laskowski, M.,
`Jr. Acta Crystallogr., Sect. A, Supp/. 1984, 40, C-50.
`(22) Bode, W.; Epp, 0.; Huber, 0.; Laskowski, M., Jr.; Ardelt, W. Eur.
`J. Biochem. 1985, 147, 387.
`(23) Weber, E.; Papamokos, E.; Bode, W.; Huber, R.; Kato, I.; Laskowski,
`M., Jr. J. Mol. Bioi. 1981, 149, 109. Papamokos, E.; Weber, E.; Bode, W.;
`Huber, R.; Empie, M. W.; Kato, I.; Laskowski, M., Jr. J. Mol. Bioi. 1982,
`158, 515.
`(24) Musil, D.; Bode, W.; Mayr, 1.; Huber, R.; Laskowski, M., Jr.; Lin,
`T.-Y.; Ardelt, W. Unpublished results.
`(25) Robertson, A.; Westler, W. M.; Markley, J. L. Biochemistry 1988,
`27, 2519.
`(26) Rhyu, G. 1.; Markley, J. L. Biochemistry 1988, 27, 2529.
`(27) Weiner, S. J.; Kollman, P. A.; Case, D. A.; Singh, U. C.; Ghio, C.;
`Alagona, G.; Profeta, S.; Weiner, P. J. J. Am. Chern. Soc. 1984, /06, 765.
`(28) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J.D.; Impey, R. W.;
`Klein, M. L. J. Chern. Phys. 1983, 70, 926.
`(29) van Gunsteren, W. F.; Berendsen, H. J. C. Mol. Phys. 1977,34, 1311.
`
`-19000
`
`I
`
`I
`I
`
`-19300 ~
`I
`>:
`t c
`
`-19600
`
`0
`
`Oll
`~
`
`5
`~
`
`-20500 i
`
`0
`
`20
`
`60
`40
`time(psJ
`
`80
`
`100
`
`Figure l. Potential energy variation during the course of the MD sim(cid:173)
`ulation.
`
`i.
`I 0 lJ
`
`0 . .\.C
`
`11.11 -r-------,----,------,---1
`211.0
`W.O
`100.0
`1111
`()()()
`~0.0
`
`Figure 2. RMS deviations between the instantaneous computed structure
`and the crystal structure for all residues as a function of simulation time.
`
`which allowed the use of a time step of 2 fs. The temperature and
`pressure were kept constant at 298 K and 1 bar (0.987 atm). A non(cid:173)
`bonded pair list was used to accelerate the calculations and was updated
`every 10 steps. This list was generated by using a residue-based cutoff
`(9 A) to avoid splitting dipoles. All the calculations utilized periodic
`boundary conditions to avoid edge effects.
`The initial coordinates for the third domain of silver pheasant ovo(cid:173)
`mucoid were obtained from the crystal structure22 deposited in the
`Brookhaven Protein Data Bank.l.3° The protein molecule, without any
`of the crystallographically located water molecules, was centered in a
`rectangular box of water obtained by periodic translations in the x, y, and
`z directions of a cube of water previously equilibrated via Monte Carlo
`calculations. Any water molecule closer than 1.5 A to any protein atom
`or farther away than 6 A from the closest protein atom in any one
`Cartesian direction was then deleted to give an initial system containing
`the solute plus 1721 water molecules (5676 total atoms) in a rectangular
`box of dimensions 43.8 X 42.0 X 34.3 A.
`The initial preparation of the system consisted of 100 steps of steepest
`descent energy minimization, followed by a short (1 ps) constant volume
`molecular dynamics run. The resulting structure was then used as the
`starting point for the MD simulation at constant temperature and pres(cid:173)
`sure. Initial atomic velocities were assigned from a Maxwellian distri(cid:173)
`bution corresponding to a temperature of 298 K. The values of the
`potential energy, RMS deviation from the crystal structure, and volume
`were monitored continuously in order to follow the equilibration of the
`system. A total time of 100 ps was covered in the simulation, during
`which the coordinates, velocities, and energies were saved every 50 time
`steps (0.1 ps) for further analysis.
`
`Results
`Convergence Behavior. The potential energy and RMS devi(cid:173)
`ations of the main-chain atoms (N, C", C, and 0) of residues 8-56
`for each instantaneous structure are plotted in Figures I and 2,
`respectively, as a function of simulation time. The first seven
`
`(30) Entry 20VO, version of Nov 8, 1985.
`
`

`
`Molecular Dynamics Simulation for Pheasant Ovomucoid
`
`J. Am. Chern. Soc., Vol. 112, No.7, /990 2775
`
`residues were excluded from the RMS deviations, since they form
`a pendant tail that is less well resolved in the crystal structure. 22
`The RMS deviations are typically 0.2 A higher when these residues
`are included. The system appears to have reached an equilibrium
`state after 30-40 ps, though the ultimate longevity of this state
`could only be determined in a much longer run. After this initial
`period, the simulation was continued to 100 ps in order to acquire
`enough data for averaging and analysis. In the 210-ps simulation
`for the other protease inhibitor, BPTI, in water, an equilibrium
`state was achieved after ca. 50 ps. 8
`Protein Structure and Dynamics. (A) Comparison with the
`Crystal Structure. The calculated mean structure of the protein
`in solution, obtained by directly averaging the Cartesian coor(cid:173)
`dinates of all the saved structures from 30 to 100 ps, shows an
`RMS deviation from the crystal structure of 1.28 A for the
`backbone atoms of residues 8-56 and 1.49 A for all residues.
`Results from previous simulations in water include 1.05 A for the
`ca of APP, averaged from 5 to 15 ps,9 and 0.77 A for all the
`backbone atoms of residues 1-56 of BPTI, averaged from I 05
`to 210 ps. 8 Deviations from simulations in vacuo are typically
`larger by at least a factor of 2. 8•11
`The value of the RMS deviation of the instantaneous structure
`at I 00 ps from the crystal structure is 1.43 A for the backbone
`atoms of residues 8-56 and 1.61 A for all residues. These values
`are in the same range as those obtained in previous, shorter sim(cid:173)
`ulations in water, e.g., 1.5 A for only theca atoms of BPTI after
`20 ps6 and 1.72 A for theca atoms of the trypsin-benzamidine
`complex at 45 psY Clearly, the AMBER/OPLS force field is
`providing a comparatively reasonable representation of the protein.
`An approximate base line for the RMS deviations can be deduced
`from the comparison of coordinates from proteins whose crystal
`structures have been solved in different crystalline forms or that
`have different molecules in the asymmetric unit. For five such
`cases, Chothia and Lesk obtained RMS deviations ranging from
`0.25 to 0.40 A, with a mean of 0.33 A.31
`A comparative plot of the backbone atoms of the crystal
`structure with the instantaneous structure at the end of the MD
`It can be seen that the tertiary structure
`run is given in Figure 3.
`is well preserved in the simulation. However, several striking
`differences with the crystal structure are revealed by more detailed
`graphical analyses:
`(I) In the crystal structure, the amino acid side chains are
`mostly folded onto the surface of the protein, while in the solution
`simulation they extend more into the solvent. This effect likely
`reflects the lower water content in the crystal and the removal
`of the interprotein interactions. It also allows polar side chains
`to be better stabilized by hydration.
`(2) The last residue of the a-helix in the crystal, Ser-44, does
`not form the required hydrogen bond with Ala-40 to be part of
`the helix in solution. This terminal hydrogen bond is also absent
`in the crystal structures of isolated Japanese quai123 and hydrolyzed
`silver pheasant24 ovomucoids and in the complexes of turkey
`ovomucoid with a-chymotrypsin21 and S. griseus protease 8. 19
`(3) As discussed in the next section, some of the residues of
`the outermost strand of the 13-sheet (Lys-29, Thr-30, Tyr-31, and
`Gly-32) are tilted out of the plane of the other two strands in the
`simulation. This displacement inhibits the formation of the in(cid:173)
`terstrand hydrogen bonds and yields greater hydrogen bonding
`between these residues and water molecules.
`(4) The hydrogen bond in the crystal structure between the side
`chains of Glu-19 and Thr-17, which surround the scissile peptide
`bond between Met- I 8 and Glu-19, is not found in the simulation
`results. Instead, it is replaced by a hydrogen bond between the
`guanidinium fragment of Arg-21 and the carboxylate in Glu-19,
`while the hydroxyl group of Thr-17 is hydrogen-bonded to solvent
`molecules. This effect appears to be related to the loss of in(cid:173)
`terprotein interactions. In the crystal lattice, the charged side
`charged side chain of Arg-21 is close to the carbonyl oxygens of
`Glu- 10 and Pro-12 from a protein molecule in a neighboring unit
`cell. Interestingly, the same interprotein contacts were found in
`
`MD, 100 ps
`
`Crystal Structure
`
`Figure 3. Comparison of backbone atoms in the instantaneous structure
`at the end of the simulation (t = 100 ps) and in the crystal structure.
`
`the crystal structure of Japanese quail ovomucoid, which has a
`very different crystalline form. 23 •32
`In addition, the Glu-19
`carboxylate group in the crystal is hydrogen-bonded to the
`side-chain ammonium group of Lys-13 in the same neighboring
`cell. Figure 4 compares the different interactions in this region
`for the crystal structure and the instantaneous structure at 80 ps.
`The overall conformation of a protein can be expressed in terms
`of the backbone torsional angles tf>; (C;_ 1-N;-Cj-CJ and !/;;
`(N;-Cj-C;-N;+ 1). The average values for these angles during the
`3Q-100-ps period are compared with the corresponding values for
`the crystal structure in Figure 5. Consistent with the statements
`above, the biggest differences are found in the region around the
`C-terminus of the a-helix (Giu-43, Ser-44, and Asn-45), where
`the middle residue is twisted away from the helix, and in the region
`near the turn connecting the central and outer strands of the
`13-sheet (Ser-26, Asp-27, and Asn-28). Some of the other dif(cid:173)
`ferences are of little consequence to the overall conformation, since
`compensation occurs when there are differences of opposing signs
`in an angle 1/1; and the subsequent tf>;+l· This is the case for Met-18
`and Glu-19 and for the Thr-47, Leu-48, and Thr-49 region. The
`
`(31) Chothia, C.; Lesk, A. M. EMBO J. 1986, 5, 823.
`
`(32) Silver pheasant ovomucoid crystallizes in the C2 space group, while
`for Japanese quail, the crystals belong in the tetragonal P42 12 space group.
`
`

`
`2776 J. Am. Chern. Soc., Vol. 112, No. 7, 1990
`
`MD, lOOps
`
`Tirado- Rives and Jorgensen
`
`I ~ll
`
`120
`
`II
`II
`
`-
`
`MDavcrag:c
`
`- - crystal
`
`0
`
`10
`
`20
`
`30
`
`40
`
`50
`
`(10
`
`fC'>idUC
`
`'
`
`Crystal Structure
`
`< N~ 3.01
`
`3·! .. :·.433 ~. Lys-!3'
`
`'" ' J .•
`of:~
`
`Pro-12'
`
`Figure 4. Hydrogen bonds near the scissile bond in the instantaneous
`structure at the end of the simulation (t = 100 ps) and in the crystal
`structure.
`
`RMS differences from the crystal structure computed for each
`backbone dihedral angle over the entire protein are 39°, 48°, and
`11° for ¢, t/;, and w, respectively. These differences are in the
`same range as the results from previous simulations, e.g., 31°,
`37°, and 8° for¢, t/;, and win trypsin 12 and 26° and 33° for¢
`and t/; for BPTI in "van der Waals water". 11"
`A more global impression of the overall conformation of the
`protein can be obtained from the plots of¢ vs t/; (Ramachandran
`maps) in Figure 6. The general trend observed in the maps is
`that residues in the a-helix stay close to the crystallographically
`observed angles, while those outside the helix are shifted toward
`values more typical of the C5 and q'l conformations. This dis(cid:173)
`placement is consistent with the results of our previous energy(cid:173)
`minimization studies on the conformations of N-acetylglycine(cid:173)
`N-methylamide and N-acetylalanine-N-methylamide, which
`showed that these conformations are lower in energy than the
`corresponding a, aR, or aL alternatives. 14
`(B) Comparison with NMR Data. The NMR data in solution
`on native and hydrolyzed turkey ovomucoid should be relevant
`to the present case, since turkey and silver pheasant ovomucoids
`differ by only one residue (18: Met/Leu). In general, the main
`effort in the NMR investigations was devoted to the complete
`assignment of resonances, and only qualitative information re(cid:173)
`garding secondary structure was obtained.25•26
`Although the principal conclusion of these studies was that the
`solution structure had to be very similar to that in the crystal,
`some details were explored further. In particular, the NOESY
`connectivities that establish the antiparallel triple-stranded ~-sheet
`were given. Figure 7 provides a comparison of the segments
`comprising the {1-sheet in the crystal and in the mean calculated
`solution structure, obtained by direct averaging of the Cartesian
`
`()
`
`10
`
`20
`
`30
`rc:-.idu~
`Figure S. Plots of the backbone angles ¢ and !J; for the crystal structure
`and the averages from the MD simulation.
`
`40
`
`50
`
`60
`
`coordinates from 30 to 100 ps in the MD simulation. The in(cid:173)
`terproton distances given with the crystal fragments are measured
`in the crystaJ33 but correspond to observed NOE pairs in the
`solution NMR. Only the distances that show significant deviations
`from the MD averages are given. As alluded to above, the largest
`differences are found in the separation between the central and
`the outermost strands of the {1-sheet. The two distances of 6-7
`A calculated in the simulation are beyond the normally accepted
`range for NOE detection. A problem with the force field or the
`preparation of the system could be indicated. However, it is also
`possible that the protein is exploring the expanded phase of a
`low-frequency "breathing mode", and the separations may decrease
`at later times in the simulation. 34
`(C) Fluctuations. The overall mobility of the different atoms
`during an MD simulation can be expressed as their RMS fluc(cid:173)
`tuations, (&2) 1/ 2 = ( (r,.- (r,.) )2) 112, computed over the averaging
`period. An examination of the accumulated average atomic
`fluctuations as the time span of their evaluation increases shows,
`as expected, a monotonic increase for about 40 ps that levels
`toward a plateau after ca. 50 ps. These fluctuations can be
`compared to the average movement observed in the crystallo(cid:173)
`graphic determination as expressed in the B thermal factors.
`Crystallographic RMS fluctuations can then be derived via the
`relation (tJ.r2 ) 112 = {3B/87r2) 112 from the B-values. 35 Figure 8
`compares the RMS fluctuations calculated during the last 25 ps
`of the simulation with the fluctuations derived from the thermal
`factors. 22
`
`(33) Any missing hydrogen atoms in the X-ray and the MD structures
`were added by using the SYBYL program from TRIPOS Associates.
`(34) Suezaki, Y.; Go, N. Int. J. Pep/. Protein Res. 1975, 7, 333.
`(35) Willis, B. T. M.; Pryor, A. W. Thermal Vibrations in Crystallogra(cid:173)
`phy; Cambridge University Press: Cambridge, England, 1975.
`
`

`
`Molecular Dynamics Simulation for Pheasant Ovomucoid
`
`J. Am. Chern. Soc., Vol. 112, No. 7, 1990 2777
`
`ISO~
`I
`
`I *
`
`I
`120-
`
`+
`
`7 "'
`
`0 -·
`
`+
`
`-60 ~
`I
`
`I
`-120 -,
`
`-lg(J
`-IRO
`
`+
`
`MD Average
`
`+
`
`-120
`
`-60
`
`60
`
`120
`
`180
`
`X
`
`><,(X :f< X X~
`X
`x~Xx
`
`loli
`
`X X
`
`X
`
`XX
`~
`X
`
`X~
`
`180
`
`120
`
`60
`
`0
`
`-60
`
`-120
`
`0
`
`c,_
`
`X
`
`~
`X
`
`Crystal
`
`-180 +-----,---,-----,----,-----,------1
`-180
`-120
`-60
`0
`60
`120
`180
`phi( 0
`Figure 6. Ramachandran maps for the MD results and the crystal
`structure.
`
`)
`
`Table I. Average RMS Fluctuations (A) for Different Atom Types
`averaging time, ps
`3G-70
`1.16
`1.24
`1.30
`1.30
`1.27
`3.96
`
`atoms
`ca
`0
`c~
`c~
`all protein
`owater
`
`3D-50
`0.73
`0.82
`0.90
`0.92
`0.86
`2.76
`
`3G-100
`1.32
`1.40
`1.50
`1.45
`1.44
`5.17
`
`It can be seen from the plot that qualitative agreement in the
`width and location of the largest peaks is obtained. The most
`noticeable difference is the comparatively high fluctuations for
`the N-terminus in the simulation, which decrease until Cys-8 is
`reached. The first seven residues form a pendant chain whose
`increased mobility in solution can be attributed to the formation
`of a four-stranded ~-channel with the N-terminal segments of
`neighboring molecules in the crystal. 22 The remaining qualitative
`differences reflect greater motion in the simulation for the side
`chains of Glu-1 0, Met-18, Asp-27, Lys-29, Glu-43, and Asn-45,
`which are more exposed to the solvent in the simulation than in
`the crystal structure.
`The RMS fluctuations for the different atoms are correlated
`to their distances from the backbone of the protein. Table I lists
`the average RMS fluctuations for different carbon atoms for three
`averaging periods. The values obtained for the C6 atoms have
`larger statistical uncertainties, since there are only 30 of these
`atoms in the protein. The fluctuations calculated in the initial
`20 ps of the averaging period are comparable to those obtained
`in the simulations of aqueous trypsin-benzamidine during 19 ps
`(C" = 0.52, 0~ = 0.59, C"Y = 0.69, C 8 = 0.70, all= 0.68), 12 BPTI
`in "van der Waals water~ for 25 ps (C" = 0.54, Cfj = 0.69, C'Y
`
`~ \ • • ~ 0
`li •
`~ 635
`39' y~
`r •
`~:;~,
`• • •
`
`I
`
`•
`
`361
`
`663
`
`MD <30-1 OOps>
`
`Crystal Structure
`
`Figure 7. Backbone atoms of the ~-sheet segments for the average MD
`structure and in the crsytal structure. Some interproton distances are
`shown.
`
`3 .0 - , - - - - - - - - - - - - - - - - - - ,
`
`-MD
`- - X-Ray
`
`-$
`tl
`·~
`a
`!!!
`r;:
`~
`
`Cl)
`
`2.5
`
`2.0
`
`1.5
`
`1.0
`
`0.5
`
`0.0
`
`0
`
`90
`
`180
`
`270
`Atom Number
`Figure 8. Comparison of RMS fluctuations from the MD calculation and
`from the crsytallographic B-factors.
`
`360
`
`450
`
`540
`
`= 0.89, C8 = 1.15, all= 0.78), 11" and aqueous APP for 10 ps (C"'
`= 0.53, Cfj = 0.62, C'Y = 0.73, C8 = 0.79, all= 0.75).9 However,
`the fluctuations calculated over the full averaging period of 70
`ps are considerably larger than those reported by Levitt and
`Sharon for the final 105 ps of residues 2-56 of BPTI (C" = 0.42,
`Cil = 0.50, C'Y = 0.54, C6 = 0.67).8 Greater experience is needed
`to ascertain if these differences are associated more with the
`proteins, the force fields, or the details of the simulation procedures,
`though further comment on this issue is made below.
`(D) Hydrogen Bonding. Analyses of the intra protein hydrogen
`bonding were carried out on the 700 coordinate sets saved during
`the last 70 ps of the MD simulation. The criteria used to define
`a hydrogen bond were purely geometric. A list of all potential
`donors and acceptors (hydrogens attached to heteroatoms and the
`heteroatoms themselves) was generated at the beginning of the
`analysis. For each coordinate set, every potential donor-acceptor
`
`

`
`2778
`
`J. Am. Chern. Soc., Vol. 112, No. 7, /990
`
`Tirado- Rives and Jorgensen
`
`Tyr-11
`Arg-21
`Leu-23
`Cys-24
`Gly-25
`Ser-26
`Asn-28
`Tyr-31
`Phe-37
`Cys-38
`Asn-39
`Ala-40
`Val-41
`Val-42
`Glu-43
`Ser-44
`Asn-45
`Gly-46
`Thr-47
`Ser-51
`His-52
`Cly-54
`
`N
`N
`N
`N
`N
`N
`N
`N
`N
`N
`N
`N
`N
`N
`N
`N
`N
`N
`N
`N
`N
`N
`
`N
`N
`N
`N
`N
`N
`
`OG
`Ser-26
`NZ
`Lys-34
`ND2 Asn-33
`ND2 Asn-33
`ND2 Asn-39
`OG
`Ser-44
`Ser-44
`OG
`OGI
`Thr-47
`NDI His-52
`
`Table II. Protein-Protein Hydrogen Bonds and Their Frequencies
`acceptor
`donor
`crystal
`residue distance, A
`atom
`atom
`residue
`Backbone to Backbone
`Cys-8
`3.4
`0
`Gly-32
`2.7
`0
`Tyr-31
`2.9
`0
`2.8
`His-52
`0
`Lys-29
`2.8
`0
`Thr-49
`2.8
`0
`Gly-25
`2.9
`0
`Leu-23
`2.8
`0
`Asn-33
`3.0
`0
`Lys-34
`2.8
`0
`Cys-35
`3.0
`0
`Asn-36
`3.1
`0
`Phe-37
`2.9
`0
`Cys-38
`3.1
`0
`Asn-39
`2.9
`0
`Ala-40
`2.9
`0
`Val-42
`2.9
`0
`Val-41
`2.8
`0
`Ser-44
`3.2
`0
`Cys-24
`2.9
`0
`Cys-24
`3.2
`0
`Pro-22
`3.1
`0
`Backbone to Side Chain
`Lys-13 ODI Asn-39
`3.1
`Lys-29 OD2 Asp-27
`3.1
`Asn-36 ODI Asn-33
`3.0
`Ser-44
`3.0
`Leu-48 OG
`Thr-49 OG
`Ser-26
`3.0
`Cys-56 OGI
`Thr-30
`3.0
`Side Chain to Backbone
`Thr-49
`3.2
`0
`Asp-7
`3.1
`0
`Thr-17
`2.8
`0
`Glu-19
`3.0
`0
`Lys-13
`3.1
`0
`Ala-40
`3.1
`0
`Val-41
`2.9
`0
`Ser-44
`2.9
`0
`Phe-53
`2.6
`0
`Side Chain to Side Chain
`ODI Asp-7
`Ser-9
`OG
`2.5
`Thr-17 OE2 Glu-19
`OGI
`3.0
`Ser-26 OGI
`Thr-49
`OG
`3.2
`OEH Tyr-31 ODI Asp-27
`2.7
`Lys-55 OEH Tyr-20
`NZ
`3.2
`• The reported percentage corresponds to the population of hydrogen
`bonds to either of the two oxygens in the carboxylate acceptor. The
`hydrogen bonds to each oxygen have occurrences of 50.4 and 32.9%.
`bThe reported percentage corresponds to the population of hydrogen
`bonds to either of the two oxygens in the carboxylate acceptor. The
`hydrogen bonds to each oxygen have occurrences of 52.7 and 32.0%.
`
`MD
`frequency, %
`
`22.6
`99.6
`
`99.4
`
`96.6
`
`94.1
`100.0
`98.9
`52.3
`70.6
`98.6
`85.3
`0.3
`
`8.0
`
`55.3
`77.7
`53.4
`
`99.0
`44.1
`99.9
`
`1.4
`
`0.1
`
`82.9
`64.7
`13.7
`
`23.7
`
`79.0"
`
`80.Jb
`
`60
`
`50
`
`40
`
`30
`
`20
`
`10
`
`0
`0.0
`
`2.0
`
`II
`•II
`llillvl II
`II II \1'1
`Iii
`II
`II
`~
`I I
`I
`I
`
`: Y
`
`-protein
`
`- - water
`
`I;
`1111
`111
`111
`I ~ ;\
`II 1,,
`~ c II\
`' II;~
`\ I tJ 1 r
`
`10.0
`
`12.0
`
`8.0
`4.0
`6.0
`RMS fluctuation (A)
`Figure 9. Distribution of RMS fluctuations for heavy atoms of the
`protein and water oxygens from the MD simulation.
`
`conserved in the MD simulation. This is also reasonable, since
`the backbone atoms show smaller RMS deviations from the crystal
`structure than do the side chains. A total of 27 backbone hydrogen
`bonds were observed during the period from 30 to 100 ps for more
`than 10% of the structures, with a mean occurrence of 64%.
`Among these, 14 of the 22 hydrogen bonds detected in the crystal
`structure are conserved at the 10% or greater level. Some of the
`previous MD simulations have reported conservation of 13 out
`of 16 crystallogr