Multisite ion model in concentrated solutions of divalent cations (MgCl2 and CaCl2): osmotic pressure calculations.

Accurate force field parameters for ions are essential for meaningful simulation studies of proteins and nucleic acids. Currently accepted models of ions, especially for divalent ions, do not necessarily reproduce the right physiological behavior of Ca(2+) and Mg(2+) ions. Saxena and Sept (J. Chem. Theor. Comput. 2013, 9, 3538-3542) described a model, called the multisite-ion model, where instead of treating the ions as an isolated sphere, the charge was split into multiple sites with partial charge. This model provided accurate inner shell coordination of the ion with biomolecules and predicted better free energies for proteins and nucleic acids. Here, we expand and refine the multisite model to describe the behavior of divalent ions in concentrated MgCl2 and CaCl2 electrolyte solutions, eliminating the unusual ion-ion pairing and clustering of ions which occurred in the original model. We calibrate and improve the parameters of the multisite model by matching the osmotic pressure of concentrated solutions of MgCl2 to the experimental values and then use these parameters to test the behavior of CaCl2 solutions. We find that the concentrated solutions of both divalent ions exhibit the experimentally observed behavior with correct osmotic pressure, the presence of solvent separated ion pairs instead of direct ion pairs, and no aggregation of ions. The improved multisite model for (Mg(2+) and Ca(2+)) can be used in classical simulations of biomolecules at physiologically relevant salt concentrations.

Introduction

Metal ions have many important roles in biology and environmental chemistry. Ions such as calcium, magnesium, sodium, potassium, and chloride take part in vital bodily functions such as the generation of action potentials, catalytic actions of enzymes, signal transduction, and structural stabilization of nucleic acids.[1,2] Furthermore, a large number of drugs contain metal ions, which allow those drugs to interact with endogenous proteins.[3] Metals interact with the biological molecule, either directly or through the water present in their first solvation shell.[4−6] When present in solution, metals either form monovalent or divalent ions. Divalent ions, especially calcium (Ca2+) and magnesium (Mg2+), play an important role in physiology as both act as cofactors in many enzymatic reactions. Mg2+ drives the folding of RNA into stable states[2,7,8] and also participates in energy metabolism by binding to adenosine triphosphate (ATP).[9] Ca2+ regulates neurotransmitter release, muscle contractions, and blood pressure.[10−12]

Most researchers study physiological activity of ion-related mechanisms by estimating the thermodynamic quantities of ions.[13] They find their estimates by performing experiments on mixtures of cations and anions in aqueous solution. Such experiments give thermodynamic quantities of the overall system. To tease out the properties of an individual ion, researchers have to make several assumptions, leaving the values of an individual ion to reside within an approximate range.[14−16] Some researchers have used molecular simulations as an alternate method. Through simulations, they have been able to theoretically calculate precise quantities for each ion species.[5,17−19] However, they observed that accurate calculations of thermodynamic quantities require careful selection of parameters so that their theoretical values would match with the experimentally derived values. In simulations, ions are represented using three basic parameters: a charge and two Lennard-Jones (LJ) potential parameters. The charge represents the valency of the ion, and the LJ parameters determine the size. The experimental measure most widely used for validation is the hydration free energy or the negative work required to insert a single ion in bulk water.[5,17−19] Despite careful selection of parameters and perfect match with these experimental quantities, researchers have faced challenges in reproducing accurate physiological behavior of biomolecules interacting with divalent ions. Historically, simulations of divalent ions with biomolecules have suffered from two main problems. First, an accurate structure of the coordination of ions with the biomolecules cannot be obtained, and second, the introduction of higher concentration of ions in a solution leads to excessive aggregations and formation of ion clusters.

The classical model of ions contains the entire charge of the ion at the center. This property prevents the ion from responding directly to the molecular environment and thus limits their usage with bigger and complex biomolecules in classical simulations. Researchers observed that when they fit the parameters after precise matching with the experimental hydration free energy, they were able to describe the structure of the first hydration shell of ions with considerable accuracy. For Mg2+ and Ca2+, the ion-oxygen distances in water have been reported as 2.1–2.5 Å and 2.4–3.0 Å,[20−31] which lie close to the experimental values of 2.1 and 2.4 Å,[32,33] respectively. The behavior of the metal cation when it interacts with a biomolecule, however, involves substantial charge-transfer effects, which are hard to capture with the localized, fixed, charged model mentioned above. An external constraint placed between a Ca2+ ion and oxygen atoms in a theoretical study of Ca2+ binding with Parvalbumin provided evidence of how the charge in the center of the ion by itself is not sufficient to keep the Ca2+ ion stable in its binding site.[34] Another theoretical study of Ca2+ binding in Calbindin reported that a specialized treatment of electrostatics was necessary to keep the Ca2+ ions stably bound to the two calcium binding sites of the protein.[4] The accurate structure of direct coordination of a cation with biomolecules becomes much more important in studies involving the description of conformational changes occurring after ion binding. Mg2+ binding to unstructured RNA molecules is known to shift the folding versus unfolding equilibrium toward folded states.[2] Similarly, the activation mechanism of the Calmodulin family of proteins also involve large conformational changes occurring after Ca2+ binding.[37−39] Theoretical studies of this phenomenon have been extremely challenging with the ion models in which the charge is located at the center of the ion.[35,36]

Various approximations for the interaction of ions with biomolecules also result in high-energy barriers for change of coordinating partners, a process that occurs frequently in biological systems. For example, during passage of ions through ion channels,[40] the ions either partially or completely dehydrate, changing their coordination partners from water oxygens to the protein oxygens. Also, in the presence of Ca2+ ions, the EF-hand family of proteins replaces the already bound Mg2+ ion with the more preferred Ca2+ ion. Researchers have proposed several solutions to treat this problem, including reparameterizing the charges of the coordinating atoms[41] or using polarizable parameters for all atoms.[27,42] These solutions are either too specific for the system, rendering the parameters nontransferable, or require the use of extensive computational resources. In a previous work, Saxena and Sept[43] offered a workable solution with the development of a new model for the ions. A multisite model was developed for divalent ions Mg2+ and Ca2+, where the localized charge of the ion was split into multiple sites. This model was able to accurately represent the structure of the first coordination shell of the ion in water, proteins, and nucleic acids. Saxena and Sept demonstrated that the model could be successfully used to find conformational changes and thermodynamical changes associated with ion binding in biomolecules, such as EF-hand proteins and RNA.

The second major problem faced by researchers in theoretical studies of ion-mediated mechanisms is ion aggregation. Ions at low concentrations behave as ideal solutions, but as the ion concentration increases, the solutions move away from ideality.[44] The electrostatic forces between the ions cause association and/or dissociation of solute particles. The study of biological functions of proteins requires the experiments to be performed at physiological ion concentrations. At such concentrations, water molecules occupy the first coordination shell of the ions. Thus, at these concentrations, the cation and anion interact with each other through a water molecule forming a solvent separated pair (SSP).[32,45,46] Simulation studies of ion solutions at finite concentrations have found untimely aggregation of ion pairs.[47−49] In-depth analysis of the ion clusters showed the that at conditions where the oppositely charged ions should have formed an SSP, they instead formed direct contact pairs.[50−54] This behavior is due to the inaccurate representation of ion–cation interactions.[47,48] The aforementioned artifact thus leads to excessive ion pairing, poor solubility, and spontaneous crystallization of ion solutions in theoretical calculations.

The multisite model developed by Saxena and Sept also suffers with the excessive aggregation problem in concentrated MgCl2 and CaCl2 solutions. In this work, we refine the parameters of the multisite model such that it is devoid of both artifacts. The multisite model is already able to predict accurate inner shell coordination with biomolecules. To eliminate the ion aggregation problem,we modify the interactions of the divalent cations with the monovalent anion to reproduce the physical behavior of concentrated ion solutions, while maintaining the parameters for the interaction of the ions with water and biomolecules. Experimentally, the physical behavior of regular ion solutions can be studied by measuring properties such as osmotic pressure and freezing point,[44] but the calculation of these properties through simulations has been difficult. Recently, a methodology was developed to calculate osmotic pressure through simulations.[48] Osmotic pressure calculations offer a method to measure the strength of solute interactions. At a similar concentration of solutes, aggregating molecules apply less pressure on the walls and results in reduced osmotic pressure, while completely dissociated solute molecules move around in water and exert more pressure on the walls, resulting in increased osmotic pressure. By matching the calculated osmotic pressure to the experimental values, the solute interaction can be improved to accurately reproduce the behavior of concentrated electrolytes. This method has been successfully used to obtain accurate behavior of aqueous solutions of NaCl and KCl[48] and to obtain accurate behavior of trimethylamine N-oxide.[55] We use this method here to calculate the osmotic pressure of MgCl2 solutions at various concentrations and match them to the experimental values to refine the parameters of the multisite model. The cation–anion radial distribution functions obtained after the parameter refinement show that the solutions do not form excess direct ion-pairs, matching the experimentally observed behavior. Subsequently, we use the refined parameters of the multisite model to study the behavior of CaCl2 solutions and find that the calculations reproduce the experimental osmotic pressure of CaCl2 with high accuracy, without further refinement. Structural characterizations confirm that the divalent ions show accurate inner shell coordination and are devoid of unusual ion–ion pairing at high salt concentrations in both solution types (MgCl2 or CaCl2). The model captures the exchange of inner shell waters and better coordination geometries with its binding partners. Further, in comparison with the work by Saxena and Sept,[43] we show that the model retains the ion solution properties seen in that study, indicating that the accuracy of the model is not compromised after the changes made in this work. Through this study, we present an improved multisite model, which can not only be used for prediction of free energies and selectivity of divalent cation (Mg2+ and Ca2+) binding sites in proteins and nucleic acids, but can now be used over a wide range of physiological salt concentrations. Additionally, the model can be easily used with standard fixed charged force fields with no extra computational cost.

Theory and Methods

Multisite Model

The multisite model consists of a central atom M and n dummy charge centers, d, where n denotes the coordination number of the ion. The dummy atoms are located along the vectors connecting the ion center and the coordinating atoms, determined from the coordination structure of the ion, octahedral for Mg2+ and pentagonal bipyramidal for Ca2+. The equilibrium bond length between the ion center and the charged sites was fixed at 0.9 Å. Each dummy atom in a given ion has the same atom type and carries identical charge and Lennard-Jones (LJ) parameters. The entire charge of the ion is distributed to the dummy atoms. While the LJ parameters of the dummy atoms remain the same among different ions, the LJ parameters of the central atom differ and are based on their respective match to the hydration free energy. The starting parameter values of the model are shown in Table 1.[43]

Table 1

Lennard-Jones and Charge Parameters for the Multisite Models of Ca2+ and Mg2+ Ionsa

ion	qM	qd	C12ij,M	C6ij,M	C12ij,d	C6ij,d
Ca2+	0	0.29	2.275 ×10–7	0.005	1.046 × 10–14	0.0
Mg2+	0	0.33	3.050 ×10–9	0.001	1.046 × 10–14	0.0

Units: C12 in kJ/mol·nm12 and C6 in kJ/mol·nm6.

Osmotic Pressure Theory

The behavior of ions in concentrated electrolytes causes the solutions to deviate from the van’t Hoff ideal solution.[44] Osmotic pressure has been widely used to estimate thermodynamics quantities, as it is relatively simple to measure.[13,44] Osmosis is the process where the solvent spontaneously passes through a semipermeable membrane separating the solute–solvent solution and a pure solvent, until equilibrium is reached. The pressure applied by the solute on the membrane is related to the concentration of the solute by Morse equation:where i is the dimensionless van’t Hoff factor, M is the molarity, R is the gas constant, and T is the thermodynamic (absolute) temperature. The calculation of of the osmotic pressure simulations has been challenging. Recently, Luo and Roux developed a method for calculating osmotic pressure in simulations.[48] The methodology involves creating virtual walls in a simulation box, which act as semipermeable membranes allowing only water to pass through. The average force required to keep the solute molecules inside the simulation box is used to calculate the osmotic pressure. The osmotic pressure also offers a measure of the strength of interactions between solutes. At a similar concentration of solutes, aggregating molecules apply less pressure on the walls and result in reduced osmotic pressure, while completely dissociated solute molecules move around in water and exert more pressure on the walls, resulting in increased osmotic pressure. Due to these characteristics, the osmotic pressure calculations provide a means to improve the parameters of the solutes in simulations by matching with the experimental values. Here, we follow the same methodology as Luo and Roux to determine accurate interaction between Mg2+ and Cl– ions and Ca2+ and Cl– ions. We systematically modify the van der Waals parameters of the divalent ions such that the experimental osmotic pressure values of these solutions are matched. In particular, we modify the interaction of the cation dummy sites with Cl– ions. All other parameters are kept the same as in Saxena and Sept, such that the interactions of the divalent cations with water and biomolecules do not change.

Simulation Protocol

We calculate the osmotic pressure of MgCl2 and CaCl2 solutions at various molalities using molecular dynamics simulations at constant temperature and pressure. The simulation system consists of a central box of the ionic solution flanked by a pure water box of equal volume. The boundary of the central box determines the location of the virtual walls. The two boxes were generated individually at the required size, equilibrated using Berendsen pressure coupling and then concatenated along the z-direction. Periodic boundary conditions ensured that another pure water box flanked the center box from the other side. The resulting box was energy minimized and equilibrated for 1 ns to allow equilibration of water and the proper mixing of the ions. The final state obtained from the previous step is then used to start three independent simulations, each 10 ns long. Simulations were conducted at constant pressure and constant cross sectional area for each desired molality. This process was followed individually for each ion concentration of the central box. The composition of the simulated systems are described in Table 2. The wall restraints were applied along the virtual walls by modifying Gromacs 4.0.7 to enforce a flat-bottom, half-harmonic potential. The osmotic pressure was calculated from the force exerted by the wall to keep the solute inside the simulation box, π = (⟨Fwall⟩/A), where A is the cross-sectional area of the virtual wall obtained from the box length of the central box for each ion concentration (Table 2). ⟨Fwall⟩ = (k/Np)∑∑(|z| – |zwall|), with |z| > |zwall|, k as the force constant (10 kcal/mol/2), Np as the total number of configurations in the production stage of each simulation (5 ns), and i as the index of ions. The force is averaged over two half-harmonic walls to get the osmotic pressure. Errors in the averages are estimated using 2.5 ns block averages over the production stage of the three independent simulations. Osmotic pressure for each concentration was compared with the experimentally obtained values.[56,57] The multisite model was used for Mg2+ and Ca2+ parameters,[43] TIP3P for water molecules,[58] and Cl– parameters were obtained from Luo and Roux’s work.[48] The parameters for the Cl– anions with the cation dummy sites were adjusted iteratively to match experimental osmotic pressure.

Table 2

System Details for MgCl2 Osmotic Pressure Calculations

molality	NMg or NCl	NWat	box length (nm)	molarity
0.08	16/32	11341	6.95	0.08
0.43	16/32	2096	4.00	0.42
1.17	16/32	810	2.95	1.03
1.59	16/32	605	2.75	1.28
1.95	16/32	503	2.60	1.51

Results and Discussions

MgCl2 Interaction before Fitting

We calculate the osmotic pressure of Mg2+ in solution with Cl– ions at different concentrations. When comparing with the experimental values, we find that at lower concentrations (<0.5 M) the calculated values of osmotic pressure match well with the experiments, but at higher concentration the calculated values drop significantly and the deviation grows larger as the concentration increased (Figure 1a). The decrease in osmotic pressure indicates that the solute (ions) might have a strong attraction between them. To verify this, we examined the interactions between the ions. The MgCl2 mixtures, at almost all concentrations, show spontaneous aggregation of ions, which occurs early in the simulations and remained unchanged for the entire length of the simulation (100 ns, Figure 1b). We found several instances where Cl– ions occupied the first coordination shell of Mg2+ ions and vice versa. The radial distribution function (RDF) between Mg2+ and Cl– ions is proportional to the probability of finding Cl– ions within a certain radius around Mg2+ ions. The RDFs show that there is a large population of Cl– ions located at a distance of 1.95 Å from the Mg2+ ions (Figure 1c), indicating that the two ions form direct contact pairs. This behavior does not match X-ray experiments of MgCl2 solutions, where the minimum distance between Mg2+ and Cl– ions, below 3 M concentration, is close to 4.5 Å.[45] Structurally, this value corresponds to a solvent-separated pair, which we do not see in our simulations. Further, the X-ray and neutron scattering experiments also reported that the total number of Cl– ions coordinating with Mg2+ ions does not exceed 2 below 3 M concentrations.[32,33,45,46,59] From the simulation results, however, we find the number sof Cl– ions coordinating with the Mg2+ ions are 7, 5, 3.2, and 3 for the four different concentrations studied (Figure 1d), all being larger than 2. Thus, from these observations it becomes clear that the strength of interactions between ions is overestimated in the simulations, leading to aggregation.

Figure 1

(A) Osmotic pressure obtained for different concentrations of MgCl2 in water. The black line is the osmotic pressure data obtained from experiments,[56] and the magenta line is the osmotic pressure obtained from simulations using the multisite Mg2+ model. (B) Snapshot of a simulation configuration of Mg2+ and Cl– ions in water showing aggregation (Mg2+ as red spheres, Cl– as green spheres, water molecules occupying the first solvation shell of Mg2+ shown in Licorice). Direct interaction of cations and anions can be seen in the dotted oval. (C) Radial distribution function between Mg2+ and Cl– ions at different concentrations of MgCl2 (m represents molality = mol/kg). (D) Total number of Cl– ions binding to Mg2+ ions at different concentrations of MgCl2.

We observe similar behavior in the simulations of CaCl2 solutions (Figure 2). Osmotic pressure values start to deviate for concentrations >0.5 M (Figure 2a). We find a direct correspondence between the occurrence of ion clusters and the increase in deviation of the calculated osmotic pressure from the experimental values. Visual inspection shows aggregate formation for all concentrations (Figure 2c). Through RDF calculations we find that the probability of finding Cl– around Ca2+ ions within a minimum radius of 1.98 Å is very high, indicating the presence of direct ion pairs (Figure 2b). As seen in MgCl2 simulations, the CaCl2 simulations show coordination numbers of Cl– with Ca2+ more than 2 for all concentrations (Figure 2d).

Figure 2

(A) Osmotic pressure obtained for different concentrations of CaCl2 in water. The black line is the osmotic pressure data obtained from experiments,[57] and the magenta line is the osmotic pressure obtained from simulations after using the multisite Ca2+ model. (B) Snapshot of a simulation of Ca2+ and Cl– ions in water showing aggregation of the ions (Ca2+ as magenta spheres, Cl– as green spheres, water molecules are not shown for clarity). Direct interaction of anions and cations can be observed. (C) Radial distribution function between Ca2+ and Cl– ions at different concentrations of CaCl2 (m represents molality = mol/kg). (D) Total number of Cl– ions binding to Ca2+ ions at different concentrations of CaCl2.

The results described above clearly indicate that the reduction in osmotic pressure at higher concentration occurred due to clustering and aggregation of ions. Aggregation of ions leads to fewer collisions of the solute particles with the membrane and thus leads to a reduction in the osmotic pressure. Based on the experimental results, however, such aggregation only starts to occur at solute concentrations nearing the crystallization conditions.[45,46] Simulation studies by Callahan and co-workers have also shown that the existence of direct ion pairing in MgCl2 solutions is not energetically favorable.[26] This suggests that the ion aggregation that we observe in solute concentrations ranging below 3 M is certainly an artifact of the model. Such an artifact can be attributed to the discrepancies in the interaction parameters of the atoms. As mentioned before, the parameters for the multisite model were developed and optimized to reproduce accurate coordination properties of a single ion in water and are insufficient to correct for the excessive cation–anion aggregation at physiological conditions. In the next section we focus on refining the parameters of the multisite model such that an accurate physical behavior of concentrated ionic solutions can be reproduced through simulations.

Refinement of the Multisite Ion Parameters

To reduce ion aggregation, we focus on changing the cation–anion interaction parameters, which affect the electrostatic and the van der Waals interactions. We keep the same charge distribution and LJ parameters for the central atoms, as these choices were found to be critical for calculation of sensitive thermodynamic parameters in our previous work.[43] In the calculations described above, the LJ parameters between the cation and anions were determined by the Lorentz–Berthelot mixing rule. Here, following the work of Luo and Roux,[48] we move away from the conventional mixing rules and allowed the two unique atoms of the cation to interact differently with the anion.

We modify the force field by scaling the repulsive component of the LJ interaction (C12, where i = dummy atom and j = Cl–) until the calculated osmotic pressure reached the experimental value. A systematic survey of C12 values was first carried out at a concentration of 1.0 M to get the upper and lower bounds for C12 values (Figure 3a). The fitting to experimental osmotic pressure is obtained by linear regression. We started with log10[C12] = −6 and continued the variation until we obtained an optimal value of C12 that could reproduce experimental osmotic pressure for a wide range of concentrations. Osmotic pressures were plotted for intermediate sets of LJ parameters to find the best fit (intermediate and final profiles in Figure 3b). A final value of log10[C12] = −5.5 is able to match the experimental osmotic pressure for all concentrations (Figure 3b).

Figure 3

(A) Linear regression fit of osmotic pressure values for 1 M MgCl2 at different LJ parameters for interaction between dummy atoms and Cl– ions (C12); (B) Osmotic pressure profiles against different concentrations of MgCl2 for 3 different van der Waals parameter sets. The results for the original (Saxena and Sept) are shown in magenta, intermediate values are shown in yellow, and the final fit is shown in blue. (C) Radial distribution function between all atom pairs in 1 M MgCl2 solution after using the optimized parameters (solid red line). The profile is compared with the data obtained from X-ray experiments (dotted red line).

Structural Properties of MgCl2 Solution

After obtaining the refined LJ parameters for the dummy atoms interaction with Cl, we test the structural characteristics of the MgCl2 solutions by running simulations with the modified parameters. We find that the new simulations of MgCl2 solutions do not show any signs of spontaneous ion aggregations (Figure 4a). The first coordination shell of both Mg2+ and Cl– ions contained water molecules, showing evidence of solvent separated ion pair formation. The first peak in the radial distribution function profile of Mg2+ and Cl– for 1 M MgCl2 solution is now observed at 4.5 Å (Figure 3c), in contrast to the previous value of 1.95 Å (Figure 1c). Cation–anion pair at a separation of 4.5 Å corresponds to an SSP and matches the results from X-ray crystallization experiments[45] (Figure 3c). We also find that the new simulations show a total of 2 Cl– ions coordinating with Mg2+ ions, for all concentrations (below 3 M), matching the observations from neutron scattering and X-ray experiments on MgCl2.[45,46] Similar behavior was observed in some recently improved theoretical models of Mg2+ ions.[26,27] Thus, the structural characteristics of Mg2+ and Cl– ions in the simulations of MgCl2 with modified parameters confirm that the improved parameters for the multisite model are able to eliminate the ion aggregation artifact and can now reproduce accurate physical behavior of the ions at physiological conditions.

Figure 4

(A) Snapshot of a simulation of Mg2+ and Cl– ions in water after refinement of the mutisite model (Mg2+ as red spheres, Cl– as green spheres, water occupying first solvation shell of Mg2+ are shown in Licorice). No aggregations are observed now. (B) Radial distribution functions between Mg2+ and Cl– ions at different MgCl2 concentrations, after parameter refinement. (C) Osmotic pressure obtained for different concentration of CaCl2 after using 2 sets of van der Waals parameters (original (magenta) and modified parameters (blue)). (D) Radial distribution function between Ca2+ and Cl– ions at different CaCl2 concentrations, after parameter refinement.

Structural Properties of CaCl2 Solution

To test the transferability of the new multisite model parameters to other ions, we calculate the osmotic pressure of the CaCl2 solutions using the new C12 values for dummy atoms and Cl– ions, obtained from the refinement of MgCl2 solutions. Keeping everything else the same, we replaced the old LJ parameters of the dummy atom with the improved values, thereby improving the interaction of the dummy atoms with Cl– ions. The osmotic pressure values, obtained from the simulations of CaCl2 solutions with the new parameters, match well with the experimental values, without the need of further refinement (Figure 4c). This suggests that the changes made in the interaction of the dummy atoms with the anions were sufficient to reproduce accurate physiological behavior of Ca2+ ions with Cl– ions. We further verified this result by calculating the structural properties of CaCl2 in solution. The RDF profile between Ca2+ and Cl– ions shows the first peak at an ion–ion separation distance of 5 Å. This value corresponds to a solvent-separated ion pair formation in CaCl2 solutions.[32,46,59] Also, the number of Cl– ions coordinating with Ca2+ ions is 2 for all simulated concentrations (below 3 M), matching the experimental data.[32,46,59] Thus, this result suggests that the new parameters for the dummy atoms are sufficient to improve the interaction of another divalent cation with the anion, without the need for further customization. The unique parameters of the dummy atoms could be universally used whether they are part of the Mg2+ or Ca2+ multisite model.

Comparison of Results Obtained with the Saxena and Sept Model

Next, we calculate how the modeling of concentrated ionic solutions affect the water coordination and exchange rate for the cations. The original parameters of the multisite model allowed accurate prediction of the structure of the coordination shell, with the ion-oxygen distances, and the coordination number matching the experiments. In addition, the model also gave accurate values for water exchange rates, indicating that they are able to accurately represent long distance interactions.[43] In this section, we used the parameters of the new Mg2+ and Ca2+ multisite model and performed simulations in water with and without Cl– ions. The calculations without Cl– ions use identical parameters as in the work of Saxena and Sept.[43]

The resulting ion-oxygen distances, coordination number, and water exchange rates are then calculated from these simulations (Table 3). A value of 2.1 Å has been reported for Mg–O interaction[4] and 2.4 Å for Ca–O interactions.[32,60] These values were observed as 2.1 and 2.4 Å, respectively, in previous work.[43] As expected, these values remain unchanged for dilute solutions (Table 3). The refined model further allows us to calculate cation-oxygen interactions in the presence of counterions, which was not possible in our earlier simulations due to the ion aggregation artifact. We observed that the cation–oxygen interactions remain unchanged in solutions with higher anion concentrations (1M). The second property we look at is the coordination number of Mg2+ and Ca2+ ions. These values remain unchanged at 6 and 7, respectively, before and after the parameter refinement. Next we calculate the water exchange rates for different cations in concentrated CaCl2 and MgCl2 aqueous solutions. The exchange rates are difficult to determine from experiments: therefore, the available data covers a wide range of values or bounds. For Mg2+ ions, the water exchange rates have been reported in the order of microseconds. For Ca2+ ions, the rates are much faster (less than a nanosecond).[61−64] Saxena and Sept[43] showed that the multisite model significantly enhanced the water exchange properties of the cations from the previous models. They obtained a value of 20 ps for Ca2+ ions, while for Mg2+ ions the simulations were not long enough to capture the microsecond exchange rates. In this work, we calculate the exchange rates from simulations of Mg2+ and Ca2+ ions in the presence of Cl– ions. We find that the exchange rates for Ca2+ ions still remain in the ps range (∼50 ps) but with a slight increase. The increase in the exchange rate can be attributed to the presence of counterions in the simulations. For Mg2+ ions, the inner shell waters remain unchanged over a 100 ns simulation, indicating that the exchange rates would be greater than 100 ns. Thus, these results suggest that the refined multisite model preserved its interaction with water. The properties of the single ions remained unchanged from the original model and the refined model acquired an additional ability to interact with the counterions at any concentration.

Table 3

Comparison of Single Ion Properties of in Dilute and Concentrated Solutions

	system	old model[43]	refined model	experiments
	Ca2+, no Cl–	2.4 Å	2.4 Å	2.4 Å[32,60]
ion–oxygen distance (Å)	1 M CaCl2	2.4 Å	2.4 Å
	Mg2+ no Cl–	2.1 Å	2.1 Å	2.1 Å[4]
	1 M MgCl2	2.1 Å	2.1
AA
	Ca2+, no Cl–	7	7	6–7[64]
coordination number	1 M CaCl2	7	7
	Mg2+ no Cl–	6	6	6[64]
	1 M MgCl2	6	6
water exchange rates	Ca2+	20 ps	20 ps	<1 ns[62,64]
	Mg2+	>100 ns	>100 ns	>1.5 μs[61−64]
	1 M CaCl2		50 ps
	1 M MgCl2		>100 ns

Conclusion

The multisite model for divalent ions (Mg2+ and Ca2+) developed previously[43] achieved better solvation, thermodynamics, and coordination properties of divalent cations in molecular simulations. However, it suffered with artifacts of cation–anion pairing and unusual ion aggregation in highly concentrated ion solutions. In this work, we refined the parameters of the multisite model such that the physical behavior of highly concentrated solutions is reproduced accurately. The final parameters produced a good fit with the experimental osmotic pressure values for various concentrations of MgCl2 solutions, and the common artifact of ion aggregation was eliminated. Further, the changes made were directly transferable to other ions and when used in CaCl2 solutions, and experimental values of osmotic pressure were successfully reproduced. The original multisite model had presented a number of advantages over standard spherical ion models in molecular simulations of cation binding biomolecules. We verified that these properties were not lost after the refinement. The refined model went one step further in extending the range of calculations from neutralizing salt concentrations to excess salt. Taken together, the multisite cation model now presents a complete model for the use of divalent ions in classical simulations. The model allows for precise predictions of conformational and thermodynamic changes occurring in divalent cation binding biomolecules (proteins and nucleic acids) at physiologically relevant concentrations, an environment that was challenging to obtain previously in simulations of divalent ions. Since we did not change the interaction parameters for the cation with water and the biomolecule, the properties calculated for biomolecules by Saxena and Sept will still hold for the modified model described here. Additionally, the model is compatible with currently available standard force fields and has been implemented in GROMACS.[65] The application of this model will enable numerous future studies involving cation binding to biomolecules without compromising accuracy.