Letter pubs.acs.org/JPCL
Improved Parametrization of Li+, Na+, K+, and Mg2+ Ions for All-Atom Molecular Dynamics Simulations of Nucleic Acid Systems Jejoong Yoo† and Aleksei Aksimentiev*,†,‡ †
Department of Physics and ‡Beckman Institute for Advanced Science and Technology, University of Illinois at Urbana−Champaign, 1110 West Green Street, Urbana, Illinois 61801, United States S Supporting Information *
ABSTRACT: Atomic-scale modeling of compacted nucleic acids has the ability to reveal the inner workings of spectacular biomolecular machines, yet the outcome of such modeling efforts sensitively depends on the accuracy of the underlying computational models. Our molecular dynamics simulations of an array of 64 parallel duplex DNA revealed considerable artifacts of cation−DNA phosphate interactions in CHARMM and AMBER parameter sets: both the DNA arrangement and the pressure inside the DNA arrays were found to be in considerable disagreement with experiment. To improve the models, we fine-tuned van der Waals interaction parameters for specific ion pairs to reproduce experimental osmotic pressure of binary electrolyte solutions of biologically relevant ions. Repeating the DNA array simulations using our parameters produced results consistent with experiment. Our improved parametrization can be directly applied to molecular dynamics simulations of various charged biomolecular systems, including nucleic acids, proteins, and lipid bilayer membranes. SECTION: Macromolecules, Soft Matter
B
of nucleic acids.26 Several groups have reported large-scale simulations of highly compacted nucleic acid systems such as virus particles27 and ribosomes.28,29 Undesirably, our test simulations of dense DNA array systems revealed considerable artifacts of the current parametrization of nucleic acid systems. Our simulations were designed to mimic experimental measurements of the structure of DNA arrays subject to osmotic stress.13,30 Figure 1A illustrates a typical simulation system. The system contained 64 parallel double-stranded DNA molecules confined in an ideal cylindrical wall of a fixed radius and surrounded by an electrolyte solution that could freely permeate the wall of the cylinder; see the Supporting Information for simulation details. Two simulations were carried out using the standard CHARMM force field. In the first simulation (system I), the radius of the cylinder was 12 nm and the molar concentration of sodium chloride in the volume outside the confining cylinderthe buffer concentrationwas 250 mM. An 11 nm cylinder was used in the second simulations (system II), and the buffer concentrations of sodium and magnesium chloride were 250 and 25 mM, respectively. In the simulation of both systems, DNA molecules were observed to form close-packing clusters (Figure 1B and Figure S5A of the Supporting Information). The presence of such clusters is in disagreement with the results of the X-ray diffraction measurements, which, under similar conditions, found DNA molecules arranged on a
ecause it carries the genetic blueprint of an entire living organism, DNA is one of the most important biological molecules. Unsurprisingly, the questions of how DNA is packaged, replicated, repaired, and transcribed in a living cell and how the above processes are regulated by various biomolecular machines are at the core of modern life science.1 Computer modeling of nucleic acids can offer unique insight into their biological function.2−4 In particular, all-atom molecular dynamics (MD) simulations5,6 have provided unique atomic-level views on the structural dynamics of nucleic acids,7 forces governing DNA organizations,8,9 and the molecular mechanisms of DNA processing machines.10,11 Being a highly charged polyanion, DNA is constantly surrounded by a cloud of counter- and co-ions, whose primary physical purpose is to neutralize the DNA charge. Over the past few decades, it has become apparent that the ions surrounding DNA are not just passive bystanders.12 On the contrary, the ions perform active regulatory roles in the processes of nucleic acid packaging,13,14 folding,15,16 and ribozyme activity17 and are essential mediators of DNA−protein interactions.18 The most commonly used atomic-scale models of DNA are based on the CHARMM19 and AMBER20 force fields that express the potential energy landscape of individual atoms through a small set of analytical functions parametrized to reproduce select properties of model compounds. Both force fields are equipped with complementary ion parameters.21−23 Despite providing only an approximate description of nucleic acid systems, such models have been successfully used to address a variety of questions regarding the mechanical properties of DNA,24,25 DNA−DNA interactions,8 and folding © 2011 American Chemical Society
Received: November 14, 2011 Accepted: December 6, 2011 Published: December 6, 2011 45
dx.doi.org/10.1021/jz201501a | J. Phys. Chem.Lett. 2012, 3, 45−50
The Journal of Physical Chemistry Letters
Letter
electrolyte concentration by tempering formation of ion pairs and clusters. For a pair of ions, the LJ σ parameter approximately equals the sum of the effective radii of the two atoms. Following the language used by Luo and Roux, we refer to ion-pair specific corrections to LJ parameter σ as NBFIX (non-bonded interaction FIX) corrections. In this study, we use the method of Luo and Roux32 to reparametrize interactions between cations and biologically relevant anions such as acetate and phosphate. Specifically, we consider combinations of Li+, Na+, K+, or Mg2+ cations with Cl−, acetate (Ac−), or dimethylphosphate (DMP−) anions. Here DMP− and Ac−, whose chemical structures are shown in Figure 2A, serve as model compounds for DNA phosphate and
Figure 1. Parametrization artifacts in the simulations of dense DNA and DNA mimetic systems. (A) Setup of a DNA array simulation. 64 DNA molecules (shown in black) are subject to a half-harmonic potential (sketched in orange) that confines the molecules to a 12 nm radius cylinder. The potential does not apply to water and ions, which are shown as a blue semitransparent surface and yellow (Na+) and green (Cl−) spheres, respectively. (B) Representative snapshot of the DNA array system obtained from an MD simulation performed using the CHARMM parameter set.19 The image shows an xy cross section of system. The location of each DNA molecule is depicted using a white 2 nm diameter circle. The Na+ density (averaged over the z coordinate) is shown. The concentration of NaCl outside the DNA array (not visible within the field of view of the image) is 250 mM. (C) Large clusters are formed in an MD simulation of ∼3 M bulk solution of sodium acetate. Na+ ions are shown as yellow spheres, Ac− ions are shown using the molecular bonds representation, and water is not shown.
hexagonal lattice and the distance between two neighboring DNA molecules to be considerably greater than the distance of direct contact.13 The 2D ion density plot in Figure 1B suggests a possible cause of the DNA clustering: the local concentration of sodium ions exceeds 3 M in the region where DNA molecules make direct contact. Furthermore, in our test simulations of a ∼3 M sodium acetate (NaAc) solution using either CHARMM or AMBER force fields, we observed formation of a large NaAc cluster, Figure 1C. As the experimental solubility limit of a NaAc solution is ∼5 M,31 we suspect that imperfect parametrization of sodium-phosphate oxygen and sodiumacetate oxygen interactions in DNA and NaAc systems, respectively, is responsible for the clustering behavior. Traditionally, force field parameters for ions have been derived to match experimental solvation energies.21,22 Several recent studies reevaluated the ability of standard parametrization to describe accurately ion−ion interactions.23,32−37 In particular, Luo and Roux have shown that clustering of NaCl and KCl solutions at ion concentrations exceeding 3 M results in a significant drop of the solutions’ osmotic pressure.32 The same study has also shown a method to remedy the problem: increasing Lennard-Jones (LJ) σ parameters for Na−Cl and K− Cl pairs by only 0.038 and 0.043 Å, respectively, recovered the experimental dependence of the osmotic pressure on the
Figure 2. MD simulations of osmotic pressure coefficients. (A) Chemical structures of acetate and dimethylphosphate. Oxygen atoms subject to NBFIX corrections are highlighted in blue. (B) Typical simulation system that contains an electrolyte solution separated from pure water by two planar half-harmonic potentials. Dashed lines illustrate the location of half-harmonic potentials that prevent salt ions from leaving the central compartment. The central compartment of the system shown contains ∼2.8 m NaAc in equilibrium. (C) Timeaveraged mass density profiles of the system shown in panel B.
amino acids aspartate and glutamate, respectively. For 11 electrolyte solutions (excluding MgDMP2 solution due to missing experimental data), we carry out osmotic pressure calculations using three sets of parameters based on the CHARMM and AMBER force fields (see the Supporting Information for details), aiming to derive a set of ion pairspecific LJ σ values optimized for simulations of highly charged biomolecular systems. Figure 2B illustrates a typical simulation system that we used to measure osmotic pressure πsim of an electrolyte solution.32 Each simulation system contained two compartments separated 46
dx.doi.org/10.1021/jz201501a | J. Phys. Chem.Lett. 2012, 3, 45−50
The Journal of Physical Chemistry Letters
Letter
by two virtual semipermeable membranes aligned with the xy plane. One compartment contained an electrolyte solution, whereas the other contained pure water. Whereas water could freely pass through the membranes separating the two compartments, the ions were subject to the force of halfharmonic planar potentials, whose action mimicked the ideal behaviors of a semipermeable membrane. The presence of the semipermeable membranes created a density difference between the electrolyte and pure water compartments (Figure 2C) maintained by the excess pressure applied by the membranes. This excess pressure equals the osmotic pressure of the electrolyte and can be measured by averaging the forces exerted by ions on the confining semipermeable membranes.32 To make direct comparison with experiment possible, we compute the molal concentrations of electrolytes, m, using the mass densities of the salt ions and pure water in the electrolyte compartment. Then, we calculate the osmotic coefficient, ϕ, defined as ϕ = πsim/πid(m), where πid(m) is the osmotic pressure of an ideal electrolyte solution at molal concentration, m. A complete description of the methods and protocols is provided in the Supporting Information. Figure 3 plots osmotic coefficients simulated using unmodified CHARMM parameters (blue lines) along with
tration. Inaccurate parametrization of cation−anion interactions is the source of these discrepancies. In general, the standard parameters were observed to overestimate cation−anion attraction, facilitating excessive ion pair formation (see Figure 1B,C for examples), which lowers the osmotic pressure. Compared with other electrolyte solutions considered in the present study, osmotic coefficients for Li/Na−Ac solutions deviate the most from the corresponding experimental values. For example, for both CHARMM and AMBER parameter sets, computed osmotic coefficients of Li/Na−Ac solutions at 3 m are below 0.2, whereas the experimental values are ∼1.2. To correct the imperfections of standard parametrization, we systematically varied the LJ σ parameters of the specific cation− anion pairs until the simulated osmotic coefficients matched the experimental values. Figure 4 provides an overview of our
Figure 4. Illustration of our reparametrization scheme. See the text for details. Oxygen is colored in red, magnesium in pink, carbon in cyan, phosphorus in tan, hydrogen in white, and salt ions in yellow.
approach. For monovalent cations, we chose to adjust the LJ σ parameter describing vdW interaction of the cation−Cl−, cation−acetate oxygen−, or cation−DMP terminal oxygen− pairs because these atoms were observed to come in direct contact with each other to form ion pairs. In the case of Mg2+, our reparametrization scheme focused on the six water molecules from the first solvation shell of the ions. Even though a Mg2+ ion and the six surrounding water molecules are not covalently bonded, the Mg−hexahydrate complex, Mg(H2O)62+, is expected to remain intact as long as the simulation time is considerably shorter than the typical exchange time for the first solvation shell water, which was estimated to be ∼10 μs.38 For Mg(H2O)62+, we adjust LJ σ between oxygens of a Mg(H2O)62+ complex and either Cl−, oxygens of Ac−, or terminal oxygens of DMP−. Following the results obtained using the AMOEBA polarizable force field,39 we increased the dipole moments of the Mg(H2O)62+ water by 1 debye. Without using the above adjustment to the water dipole moment, matching the simulated osmotic pressure to the experimental value was not possible because the simulations considerably underestimated the osmotic pressure for any reasonable value of LJ σ. For cation-DMP− solutions, experimental osmotic pressure is known only at low salt concentration (3
Figure 3. Simulated and experimentally measured osmotic coefficients of NaAc (A,B), NaDMP (C), and MgAc2 (D) solutions. In all panels, the black lines depict the experimental dependence, blue lines show the dependence simulated using the standard CHARMM parameter set (Δσ = 0), and the red lines show the data obtained using our optimized values of LJ σ. In panel A, black diamonds show the osmotic coefficient and the molal concentration of the same system for several values of Δσ. Panel B shows a close-up view of the dependence shown in panel A.
the experimental values (black lines). Disagreement between the simulated and experimental data is apparent for all solutions tested. In general, simulated dependences deviate very little from experiment at low ion concentrations (
