Search

We present a systematic application of reaction path sampling to computer simulations of the folding of peptides and small proteins at atomic resolution in the presence of solvent. We use a simulated annealing protocol to generate an ensemble of room temperature folding trajectories of fixed length, which connect predetermined initial and final states. The trajectories are distributed according to a discretized version of the Onsager–Machlup action functional. We show that, despite the enormous practical restrictions placed on the number of time slices which can be explored, some of the basic kinetic features found experimentally for the folding of peptides and small proteins are exhibited in the nature of the reaction paths sampled. We test the method on three systems: A 12 residue α-helical peptide, a 16 residue β-hairpin peptide, and the 36 residue avian Pancreatic Polypeptide (aPP). All systems are represented at atomic resolution, and include explicit water molecules. For the 12 residue α-helix, we find that (i,i+3) hydrogen bonds can play a significant role in the folding pathway, with specific (i,i+3) bonds appearing, then transforming to the corresponding (i,i+4) hydrogen bond for some, but not all of the native hydrogen bonds. For the β-hairpin and aPP, hydrophobic interactions play a dominant role, with nonbonded interactions consistently appearing before hydrogen bonds. This is true both at the level of tertiary structure, and at the level of individual hydrogen bonds which tend to form only after stabilizing nonbonded interactions have already formed between the residues involved.

We describe a technique for generating potentials of mean force (PMF) between solutes in an aqueous solution. We first generate solute–solvent correlation functions (CF) using Monte Carlo (MC) simulations in which we place a single atom solute in a periodic boundary box containing a few hundred water molecules. We then make use of the Kirkwood superposition approximation, where the 3‐body correlation function is approximated as the product of 2‐body CFs, to describe the mean water density around two solutes. Computing the force generated on the solutes by this average water density allows us to compute potentials of mean force between the two solutes. For charged solutes an additional approximation involving dielectric screening is made, by setting the dielectric constant of water to e=80. These potentials account, in an approximate manner, for the average effect of water on the atoms. Following the work of Pettitt and Karplus [Chem. Phys. Lett. 121, 194 (1985)], we approximate the n‐body potential of mea...

The standard free energy perturbation (FEP) techniques for the calculation of conformational free energy changes of a solvated biomolecule involve long molecular dynamics (MD) simulations. We have developed a method for performing the same calculations many orders of magnitude faster. We model the average solvent density around a solute as the product of the relevant solute–solvent correlation functions (CF), following the work of Garcia, Hummer, and Soumpasis. We calculate the CF’s by running Monte Carlo simulations of a single solute atom in a box of explicit water molecules and also angular dependent CF’s for selected pairs of solute atoms. We then build the water shell around a larger solute (e.g., alanine dipeptide) by taking the product of the appropriate CF’s. Using FEP techniques we are able to calculate free energy changes as we rotate the dihedral angles of the alanine dipeptide and we find they are in close agreement with the MD results. We also compute the potential of mean force as a function...

Bond constraint algorithms for molecular dynamics typically take, as the target constraint lengths, the values of the equilibrium bond lengths defined in the potential. In Langevin form, the equations of motion are temperature dependent, which gives the average value for the individual bond lengths a temperature dependence. In addition to this, locally constant force fields can shift the local equilibrium bond lengths. To restore the average bond lengths in constrained integration to their unconstrained values, we suggest changing the constraint length used by popular constraint methods such as RATTLE [H. C. Andersen, J. Comput. Phys. 52, 23 (1983)] at each step. This allows us to more accurately capture the equilibrium bond length changes (with respect to the potential) due to the local equilibration and temperature effects. In addition, the approximations to the unconstrained nonbonded energies are closer using the dynamic constraint method than a traditional fixed constraint algorithm. The mechanism for finding the new constrained lengths involves one extra calculation of the bonded components of the force, and therefore adds O(N) time to the constraint algorithm. Since most molecular dynamics calculations are dominated by the O(N2) nonbonded forces, this new method does not take significantly more time than a fixed constraint algorithm.

We present an adaptive time stepping scheme based on the extrapolative method of Barth and Schlick [LN, J. Chem. Phys. 109, 1633 (1998)] to numerically integrate the Langevin equation with a molecular-dynamics potential. This approach allows us to use (on average) a time step for the strong nonbonded force integration corresponding to half the period of the fastest bond oscillation, without compromising the slow degrees of freedom in the problem. We show with simple examples how the dynamic step size stabilizes integration operators, and discuss some of the limitations of such stability. The method introduced uses a slightly more accurate inner integrator than LN to accommodate the larger steps. The adaptive time step approach reproduces temporal features of the bovine pancreatic trypsin inhibitor (BPTI) test system (similar to the one used in the original introduction of LN) compared to short-time integrators, but with energies that are shifted with respect to both LN, and traditional stochastic versions of Verlet. Although the introduction of longer steps has the effect of systematically heating the bonded components of the potential, the temporal fluctuations of the slow degrees of freedom are reproduced accurately. The purpose of this paper is to display a mechanism by which the resonance traditionally associated with using time steps corresponding to half the period of oscillations in molecular dynamics can be avoided. This has theoretical utility in terms of designing numerical integration schemes--the key point is that by factoring a propagator so that time steps are not constant one can recover stability with an overall (average) time step at a resonance frequency. There are, of course, limitations to this approach associated with the complicated, nonlinear nature of the molecular-dynamics (MD) potential (i.e., it is not as straightforward as the linear test problem we use to motivate the method). While the basic notion remains in the full Newtonian problem, it is easier to see the effects when damping is considered to be physical--that is, we do not view our method as a perturbation of Newtonian dynamics, we associate the damping with the environment, for example, a water bath (with gamma approximately 90 ps(-1)) [Zagrovic and Pande, J. Comp. Chem. 24, 1432 (2003)]. All stochastic approaches to MD are stabilized by large physical damping, but here, we are really using it only to show that the resonance frequency can be obtained. Another simplifying assumption used in this paper is "heavy" hydrogen (we take the hydrogen mass to be 10 amu)--the view here is that we are interested primarily in the slowest degrees of freedom, and this approach has effects similar to bond freezing and united atom treatments of hydrogen. So from the point of view of biomolecular applications, the method described here is best suited to studies in which water is not explicit (so that damping in the problem can really be viewed as environmental interaction), and the interest is in slow dynamics where the effects of hydrogen are neglectable. There are a number of parameters in the LN method and the one derived here, and we cannot in a short paper address all adjustments, so our primary goal as a first pass is to show that stability can be recovered for a set of numerically forced (and hence artificial) bond oscillations, and compare stability to fixed-step methods.

A simplified statistical mechanical model of the solid–fluid melting transition of wet lipid bilayers is developed which allows for calculation of thermodynamic fluctuation effects. The model is used to estimate the lowering of the free energy barrier for transbilayer ion permeability due to enhanced lateral compressibility near the melting point. The calculation is used to fit permeability data for dipalmitoyl phosphatidyl choline vesicles and the fit provides evidence that the bilayer is quite close to a critical temperature at which the melting transition would go from first to second order.

We report the results of multiple scattered wave SCF X‐alpha calculations of the one‐electron cross section for K‐shell photoabsorption in the molecular complexes MoO4−−, CrO4−−, and MoS4−−. We show that the method can successfully account for energy separations and relative cross sections of spectral features both below and above the K‐shell ionization threshold. Furthermore, we show: (a) that the first fairly intense peak on the low energy side of the rising edge for molybdate and chromate is due to a dipole allowed transition to a bound antibonding state of mainly nd character on the metal ion; this transition is possible because of the mixing with the ligand p orbitals having the proper T2 symmetry induced by the tetrahedral molecular potential; (b) the shoulder on the rising absorption edge can be explained by the beginning of the steplike continuum absorption when convolved with a Lorentzian function of frequency to imitate lifetime and monochromator broadening: (c) the main absorption peak is due t...