Simulation from molecular dynamics trajectories

Methods for simulating restrained x>ray refinement data from molecular dynamics trajectories. 

The first computations of ionization constants of residues in proteins for structures derived from molecular dynamics trajectories were described by Wendoloski and Matthew for tuna cytochrome c. In that study, conformers were generated using molecular dynamics simulations with a range of solvents, simulating macroscopic dielectric formalisms, and one solvent model that explicitly included solvent water molecules. The authors calculated individual pR values, overall titration curves, and electrostatic potential surfaces for average structures and structures along each simulation trajectory. However, the computational scheme for predicting electrostatic interactions in proteins used by Wendoloski and Matthew was not based on a FDPB model but on the modified Tanford-Kirkwood approach, which is not discussed in this chapter. 

In molecular mechanics and molecular dynamics studies of proteins, assig-ment of standard, non-dynamical ionization states of protein titratable groups is a common practice. This assumption seems to be well justified because proton exchange times between protein and solution usually far exceed the time range of the MD simulations. We investigated to what extent the assumed protonation state of a protein influences its molecular dynamics trajectory, and how often our titration algorithm predicted ionization states identical to those imposed on the groups, when applied to a set of structures derived from a molecular dynamics trajectory [34]. As a model we took the bovine... 

It is appropriate to consider first the question of what kind of accuracy is expected from a simulation. In molecular dynamics (MD) very small perturbations to initial conditions grow exponentially in time until they completely overwhelm the trajectory itself. Hence, it is inappropriate to expect that accurate trajectories be computed for more than a short time interval. Rather it is expected only that the trajectories have the correct statistical properties, which is sensible if, for example, the initial velocities are randomly generated from a Maxwell distribution. 

An important issue, the significance of which is sometime underestimated, is the analysis of the resulting molecular dynamics trajectories. Clearly, the value of any computer simulation lies in the quality of the information extracted from it. In fact, it is good practice to plan the analysis procedure before starting the simulation, as the goals of the analysis will often detennine the character of the simulation to be performed. 

Semiclassical techniques like the instanton approach [211] can be applied to tunneling splittings. Finally, one can exploit the close correspondence between the classical and the quantum treatment of a harmonic oscillator and treat the nuclear dynamics classically. From the classical trajectories, correlation functions can be extracted and transformed into spectra. The particular charm of this method rests in the option to carry out the dynamics on the fly, using Born Oppenheimer or fictitious Car Parrinello dynamics [212]. Furthermore, multiple minima on the hypersurface can be treated together as they are accessed by thermal excitation. This makes these methods particularly useful for liquid state or other thermally excited system simulations. Nevertheless, molecular dynamics and Monte Carlo simulations can also provide insights into cold gas-phase cluster formation [213], if a reliable force field is available [189]. 

Of great interest to the molecular biologist is the relationship of protein form to function. Recent years have shown that although structural information is necessary, some appreciation of the molecular flexibility and dynamics is essential. Classically this information has been derived from the crystallographic atomic thermal parameters and more recently from molecular dynamics simulations (see for example McCammon 1984) which yield independent atomic trajectories. A diaracteristic feature of protein crystals, however, is that their diffraction patterns extend to quite limited resolution even employing SR. This lack of resolution is especially apparent in medium to large proteins where diffraction data may extend to only 2 A or worse, thus limiting any analysis of the protein conformational flexibility from refined atomic thermal parameters. It is precisely these crystals where flexibility is likely to be important in the protein function. 

From the sphere-dimer studies, two major conclusions emerge. The first is that the trajectory method can be extended to structured reactants with anisotropic reactivity and anisotropic direct forces and hydrodynamic interactions. The second major conclusion is that complicated electrostatic interactions between species with anisotropic reactivity can "steer" the approaching particles into favorable orientations and enhance the reaction rate. For these model studies, rate enhancements up to 20% have been obtained. The second conclusion is likely to be of considerable relevance to molecular biology. In the third and final series of simulations, the Brownian dynamics trajectory method is applied to a particular biological system. 

From molecular dynamics simulations, information on water exchange is available for the three regions of coordination. The simulation qualitatively confirmed the lability maximum of the first hydration shell found in the middle of the series [70]. From the trajectories obtained from the simulations, all transitions of water molecules between the first hydration shell and the bulk were identified and grouped into pairs of water molecules that together form a coupled exchange [70]. 

A number of MD studies on various unimolecular reactions over the years have shown that there can sometimes be large discrepancies (an order of magnitude or more) between reaction rates obtained from molecular dynamics simulations and those predicted by classical RRKM theory. RRKM theory contains certain assumptions about the nature of prereactive and postreactive molecular dynamics it assumes that all prereactive motion is statistical, that all trajectories will eventually react, and that no trajectory will ever recross the transition state to reform reactants. These assumptions are apparently not always valid otherwise, why would there be discrepancies between trajectory studies and RRKM theory Understanding the reasons for the discrepancies may therefore help us learn something new and interesting about reaction dynamics. 

Figure 11.1. Conolly surface plot of the oxygen atoms obtained from molecular dynamic simulation of the water/l,2-dichloroethane (DCE) junction (a). A flexible SPC model was employed for the 343 water molecules, while the 108 DCE molecules were described by a four-centre, simple charge, flexible model. A spherical probe was used with a radius of SA and the Gibbs dividing surface is located at z = 4A. Density variations of the water/vapour, water/DCE and DCE/vapour junction as obtained from the average of trajectories over 200 ps at 300 K (b). Reprinted with permission from ref.[6]. Copyright (1996) American Chemical Society.

A model of immiscible Lennard-Jones atomic solvents has been used to study the adsorption of a diatomic solute [71]. Subsequently, studies of solute transfer have been performed for atoms interacting through Lennard-Jones potentials [69] and an ion crossing an interface between a polar and a nonpolar liquid [72]. In both cases the potential of mean force experienced by the solute was computed the results of the simulation were compared with the result from the transition state theory (TST) in the first case, and with the result from a diffusion equation in the second case. The latter comparison has led to the conclusion that the rate calculated from the molecular dynamics trajectories agreed with the rate calculated using the diffusion equation, provided the mean-force potential and the diffusion coefficient were obtained from the microscopic model. 

---