POSSIM: Parameterizing Complete Second-Order Polarizable Force Field for Proteins.

Previously, we reported development of a fast polarizable force field and software named POSSIM (POlarizable Simulations with Second order Interaction Model). The second-order approximation permits the speed up of the polarizable component of the calculations by ca. an order of magnitude. We have now expanded the POSSIM framework to include a complete polarizable force field for proteins. Most of the parameter fitting was done to high-level quantum mechanical data. Conformational geometries and energies for dipeptides have been reproduced within average errors of ca. 0.5 kcal/mol for energies of the conformers (for the electrostatically neutral residues) and 9.7° for key dihedral angles. We have also validated this force field by running Monte Carlo simulations of collagen-like proteins in water. The resulting geometries were within 0.94 Å root-mean-square deviation (RMSD) from the experimental data. We have performed additional validation by studying conformational properties of three oligopeptides relevant in the context of N-glycoprotein secondary structure. These systems have been previously studied with combined experimental and computational methods, and both POSSIM and benchmark OPLS-AA simulations that we carried out produced geometries within ca. 0.9 Å RMSD of the literature structures. Thus, the performance of POSSIM in reproducing the structures is comparable with that of the widely used OPLS-AA force field. Furthermore, our fitting of the force field parameters for peptides and proteins has been streamlined compared with the previous generation of the complete polarizable force field and relied more on transferability of parameters for nonbonded interactions (including the electrostatic component). The resulting deviations from the quantum mechanical data are similar to those achieved with the previous generation; thus, the technique is robust, and the parameters are transferable. At the same time, the number of parameters used in this work was noticeably smaller than that of the previous generation of our complete polarizable force field for proteins; thus, the transferability of this set can be expected to be greater, and the danger of force field fitting artifacts is lower. Therefore, we believe that this force field can be successfully applied in a wide variety of applications to proteins and protein-ligand complexes.

Polarizable simulations with second order interaction model (POSSIM) force field: developing parameters for protein side-chain analogues.

A previously introduced polarizable simulations with second-order interaction model (POSSIM) 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 used in further development of the POSSIM fast polarizable force field for proteins.

Importance of electrostatic polarizability in calculating cysteine acidity constants and copper(I) binding energy of Bacillus subtilis CopZ.

CopZ is a copper chaperone from Bacillus subtilis. It is an important part of Cu(I) trafficking. We have calculated pK(a) values for the CXXC motif of this protein, which is responsible for the Cu(I) binding, and the Cu(I) binding constants. Polarizable and fixed-charges formalisms were used, and solvation parameters for the both models have been refitted. We had to partially redevelop parameters for the protonated and deprotonated cysteine residues. We have discovered that the polarizable force field (PFF) is qualitatively superior and allows a uniformly better level of energetic results. The PFF pK(a) values for cysteine are within about 0.8-2.8 pH units of the experimental data, while the fixed-charges OPLS formalism yields errors of up to tens of units. The PFF magnitude of the copper binding energy is about 10 kcal/mol or 50% higher than the experimental value, while the using the refitted OPLS parameters leads to an overall positive binding energy, thus predicting no thermodynamically stable complex. At the same time, the agreement of the polarizable S···Cu(I) distances with the experimental results is within 0.08 Å range, and the nonpolarizable calculations lead to an error of about 0.4 Å. Moreover, the accuracy of the PFF has been achieved without any explicit fitting to either pK(a) or CopZ···Cu(I) binding energies. We believe that this makes our polarizable technique a choice method in reproducing protein-copper binding and further supports the notion that explicit treatment of electrostatic polarization is crucial in many biologically relevant studies, especially ion binding and transport.

Effects of lysine substitution on stability of polyalanine alpha-helix.

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.

Polarizable Simulations with Second order Interaction Model (POSSIM) force field: Developing parameters for alanine peptides and protein backbone.

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.

Electrostatic polarization is crucial in reproducing Cu(I) interaction energies and hydration.

We have explored the suitability of fixed-charges and polarizable force fields for modeling interactions of the monovalent Cu(I) ion. Parameters for this ion have been tested and refitted within the fixed-charges OPLS-AA and polarizable force field (PFF) frameworks. While this ion plays an important role in many protein interactions, the attention to it in developing empirical force fields is limited. Our PFF parameters for the copper ion worked very well for the Cu(I) interactions with water, while both the original OPLS2005 and our refitted OPLS versions moderately underestimated the copper-water interaction energy. However, the greatest problem in using the nonpolarizable fixed-charges OPLS force field was observed while calculating interaction energies and distances for Cu(I)-benzene complexes. The OPLS2005 model underestimates the interaction energy by a factor of 4. Refitting the OPLS parameters reduced this underestimation to a factor of 2.2-2.4, but only at a cost of distorting the complex geometry. At the same time, the polarizable calculations had an error of about 4%. Moreover, we then used the PFF and nonpolarizable refitted OPLS models for finding free energy of hydration for copper ion via molecular dynamics simulations. While the OPLS calculations lead to a 22% error in the solvation energy, the PFF result was off by only 1.8%. This was achieved with no refitting of the parameters but simply by employing the model developed for the Cu(I) interaction with a single water molecule. We believe that the presented results not only lead to a conclusion about a qualitatively greater suitability of polarizable force fields for simulating molecular interactions with ions but also attest to the excellent level of transferability of PFF parameters.

Computational prediction and analysis of the DR6-NAPP interaction.

Alzheimer's disease (AD) is a progressive neurodegenerative disorder that involves a devastating clinical course and that lacks an effective treatment. A biochemical model for neuronal development, recently proposed by Nikolaev et al., that may also have implications for AD, hinges on a novel protein­protein interaction between the death cell receptor 6 (DR6) ectodomain and an Nterminal fragment of amyloid precursor protein (NAPP), specifically, the growth factor-like domain of NAPP (GFD NAPP). Given all of this, we used a pure computational work-flow to dock a binding competent homology model of the DR6 ectodomain to a binding competent crystal structure of GFD NAPP. The DR6 homology model was built according to a template supplied by the neurotrophin p75 receptor. The best docked model was selected according to an empirical estimate of the binding affinity and represents a high quality model of probable structural accuracy, especially with respect to the residue-level contribution of GFD NAPP. The final model was tested and verified against a variety of biophysical and theoretical data sets. Particularly, worth noting is the excellent observed agreement between the theoretically calculated DR6­GFD NAPP binding free energy and the experimental quantity. The model is used to provide a satisfying structural and energetic interpretation of DR6­GFD NAPP binding and to suggest the possibility of and a mechanism for spontaneous apoptosis. The evidence suggests that the DR6­NAPP model proposed here is of probable accuracy and that it will prove useful in future studies, modeling work, and structure-based AD drug design.

Dynamics of water and ions near DNA: comparison of simulation to time-resolved stokes-shift experiments.

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 [ Andreatta, D., et al. J. Am. Chem. Soc. 2005, 127, 7270 ]. 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 rate heterogeneity.

Polarizable Simulations with Second order Interaction Model - force field and software for fast polarizable calculations: Parameters for small model systems and free energy calculations.

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.

Relaxation dynamics of nucleosomal DNA.

Recent experimental and theoretical evidence demonstrates that proteins and water in the hydration layer can follow complex stretched exponential or power law relaxation dynamics. Here, we report on a 50 ns all atom molecular dynamics (MD) simulation of the yeast nucleosome, where the interactions between DNA, histones, surrounding water and ions are explicitly included. DNA interacts with the histone core in 14 locations, approximately every 10.4 base pairs. We demonstrate that all sites of interaction exhibit anomalously slow power law relaxation, extending up to 10 ns, while fast exponential relaxation dynamics of hundreds of picoseconds applies to DNA regions outside these locations. The appearance of 1/f(alpha) noise or pink noise in DNA dynamics is ubiquitous. For histone-bound nucleotide dynamics alpha --> 1 and is a signature of complexity of the protein-DNA interactions. For control purposes two additional DNA simulations free of protein are conducted. Both utilize the same sequence of DNA, as found the in the nucleosome. In one simulation the initial conformation of the double helix is a straight B-form. In the other, the initial conformation is super helical. Neither of these simulations exhibits the variation of alpha as a function of position, the measure of power law for dynamical behavior, which we observe in the nucleosome simulation. The unique correspondence (high alpha to DNA-histone interaction sites, low alpha to free DNA sites), suggests that alpha may be an important and new quantification of protein-DNA interactions for future experiments.

Spectroscopic and molecular dynamics evidence for a sequential mechanism for the A-to-B transition in DNA.

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.

Root mean square deviation probability analysis of molecular dynamics trajectories on DNA.

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.

Ion motions in molecular dynamics simulations on DNA.

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.