On the truncation of long-range electrostatic interactions in DNA

Long-range interactions are known to play an important role in highly polar biomolecules like DNA. In molecular dynamics simulations of nucleic acids and proteins, an accurate treatment of the long-range interactions are crucial for achieving stable nanosecond trajectories. In this report, we evaluate the structural and dynamic effects on a highly charged oligonucleotide in aqueous solution from different long-range truncation methods. Two group-based truncation methods, one with a switching function and one with a force-switching function were found to fail to give accurate stable trajectories close to the crystal structure. For these group-based truncation methods, large root mean square (rms) deviations from the initial structure were obtained and severe distortions of the oligonucleotide were observed. Another group-based truncation scheme, which used an abrupt truncation at 8. 0 A or at 12.0 A was also investigated. For the short cutoff distance, the conformations deviated far away from the initial structure and were significantly distorted. However, for the longer cutoff, where all necessary electrostatic interactions were included, the trajectory was quite stable. For the particle mesh Ewald (PME) truncation method, a stable DNA simulation with a heavy atom rms deviation of 1.5 A was obtained. The atom-based truncation methods also resulted in stable trajectories, according to the rms deviation from the initial B-DNA structure, of between 1.5 and 1.7 A for the heavy atoms. In these stable simulations, the heavy atom rms deviations were approximately 0.6-1.0 A lower for the bases than for the backbone. An increase of the cutoff radius from 8 to 12 A decreased the rms deviation by approximately 0.2 A for the atom-based truncation method with a force-shifting function, but increased the computational time by a factor of 2. Increasing the cutoff from 12 to 18 A for the atom-based truncation method with a force-shifting function requires 2-3 times more computational time, but did not significantly change the rms deviation. Similar rms deviations from the initial structure were found for the atom-based method with a force-shifting function and for the PME method. The computational cost was longer for the PME method with a cutoff of 12. 0 A for the direct space nonbonded calculations than for the atom-based truncation method with a force-shifting function and a cutoff of 12.0 A. If a nonperiodic boundary, e.g., a spherical boundary, was used, a considerable speedup could be achieved. From the rms fluctuations, the terminal nucleotides and especially the cytidines were found to be more flexible than the nonterminal nucleotides. The B-DNA form of the oligonucleotide was maintained throughout the simulations and is judged to depend on the parameters of the energy function and not on the truncation method used to handle the long-range electrostatic interactions. To perform accurate and stable simulations of highly charged biological macromolecules, we recommend that the atom-based force-shift method or the PME method should be used for the long-range electrostatics interactions.     

The Application of the Reaction-Field Method to the Calculation of Dielectric Constants

Abstract

The behaviour of the popular TIP3P water model has been investigated using both molecular dynamics and Monte Carlo simulation procedures. Long-range electrostatic interactions were included through a reaction-field treatment, and the nonbonded interactions were either truncated at the cutoff distance, or smoothly scaled to zero using a switching function. The thermodynamic observables, and in particular the dipole-dipole correlation functions, are found to differ between the two simulation techniques if a rigid nonbonded cutoff is applied. However, use of a switching function gives exact agreement between the simulation methodologies. This difference is ascribed to the effect of energy pumping in the molecular dynamics simulations, and suggests that dielectric constants calculated using this simulation method with the fluctuation procedure in conjunction with a reaction field should be reappraised. Thus the Monte Carlo simulation procedure offers a number of intrinsic advantages over molecular dynamics for the calculation of dielectric constants with a reaction field. The most precise value for the dielectric constant of TIP3P is calculated to be 102 ± 3 at 298 K.     

An evaluation of implicit and explicit solvent model systems for the molecular dynamics simulation of bacteriophage T4 lysozyme

In this report we examine several solvent models for use in molecular dynamics simulations of protein molecules with the Discover program from Biosym Technologies. Our goal was to find a solvent system which strikes a reasonable balance among theoretical rigor, computational efficiency, and experimental reality. We chose phage T4 lysozyme as our model protein and analyzed 14 simulations using different solvent models. We tested both implicit and explicit solvent models using either a linear distance-dependent dielectric or a constant dielectric. Use of a linear distance-dependent dielectric with implicit solvent significantly diminished atomic fluctuations in the protein and kept the protein close to the starting crystal structure. In systems using a constant dielectric and explicit solvent, atomic fluctuations were much greater and the protein was able to sample a larger portion of conformational space. A series of nonbonded cutoff distances (9.0, 11.5, 15.0, 20.0 Å) using both abrupt and smooth truncation of the nonbonded cutoff distances were tested. The method of dual cutoffs was also tested. We found that a minimum nonbonded cutoff distance of 15.0 Å was needed in order to properly couple solvent and solute. Distances shorter than 15.0 Å resulted in a significant temperature gradient between the solvent and solute. In all trajectories using the proprietary Discover switching function, we found significant denaturation in the protein backbone; we were able to run successful trajectories only in those simulations that used no switching function. We were able to significantly reduce the computational burden by using dual cutoffs and still calculate a quality trajectory. In this method, we found that an outer cutoff distance of 15.0 Å and an inner cutoff distance of 11.5 worked well. While a 10 Å shell of explicit water yielded the best results, a 6 A shell of water yielded satisfactory results with nearly a 40% reduction in computational cost. © 1994 John Wiley & Sons, Inc.     

New nonbonded interactions calculation strategy for rectangular systems. I. Preliminary molecular dynamics study of solvated Na(+) ion

A new molecular nonbonded interactions treatment strategy is proposed in the context of rectangular periodic boundary conditions simulations. Several molecular dynamics simulations are performed on a sodium ion in aqueous solution. Box sizes are modified from a cubic to a rectangular shape. The results are compared with those found using a classical spherical cutoff. This new method yields ion-oxygen radial distribution functions in good agreement with experimental results, thus showing its reliability. Severe perturbations in the structural orientation of water molecules in the first shell with the increase of the box length are observed under the classical cutoff method. However, these distorting effects are reduced with the present nonbonded interactions treatment.     

Long-range electrostatic interactions in molecular dynamics: an endothelin-1 case study

An extensive conformational search in explicit solvent was performed in order to compare the influence of different long-range electrostatic interaction treatments in molecular dynamics. The short peptide endothelin-1 was selected as the subject of molecular dynamics studies that started from both X-ray and NMR obtained structures. Electrostatic interactions were treated using two of the most common methods--residue-based cutoff and particle mesh Ewald (PME). Analyses of free energy calculations (MM-PBSA method used), secondary structure elements and hydrogen bonds were performed, and there suggested that there is no unambiguous conclusion about which of the two methods of long-range electrostatics treatment should be used in MD simulations in this case. The most reliable data was provided by a trajectory that started with the NMR structure and used the cutoff method to treat electrostatic interactions. This leads to a recommendation that the choice of electrostatics treatment should be made carefully and not automatically by choosing the PME method simply because it is the most widely used.     

Long-range Electrostatic Interactions in Molecular Dynamics: An Endothelin-1 Case Study

Abstract

An extensive conformational search in explicit solvent was performed in order to compare the influence of different long-range electrostatic interaction treatments in molecular dynamics. The short peptide endothelin-1 was selected as the subject of molecular dynamics studies that started from both X-ray and NMR obtained structures. Electrostatic interactions were treated using two of the most common methods—residue-based cutoff and particle mesh Ewald (PME). Analyses of free energy calculations (MM-PBSA method used), secondary structure elements and hydrogen bonds were performed, and there suggested that there is no unambiguous conclusion about which of the two methods of long-range electrostatics treatment should be used in MD simulations in this case. The most reliable data was provided by a trajectory that started with the NMR structure and used the cutoff method to treat electrostatic interactions. This leads to a recommendation that the choice of electrostatics treatment should be made carefully and not automatically by choosing the PME method simply because it is the most widely used.     

The effects of truncating long-range forces on protein dynamics

This paper considers the effects of truncating long-range forces on protein dynamics. Six methods of truncation that we investigate as a function of cutoff criterion of the long-range potentials are (1) a shifted potential; (2) a switching function; (3) simple atom-atom truncation based on distance; (4) simple atom-atom truncation based on a list which is updated periodically (every 25 steps); (5) simple group-group truncation based on distance; and (6) simple group-group truncation based on a list which is updated periodically (every 25 steps). Based on 70 calculations of carboxymyoglobin we show that the method and distance of long range cutoff have a dramatic effect on overall protein behavior. Evaluation of the different methods is based on comparison of a simulation's rms fluctuation about the average coordinates, the rms deviation from the average coordinates of a no cutoff simulation and from the X-ray structure of the protein. The simulations in which long-range forces are truncated by a shifted potential shows large rms deviations for cutoff criteria less than 14 A, and reasonable deviations and fluctuations at this cutoff distance or larger. Simulations using a switching function are investigated by varying the range over which electrostatic interactions are switched off. Results using a short switching function that switches off the potential over a short range of distances are poor for all cutoff distances. A switching function over a 5-9 A range gives reasonable results for a distance-dependent dielectric, but not using a constant dielectric. Both the atom-atom and group-group truncation methods based on distance shows large rms deviation and fluctuation for short cutoff distances, while for cutoff distances of 11 A or greater, reasonable results are achieved. Although comparison of these to distance-based truncation methods show surprisingly larger rms deviations for the group-group truncation, contrary to simulation studies of aqueous ionic solutions. The results of atom-atom or group-group list-based simulations generally appear to be less stable than the distance-based simulations, and require more frequent velocity scaling or stronger coupling to a heat bath.     

Structure and dynamics of calcium-activated calmodulin in solution.

Two 4-ns molecular dynamics simulations of calcium loaded calmodulin in solution have been performed, using both standard nonbonded cutoffs and Ewald summation to treat electrostatic interactions. Our simulation results are generally consistent with solution experimental studies of calmodulin structure and dynamics, including NMR, cross-linking, fluorescence and x-ray scattering. The most interesting result of the molecular dynamics simulations is the detection of large-scale structural fluctuations of calmodulin in solution. The globular N- and C-terminal domains tend to move approximately like rigid bodies, with fluctuations of interdomain distances within a 7 A range and of interdomain angles by up to 60 deg. Essential dynamics analysis indicates that the three dominant types of motion involve bending of the central helix in two perpendicular planes and a twist in which the domains rotate in opposite directions around the central helix. In the more realistic Ewald trajectory the protein backbone remains mostly within a 2-3 A root-mean-square distance from the crystal structure, the secondary structure within the domains is conserved and middle part of the central helix becomes disordered. The central helix itself exhibits limited fluctuations, with its bend angle exploring the 0-50 degrees range and the end-to-end distance falling in 39-43 A. The results of the two simulations were similar in many respects. However, the cutoff trajectory exhibited a larger deviation from the crystal, loss of several helical hydrogen bonds in the N-terminal domain and lack of structural disorder in the central helix.     

Molecular dynamics simulations of lipid bilayers: major artifacts due to truncating electrostatic interactions

We study the influence of truncating the electrostatic interactions in a fully hydrated pure dipalmitoylphosphatidylcholine (DPPC) bilayer through 20 ns molecular dynamics simulations. The computations in which the electrostatic interactions were truncated are compared to similar simulations using the particle-mesh Ewald (PME) technique. All examined truncation distances (1.8-2.5 nm) lead to major effects on the bilayer properties, such as enhanced order of acyl chains together with decreased areas per lipid. The results obtained using PME, on the other hand, are consistent with experiments. These artifacts are interpreted in terms of radial distribution functions g(r) of molecules and molecular groups in the bilayer plane. Pronounced maxima or minima in g(r) appear exactly at the cutoff distance indicating that the truncation gives rise to artificial ordering between the polar phosphatidyl and choline groups of the DPPC molecules. In systems described using PME, such artificial ordering is not present.     

Explicit-solvent molecular dynamics simulations of a DNA tetradecanucleotide duplex: lattice-sum versus reaction-field electrostatics

Five long-timescale (10 ns) explicit-solvent molecular dynamics simulations of a DNA tetradecanucleotide dimer are performed using the GROMOS 45A4 force field and the simple-point-charge water model, in order to investigate the effect of the treatment of long-range electrostatic interactions as well as of the box shape and size on the structure and dynamics of the molecule (starting from an idealised B-DNA conformation). Long-range electrostatic interactions are handled using either a lattice-sum (LS) method (particle–particle–particle–mesh; one simulation performed within a cubic box) or a cutoff-based reaction-field (RF) method (four simulations, with long-range cutoff distances of 1.4 or 2.0 nm and performed within cubic or truncated octahedral periodic boxes). The overall double-helical structure, including Watson–Crick (WC) base-pairing, is well conserved in the simulation employing the LS scheme. In contrast, the WC base-pairing is nearly completely disrupted in the four simulations employing the RF scheme. These four simulations result in highly distorted compact (cutoff distance of 1.4 nm) or extended (cutoff distance of 2 nm) structures, irrespective of the shape and size of the computational box. These differences observed between the two schemes seem correlated with large differences in the radial distribution function between charged entities (backbone phosphate groups and sodium counterions) within the system.     

Insights into the mode of flavin mononucleotide binding and catalytic mechanism of bacterial chromate reductases: A molecular dynamics simulation study

Enzymes from natural sources protect the environment via complex biological mechanisms, which aid in reductive immobilization of toxic metals including chromium. Nevertheless, progress was being made in elucidating high-resolution crystal structures of reductases and their binding with flavin mononucleotide (FMN) to understand the underlying mechanism of chromate reduction. Therefore, herein, we employed molecular dynamics (MD) simulations, principal component analysis (PCA), and binding free energy calculations to understand the dynamics behavior of these enzymes with FMN. Six representative chromate reductases in monomeric and dimeric forms were selected to study the mode, dynamics, and energetic component that drive the FMN binding process. As evidenced by MD simulation, FMN prefers to bind the cervix formed between the catalytic domain surrounded by strong conserved hydrogen bonding, electrostatic, and hydrophobic contacts. The slight movement and reorientation of FMN resulted in breakage of some crucial H-bonds and other nonbonded contacts, which were well compensated with newly formed H-bonds, electrostatic, and hydrophobic interactions. The critical residues aiding in tight anchoring of FMN within dimer were found to be strongly conserved in the bacterial system. The molecular mechanics combined with the Poisson-Boltzmann surface area binding free energy of the monomer portrayed that the van der Waals and electrostatic energy contribute significantly to the total free energy, where, the polar solvation energy opposes the binding of FMN. The proposed proximity relationships between enzyme and FMN binding site presented in this study will open up better avenues to engineer enzymes with optimized chromate reductase activity for sustainable bioremediation of heavy metals.     

Molecular Dynamics Simulations of Double-Stranded DNA in an Explicit Solvent Model with the Zero-Dipole Summation Method

Molecular dynamics (MD) simulations of a double-stranded DNA with explicit water and small ions were performed with the zero-dipole summation (ZD) method, which was recently developed as one of the non-Ewald methods. Double-stranded DNA is highly charged and polar, with phosphate groups in its backbone and their counterions, and thus precise treatment for the long-range electrostatic interactions is always required to maintain the stable and native double-stranded form. A simple truncation method deforms it profoundly. On the contrary, the ZD method, which considers the neutralities of charges and dipoles in a truncated subset, well reproduced the electrostatic energies of the DNA system calculated by the Ewald method. The MD simulations using the ZD method provided a stable DNA system, with similar structures and dynamic properties to those produced by the conventional Particle mesh Ewald method.     

Determination of the protonation state of the Asp dyad: conventional molecular dynamics versus thermodynamic integration

The protonation state of the Asp dyad is important as it can reveal enzymatic mechanisms, and the information this provides can be used in the development of drugs for proteins such as memapsin 2 (BACE-1), HIV-1 protease, and rennin. Conventional molecular dynamics (MD) simulations have been successfully used to determine the preferred protonation state of the Asp dyad. In the present work, we demonstrate that the results obtained from conventional MD simulations can be greatly influenced by the particular force field applied or the values used for control parameters. In principle, free-energy changes between possible protonation states can be used to determine the protonation state. We show that protonation state prediction by the thermodynamic integration (TI) method is insensitive to force field version or to the cutoff for calculating nonbonded interactions (a control parameter). In the present study, the protonation state of the Asp dyad predicted by TI calculations was the same regardless of the force field and cutoff value applied. Contrary to the intuition that conventional MD is more efficient, our results clearly show that the TI method is actually more efficient and more reliable for determining the protonation state of the Asp dyad.     

Developing polarized protein-specific charges for protein dynamics: MD free energy calculation of pKa shifts for Asp26/Asp20 in thioredoxin

Ab initio quantum mechanical calculation of protein in solution is carried out to generate polarized protein-specific charge(s) (PPC) for molecular dynamics (MD) stimulation of protein. The quantum calculation of protein is made possible by developing a fragment-based quantum chemistry approach in combination with the implicit continuum solvent model. The computed electron density of protein is utilized to derive PPCs that represent the polarized electrostatic state of protein near the native structure. These PPCs are atom-centered like those in the standard force fields and are thus computationally attractive for molecular dynamics simulation of protein. Extensive MD simulations have been carried out to investigate the effect of electronic polarization on the structure and dynamics of thioredoxin. Our study shows that the dynamics of thioredoxin is stabilized by electronic polarization through explicit comparison between MD results using PPC and AMBER charges. In particular, MD free-energy calculation using PPCs accurately reproduced the experimental value of pK_a shift for ionizable residue Asp²⁶ buried inside thioredoxin, whereas previous calculations using standard force fields overestimated pK_a shift by twice as much. Accurate prediction of pK_a shifts by rigorous MD free energy simulation for ionizable residues buried inside protein has been a significant challenge in computational biology for decades. This study presented strong evidence that electronic polarization of protein plays an important role in protein dynamics.     

Cutoff size does strongly influence molecular dynamics results on solvated polypeptides.

The behavior of a 17-residue model peptide is analyzed by means of molecular dynamics simulations including explicitly more than a thousand water molecules. On the basis of the charge-group concept, Coulomb interactions are truncated for three values of the cutoff radius: 0.6, 1.0, and 1.4 nm. It is found that the stability of an alpha-helix, which acts as a common starting configuration, is a function of the cutoff size. While the overall stability of the helix is conserved in a simulation using a cutoff of 1.0 nm, it is lost within a very short period of 100 ps when the cutoff is increased to 1.4 nm. This demonstrates that the commonly used cutoff size of 1.0 nm is inappropriate because it does not ensure the convergence of Coulomb interactions. In order to permit an independent judgment, we have performed a 225-ps simulation using the Ewald summation technique, which is more elaborate but circumvents the problem to find an appropriate cutoff value. In contrast to the 1.4-nm cutoff trajectory, the Ewald technique simulation conserves the helical character of the peptide conformation. This demonstrates that even 1.4 nm is too short a cutoff. Due to the fundamental uncertainty introduced by the use of a simple cutoff, this truncation scheme seems questionable for molecular dynamics simulations of solvated biomolecules.     

Electrostatic interactions between α-helix dipoles in crystals of an uncharged helical undecapeptide

The electrostatic interactions between α-helix dipoles in the crystals of an uncharged helical undecapeptide have been studied in detail. The electrostatic interaction energy between one helix dipole and its 26 nearest neighbors is approximately ?23 kcal mol^?1. A very similar result is obtained when calculating the interactions between one helix dipole and all 988 helix dipoles occurring within a distance of 75 Å. It therefore appears that in these crystals of completely uncharged molecules large, favorable electrostatic interactions occur.     

Analysis of the effect of electrostatic energy truncation in molecular dynamics simulations of immunoglobulin G light chain dimer

Molecular dynamics (MD) simulations of immunoglobulin G (IgG) light chain dimer using particle mesh Ewald (PME) and cutoff methods of treating electrostatic interactions were performed. The results indicate that structural parameters (RMSD, radius of gyration, solvent accessible surface) are very similar for both schemes; however, PME simulation shows increased mobility of side chains. This leads to larger fluctuations in the distance between the monomers in the dimer molecule, and, as a consequence, results in decreased number of interactions across the dimer interface. The wall clock time of the simulations was also compared. It was shown that the PME method is approximately 30% faster than the cutoff method for the system studied on a single processor.Figure Backbone order parameters for PME (red) and cutoff (green) calculations. Thick, horizontal lines show stable secondary structures