Improved model of hydrated calcium ion for molecular dynamics simulations using classical biomolecular force fields

Abstract

Calcium ions (Ca(2+) ) play key roles in various fundamental biological processes such as cell signaling and brain function. Molecular dynamics (MD) simulations have been used to study such interactions, however, the accuracy of the Ca(2+) models provided by the standard MD force fields has not been rigorously tested. Here, we assess the performance of the Ca(2+) models from the most popular classical force fields AMBER and CHARMM by computing the osmotic pressure of model compounds and the free energy of DNA-DNA interactions. In the simulations performed using the two standard models, Ca(2+) ions are seen to form artificial clusters with chloride, acetate, and phosphate species; the osmotic pressure of CaAc2 and CaCl2 solutions is a small fraction of the experimental values for both force fields. Using the standard parameterization of Ca(2+) ions in the simulations of Ca(2+) -mediated DNA-DNA interactions leads to qualitatively wrong outcomes: both AMBER and CHARMM simulations suggest strong inter-DNA attraction whereas, in experiment, DNA molecules repel one another. The artificial attraction of Ca(2+) to DNA phosphate is strong enough to affect the direction of the electric field-driven translocation of DNA through a solid-state nanopore. To address these shortcomings of the standard Ca(2+) model, we introduce a custom model of a hydrated Ca(2+) ion and show that using our model brings the results of the above MD simulations in quantitative agreement with experiment. Our improved model of Ca(2+) can be readily applied to MD simulations of various biomolecular systems, including nucleic acids, proteins and lipid bilayer membranes. © 2016 Wiley Periodicals, Inc. Biopolymers 105: 752-763, 2016.

Keywords: AMBER; CHARMM; calcium; force field; molecular dynamics; nucleic acid.

Figure 1

Standard molecular force fields overestimate association of Ca²⁺ with anionic species. (A,B) Representative configuration of solutes in MD simulations of 2 m CaAc₂ (A) and 3 m CaCl₂ (B) solutions performed using the standard CHARMM force field. Only a 1 nm wide slice through the system is shown for clarity. In each panel, the blue semi-transparent surface illustrates the dimensions of a unit simulation cell; dashed lines indicate the presence of ideal semi-permeable membranes that confine solutes to one of the two compartments; the other compartment contains pure water. (C,D) Comparison of experimental (black square or black line) and simulated (red or blue) osmotic pressure values for the CaAc₂ (C) and CaCl₂ (D) solutions as a function of the solutes’ concentration. The data were obtained using the standard AMBER (blue) and CHARMM (red) parameter sets. The osmotic pressure of an ideal 2:1 solution (osmotic coefficient = 1) is shown using dotted lines. The standard error of 1 ns block averages of 2 ps sampled osmotic pressure data are smaller than the symbols.

Figure 2

MD simulation of Ca²⁺-mediated DNA–DNA interaction free energy using the standard CHARMM and AMBER force fields. (A) Representative configuration of DNA and ions in umbrella sampling simulations of the DNA–DNA potential of mean force. The DNA molecules are shown in gray. A blue semitransparent surface depicts the volume occupied by the solvent in a unit simulation cell. Colored spheres indicate the locations of calcium (blue), phosphate oxygen (red) and chloride (green). This particular image shows the DNA system at an inter-DNA distance of 24 Å; the bulk concentration of CaCl₂ is 25 mM. The inset shows a pair of phosphate groups stabilized by multiple Ca²⁺ ions directly bound to the phosphates forming an ion bridge between the two DNA molecules. (B) The interaction free energy of two parallel dsDNA molecules versus the DNA–DNA distance at 25 mM CaCl₂ computed using the standard parameterization of the CHARMM (red) and AMBER (blue) force fields. The experimental dependence is shown as a black line.

Figure 3

Reparameterization of a hydrated calcium ion using a heptahydrate model. In each heptahydrate complex, Ca(H₂O)₇, seven water oxygens are harmonically restrained to a calcium atom. To account for Ca²⁺-induced polarization of water molecules,^, the dipole moment of each of the seven water molecules in the heptahydrate complex is increased by 0.5 debye by adjusting the partial charges of water oxygen and water hydrogen atoms. Following that, the LJ σ parameters describing specific interactions between water oxygens of Ca²⁺ (H₂O)₇ and the target anionic species (phosphate oxygen, acetate oxygen, and Cl⁻) are refined to match the experimental osmotic pressure data.

Figure 4

MD simulations of CaAc₂ and CaCl₂ solutions using the refined force field models. (A,B) Representative configuration of solutes in MD simulations of 2 m CaAc₂ (A) and 3 m CaCl₂ (B) solutions performed using our improved parameterization of Ca²⁺ ions. In each panel, the blue semi-transparent surface illustrates the dimensions of a unit simulation cell; dashed lines indicate the presence of ideal semi-permeable membranes that confine the solutes to one of the two compartments; the other compartment contains pure water. Only a 1 nm-wide slice of the system is shown for clarity. (C,D) Comparison of experimental (black square or black line) and simulated (red or blue) osmotic pressure values for the CaAc₂ (C) and CaCl₂ (D) solutions as a function of the solute concentration. The data were obtained using our calcium heptahydrate parameterization of Ca²⁺ ions in AMBER (blue) and CHARMM (red) force fields. The osmotic pressure of an ideal 2:1 solution (osmotic coefficient = 1) is shown using dotted lines. The standard error of 1 ns block averages of 2 ps sampled osmotic pressure data are smaller than the symbols.

Figure 5

Simulations of DNA–DNA potential of mean force using our refinement of the CHARMM and AMBER force fields. (A,B) The interaction free energy of two parallel dsDNA molecules versus the DNA–DNA distance at 25 mM CaCl₂ computed using our refined parameterization of the CHARMM (A) and AMBER (B) force fields. The experimental dependence is shown as a black line. The red lines show the results obtained using the same adjustments of the σ parameter for the Ca²⁺ (H₂O)₇ oxygen-DNA phosphate oxygen pair as for the Ca²⁺ (H₂O)₇ oxygen-Ac⁻ oxygen pair (Δσ = 0.07 Å for both force fields). Blue lines show the results obtained using the σ parameters adjusted to reproduce experimental DNA–DNA interaction free energy (Δσ = 0.02 and 0.04 Å for CHARMM and AMBER, respectively). (C) Representative configuration of DNA and ions in umbrella sampling simulations of the DNA–DNA PMF carried out using our refined version of the CHARMM force field. The DNA molecules are shown in gray, blue semi-transparent surface depicts the volume occupied by the solvent in a unit simulation cell. Colored spheres indicate the locations of calcium (blue), phosphate oxygen (red) and chloride (green). This particular image shows the DNA system at the inter-DNA distance of 24 Å; the bulk concentration of CaCl₂ is 25 mM.

Figure 6

Simulations of electric field-driven nanopore transport of DNA in the presence of Ca²⁺ ions. (A) The microscopic configuration of the nanopore system at the beginning of the DNA translocation simulations. DNA is shown using tan spheres; calcium ions are shown as blue spheres. The volume occupied by water is indicated by a semitransparent surface. The Si₃N₄ membrane is shown cut away, with surface atoms shown as black spheres and the interior atoms shown as grey spheres. In addition to calcium, the solution contains 0.4 M of KCl (not shown). The arrow illustrates the direction of applied electric field producing a positive transmembrane bias. (B–D) The z coordinate of the DNA’s center of mass (CoM) during MD simulations at a transmembrane bias of 1 V (B), 200 mV (C), and 100 mV (D). Data from simulations performed using the standard parameterization of the CHARMM force field are shown in blue; orange lines illustrate the results of simulations performed using our parameterization of Ca²⁺ ions for the CHARMM force field. The z coordinate is defined in panel A. (E–G) The number of calcium ions directly bound to DNA (black) and the number of chloride ions directly bound to calcium ions that are already bound to DNA (blue) during the simulations of DNA translocation performed at a transmembrane bias of 1 V (E), 200 mV (F), and 100 mV (G). These data characterize the simulations carried out using the standard CHARMM force field. No Cl⁻-Ca²⁺-DNA binding (defined as a Cl⁻ ion remaining within 0.4 nm of a Ca²⁺ ion for more than 30 ps) was observed in MD simulations employing our parameterization of Ca²⁺ ions.