International
Tables for Crystallography Volume F Crystallography of biological macromolecules Edited by E. Arnold, D. M. Himmel and M. G. Rossmann © International Union of Crystallography 2012 
International Tables for Crystallography (2012). Vol. F, ch. 20.2, pp. 642648
https://doi.org/10.1107/97809553602060000878 Chapter 20.2. Moleculardynamics simulations of biological macromolecules ^{a}Department of Medicinal Chemistry and Molecular Pharmacology, Purdue University, West Lafayette, Indiana 47907–1333, USA Advances in the theory of atomic interactions and the increasing availability of highpower computers have led to rapid development of the moleculardynamics field and greater understanding of macromolecular motions. It is now practical to use molecular dynamics, in combination with crystallographic and NMR data, to search the large conformational space of proteins and nucleic acids to find structures consistent with the data and to improve the agreement with the data. The topics covered in this chapter include: the simulation method; potentialenergy functions; empirical parameterization of the force field; modifications in the force field for structure determination; internal dynamics and average structures; assessment of the simulation procedure; and effect of crystallographic atomic resolution on structural stability during moleculardynamics simulation. 
Molecular dynamics (MD) is the simulation of motion for a system of particles. Advances in the theory of atomic interactions and the increasing availability of highpower computers have led to rapid development of this field and greater understanding of macromolecular motions. In the earliest moleculardynamics simulations of protein molecules (McCammon et al., 1977; McCammon & Harvey, 1987), the systems were greatly simplified in order to fit within the computing capabilities of that time. Simplifications included the exclusion of water molecules and even of explicit hydrogen atoms; the effect of hydrogen atoms was built into the heavyatom properties using socalled extendedatom parameters. Simulation time periods were limited to tens of picoseconds for systems of less than 10^{3} atoms. Modern simulations, by contrast, are based on improved force fields (MacKerell et al., 1998) and benefit from considerable development in algorithms. In addition, the possible size and time period of simulations have increased by orders of magnitude; large systems of the order of 10^{4} atoms (including explicit solvent molecules) and nanosecond time periods are accessible. With dedicated computer time, the microsecond regime is possible (Duan & Kollman, 1998). Interestingly, the first 100 ps simulation of an enzyme complex was of hen eggwhite lysozyme (Post et al., 1986), the first enzyme whose structure was solved by Xray crystallography. Then the simulation required several months of dedicated time on a Cray supercomputer, but now it can be accomplished in less than a week on a common workstation.
A consequence of this enormous growth in computing power has been the particularly successful application of molecular dynamics of biological molecules to threedimensional structure determination and refinement. It is now practical to use molecular dynamics, in combination with crystallographic and NMR data, to search the large conformational space of proteins and nucleic acids to find structures consistent with the data and to improve the agreement with the data. The advantages of molecular dynamics over manual rebuilding and leastsquares refinement are the abilities to overcome the local minimum problem in an automated fashion and to search the complex conformational space of a macromolecule more extensively (Brünger et al., 1987).
Molecular mechanics, whereby the energy of the system is expressed in classical terms as a function of atomic coordinates, is well established as a useful approach for describing atomic interactions (Brooks et al., 1988; Goodfellow & Levy, 1998). Owing to the size of proteins and nucleic acids, the potentialenergy function for large biomolecules is empirically based rather than derived from quantummechanical calculations. The total force on each atom, , is calculated from the gradient, or the first derivative of this potential energy, with respect to the atomic coordinates. The motion of the atom resulting from the net force is described by Newton's equation of motion, where is the mass of atom i and is the acceleration. Integration of equation (20.2.2.1) gives where is the time step in the integration, is the atomic position at time , given the position at time t, and is the velocity. The forces on the particles change continuously so that a numerical solution of the equation is required. The Verlet algorithm (Verlet, 1967) or a variation, Leapfrog, is commonly used.
The central element of simulations is the interaction potential between atoms as a function of atomic position, r. The success of simulations in describing the average structure of proteins and other biological features suggests that such relatively simple potential functions adequately represent proteins and nucleic acids. The empirically based components of the energy function, , include geometric terms for bond lengths, bond angles and torsion angles, and nonbonding terms for steric van der Waals interactions and electrostatic interactions. A commonly used energy function is where and are force constants, and are equilibrium values for bond lengths, b, and angles, θ, respectively, and ϕ is the torsion angle of periodicity n and phase δ. The nonbonded terms depend on the interatomic distance , the dielectric constant D, the partial atomic charge , and the van der Waals parameters and . The bondstretching and anglebending contributions are represented by harmonic potentials, while the energy associated with rotation about a bond, the torsional potential, is modelled by a cosine function [equation (20.2.3.2)]. The electrostatic component of the nonbonded interactions [first term of equation (20.2.3.3)] follows Coulomb's Law, and a Lennard–Jones 6–12 potential function [second term of equation (20.2.3.3)] is used to model steric repulsion and attractive dispersion interactions. , as a sum over pairs of atoms not involved in either a bond or bond angle, requires the use of a pairwise list between atoms. The small contribution from pairs separated by a large distance allows the use of cutoff limits for this list, but at some cost in accuracy.
Initial values for the atomic coordinates and velocities are required to begin the moleculardynamics simulation. While initial coordinates are obtained from the model built into the electrondensity map, it is necessary to generate initial velocities computationally. The most common approach is to assign random values for each atom, i, consistent with the temperature chosen for the system: .
Integration of the equation of motion [equation (20.2.2.1)] also requires specification of the time step Δt. In the case of structure determination, this choice is limited only by the numerical stability of the calculation. Too large a value for Δt results in errors in the integration, manifested by a rapid and unacceptable increase in energy. Whereas a value for Δt of 1 to 2 fs is required for accurate trajectories and strict conservation of the energy, structuredetermination protocols can employ larger values and are limited only by the need for numerical stability.
Both the temperature and Δt influence the sampling rate of conformational space. Enhanced sampling increases the rate of convergence to the structure solution. Issues related to sampling are detailed in Chapter 18.2 .
In an MD simulation, the accurate and rapid calculation of longrange electrostatic interactions is a central issue for the correct physical representation of the system. The Coulombic potential [first term of equation (20.2.3.3)] has been used in most cases, but the Coulombic interactions must be limited for practical reasons to a subset of pair interactions so that longdistance interactions are truncated. More recently, the Ewald method has been implemented to avoid the need for truncating the nonbonded pair list while maintaining computational efficiency. The electrostatic energy is calculated using periodic boundary conditions^{1} and requires a double summation over all atoms in the central unit cell, as well as their infinite number of periodic images.
The total electrostatic energy of a periodic system with a neutral unit cell containing N point charges located at positions is The sum over takes into account all the periodic unit cells, and n reflects the shape of the unit cell. For a simple periodic cubic lattice, , where L is the length of the unit cell, and and are integers. The prime indicates that in the central unit cell the terms are not included. The sum is convergent for very large .
An efficient way of evaluating the triple sum above is represented by the Ewald method. In this method, each point charge is surrounded by a Gaussian charge distribution of opposite sign and same integrated magnitude as the original charge: β is inversely proportional to the width of the charge distribution and r is the position relative to the centre of the distribution. The result of adding this charge distribution is a screening of the electrostatic interactions, so that the interaction between neighbouring charges is effectively shortranged. A total screened potential is calculated by summation in the realspace lattice over the initial pointcharge distribution together with the screeningcharge distribution. To end up with the original pointcharge distribution, a cancelling distribution having the opposite sign and same functional form as the Gaussian screening distribution must be added. Together, these Gaussian distributions form a smooth varying function of r, which can be approximated by a superposition of continuous functions. The contribution of the cancelling distribution is calculated by adding the Fourier transforms of the distributions in reciprocal space at each point charge and transforming the total back into real space.
The final Ewald expression for the total electrostatic energy contains a realspace sum plus a reciprocalspace sum, minus a selfterm introduced by the interaction of the cancelling distribution with itself: (de Leeuw et al., 1980) depends on the total dipole moment of the unit cell, the shape of the macroscopic crystal lattice and the dielectric constant of the surrounding medium. The term is the complementary error function which falls to zero as increases. The term is the reciprocal Ewald sum: V is the volume and m is a vector in the reciprocal space. The Ewald sum in reciprocal space is the solution to Poisson's equation, with Gaussian charge densities as sources and with periodic boundary conditions. For large β, the only contribution to the sum in the real space comes from the minimum image terms. When β is large, the charge distributions are sharp and the summation in reciprocal space must be done over a large number of reciprocalspace vectors to achieve convergence.
Calculation of the Ewald sum in the original form is slow for large macromolecular solution simulations. The usual implementation of the Ewald summation for calculating electrostatic interaction is an order algorithm, or at best, (Christiansen et al., 1993) by adjusting β to optimize computational effort. The particle mesh Ewald (PME) method (Essmann et al., 1995; Darden et al., 1993) is an approximation of the Ewald sum that allows an order N log N algorithm for the calculation of the total electrostatic energy. This reduction in the number of steps is accomplished by choosing β large enough that atom pairs separated by distances larger than a specified cutoff (9–12 Å) make a negligible contribution to the directspace sum, which thus becomes an order N computation. Moreover, the reciprocalspace sum is expressed in terms of the electrostatic structure factors. Each atomic charge distribution is approximated by a gridded charge distribution. The resulting approximate structure factors are calculated by a threedimensional fast Fourier transform applied to the grid. Using an optimized β and grid density for each simulated system, the PME method can compute the electrostatic energy for a periodic boundary system at the same level of accuracy as the classical Ewald summation in a much shorter time (Darden et al., 1993).
For the purpose of structure determination, the potentialenergy function used for moleculardynamics calculation incorporates the information from experimental data in the form of nonphysical restraint terms. These restraint terms, introduced to bias the conformational sampling toward structures consistent with the experiment, are used in addition to the total potential function and, in some sense, can fulfil the requirements of a physical term in equation (20.2.2.2) (see below). The experimentally based restraint terms are added to the potentialenergy function to give a total effective potential, . Whereas structuredetermination protocols based on NMR data employ a number of types of restraint terms, data from Xray crystallography provide a single restraint term, , the residual between the observed and calculated structurefactor amplitudes; where where w and k are scale factors.
Considerable effort has gone into the development of a number of force fields for use in moleculardynamics simulations of biomolecules (Jorgensen & TiradoRives, 1988; MacKerell et al., 1995, 1998; van Gunsteren et al., 1996). The parameters described here are those of the CHARMM22 force field, the force field used in the XPLOR program. Estimation of the force constants, equilibrium values and nonbonding parameters in equations (20.2.3.2) and (20.2.3.3) involves a selfconsistent approach that balances the bonding and nonbonding interaction terms among the macromolecule and solvent molecules (MacKerell et al., 1998). A wide range of data are taken into account during an interactive process of optimization in order to adequately account for the extensive and correlated nature of the parameters in a consistent fashion. Smallmolecule model compounds representative of proteins or nucleic acids are considered in detail, and a hierarchical approach is applied to extend the parameters to larger molecules with minimal adjustment at the points of connection.
The empirical basis of the parameters is broad. Gasphase geometries and crystal structures are used to determine equilibrium bond lengths, bond angles, and dihedral phase and periodicity. Vibrational spectra, primarily from gasphase infrared and Raman spectroscopy, are used to fit values for the force constants. Torsionangle terms are estimated from relative energies of different conformers of model compounds, such as 4ethylimidazole and ethylbenzene, based on gasphase data. In cases where no satisfactory experimental data are available, ab initio calculations are used to obtain the required energy surfaces. Adjustments are made to describe the energy barriers and positions of saddle points, as well as the minimumenergy structures.
Optimization of the nonbonded parameters includes fitting the van der Waals and electrostatic terms of equation (20.2.3.3), while maintaining a balance among the protein–protein, water–water and protein–water interactions. The parameterization of the CHARMM22 force field is based on the water model and water–water interactions of the TIP3P model (Jorgensen et al., 1983). As such, use of this parameter set with another water model will lead to inconsistencies in the balance of intermolecular interactions. Data from dipole moments, heats and free energies of vaporization, solvation and sublimation, and molecular volumes, as well as ab initio calculations of interaction energies and geometries are used to optimize intermolecular interactions. Partial charges of atoms are determined by fitting ab initio interaction energies and geometries of smallmolecule compounds that model the peptide backbone and aminoacid side chains. Magnitudes and directions of dipolemoment values are also used to optimize partial charges. Experimental gasphase dipolemoment values are used when available, while ab initio calculated values are adopted otherwise. The van der Waals parameters are then refined by comparing results of condensedphase simulations on pure solvents with heats of vaporization and molecular volumes.
The crystallographic restraint term in the potentialenergy function, , must also be parameterized to optimize the agreement with the experimental structurefactor amplitudes while simultaneously retaining good geometry and nonbonding interactions. Optimization of involves only the estimation of w. Unlike the parameters in , w has no physical basis and is usually chosen to make the force due to balance the total force contributed by all terms in . As refinement of the structure progresses, these forces, and hence w, necessarily change since the quality, in terms of geometry and nonbonding interactions, of the structure improves and the crystallographic residual is reduced.
Simulatedannealing protocols require modification of the parameters to maintain the correct geometry and local structural integrity of the molecule in order to allow heating to very high, nonphysical temperatures for several thousand integration steps. Such modifications are acceptable in the case of structure determination since the primary goal is to define the optimum equilibrium structure in best agreement with the crystallographic or NMR data. Simulations intended to reproduce the fluctuations or dynamic properties of the system must employ carefully defined parameters without such modifications. These modifications include substantial increases in the force constants for bond lengths and angles, e.g., factors of two to ten are used in the parameters specified in the XPLOR file parallhdg.pro. A number of improper torsional terms are added to maintain proper chirality.
The specific terms in are also modified for the purpose of structure determination. In this methodology, the goal is to converge efficiently to a model that satisfies the experimental data, rather than to obtain an accurate description of the conformational surface, such as estimating fluctuations and equilibrium distributions. Alterations in include the replacement of the computationally expensive by a quartic or harmonic repulsive term, which prevents steric conflict among atoms, but ignores dispersive attraction. The electrostatic term, , is frequently excluded altogether, since the 1/r dependence of the Coulombic potential allows charge interactions to dominate the interatomic forces far from the global minimum in a fashion that hinders movement toward the global minimum. Exclusion of this important physical property of biological systems is possible, because the crystallographic structure factors contain sufficient information to reflect adequately the imprint of electrostatics on the average structure.
It is most often the goal of the structural biologist to define a single average structure of a macromolecule. The well recognized internal motions arising from thermal fluctuations of a macromolecule may be necessary for function, but, nonetheless, the methods of structure determination generally aim to model a single average structure. Internal motions range from the high frequency, small amplitude motions (i.e. those modelled by crystallographic B values) to low frequency, larger amplitude motions of loops and whole domains. Some studies (Kuriyan et al., 1986; Post, 1992) have examined the validity of the assumptions about fast timescale motions made by the methods of structure determination. It is reasonable that some of the differences between the structure solutions of a protein obtained by NMR spectroscopy and Xray crystallography are due to differences in the effects of internal motions. The application of moleculardynamics algorithms for structure determination has allowed the use of protocols that account for effects of internal motions by employing timeaveraged restraints (Schiffer et al., 1995).
Equation (20.2.3.1) is a reasonable representation of the energy function of proteins. This point is illustrated here with results from 0.8 ns moleculardynamics simulations of hen eggwhite lysozyme (PDB entry 1lzt, 1.97 Å resolution), bovine pancreas ribonuclease A (5rsa, 2.0 Å resolution), bovine αlactalbumin (1hfz, 2.3 Å resolution) and trypsin (2ptn, 1.55 Å resolution). [The coordinates for bovine αlactalbumin were kindly provided by K. R. Acharya prior to their publication (Pike et al., 1996).]
The proteins were fully hydrated, and the simulations were calculated with the CHARMM22 force field, using truncated octahedron periodic boundary conditions. The four proteins were overlaid with bulk water from an equilibrated simulation, and a 10 ps trajectory was calculated for the rearrangement of the water molecules around the protein with fixed protein atoms. The number of water molecules (4000 to 6000) required to hydrate the protein varied with the protein size (123–223 residues) and shape, with at least four layers of water molecules between the peripheral protein atoms and the walls of the boxes. The simulations were performed at constant pressure and temperature (T = 300 K and p = 1 atm) using the extendedsystem algorithms (Hoover, 1985; Nose, 1984) implemented in CHARMM. The 300 K constant temperature was maintained by coupling to an external bath, with a coupling constant of 25 ps. The SHAKE algorithm (Ryckaert et al., 1977) was used to constrain bond lengths between hydrogen atoms and heavy atoms, allowing for a time step of 2 fs in the integration of the equations of motion. A nonbonded cutoff of 12 Å was used for the Lennard–Jones potential calculation. The electrostatic forces and energies were computed using the PME method (Darden et al., 1993; Essmann et al., 1995). The PME charge grid spacing was 0.7 Å, and the charge grid was interpolated with the direct sum tolerance set to 4.0 × 10^{−6}. The nonbonded pair lists were updated every 10 steps. Structures for analysis were saved every 0.1 ps.
Fig. 20.2.7.1 shows the deviation during the simulation period of the mainchain coordinates between the simulation structures and the crystallographic starting structure. The r.m.s. coordinate deviations in the case of lysozyme are particularly stable and small: approximately 1.0 Å over the course of the trajectory. Other r.m.s. values are more typical and range from 1.0 to 2.0 Å. The time dependence of other properties, such as the radius of gyration, can also be used to follow the stability and behaviour of a trajectory. These time series, also shown in Fig. 20.2.7.1, are constant. Jumps in such a time series can be used to detect conformational transitions (Post et al., 1989).
Other experimental properties have been compared in the literature with those calculated from moleculardynamics trajectories. Of particular interest is comparing timedependent properties measured by NMR spectroscopy. An approach to calculating NMR relaxation rates was recognized early on when development of both the moleculardynamics simulations of proteins and a modelindependent theory for NMR relaxation was started (Levy, Karplus & McCammon, 1981; Levy, Karplus & Wolynes, 1981; Lipari & Szabo, 1982). Since then, the common practice of isotopic labelling of proteins for NMR structure determination has allowed the measurement of numerous NMR relaxation rates, particularly rates that characterize the motion of backbone atoms. Long simulations have been conducted to compare the calculated and experimental values (Abseher et al., 1995; Chatfield et al., 1998; Smith et al., 1995). In a particularly long simulation, an 11 ns trajectory period was used to estimate relaxation rates associated with the motions of the vectors N—H, Cα—H and C—H methyl groups from alanine and leucine (Chatfield et al., 1998). Trends in the general order of mobility of these vectors are reproduced, although a residuebyresidue comparison shows some differences.
20.2.8. Effect of crystallographic atomic resolution on structural stability during molecular dynamics
The variation in r.m.s. deviation between the initial crystallographic structure and the simulation coordinates for different protein trajectories (Fig. 20.2.7.1) raises the question of whether the atomic resolution of the starting Xray structure influences the magnitude of this deviation. In order to investigate this issue, we calculated trajectories for bovine pancreatic trypsin inhibitor (BPTI), starting with crystallographic structures determined from data at three different atomic resolutions: 1bpi at 1.1 Å resolution (Parkin et al., 1999), 6pti at 1.7 Å resolution (Wlodawer et al., 1987) and 1bhc at 2.7 Å resolution (Hamiaux et al., 1999). The errors in the atomic coordinates estimated from the Luzzati plots are 0.06 Å for the 1.1 Å resolution structure and 0.41 Å for the 2.7 Å resolution structure. The protocol described in the previous section was followed for simulations starting with each of the three crystallographic structures over a 500 ps simulation time. The net charge of +6 e on BPTI was neutralized by adding six chloride anions to the solvated protein system, thus accomplishing the ideal conditions for a PME calculation for the electrostatic interaction. The truncated octahedra contain approximately 3700 water molecules, and the total number of atoms in the simulations is over 12 000. The simulations were carried out on an eightnode IBM/SP2 and required 4.5 h of CPU time per 10 ps of dynamics run.
Rootmeansquare differences (r.m.s.d.'s) in atomic coordinates were calculated between all pairs of coordinates from the Xray structures, the energyminimized Xray structures and the 10 ps average MD structure obtained near 300 ps of the simulation period. In Table 20.2.8.1, the upper diagonal r.m.s.d. values are the mainchainatom differences, while the lower diagonal ones are the sidechainatom differences. The r.m.s.d.'s between the three Xray structures range from 0.4–0.5 Å for the mainchain atoms and 1.4–1.6 Å for the sidechain atoms. The larger r.m.s. values for averages over sidechain atoms imply alternative sidechain orientations. This degree of structural deviation is visualized with an overlay of the three Xray structures of BPTI shown on the lefthand side of Fig. 20.2.8.1. Energy minimization with respect to the CHARMM force field alters the mainchain atoms of the Xray structures by approximately 0.4 Å and increases the differences between two Xray structures to 0.6–0.7 Å. The differences in the mainchain atoms between an MD average structure and an energyminimized Xray structure are only somewhat larger: 1.0–1.1 Å. Interestingly, these values are comparable to the values obtained when comparing the three MD average structures: 0.9–1.1 Å. The mainchain structural differences among the three 10 ps average MD structures are shown on the righthand side of Fig. 20.2.8.1. The general trends observed for mainchain atoms are also found for sidechain atoms. Thus, the differences between the Xray structures increase somewhat as a result of energy minimization, and the differences between MD average structures and Xray structures (1.9–2.9 Å) are similar to those between two Xray structures (1.5–2.8 Å) or two MD average structures (approximately 2.2 Å).

Similar r.m.s. values are found if the starting velocites for a simulation are varied while maintaining the same starting coordinates (Caves et al., 1998); the r.m.s.d.'s obtained from 120 ps MD simulations were 0.7–1.1 Å for the mainchain atoms with respect to the reference Xray structure and 0.8–1.5 Å between MD individual trajectory averages. The results given in Table 20.2.8.1, together with those of Caves et al. (1998), suggest that sampling on the nanosecond timescale largely reflects the conformational variations due to thermal fluctuations that result from a potentialenergy surface with multiple minima separated by low barriers (Cooper, 1976). In this context, MD simulations starting with different Xray structures offer a more extensive sampling of the conformational space of the specific protein than simulations carried out from a single Xray structure, although this conclusion remains to be demonstrated by a more thorough analysis. Our results do not support the conclusion that overall r.m.s.d.'s between MD average structures and the starting Xray structures correlate with atomic resolution.
The r.m.s.d.'s between mainchain atoms in the starting Xray structures and simulation snapshots as a function of time are presented in Fig. 20.2.8.2. The 1.1 Å resolution structure has the most stable trajectory during the 500 ps trajectory, with an average r.m.s. value of 1.01 (9) Å. The 1.7 Å resolution structure has an r.m.s. value of 0.98 (22) Å. In this simulation, the r.m.s.d.'s fluctuate more widely from the average value, with small differences in the first 200 ps, larger ones between 200 and 400 ps, and again smaller ones in the last 100 ps. For the 2.7 Å resolution structure, the average over the simulation is 1.28 (21) Å. From the results presented here, it is concluded that the higherresolution structures are more stable during MD simulations and have a shorter equilibration period (50 ps for 1.1 Å resolution and over 300 ps for 2.7 Å resolution). This conclusion is consistent with larger errors in the atomic coordinates of Xray structures determined from lowerresolution data.
Acknowledgements
This work was supported by grants to CBP from the NIH (R01GM39478, AI39639). CBP was supported by a Research Career Development Award from the NIH (K04GM00661) and VMD is a DOE/SLOAN Postdoctoral Fellow in computational biology. The computing facilities shared by the Structural Biology group were supported by grants from the Lucille P. Markey Foundation and the Purdue University Academic Reinvestment Program.
References
Abseher, R., Ludemann, S., Schreiber, H. & Steinhauser, O. (1995). NMR crossrelaxation investigated by molecular dynamics simulation: a case study of ubiquitin in solution. J. Mol. Biol. 249, 604–624.Brooks, C. L. III, Karplus, M. & Pettitt, B. M. (1988). Proteins: a theoretical perspective of dynamics, structure and thermodynamics. Adv. Chem. Phys. 71, 1–259.
Brünger, A. T., Kuriyan, J. & Karplus, M. (1987). Crystallographic R factor refinement by molecular dynamics. Science, 235, 458–460.
Caves, L. S. D., Evanseck, J. D. & Karplus, M. (1998). Locally accessible conformations of proteins: multiple molecular dynamics simulations of crambin. Protein Sci. 7, 649–666.
Chatfield, D. C., Szabo, A. & Brooks, B. R. (1998). Molecular dynamics of staphylcoccal nuclease: comparison of simulation with 15N and 13C NMR relaxation data. J. Am. Chem. Soc. 120, 5301–5311.
Christiansen, D., Peram, J. W. & Petersen, H. G. (1993). On the fast multipole method for computing the energy of periodic assemblies of charged and dipolar particles. J. Comput. Phys. 107, 403–405.
Cooper, A. (1976). Thermodynamic fluctuations in protein molecules. Proc. Natl Acad. Sci. USA, 92, 2740–2741.
Darden, T., York, D. & Pedersen, L. (1993). Particle mesh Ewald (PME): a N log(N) method for Ewald sums in large systems. J. Chem. Phys. 98, 10089–10092.
Duan, Y. & Kollman, P. A. (1998). Pathways to a protein folding intermediate observed in a 1microsecond simulation in aqueous solution. Science, 282, 740–744.
Essmann, U., Perrera, L., Berkovitz, M. L., Darden, T., Lee, H. & Pedersen, L. G. (1995). A smooth particle mesh Ewald method. J. Chem. Phys. 103, 8577–8593.
Goodfellow, J. M. & Levy, R. M. (1998). Theory and simulation. Curr. Opin. Struct. Biol. 8, 209–210.
Gunsteren, W. F. van, Billeter, S. R., Eising, A. A., Hünenberger, P. H., Mark, A. E., Scott, W. R. P. & Tironi, I. G. (1996). Biomolecular Simulation: the GROMOS96 Manual and User Guide. Vdf Hochschulverlag, Zurich, Switzerland.
Hamiaux, C., Prangé, T., RièsKautt, M., Ducruix, A., Lafont, S., Astier, J. P. & Veesler, S. (1999). The decameric structure of bovine pancreatic trypsin inhibitor (BPTI) crystallized from thiocyanate at 2.7 Å resolution. Acta Cryst. D55, 103–113.
Hoover, W. G. (1985). Canonical dynamics: equilibrium phasespace distributions. Phys. Rev. A, 31, 1695–1697.
Jorgensen, W. L., Chandrasekhar, J., Madura, J. D., Impey, R. W. & Klein, M. L. (1983). Comparison of simple potential functions for simulating liquid water. J. Chem. Phys. 79, 926–935.
Jorgensen, W. L. & TiradoRives, J. (1988). The OPLS potential functions for proteins. Energy minimizations for crystals of cyclic peptides and crambin. J. Am. Chem. Soc. 110, 1657–1666.
Kuriyan, J., Petsko, G. A., Levy, R. M. & Karplus, M. (1986). Effect of anisotropy and anharmonicity on protein crystallographic refinement: an evaluation by molecular dynamics. J. Mol. Biol. 190, 455–479.
Leeuw, S. W. de, Perram, J. W. & Smith, E. R. (1980). Simulation of electrostatic systems in periodic boundary conditions. I. Lattice sums and dielectric constants. Proc. R. Soc. London Ser. A, 373, 27–56.
Levy, R. M., Karplus, M. & McCammon, J. A. (1981). Increase of ^{13}C NMR relaxation times in proteins due to picosecond motional averaging. J. Am. Chem. Soc. 103, 994–996.
Levy, R. M., Karplus, M. & Wolynes, P. G. (1981). NMR relaxtion parameters in molecules with internal motion: exact Langevin trajectory results compared with simplified relaxation models. J. Am. Chem. Soc. 103, 5998–6011.
Lipari, G. & Szabo, A. (1982). Modelfree approach to the interpretation of nuclear magnetic resonance relaxation in macromolecules. 1. Theory and range of validity. J. Am. Chem. Soc. 104, 4546–4559.
McCammon, J. A., Gelin, B. R. & Karplus, M. (1977). Dynamics of folded proteins. Nature (London), 267, 585–590.
McCammon, J. A. & Harvey, S. C. (1987). Dynamics of Proteins and Nucleic Acids. Cambridge University Press.
MacKerell, A., WiorkiewiczKuczera, J. & Karplus, M. (1995). An all atom empirical energy function for the simulation of nucleic acids. J. Am. Chem. Soc. 117, 11946–11975.
MacKerell, A. D. Jr, Bashford, D., Bellott, M., Dunbrack, R. L. Jr, Evanseck, J. D., Field, M. J., Fischer, S., Gao, J., Guo, H., Ha, S., JosephMcCarthy, D., Kuchnir, L., Kuczera, K., Lau, F. T. K., Mattos, C., Michnick, S., Ngo, T., Nguyen, D. T., Prodhom, B., Reiher, W. E. III, Roux, B., Schlenkrich, M., Smith, J. C., Stote, R., Straub, J. & Karplus, M. (1998). Allatom empirical potential for molecular modeling and dynamics studies of proteins. J. Phys. Chem. B, 102, 3586–3616.
Nose, S. (1984). A unified formulation of the constant temperature molecular dynamics methods. J. Chem. Phys. 81, 511–519.
Parkin, S., Rupp, B. & Hope, H. (1999). The structure of bovine pancreatic trypsin inhibitor at 125 K: definition of carboxylterminal residues glycine57 and alanine 58. In preparation.
Pike, A. C., Brew, K. & Acharya, K. R. (1996). Crystal structures of guineapig, goat and bovine αlactalbumin highlight the enhanced conformational flexibility of regions that are significant for its action in lactose synthase. Structure, 4, 691–703.
Post, C. B. (1992). Internal motional averaging and threedimensional structure determination by NMR. J. Mol. Biol. 224, 1087–1101.
Post, C. B., Brooks, B. R., Karplus, M., Dobson, C. M., Artymiuk, P. J., Cheetham, J. C. & Phillips, D. C. (1986). Molecular dynamics simulations of native and substratebound lysozyme. J. Mol. Biol. 190, 455–479.
Post, C. B., Dobson, C. M. & Karplus, M. (1989). A molecular dynamics analysis of protein structural elements. Proteins, 5, 337–354.
Ryckaert, J.P., Ciccotti, G. & Berendsen, H. J. C. (1977). Numerical integration of the Cartesian equations of motion of a system with constraints: molecular dynamics of nalkanes. J. Comput. Phys. 23, 327–341.
Schiffer, C. A., Gros, P. & van Gunsteren, W. F. (1995). Timeaveraging crystallographic refinement: possibilities and limitations using αcyclodextrin as a test system. Acta Cryst. D51, 85–92.
Smith, P. E., van Schaik, R. C., Szyperski, T., Wuthrich, K. & van Gunsteren, W. F. (1995). Internal mobility of the basic pancreatic trypsin inhibitor in solution: a comparison of NMR spin relaxation measurements and molecular dynamics simulation. J. Mol. Biol. 246, 356–365.
Verlet, L. (1967). Computer experiments on classical fluids. I. Thermodynamical properties of Lennard–Jones molecules. Phys. Rev. 159, 98–103.
Wlodawer, A., Nachman, J., Gilliland, G. L., Gallagher, W. & Woodward, C. (1987). Structure of form III crystals of bovine pancreatic trypsin inhibitor. J. Mol. Biol. 198, 469–480.