Application of isotropic periodic sum method for 4-pentyl-4′-cyanobiphenyl liquid crystal

In future large-scale molecular dynamics (MD) simulations that will use parallel computing, the isotropic periodic sum (IPS) method is expected to effectively reduce the cost of interaction calculations while maintaining adequate accuracy. To assess the accuracy of this method in estimating low-charge-density polymer systems, we performed atomistic MD simulations of the bulk state of liquid crystal systems based on 4-pentyl-4′-cyanobiphenyl (5CB). In conditions of 270 K ≤ T ≤ 320 K and a normal pressure, the temperature dependence of the density, potential energy and order parameter was estimated using the IPS and Ewald sum method. The results of the IPS method and Ewald sum were consistent within the range of error. In conditions close to the phase transition point, however, the averaged values of potential energy and order parameter had a small difference. We concluded that the fundamental physical properties for the bulk state of 5CB systems are determined reasonably by using the IPS method, at least in conditions that are not close to the phase transition point.     

Cutoff Errors in the Ewald Summation Formulae for Point Charge Systems

Abstract

Closed formulae for both real and reciprocal space parts of cutoff errors in the Ewald summation method in cubic periodic boundary conditions are derived. Such estimates are useful in tuning parameters in molecular simulations. Errors in both the electrostatic energy and forces are considered. The estimates apply to a disordered configuration of point charges and, with some limitations, also to point-charge molecular models. The accuracy of our estimates is tested and confirmed using simulated configurations of two systems (molten salt and diethylether) under a variety of conditions.     

A determination of liquid–vapour interfacial properties for methanol using a linear-combination-based isotropic periodic sum

We performed molecular dynamics simulations of the methanol liquid–vapour interfacial systems for a wide temperature range of 200–350 K. The linear-combination-based isotropic periodic sum (LIPS) method was used for estimating Coulombic interactions of the systems. The temperature dependence of the liquid and vapour phase densities, surface thickness and surface tension using the LIPS method was almost equal to that of the Ewald sum method. The temperature dependence of the electrostatic potential profile using the LIPS method was very similar to that of the Ewald sum method. We conclude that the LIPS method provides a reasonable accuracy in determining the liquid–vapour interfacial properties of polar molecules for many temperature conditions.     

Energetics of water proton configurations in gas hydrates: comparison of various water models

The total interaction energies for a large number of water proton configurations in the unit cell of hydrate structure I consisting of 46 molecules are compared for qualitatively different water models, such as SPC/E, TIP4P, TIP5P, TIP 3f and AMOEBA. All calculations were carried out using the TINKER molecular modelling package. The Ewald summation method with metallic tin-foil boundary conditions is used to account for long-range electrostatic interactions. It was established that there is a high correlation between the energies calculated using the five water models (interaction potentials). The average correlation coefficient for all pairs of potentials is equal to 0.91. Analogous calculations were carried out to evaluate the consistency of the different water models with respect to a new property of the ice-like system: the hydrogen-bond-reversal asymmetry. It was established that, for all water models, there is relatively high correlation between the energy differences for proton configurations with opposite direction of all hydrogen bonds. In this case, the average correlation coefficient is 0.77. Data for the TIP4P potential differ noticeably from the others, especially owing to the variation in the total interaction energy. The validity and usefulness of simple discrete models of inter-molecular interactions are discussed.     

On Multiple Time-step Algorithms and the Ewald Sum

Abstract

We show how standard multiple time-step algorithms devised for systems with short-range potentials can be used successfully in simulations of periodic systems with long-range (Coulombic) potentials. Three strategies for incorporating the Ewald sum into a multiple time-step algorithm are considered. These are (i) evaluation of reciprocal space terms every time-step (ii) evaluating reciprocal space terms once every n time-steps and placing these terms in with the slowly varying forces and energies (iii) a modified form of the second strategy in which primary shell (close) electrostatic interactions are evaluated directly and the more distant interactions handled by the Ewald sum (once every n time-steps). Only the first and third approaches give satisfactory thermodynamic results. The third strategy is much more efficient than the first. With the third strategy substantial savings in cpu time are acheived in both the real space and, most importantly, the reciprocal space terms of the Ewald sum. This is achieved without significant loss of accuracy or stability. Overall execution time is decreased by a factor of between 2 and 3.     

Combined use of periodic reaction field and coarse-grained molecular dynamics simulations. I. phospholipid monolayer systems

Abstract

A periodic reaction field based on a linear-combination-based isotropic periodic sum (LIPS) method was applied for coarse-grained molecular dynamics simulations of zwitterionic lipid systems. In phospholipid monolayer systems with various number of lipid molecules, the density profile, lipid orientation and surface tension were mainly calculated using the periodic reaction field and Ewald sum. The results from the periodic reaction field were almost equal to that from the Ewald sum. It is concluded that the periodic reaction field method has a great possibility to provide a high accuracy in determining coarse-grained zwitterionic lipid systems.     

Importance of explicit salt ions for protein stability in molecular dynamics simulation.

The accurate and efficient treatment of electrostatic interactions is one of the challenging problems of molecular dynamics simulation. Truncation procedures such as switching or shifting energies or forces lead to artifacts and significantly reduced accuracy. The particle mesh Ewald (PME) method is one approach to overcome these problems by providing a computationally efficient means of calculating all long-range electrostatic interactions in a periodic simulation box by use of fast Fourier transformation techniques. For the application of the PME method to the simulation of a protein with a net charge in aqueous solution, counterions are added to neutralize the system. The usual procedure is to add charge-balancing counterions close to charged residues to neutralize the protein surface. In the present article, we show that for MD simulation of a small protein of marginal stability, the YAP-WW domain, explicit modeling of 0.2 M ionic strength (in addition to the charge-balancing counterions) is necessary to maintain a stable protein structure. Without explicit ions throughout the periodic simulation box, the charge-balancing counterions on the protein surface diffuse away from the protein, resulting in destruction of the beta-sheet secondary structure of the WW domain.     

Monte-Carlo simulations of strongly interacting dipolar fluids between two conducting walls

We report Monte-Carlo simulation results for a strongly interacting dipolar soft sphere (DSS) fluid confined between two conducting, planar walls. The long-range dipolar interactions, including contributions from the “image dipoles” in the metal, are handled by mapping onto a problem with three-dimensional (3d) periodicity which can be treated by conventional Ewald summation methods. Considering two different wall separations our results indicate the occurence of wall-induced local and long-range ordering very similar to the related case of insulating walls. To understand this behavior we present ground-state (lattice) calculations on the basis of the Ewald sums, as well as simple macroscopic arguments appropriate for dipolar systems between conducting walls.     

On the Ewald artifacts in computer simulations. The test-case of the octaalanine peptide with charged termini

The treatment of electrostatic interactions in molecular simulations is of fundamental importance. Ewald and related methods are being increasingly used to the detriment of cutoff schemes, which are known to produce several artifacts. A potential drawback of the Ewald method is the spatial periodicity that is imposed to the system, which could produce artifacts when applied in the simulation of liquids. In this work we analyze the octaalanine peptide with charged termini in explicit solvent, for which severe effects due to the use of Ewald sums were predicted using continuum electrostatics. Molecular Dynamics simulations for a total of 158 nanoseconds were performed in cells of different sizes. From the comparison of the results of different system sizes, no significant periodicity-induced artifacts were observed. It is argued that in current biomolecular simulations, the incomplete sampling is likely to affect the results to a larger extent than the artifacts induced by the use of Ewald sums.     

On the Ewald Artifacts in Computer Simulations. The Test-Case of the Octaalanine Peptide With Charged Termini

Abstract

The treatment of electrostatic interactions in molecular simulations is of fundamental importance. Ewald and related methods are being increasingly used to the detriment of cutoff schemes, which are known to produce several artifacts. A potential drawback of the Ewald method is the spatial periodicity that is imposed to the system, which could produce artifacts when applied in the simulation of liquids. In this work we analyze the octaalanine peptide with charged termini in explicit solvent, for which severe effects due to the use of Ewald sums were predicted using continuum electrostatics. Molecular Dynamics simulations for a total of 158 nanoseconds were performed in cells of different sizes. From the comparison of the results of different system sizes, no significant periodicity-induced artifacts were observed. It is argued that in current biomolecular simulations, the incomplete sampling is likely to affect the results to a larger extent than the artifacts induced by the use of Ewald sums.     

Computational characterization of the sequence landscape in simple protein alphabets

We characterize the "sequence landscapes" in several simple, heteropolymer models of proteins by examining their mutation properties. Using an efficient flat-histogram Monte Carlo search method, our approach involves determining the distribution in energy of all sequences of a given length when threaded through a common backbone. These calculations are performed for a number of Protein Data Bank structures using two variants of the 20-letter contact potential developed by Miyazawa and Jernigan [Miyazawa S, Jernigan WL. Macromolecules 1985;18:534], and the 2-monomer HP model of Lau and Dill [Lau KF, Dill KA. Macromolecules 1989;22:3986]. Our results indicate significant differences among the energy functions in terms of the "smoothness" of their landscapes. In particular, one of the Miyazawa-Jernigan contact potentials reveals unusual cooperative behavior among its species' interactions, resulting in what is essentially a set of phase transitions in sequence space. Our calculations suggest that model-specific features can have a profound effect on protein design algorithms, and our methods offer a number of ways by which sequence landscapes can be quantified.     

Molecular Dynamics as a Mathematical Mapping. I. Differentiable Force Functions

Abstract

The molecular dynamics technique can be viewed as a deterministic mathematical mapping between, on one side, the force field parameters that describe the potential energy interactions and the input macroscopic conditions, and, on the other, the calculated macroscopic properties of the bulk molecular system.

The differentiability of such a mapping in the conventional molecular dynamics calculations is affected by the discontinuities in particle positions introduced by the periodic boundary conditions and the discontinuities introduced by the minimum image convention and other methods commonly employed to approximate the calculation of interparticle potential and force.

This paper proposes an alternative molecular dynamics framework based on modified force functions which are almost everywhere continuous and differentiable, and exhibit a natural periodicity. These characteristics obviate the need for both the periodic boundary conditions and the minimum image convention, as well as for any corrections for long-range interactions. They also make it possible to apply standard methods of variational calculus for the computation of partial derivatives of the molecular dynamics mapping.

The modified framework is first introduced for the case of simple monoatomic fluids where the nature of the forces exerted between any pair of two particles is identical. A more general model describing the interactions of flexible molecules is then developed. We describe the application of this approach to mixtures of alkane molecules interacting via the NERD force field.     

Monte Carlo studies on water in the dCpG/proflavin crystal hydrate

The extensive water network identified in the crystallographic studies of the dCpG/Proflavin hydrate by Neidle, Berman and Shieh (Nature 288, 129, 1980) forms an ideal test case for a) assessing the accuracy of theoretical calculations on nucleic acid--water systems based on statistical thermodynamic computer simulation, and b) the possible use of computer simulation in predicting the water positions in crystal hydrates for use in the further refinement and interpretation of diffraction data. Monte Carlo studies have been carried out on water molecules in the unit cell of dCpG/proflavin, with the nucleic acid complex fixed and the condensed phase environment of the system treated by means of periodic boundary conditions. Intermolecular interactions are described by potential functions representative of quantum mechanical calculations developed by Clementi and coworkers, and widely used in recent studies of the aqueous hydration of various forms of DNA fragments. The results are analyzed in terms of hydrogen bond topology, hydrogen bond distances and energies, mean water positions, and water crystal probability density maps. Detailed comparison of calculated and experimentally observed results are given, and the sensitivity of results to choice of potential is determined by comparison with simulation results based on a set of empirical potentials.     

Effective knowledge‐based potentials

Empirical or knowledge‐based potentials have many applications in structural biology such as the prediction of protein structure, protein–protein, and protein–ligand interactions and in the evaluation of stability for mutant proteins, the assessment of errors in experimentally solved structures, and the design of new proteins. Here, we describe a simple procedure to derive and use pairwise distance‐dependent potentials that rely on the definition of effective atomic interactions, which attempt to capture interactions that are more likely to be physically relevant. Based on a difficult benchmark test composed of proteins with different secondary structure composition and representing many different folds, we show that the use of effective atomic interactions significantly improves the performance of potentials at discriminating between native and near‐native conformations. We also found that, in agreement with previous reports, the potentials derived from the observed effective atomic interactions in native protein structures contain a larger amount of mutual information. A detailed analysis of the effective energy functions shows that atom connectivity effects, which mostly arise when deriving the potential by the incorporation of those indirect atomic interactions occurring beyond the first atomic shell, are clearly filtered out. The shape of the energy functions for direct atomic interactions representing hydrogen bonding and disulfide and salt bridges formation is almost unaffected when effective interactions are taken into account. On the contrary, the shape of the energy functions for indirect atom interactions (i.e., those describing the interaction between two atoms bound to a direct interacting pair) is clearly different when effective interactions are considered. Effective energy functions for indirect interacting atom pairs are not influenced by the shape or the energy minimum observed for the corresponding direct interacting atom pair. Our results suggest that the dependency between the signals in different energy functions is a key aspect that need to be addressed when empirical energy functions are derived and used, and also highlight the importance of additivity assumptions in the use of potential energy functions.     

Structure-guided forcefield optimization

Accurate modeling of biomolecular systems requires accurate forcefields. Widely used molecular mechanics (MM) forcefields obtain parameters from experimental data and quantum chemistry calculations on small molecules but do not have a clear way to take advantage of the information in high-resolution macromolecular structures. In contrast, knowledge-based methods largely ignore the physical chemistry of interatomic interactions, and instead derive parameters almost exclusively from macromolecular structures. This can involve considerable double counting of the same physical interactions. Here, we describe a method for forcefield improvement that combines the strengths of the two approaches. We use this method to improve the Rosetta all-atom forcefield, in which the total energy is expressed as the sum of terms representing different physical interactions as in MM forcefields and the parameters are tuned to reproduce the properties of macromolecular structures. To resolve inaccuracies resulting from possible double counting of interactions, we compare distribution functions from low-energy modeled structures to those from crystal structures. The structural and physical bases of the deviations between the modeled and reference structures are identified and used to guide forcefield improvements. We describe improvements resolving double counting between backbone hydrogen bond interactions and Lennard-Jones interactions in helices; between sidechain-backbone hydrogen bonds and the backbone torsion potential; and between the sidechain torsion potential and Lennard-Jones interactions. Discrepancies between computed and observed distributions are also used to guide the incorporation of an explicit Cα-hydrogen bond in β sheets. The method can be used generally to integrate different sources of information for forcefield improvement.     

Monte Carlo Studies on Water in the dCpG/Proflavin Crystal Hydrate

Abstract

The extensive water network identified in the crystallographic studies of the dCpG/Proflavin hydrate by Neidle, Berman and Shieh (Nature 288, 129, 1980) forms an ideal test case for a) assessing the accuracy of theoretical calculations on nucleic acid—water systems based on statistical thermodynamic computer simulation, and b) the possible use of computer simulation in predicting the water positions in crystal hydrates for use in the further refinement and interpretation of diffraction data. Monte Carlo studies have been carried out on water molecules in the unit cell of dCpG/proflavin, with the nucleic acid complex fixed and the condensed phase environment of the system treated by means of periodic boundary conditions. Intermolecular interactions are described by potential functions representative of quantum mechanical calculations developed by Clementi and coworkers, and widely used in recent studies of the aqueous hydration of various forms of DNA fragments. The results are analyzed in terms of hydrogen bond topology, hydrogen bond distances and energies, mean water positions, and water crystal probability density maps. Detailed comparison of calculated and experimentally observed results are given, and the sensitivity of results to choice of potential is determined by comparison with simulation results based on a set of empirical potentials.     

One-Dimensional Harmonic Model for Biomolecules

Following in spirit a paper by Rosen, we propose a one-dimensional harmonic model for biomolecules. Energy bands with gaps of the order of semi-conductor gaps are found. The method is discussed for general symmetric and periodic potential functions.     

Absolute binding free energy calculations using molecular dynamics simulations with restraining potentials

The absolute (standard) binding free energy of eight FK506-related ligands to FKBP12 is calculated using free energy perturbation molecular dynamics (FEP/MD) simulations with explicit solvent. A number of features are implemented to improve the accuracy and enhance the convergence of the calculations. First, the absolute binding free energy is decomposed into sequential steps during which the ligand-surrounding interactions as well as various biasing potentials restraining the translation, orientation, and conformation of the ligand are turned "on" and "off." Second, sampling of the ligand conformation is enforced by a restraining potential based on the root mean-square deviation relative to the bound state conformation. The effect of all the restraining potentials is rigorously unbiased, and it is shown explicitly that the final results are independent of all artificial restraints. Third, the repulsive and dispersive free energy contribution arising from the Lennard-Jones interactions of the ligand with its surrounding (protein and solvent) is calculated using the Weeks-Chandler-Andersen separation. This separation also improves convergence of the FEP/MD calculations. Fourth, to decrease the computational cost, only a small number of atoms in the vicinity of the binding site are simulated explicitly, while all the influence of the remaining atoms is incorporated implicitly using the generalized solvent boundary potential (GSBP) method. With GSBP, the size of the simulated FKBP12/ligand systems is significantly reduced, from approximately 25,000 to 2500. The computations are very efficient and the statistical error is small ( approximately 1 kcal/mol). The calculated binding free energies are generally in good agreement with available experimental data and previous calculations (within approximately 2 kcal/mol). The present results indicate that a strategy based on FEP/MD simulations of a reduced GSBP atomic model sampled with conformational, translational, and orientational restraining potentials can be computationally inexpensive and accurate.     

Effect of periodic box size on aqueous molecular dynamics simulation of a DNA dodecamer with particle-mesh Ewald method

The particle-mesh Ewald (PME) method is considered to be both efficient and accurate for the evaluation of long-range electrostatic interactions in large macromolecular systems being studied by molecular dynamics simulations. This method assumes "infinite" periodic boundary conditions resembling the symmetry of a crystal environment. Can such a "solid-state" method accurately portray a macromolecular solute such as DNA in solution? To address this issue, we have performed three 1500-ps PME molecular dynamics (MD) simulations, each with a different box size, on the d(CGCGA6CG)-(CGT6CGCG) DNA dodecamer. The smallest box had the DNA solvated by a layer of water molecules of at least 5 A along each orthogonal direction. The intermediate size box and the largest box had the DNA solvated by a layer of water molecules of at least 10 A and 15 A, respectively, along each orthogonal direction. The intermediate size box in the present study is similar to the box size currently chosen by most workers in the field. Based on a comparison of RMSDs and curvature for this single DNA dodecamer sequence, the larger two box sizes do not appear to afford any extra benefit over the smallest box. The implications of this finding are briefly discussed.