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.     

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.     

Molecular dynamics simulations of charged and neutral lipid bilayers: treatment of electrostatic interactions

Molecular dynamics (MD) simulations complement experimental methods in studies of the structure and dynamics of lipid bilayers. The choice of algorithms employed in this computational method represents a trade-off between the accuracy and real calculation time. The largest portion of the simulation time is devoted to calculation of long-range electrostatic interactions. To speed-up evaluation of these interactions, various approximations have been used. The most common ones are the truncation of long-range interactions with the use of cut-offs, and the particle-mesh Ewald (PME) method. In this study, several multi-nanosecond cut-off and PME simulations were performed to establish the influence of the simulation protocol on the bilayer properties. Two bilayers were used. One consisted of neutral phosphatidylcholine molecules. The other was a mixed lipid bilayer consisting of neutral phosphatidylethanolamine and negatively charged phosphatidylglycerol molecules. The study shows that the cut-off simulation of a bilayer containing charge molecules generates artefacts; in particular the mobility and order of the charged molecules are vastly different from those determined experimentally. In the PME simulation, the bilayer properties are in general agreement with experimental data. The cut-off simulation of bilayers containing only uncharged molecules does not generate artefacts, nevertheless, the PME simulation gives generally better agreement with experimental data.     

Molecular dynamics simulations on HIV-1 Tat

Molecular dynamics simulations are used to investigate dynamics and intramolecular interactions of the HIV-1 transactivator (Tat) in aqueous solution. The calculations are based on the AMBER force field with particle mesh Ewald treatment for long-range electrostatics. The Tat structure exhibits a large flexibility, consistent with its absence of secondary structure elements. From an analysis of the correlation matrix and of electrostatic interactions we suggest that segments expressed by the two exons (amino acids 1-72 and 73-86, respectively) exhibit rather separated dynamic and energetic properties. We also identify intramolecular interactions of importance for structure stabilization. In particular, significant electrostatic interactions are recognized between the N-terminus and the basic domain of the protein, consistent with site-directed mutagenesis performed in this work.     

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.     

Molecular dynamics simulations of a protein-protein dimer: particle-mesh Ewald electrostatic model yields far superior results to standard cutoff model.

In this article we present two 1000 ps molecular dynamics simulations on the rat micro-glutathione S-transferase dimeric enzyme in complex with the product 1-(S-glutathionyl)-2,4-dinitrobenzene, in a periodic box with explicit solvent molecules, and investigate the effect of long-range electrostatics models on the structure and dynamics of the dimer and its components. One simulation used the standard cutoff method (10A), whilst the other used the particle-mesh Ewald (PME) method. We monitored the root mean-square atomic deviation (RMSD) from the initial crystal structure to examine the convergence of both simulations, as well as several other structural parameters such as the distance between active sites, rigid body rotation between domains in subunits, radius of gyration, B-factors, number of hydrogen bonds and salt bridges and solvent-accessible surface area. For example, with the PME method, the dimer structure remains much closer to the initial crystallographic structure with an average RMSD of 1.3A +/- 0.1A and 1.0A +/- 0.1A for all heavy and backbone atoms, respectively, in the last 200 ps; the respective values for the cutoff simulation are 4.7A +/- 0.3A and 4.2A +/- 0.3A. The large deviations observed in the cutoff simulation severely affected the stability of the enzyme dimer and its complex with the bound product. This finding is contrary to that found in a similar study of the monomeric protein ubiquitin [Fox, T. & Kollman, P. A. Proteins Struct. Func. Genet. 25, 315-334 (1996)]. Unlike the earlier published work, the present study provides evidence that the standard cutoff method is not generally valid for the study of protein complexes, or their subunits.     

Comparison of different schemes to treat long-range electrostatic interactions in molecular dynamics simulations of a protein crystal

Eight molecular dynamics simulations of a ubiquitin crystal unit cell were performed to investigate the effect of different schemes to treat the long-range electrostatic interactions as well as the need to include counter ions. A crystal system was chosen as the test system, because the higher charge density compared with a protein in solution makes it more sensitive to the way of treating the electrostatic interactions. Three different schemes of treating the long-range interactions were compared: straight cutoff, reaction-field approximation, and a lattice-sum method (P3M). For each of these schemes, two simulations were performed, one with and one without the counter ions. Two additional simulations with a reaction-field force and different initial placements of the counter ions were performed to examine the effect of the initial positions of the ions. The inclusion of long-range electrostatic interactions using either a reaction-field or a lattice-sum method proved to be necessary for the simulation of crystals. These two schemes did not differ much in their ability to reproduce the crystallographic structure. The inclusion of counter ions, on the other hand, seems not necessary for obtaining a stable simulation. The initial positions of the ions have a visible but small effect on the simulation.     

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.     

Long-range electrostatic effects on peptide folding.

The reversible folding/unfolding of a short peptide in solution is studied by molecular dynamics simulations. The effects of long-range electrostatic interactions are examined and found to be important both for the equilibrium between folded and unfolded states and the dynamics of the folding process. The neglect of long-range electrostatics leads to an increased population of unfolded states and increased structural fluctuations. When such interactions are taken into account, the peptide unfolds and folds to the experimentally determined structure several times during a 25 ns trajectory, with approximately equal populations of folded and unfolded states in the neighborhood of its proposed melting temperature. The effect of using spherical boundary conditions rather than periodic ones does not appear to have any major effect on the folding dynamics.     

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     

Charge grouping approaches to calculation of electrostatic forces in molecular dynamics of macromolecules.

Calculation of long-range electrostatic interactions is the most time-consuming step in theoretical simulation of the structure and dynamics of macromolecules. In practice very short cutoff distances are used, which may distort the behavior of the model system. We describe two accurate approaches to calculation of electrostatic forces based on hierarchical grouping of charges into cubes. The first is similar to the O(NlogN) algorithm developed by Barnes, J. and Hut, P., Nature (London) 324, 446-449 (1986), for simulation of a gravitational motion of N bodies. The second approach we formulate for a system with periodic boundary conditions in the nearest image approximation. The calculation of electrostatic interactions and a charge grouping procedure are faster than O(N2). The average inaccuracy in the force introduced by the grouping does not exceed 1%. We describe a small modification of the same approach which makes it suitable for long strongly charged polymers as well. This accurate approach to calculation of electrostatic interactions is illustrated with an example of the dynamics of ions near DNA. Quick equilibration of the ionic distribution is observed during molecular dynamics simulation if electrostatic forces are properly calculated, while the behavior and distribution of ions are less realistic when the conventional cutoff distances are used.     

Restrained molecular dynamics of solvated duplex DNA using the particle mesh Ewald method

Restrained and unrestrained aqueous solution molecular dynamics simulations applying the particle mesh Ewald (PME) method to DNA duplex structures previously determined via in vacuo restrained molecular dynamics with NMR-derived restraints are reported. Without experimental restraints, the DNA decamer, d(CATTTGCATC)d(GATGCAAATG) and trisdecamer, d(AGCTTGCCTTGAG)d(CTCAAGGCAAGCT), structures are stable on the nanosecond time scale and adopt conformations in the B-DNA family. These free DNA simulations exhibit behavior characteristic of PME simulations previously performed on DNA sequences, including a low helical twist, frequent sugar pucker transitions, BI- BII(–) transitions and coupled crankshaft (–) motion. Refinement protocols similar to the original in vacuo restrained molecular dynamics (RMD) refinements but in aqueous solution using the Cornell et al. force field [Cornell et al. (1995) J. Am. Chem. Soc., 117, 5179–5197] and a particle mesh Ewald treatment produce structures which fit the restraints very well and are very similar to the original in vacuo NMR structure, except for a significant difference in the average helical twist. Figures of merit for the average structure found in the RMD PME decamer simulations in solution are equivalent to the original in vacuo NMR structure while the figures of merit for the free MD simulations are significantly higher. The free MD simulations with the PME method, however, lead to some sequence-dependent structural features in common with the NMR structures, unlike free MD calculations with earlier force fields and protocols. There is some suggestion that the improved handling of electrostatics by PME improves long-range structural aspects which are not well defined by the short-range nature of NMR restraints.     

Acidic range titration of HEWL using a constant-pH molecular dynamics method

In this work, we present the first application to a protein of the stochastic constant-pH molecular dynamics (MD) method with the inclusion of proton tautomerism. The acidic titration of HEWL was performed under different conditions. Both generalized reaction field (GRF) and particle mesh Ewald (PME) methods were used in the treatment of the long range electrostatics and, even though the PME simulations revealed to be more stable, the better results were obtained using GRF (pK(a) RMSD of 0.82 for GRF and 1.13 for PME). The results using PME at different dielectric constants (2, 4, and 8) also revealed that there was no significant improvement in pK(a)'s prediction upon increasing the dielectric constant. The secondary structure analysis of HEWL revealed a remarkably stable protein in the acidic pH range. The beta-sheet strands (unlike the alpha-helices) seem to be destabilized upon pH decrease, suggesting that the beta-domain is less stable than the alpha-domain. The four principal alpha-helices were also ordered according to their stability in the acidic pH range and the results (4 < 1 < 2 approximately = 3) were consistent with the ones obtained in thermal denaturation studies.     

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.     

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.     

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.     

The range and shielding of dipole-dipole interactions in phospholipid bilayers

In molecular dynamics simulations of lipid bilayers, the structure is sensitive to the precise treatment of electrostatics. The dipole-dipole interactions between headgroup dipoles are not long-ranged, but the area per lipid and, through it, other properties of the bilayer are very sensitive to the detailed balance between the perpendicular and in-plane components of the headgroup dipoles. This is affected by the detailed properties of the cutoff scheme or if long-range interactions are included by Ewald or particle-mesh Ewald techniques. Interaction between the in-plane components of the headgroup dipoles is attractive and decays as the inverse sixth power of distance. The interaction is screened by the square of a dielectric permittivity close to the value for water. Interaction between the components perpendicular to the membrane plane is repulsive and decays as the inverse third power of distance. These interactions are screened by a dielectric permittivity of the order 10. Thus, despite the perpendicular components being much smaller in magnitude than the in-plane components, they will dominate the interaction energies at large distances.     

Cutoff size need not strongly influence molecular dynamics results for solvated polypeptides

The correct treatment of van der Waals and electrostatic nonbonded interactions in molecular force fields is essential for performing realistic molecular dynamics (MD) simulations of solvated polypeptides. The most computationally tractable treatment of nonbonded interactions in MD utilizes a spherical distance cutoff (typically, 8-12 A) to reduce the number of pairwise interactions. In this work, we assess three spherical atom-based cutoff approaches for use with all-atom explicit solvent MD: abrupt truncation, a CHARMM-style electrostatic shift truncation, and our own force-shifted truncation. The chosen system for this study is an end-capped 17-residue alanine-based alpha-helical peptide, selected because of its use in previous computational and experimental studies. We compare the time-averaged helical content calculated from these MD trajectories with experiment. We also examine the effect of varying the cutoff treatment and distance on energy conservation. We find that the abrupt truncation approach is pathological in its inability to conserve energy. The CHARMM-style shift truncation performs quite well but suffers from energetic instability. On the other hand, the force-shifted spherical cutoff method conserves energy, correctly predicts the experimental helical content, and shows convergence in simulation statistics as the cutoff is increased. This work demonstrates that by using proper and rigorous techniques, it is possible to correctly model polypeptide dynamics in solution with a spherical cutoff. The inherent computational advantage of spherical cutoffs over Ewald summation (and related) techniques is essential in accessing longer MD time scales.     

Molecular dynamics simulations of human prion protein: importance of correct treatment of electrostatic interactions.

Molecular dynamics simulations have been used to investigate the dynamical and structural behavior of a homology model of human prion protein HuPrP(90-230) generated from the NMR structure of the Syrian hamster prion protein ShPrP(90-231) and of ShPrP(<90-231) itself. These PrPs have a large number of charged residues on the protein surface. At the simulation pH 7, HuPrP(90-230) has a net charge of -1 eu from 15 positively and 14 negatively charged residues. Simulations for both PrPs, using the AMBER94 force field in a periodic box model with explicit water molecules, showed high sensitivity to the correct treatment of the electrostatic interactions. Highly unstable behavior of the structured region of the PrPs (127-230) was found using the truncation method, and stable trajectories could be achieved only by including all the long-range electrostatic interactions using the particle mesh Ewald (PME) method. The instability using the truncation method could not be reduced by adding sodium and chloride ions nor by replacing some of the sodium ions with calcium ions. The PME simulations showed, in accordance with NMR experiments with ShPrP and mouse PrP, a flexibly disordered N-terminal part, PrP(90-126), and a structured C-terminal part, PrP(127-230), which includes three alpha-helices and a short antiparallel beta-strand. The simulations showed some tendency for the highly conserved hydrophobic segment PrP(112-131) to adopt an alpha-helical conformation and for helix C to split at residues 212-213, a known disease-associated mutation site (Q212P). Three highly occupied salt bridges could be identified (E146/D144<-->R208, R164<-->D178, and R156<-->E196) which appear to be important for the stability of PrP by linking the stable main structured core (helices B and C) with the more flexible structured part (helix A and strands A and B). Two of these salt bridges involve disease-associated mutations (R208H and D178N). Decreased PrP stability shown by protein unfolding experiments on mutants of these residues and guanidinium chloride or temperature-induced unfolding studies indicating reduced stability at low pH are consistent with stabilization by salt bridges. The fact that electrostatic interactions, in general, and salt bridges, in particular, appear to play an important role in PrP stability has implications for PrP structure and stability at different pHs it may encounter physiologically during normal or abnormal recycling from the pH neutral membrane surface into endosomes or lysomes (acidic pHs) or in NMR experiments (5.2 for ShPrP and 4.5 for mouse PrP).