Counterions play a significant role in DNA structure and function, and molecular dynamics (MD) simulations offer the prospect of detailed description of the dynamical structure of ions at the molecular level. However, the motions of mobile counterions are notably slow to converge in MD on DNA. Obtaining accurate and reliable MD simulations requires knowing just how much sampling is required for convergence of each of the properties of interest. To address this issue, MD on a d(CGCGAATTCGCG) duplex in a dilute aqueous solution of water and 22 Na ؉ counterions was performed until convergence was achieved. The calculated first shell ion occupancies and DNA-Na ؉ radial distribution functions were computed as a function of time to assess convergence, and compared with relaxation times of the DNA internal parameters shift, slide, rise, tilt, roll, and twist. The sequence dependence of fractional occupancies of ions in the major and minor grooves of the DNA is examined, and the possibility of correlation between ion proximity and DNA minor groove widths is investigated. C ounterion structures and motions have been implicated in recent ideas about sequence effects on DNA structure axis curvature and ligand-induced bending (1-7). In each of these theories, the local interactions of ions with polyionic DNA complement the effects described by counterion condensation (CC) theory (8) and the Poisson Boltzmann (PB) equation (9,10). Although the effect of ions and ionic strength on DNA structural and thermodynamics properties are implied from diverse experiments, obtaining a fully detailed molecular model of the ions interacting with DNA based on these results is usually not possible. Oligonucleotide crystal structures reveal only the few ions that are ordered and can be unequivocally assigned. The positions of ions around DNA are generally underdetermined in experiments based on biophysical methods. Recent studies requiring this information have turned to large-scale molecular dynamics (MD) simulations to obtain computational models (11). However, in MD modeling, the complex aggregate of oligonucleotide, water, counterions, and coions is slow to fully stabilize, and ion motions are a rate-determining step in total convergence (12). Current reviews of MD on DNA (12)(13)(14) indicate that all simulations to date are based on considerably shorter trajectories. Therefore, we initiated a project aimed at obtaining demonstrably converged results on ion structures and motions. A related question is the sensitivity of the fast internal motions of the DNA to ion convergence. Our analysis is based on an MD simulation on the prototype B-form duplex d(CGC-GAATTCGCG) in a dilute aqueous solution of water and Na ϩ counterions carried out on a supercomputer until the point of convergence could be reliably determined. The trajectories are used to study the sequence dependence of ion distributions, the DNA-Na ϩ radial distribution functions, and the sensitivity of groove widths to ion proximity. The simulations are used as a basis for a compa...

Time-resolved Stokes-shift experiments measure the dynamics of biomolecules and of the perturbed solvent near them on subnanosecond time scales, but molecular dynamics simulations are needed to provide a clear interpretation of the results. Here we show that simulations using standard methods quantitatively reproduce the main features of TRSS experiments in DNA and provide a molecular assignment for the dynamics. The simulations reproduce the magnitude and unusual power-law dynamics of the Stokes shift seen in recent experiments [D. Andreatta, et al., J. Am. Chem. Soc. 127, 7270 (2005)]. A polarization model is introduced to eliminate cross correlations between the different components contributing to the signal. Using this model, well-defined contributions of the DNA, water and counterion to the experimental signal are extracted. Water is found to have the largest contribution and to be responsible for the power-law dynamics. The counterions have a smaller, but non-negligible contribution with a time constant of 220 ps. The contribution to the signal of the DNA itself is minor and fits a 30 ps stretched exponential. Both time-averaged and dynamic distributions are calculated. They show a small subset of ions with a different coupling, but no other evidence of substates or dynamic heterogeneity.

We are presenting POSSIM (POlarizable Simulations with Second order Interaction Model) – a software package and a set of parameters designed for molecular simulations. The key feature of POSSIM is that the electrostatic polarization is taken into account using a previously introduced fast formalism. This permits cutting computational cost of using the explicit polarization by about an order of magnitude. In this article, parameters for water, methane, ethane, propane, butane, methanol and NMA are introduced. These molecules are viewed as model systems for protein simulations. We have achieved our goal of ca. 0.5 kcal/mol accuracy for gas-phase dimerization energies and no more than 2% deviations in liquid state heats of vaporization and densities. Moreover, free energies of hydration of the polarizable methane, ethane and methanol have been calculated using the statistical perturbation theory. These calculations are a model for calculating protein pKa shifts and ligand binding affinities. The free energies of hydration were found to be 2.12 kcal/mol, 1.80 kcal/mol and −4.95 kcal/mol for methane, ethane and methanol, respectively. The experimentally determined literature values are 1.91 kcal/mol, 1.83 kcal/mol and −5.11 kcal/mol. The POSSIM average error in these absolute free energies of hydration is only about 0.13 kcal/mol. Using the statistical perturbation theory with polarizable force fields is not widespread, and we believe that this work opens road to further development of the POSSIM force field and its applications for obtaining accurate energies in protein-related computer modeling.

The A-to-B form transition has been examined in three DNA duplexes, d(CGCGAATTCGCG)(2), d(CGCGAATTGCGC), and d(CGCAAATTTCGC), using circular dichroism spectroscopy, ultraviolet resonance Raman (UVRR) spectroscopy, and molecular dynamics (MD) simulation. Circular dichroism spectra confirm that these molecules adopt the A form under conditions of reduced water activity. UVRR results, obtained under similar conditions, suggest that the transition involves a series of intermediate forms between A and B. Cooperative and distinct transitions were observed for the bases and the sugars. Independent MD simulations on d(CGCGAATTCGCG)(2) show a spontaneous change from the A to B form in aqueous solution and describe a kinetic model that agrees well with UVRR results. Based on these observations, we predict that the mechanism of the transition involves a series of A/B hybrid forms and is sequential in nature, similar to previous crystallographic studies of derivatized duplexes. A simulation in which waters were restrained in the major groove of B DNA shows a rapid, spontaneous change from B to A at reduced water activity. These results indicate that a quasiergodic sampling of the solvent distribution may be a problem in going from B to A at reduced water activity in the course of an MD simulation.

A previously introduced POSSIM (POlarizable Simulations with Second order Interaction Model) force field has been extended to include parameters for alanine peptides and protein backbones. New features were introduced into the fitting protocol, as compared to the previous generation of the polarizable force field for proteins. A reduced amount of quantum mechanical data was employed in fitting the electrostatic parameters. Transferability of the electrostatics between our recently developed NMA model and the protein backbone was confirmed. Binding energy and geometry for complexes of alanine dipeptide with a water molecule were estimated and found in a good agreement with high-level quantum mechanical results (for example, the intermolecular distances agreeing within ca. 0.06Å). Following the previously devised procedure, we calculated average errors in alanine di- and tetra-peptide conformational energies and backbone angles and found the agreement to be adequate (for example, the alanine tetrapeptide extended-globular conformational energy gap was calculated to be 3.09 kcal/mol quantim mechanically and 3.14 kcal/mol with the POSSIM force field). However, we have now also included simulation of a simple alpha-helix in both gas-phase and water as the ultimate test of the backbone conformational behavior. The resulting alanine and protein backbone force field is currently being employed in further development of the POSSIM fast polarizable force field for proteins.

A previously introduced POSSIM (POlarizable Simulations with Second order Interaction Model) force field has been extended to include parameters for small molecules serving as models for peptide and protein side-chains. Parameters have been fitted to permit reproducing many-body energies, gas-phase dimerization energies and geometries and liquid-phase heats of vaporization and densities. Quantum mechanical and experimental data have been used as the target for the fitting. The POSSIM framework combines accuracy of a polarizable force field and computational efficiency of the second-order approximation of the full-scale induced point dipole polarization formalism. The resulting parameters can be used for simulations of the parameterized molecules themselves or their analogues. In addition to this, these force field parameters are currently being employed in further development of the POSSIM fast polarizable force field for proteins.

We have studied stability of polyalanine alpha-helices with lysine residues added at C-and N-termini in gas-phase and aqueous solution. Monte Carlo simulations with the fixed-charges OPLS-AA and our polarizable POSSIM force fields were carried out. The results of the simulations confirm previously observed phenomena of the helix being stable with the LYS residue on the C-terminus and losing its helical structure if the charged LYS residue is located at the N-terminus of the polypeptide in gas-hase. Both OPLS-AA and POSSIM force fields performed essentially similarly, thus validity of the both for reproducing and predicting structures of such polypeptides has been confirmed. We have also studied the effect of replacing the normal N- and C-termini with methyl capping (this approach is often used in computational studies). Our results have demonstrated that the structure and stability of the polypeptides do not depend significantly on such a substitution although details of the resulting structure may change. The liquid-state simulations produced stable alpha-helixes regardless of the position of the protonated lysine residue. Overall, we have validated our polarizable POSSIM force field and the techniques used in the simulations, since the change of the helix structure as a function of the position of the LYS residue depends on a fine balance of energy contributions, and our methodology reproduced this balance well.

The comparison and detection of the commonalities and differences in multiple structural ensembles is an important step in the use of molecular simulations to gain insight into the conformation and dynamics of complex biomacromolecules. While the average structure is often employed as the representative of an ensemble of structures in such comparisons, dynamic molecular systems with multiple conformational substates call for a more accurate representation that captures the complete dynamical range of the ensemble. We present a probability analysis procedure based on the root-mean-square differences among the structural ensembles that efficiently and accurately performs the relevant comparison.