Optimized Magnesium Force Field Parameters for Biomolecular Simulations with Accurate Solvation, Ion-Binding, and Water-Exchange Properties in SPC/E, TIP3P-fb, TIP4P/2005, TIP4P-Ew, and TIP4P-D.

Magnesium is essential in many vital processes. To correctly describe Mg2+ in physiological processes by molecular dynamics simulations, accurate force fields are fundamental. Despite the importance, force fields based on the commonly used 12-6 Lennard-Jones potential showed significant shortcomings. Recently progress was made by an optimization procedure that implicitly accounts for polarizability. The resulting microMg and nanoMg force fields (J. Chem. Theory Comput. 2021, 17, 2530-2540) accurately reproduce a broad range of experimental solution properties and the binding affinity to nucleic acids in TIP3P water. Since countless simulation studies rely on available water models and ion force fields, we here extend the optimization and provide Mg2+ parameters in combination with the SPC/E, TIP3P-fb, TIP4P/2005, TIP4P-Ew, and TIP4P-D water models. For each water model, the Mg2+ force fields reproduce the solvation free energy, the distance to oxygens in the first hydration shell, the hydration number, the activity coefficient derivative in MgCl2 solutions, and the binding affinity and distance to the phosphate oxygens on nucleic acids. We present two parameter sets: MicroMg yields water exchange on the microsecond time scale and matches the experimental exchange rate. Depending on the water model, nanoMg yields accelerated water exchange in the range of 106 to 108 exchanges per second. The nanoMg parameters can be used to enhance the sampling of binding events, to obtain converged distributions of Mg2+, or to predict ion binding sites in biomolecular simulations. The parameter files are freely available at https://github.com/bio-phys/optimizedMgFFs.

Introduction

Magnesium is the second most abundant intracellular cation.[1] It is involved in more than 600 enzymatic reactions[2] and plays a crucial role in vital processes such as ATP hydrolysis,[3] cellular signaling,[2] or the catalytic activity of ribozymes.[4] The specific requirement for Mg2+ is prevalent in nucleic acid systems where Mg2+ is crucial for the stability, folding, and biological function.[5−9] In addition to its physiological relevance, Mg2+ is important in DNA nanotechnology[10−12] and is a promising candidate for divalent batteries.[13,14]

Due to the distinct role of Mg2+, the modeling of Mg2+ by all-atom molecular dynamics (MD) simulations has received significant attention.[15−26] Despite tremendous efforts, Mg2+ force fields based on the commonly used 12-6 Lennard-Jones interaction potential failed to provide a quantitative description (see Table S1) for three reasons. (i) No parameter combination existed that simultaneously reproduced the solvation free energy and the size of the first hydration shell.[20,21,23,27] (ii) The parameters yielded too low water exchange rates leading to unrealistically slow exchange kinetics such that transitions from outer to inner sphere binding and back could never be observed on the typical time scale of MD simulations.[28] (iii) The binding affinity of Mg2+ to ion binding sites on biomolecules was overrated significantly.[19,25,29,30] Regarding these shortcomings, the immediate question arises why classical, nonpolarizable Mg2+ force fields fail to provide an accurate description. Clearly, Mg2+ ions strongly polarize their environment, and the lack of charge transfer and polarization effects in classical simulations likely leads to the observed deviations between experiments and simulations. Possibilities to provide improvement include the use of computationally more demanding polarizable force fields,[31] to include additional parameters in the interaction potential,[22,30,32] or to scale the charge of the Mg2+ ion.[18,33] Another possibility is to include the effect of polarizability implicitly without introducing additional terms in the interaction potential. As shown in our previous work,[26] polarizability can be taken into account by two measures. First, optimizing the parameters of the 12-6 Lennard-Jones parameters in an enlarged parameter space captures the more attractive and long-ranged Mg2+–water interactions, thereby facilitating the simultaneous matching of the experimental solvation free energy, size of the first hydration shell, and water exchange rate. Second, by introducing scaling factors[34,35] in the standard combination rules for the Mg2+–Cl– and Mg2+–RNA interactions, deviations due to polarization are taken into account, and the experimental activity derivative and the binding affinity to the phosphate oxygen on nucleic acids is reproduced. The resulting force field parameters provide an efficient and highly accurate model for Mg2+ in biomolecular simulations in combination with the TIP3P water.

The TIP3P water model is frequently used in biomolecular simulations since protein and nucleic acid force fields were frequently optimized in the presence of TIP3P. However, to date, a large variety of different water models exist, some of which reproduce the physical properties of water better compared to TIP3P.[36] In principle, it is possible to combine more advanced water models and ion parameters with the force fields for biomolecules. In some cases, such a combination can leverage the strengths of the respective parameters and improve the agreement between simulated and experimental structures.[25,37−42] However, the combination of force field parameters does not guarantee that the physical properties, which were targeted in the first place, are reproduced. For metal cations, previous studies indicate limited transferability and different water models can have significant effects on the solvation free energy, the exchange kinetics, and even the reaction mechanism.[43−45] It is therefore crucial to assess the transferability of Mg2+ parameters to different water models and optimize the force field parameters if necessary in order to ensure physically meaningful results.

The aim of our current work is to provide optimized Mg2+ force field parameters in combination with five different water models. We focus on some of the most popular rigid nonpolarizable water models, namely, SPC/E,[46] TIP3P-fb,[47] TIP4P/2005,[48] TIP4P-Ew,[49] and TIP4P-D.[50] Our results illustrate that the transferability of the Mg2+ parameters[26] developed in combination with TIP3P is limited. In particular, for the 4-site water models, significant deviations from the experimental properties are observed, in agreement with similar studies.[44] In order to provide accurate parameters for the different water models, we systematically derive Mg2+ parameters that reproduce the solvation free energy, the distance to oxygens in the first hydration shell, the hydration number, the activity coefficient derivative in MgCl2 solutions, and the binding affinity and distance to nonbridging phosphate oxygens on nucleic acids, using our previously developed optimization procedure.[26] In order to provide consistent and robust results, Cl– was chosen as the reference ion, and its parameters were optimized in a preceding step. For each water model, we present two parameter sets: MicroMg yields water exchange on the microsecond time scale and matches the experimental exchange rate. NanoMg yields accelerated water exchange in the range of 106 to 108 exchanges per second depending on the water model used. The nanoMg parameters are suited to accelerate the binding kinetics in biomolecular simulations and improve the sampling of ionic distributions. Our results reveal that the largest speed-up is obtained in combination with TIP3P or SPC/E.

Methods

Molecular Dynamics Simulations

The ions are modeled as point charges, and the electrostatic, dispersion, and excluded volume interactions are taken into account by a pairwise interaction potential. We use the most common form of the Lennard-Jones (LJ) potential with a repulsive r–12 and an attractive r–6 termwhere q and q are the charges of atoms i and j, r is the distance between them, and ϵ0 is the dielectric constant of the vacuum. The parameters σ and ε describe the LJ diameter and interaction strength, respectively. Following the procedure to optimize the force field parameters for Mg2+ in TIP3P water from our previous work,[26] we systematically optimize the two parameters of the LJ potential and the combination rules for 5 different nonpolarizable rigid water models. From the different three site models, we chose the most commonly used SPC/E[46] and the more recent TIP3P-fb[47] water model. From the different four site models, we chose the TIP4P/2005,[48] TIP4P-Ew,[49] and TIP4P-D.[50] TIP4P/2005 has gained popularity and is often quoted as the best nonpolarizable general-purpose model.[36] In addition to TIP3P, TIP4P-Ew is frequently used in the AMBER and CHARMM force fields.[30,51,52] TIP4P-D is one of the newer offsprings in the TIP4P family and was designed to improve water dispersion interactions.

Partial charges, Lennard-Jones parameters, bond lengths, and bond angles of the different waters are shown in Section S1.2.

Since the standard combination rules do not take polarizability and charge transfer effects into account, we use scaling parameters in the Lorentz–Berthelot combination rules to describe the Mg2+–Cl– and the Mg2+–RNA interactions.[26] Following the work by Fyta and Netz[35] and previous similar studies,[29,34,53] the modified Lorentz–Berthelot combination rules have the following fromwhere X denotes Cl– or the atoms of the RNA. Note that the solvation free energy, the structural properties of the first hydration shell, and the rate of water exchange remain unchanged upon changing the Mg2+–Cl– and the Mg2+–RNA combination rule. In this work four scaling parameters are introduced, λσRNA, λϵRNA, λσCl, and λϵCl, to optimize ion–ion interactions and Mg2+–DMP interactions. Note that, in the RNA system, the same scaling factors λσRNA and λϵRNA are applied to all atoms of the RNA.

To compare the optimized parameters (Table ) to the available force fields from the literature, we performed simulations using the 12-6 based parameters by Åquist,[15] Mamatkulov–Netz,[20] and Li–Merz[21] (HFE set) and the 12-6-4 based parameters by Li–-Merz[22] for SPC/E water. In addition, the parameters by Li–Merz[21,22,32] (12-6 and 12-6-4) for TIP4P-Ew and TIP3P-fb were used (data not shown). For SPC/E water, the Cl– parameters were taken from Smith–Dang.[54] For all other water models, the Cl– parameters were adjusted as described below.

Table 1

Optimized Force Field Parameters for Mg2+ for Simulations with Different Water Modelsa

	σii [nm]	εii [kJ/mol]	σio [nm]	εio [kJ/mol]	λσCl	λεCl	σMgCl [nm]	εMgCl [kJ/mol]	λσRNA	λεRNA	σMgOP [nm]	εMgOP [kJ/mol]
TIP3P
microMg(26)	0.1019	235.80	0.2085	12.25	1.80	0.1	0.4878	0.8181	1.1375	0.3200	0.2262	4.6061
nanoMg(26)	0.1025	389.80	0.2088	15.75	1.80	0.1	0.4884	1.0518	1.1435	0.2500	0.2277	4.6266
SPC/E
microMg	0.1036	290.58	0.2101	13.75	1.59	0.1	0.4318	1.0989	1.1019	0.4856	0.2202	7.7589
nanoMg	0.1046	470.70	0.2106	17.50	1.59	0.1	0.4325	1.3986	1.1107	0.3300	0.2225	6.7111
TIP3P-fb
microMg	0.1032	311.38	0.2105	14.25	1.59	0.1	0.4743	0.3983	1.0957	0.4913	0.2187	8.1266
nanoMg	0.1034	380.38	0.2106	15.75	1.59	0.1	0.4744	0.4403	1.1002	0.4172	0.2197	7.6266
TIP4P/2005
microMg	0.0901	712.67	0.2030	23.50	1.59	0.1	0.4719	0.5529	1.1484	0.2648	0.2217	6.6266
nanoMg	0.0913	774.62	0.2036	24.50	1.59	0.1	0.4728	0.5764	1.1345	0.2923	0.2217	6.6266
TIP4P-Ew
microMg	0.0910	647.63	0.2037	21.00	1.59	0.1	0.4684	0.5622	1.1456	0.2778	0.2217	6.6266
nanoMg	0.0926	760.06	0.2045	22.75	1.59	0.1	0.4697	0.6096	1.1489	0.2371	0.2232	6.1266
TIP4P-D
microMg	0.0960	621.50	0.2057	20.08	1.59	0.1	0.4811	0.4697	1.3255	0.3264	0.2197	7.6266
nanoMg	0.0970	680.01	0.2062	21.00	1.59	0.1	0.4819	0.4913	1.1231	0.2834	0.2207	6.9266

σii, εii, σio, and εio are the ion–ion and ion–water LJ parameters. λσX and λεX are the scaling factors for the Lorentz–Berthelot combination rules (eq ) for the interaction with Cl– or the RNA atoms, shown exemplary for the interaction between Mg2+ and OP. Note that the scaling factors are only valid in combination with the Cl– parameters from Smith–Dang[54] for SPC/E water and those listed in Table for the other water models, and the parmBSC0χOL3 RNA parameters.[84−86] Values for parameters in TIP3P are taken from ref (26).

Table 2

Optimized Force Field Parameters for Cl– for Simulations with Different Water Modelsa

	σii [nm]	εii [kJ/mol]	σio [nm]	εio [kJ/mol]	ΔGsolvCl– [kJ/mol]	R1Cl– [nm]
Cl(TIP3P-fb)	0.493 358	0.050 960	0.4056	0.1823	–304.2 ± 1	0.320
Cl(TIP4P/2005)	0.503 464	0.042 887	0.4097	0.1823	–304.1 ± 1	0.321
Cl(TIP4P-Ew)	0.498 225	0.048 805	0.4073	0.182 3	–304.2 ± 1	0.324
Cl(TIP4P-D)	0.509 112	0.035 496	0.4128	0.1823	–304.2 ± 1	0.319
exp.					–304.2[63]	0.319 ± 0.007[87]

σii and εii are the ion–ion and σio and εio the ion–water LJ parameters. The parameters were obtained starting from Joung–Cheatham parameters[70] for TIP4P-Ew and modifying the LJ parameter σio until the experimental single-ion properties ΔGsolvCl– and R1Cl– are matched. The experimental value for R1Cl– uses Tissandier’s estimate for the proton solvation free energy.[63]

Simulations with force fields of the 12-6 type were performed with GROMACS[55] (versions 2018, 2020). Simulations with force fields of the 12-6-4 type were done with AMBER[52] (version 2018) and PLUMED[56] (version 2.6.4) since GROMACS does not support 12-6-4 interaction potentials. Section S1.3 (Tables S3 and S4) lists all simulation setups. Similar to previous work,[19,26,30,57,58] we used dimethylphosphate (DMP) to optimize the ion–RNA interactions. The force field parameters for the DMP molecule are based on a parametrization with GAFF[59] and can be found in ref (26). The analysis was performed with built-in GROMACS[55] code and by using the MDAnalysis package[60,61] for python.

Optimization Procedure

The optimization procedure follows the same three step strategy described in ref (26), starting with an optimization of the ion–water interaction by performing a grid search in σio–εio space. Here, in a first step, all σio–εio parameter combinations matching the experimental solvation free energy ΔGsolv, the Mg2+-oxygen distance in the first hydration shell R1, and the coordination number n1 are selected.

Second, by calculating the rate of water exchange for the above-mentioned parameter combinations, we optimize the water exchange dynamics. For each water model two parameter sets were chosen: The microMg parameter sets yield water exchange on the microsecond time scale and reproduce the experimental rate within errors. The second parameter sets yield faster water exchange time scales providing exchange dynamics that range between 106 to 108 exchanges per second while simultaneously reproducing all thermodynamic and structural properties. For consistency with our previous work,[26] we refer to the second set as nanoMg, even though the nanosecond time scale could not be reached for all water models.

In the final step, we optimize the ion–ion and ion–RNA interactions by calculating activity coefficient derivatives and ion binding affinities. We performed a grid search in λσX and λεX parameter space (eq ). The activity coefficient derivatives acc are collected and the scaling factors λσCl and λεCl are selected by reproducing the experimental activity derivative over a broad range of MgCl2 concentrations, using the Kirkwood–Buff theory.[62] To calculate the binding affinity of Mg2+ toward one of the nonbridging phosphate oxygens of DMP, we used alchemical transformation calculations. Reproducing the experimental binding affinity ΔGb0 and binding distance Rb toward the phosphate oxygen, we select the scaling factors λσRNA, λεRNA.

Solvation Free Energy

The solvation free energy is considered the most important thermodynamic property in the development of accurate force field parameters.[27] Yet, the proton solvation free energy can vary depending on the experimental sources. In most cases, the solvation free energy of the proton by Tissandier et al.,[63] −1104 kJ/mol, or Marcus,[64] −1056 kJ/mol, are used. Since, depending on the exact choice, the estimates can vary by more than 50 kJ/mol, we take the more robust solvation free energy of neutral ion pairs into account. In the following, we use chloride as the reference ion as in previous works.[20,23,27] In order to obtain consistent results, appropriate corrections for finite size effects, compression, and interfacial crossing have to be applied. Finite size effects are taken into account via[65]where z is the valency, NA Avogardo’s number, e the unit charge, ϵ0 the vacuum permittivity, and R1 the first peak of the ion–water radial distribution function. Hence, the effective ion radius, ζew = −2.837 279/L, is the Wigner potential, with L being the edge length of the simulation box in nanometers.[65,66] ϵr is the relative dielectric constant of the different water models. We used ϵr(SPC/E) = 68, ϵr(TIP3P-fb) = 81.3, ϵr(TIP4P/2005) = 58, ϵr(TIP4P-Ew) = 63, and ϵr(TIP4P-D) = 68. The values for SPC/E and TIP4P/2005 are taken from ref (36), and all other values are from their original publications.[47,49,50]

The correction term related to the compression of the gas is given by[27]and is independent from the choice of the water model, with p0 = 1 atm being the pressure of ideal gas and p1 = 24.6 atm corresponding to an ideal solution under pressure at a density of 1 mol/L. In the experiments,[63] the ions have to pass the air–water interface in order to enter into the aqueous phase. The corresponding free energy correction term is[27]Here, we chose the surface potential as ϕsurf = −0.527 V obtained for polarizable TIP4P-FQ water.[67,68] The reason is that this choice closely matches with the experimentally based estimation of −0.50 V.[69] Even more importantly, the interfacial crossing term almost exactly cancels the shift between the solvation free energies by Marcus[64] and by Tissandier[63] for different anions and cations. Consequently, we obtain consistent results if we compare the simulation results with the interfacial crossing term to the free energies by Tissandier (since those values include a full contribution from the bulk water surface potential[69]) or if we compare the simulation results without interfacial crossing term to the bulk free energies by Marcus (which omit the surface potential[69]).

Using the former approach, the solvation free energy is given byThe solvation free energy of the neutral MgCl2 ion pairs is given by

Parametrization of Cl– and Mg2+

Initially, we optimize the parameters of the reference chloride ion. For SPC/E water, we use the well established Smith–Dang parameters.[54] In particular, the calculated solvation free energy is ΔGsolvCl– = −306 kJ/mol[27] and closely matches the experimental results by Tissandier[63] (−304.2 kJ/mol). For the optimization of the parameters for Cl– in combination with TIP3P-fb, TIP4P/2005, TIP4P-Ew, and TIP4P-D water, we used two different approaches: In the first approach, we modified εio starting from the Smith–Dang parameters[54] (σio = 0.378 nm, εio = 0.52 kJ/mol). We selected the values for εio that reproduced the experimental results while keeping σio fixed. In the second approach, we modified σio starting from the Joung–Cheatham chloride parameters[70] for TIP4P-Ew water (σio = 0.404 105 546 nm, εio = 0.182 336 936 kJ/mol). We selected the values for σio that reproduced the experimental results while keeping εio fixed. As shown in Figure , the second approach allowed us to simultaneously reproduce ΔGsolvCl– and R1, the radius of the first hydration shell, for all water models. The parameters from the second approach were therefore used in the optimization of Mg2+. The resulting LJ parameters are listed in Table .

Figure 1

Solvation free energy ΔGsolvCl– of Cl– in correlation with the inverse of the Cl––oxygen distance of the first hydration shell 1/R1Cl–. The gray area indicates the experimental results from refs (64 and 87). One set of parameters is obtained by starting from Smith–Dang parameters[54] for SPC/E water and modifying the LJ parameter εio (orange, Cl_set1). The second set is obtained by starting from Joung–Cheatham parameters[70] for TIP4P-Ew water and modifying the LJ parameter σio (purple, Cl_set2).

Subsequently, the solvation free energy for Mg2+ was calculated on a grid with σio = 0.195–0.22 nm and εio = 0.02–28 kJ/mol. In addition, R1, n1, and the self-diffusion coefficient D0 were calculated for all waters as described in ref (26).

Free Energy Profiles, Water Exchange Rate, Activity Derivative, and Binding Affinity

We used umbrella sampling[71,72] to calculate one- and two-dimensional free energy profiles, straightforward simulations to calculate the water exchange rate (see Section S1.5 for more details), Kirkwood–Buff theory[62] to calculate the activity derivative with respect to the natural logarithm of the number density of the Mg2+ ions (see Section S1.6 for more details), and alchemical transformation to calculate the binding affinities. Additional details regarding the calculation of free energy profiles (Section S1.7) and binding affinities (Section S1.8) can be found in the Supporting Information of this work.

Results and Discussion

In the following, we present the results from our parameter optimization of Mg2+ in combination with the SPC/E, TIP3P-fb, TIP4P/2005, TIP4P-Ew, and TIP4P-D water models. Our optimization procedure is done in three sequential steps[26] and allows us to reproduce a broad range of thermodynamic properties including the solvation free energy, the distance to oxygens in the first hydration shell, the hydration number, the activity coefficient derivative in MgCl2 solutions, and the binding affinity and distance to the phosphate oxygens of RNA (Tables and 5). For each water model, we present two sets of optimal parameters: MicroMg yields water exchange on the microsecond time scale and matches the experimental exchange rate. NanoMg yields accelerated water exchange in the range of 106 to 108 exchanges per second depending on the water model used (Table ).

Table 3

Results for Single-Ion and Ion–Ion Properties for the Optimized Parameters in Direct Comparison with Experimental Resultsa

	ΔGsolv [kJ/mol]	R1 [nm]	n1	D0 [10–5 cm2/s]	acc
microMg(TIP3P)[26]	–2532.9 ± 1	0.207 ± 0.004	6	0.754 ± 0.006	0.93 ± 0.01
nanoMg(TIP3P)[26]	–2532.0 ± 1	0.209 ± 0.004	6	0.750 ± 0.004	0.97 ± 0.01
microMg(SPC/E)	–2530.5 ± 1	0.209 ± 0.004	6	0.475 ± 0.004	0.76 ± 0.02
nanoMg(SPC/E)	–2530.5 ± 1	0.213 ± 0.004	6	0.497 ± 0.009	0.86 ± 0.03
microMg(TIP3P-fb)	–2531.3 ± 1	0.211 ± 0.004	6	0.509 ± 0.014	0.96 ± 0.05
nanoMg(TIP3P-fb)	–2532.0 ± 1	0.211 ± 0.004	6	0.463 ± 0.017	0.95 ± 0.02
microMg(TIP4P/2005)	–2531.8 ± 1	0.210 ± 0.004	6	0.548 ± 0.007	0.87 ± 0.04
nanoMg(TIP4P/2005)	–2530.8 ± 1	0.211 ± 0.004	6	0.495 ± 0.013	0.84 ± 0.05
microMg(TIP4P-Ew)	–2530.9 ± 1	0.209 ± 0.004	6	0.475 ± 0.011	0.85 ± 0.03
nanoMg(TIP4P-Ew)	–2530.8 ± 1	0.211 ± 0.004	6	0.469 ± 0.008	0.87 ± 0.05
microMg(TIP4P-D)	–2530.8 ± 1	0.212 ± 0.004	6	0.531 ± 0.016	1.01 ± 0.02
nanoMg(TIP4P-D)	–2530.6 ± 1	0.213 ± 0.004	6	0.530 ± 0.013	0.89 ± 0.03
exp.	–2532[64]	0.209 ± 0.004[87]	6[87]	0.706[64]	0.93[88]

Solvation free energy of neutral MgCl2 ion pairs ΔGsolv, Mg2+–oxygen distance in the first hydration shell R1, coordination number of the first hydration shell n1, self-diffusion coefficient D0, and acc the activity derivative of a MgCl2 solution at 0.25 M concentration. Values for parameters in TIP3P are taken from ref (26).

Table 5

Results for Ion–RNA Properties for the Optimized Parameters in Direct Comparison with Experimental Binding Affinity Towards the Phosphate Oxygen of DMP, ΔGb0, and Mg2+–Phosphate Oxygen Distance in Inner-Sphere Conformation, Rba

	ΔGb0 [kBT]	Rb [nm]	ΔF* [kBT]
microMg(TIP3P)[26]	–0.633 ± 0.6	0.207 ± 0.004	17.2
nanoMg(TIP3P)[26]	–0.375 ± 0.1	0.207 ± 0.004	13.6
microMg(SPC/E)	–1.26 ± 0.62	0.209 ± 0.004	16.5
nanoMg(SPC/E)	–1.19 ± 0.50	0.208 ± 0.004	13.5
microMg(TIP3P-fb)	–1.13 ± 0.90	0.208 ± 0.004	17.6
nanoMg(TIP3P-fb)	–1.15 ± 0.24	0.208 ± 0.004	16.5
microMg(TIP4P/2005)	–0.97 ± 0.24	0.208 ± 0.004	15.5
nanoMg(TIP4P/2005)	–0.88 ± 0.48	0.207 ± 0.004	14.5
microMg(TIP4P-Ew)	–0.95 ± 0.46	0.208 ± 0.004	15.3
nanoMg(TIP4P-Ew)	–0.91 ± 0.32	0.208 ± 0.004	13.7
microMg(TIP4P-D)	–1.02 ± 0.27	0.207 ± 0.004	16.4
nanoMg(TIP4P-D)	–1.15 ± 0.23	0.207 ± 0.004	15.6
exp.	–1.036[77]	0.206–0.208[78]	n.a.

ΔGb0 is derived from the log stability constant of the DMP (log K = 0.45) given in ref (77). Barrier heights ΔF* of microMg and nanoMg separate the first and second minima of the free energy profile along the distance between one of the nonbridging phosphate oxygens of the DMP and Mg2+. ΔF* is the free energy difference between the top and the first minimum (Figure B,C). Values for usage in TIP3P are taken from ref (26).

Table 4

Properties of Water Exchange from Simulations and Experimentsa

	N	k [s–1]	ΔF* [kBT]
microMg(TIP3P)[26]	376 ± 56	(8.04 ± 1) × 105	15.9
nanoMg(TIP3P)[26]	52 086 ± 120	(1.11 ± 0.003) × 108	11.5
microMg(SPC/E)	452 ± 52	(9.62 ± 1) × 105	15.5
nanoMg(SPC/E)	47 472 ± 3620	(1.01 ± 0.08) × 108	11.2
microMg(TIP3P-fb)	184 ± 4	(3.94 ± 0.09) × 105	16.1
nanoMg(TIP3P-fb)	1344 ± 20	(2.88 ± 0.05) × 106	14.1
microMg(TIP4P/2005)	308 ± 16	(6.56 ± 0.4) × 105	15.5
nanoMg(TIP4P/2005)	1554 ± 90	(3.32 ± 0.2) × 106	14.1
microMg(TIP4P-Ew)	660 ± 16	(1.41 ± 0.03) × 106	15.0
nanoMg(TIP4P-Ew)	8618 ± 262	(1.84 ± 0.06) × 107	12.6
microMg(TIP4P-D)	312 ± 8	(6.65 ± 0.2) × 105	15.4
nanoMg(TIP4P-D)	1780 ± 20	(3.80 ± 0.05) × 106	13.7
exp.	248,[74] 314[73]	5.3 × 105 from ref (73), 6.7 × 105 from ref (74)	n.a.

Number of transitions N in 1 μs for different Mg2+ parameters in 1 M MgCl2 solutions. The experimental value[73,74] is obtained from eq S3. The rate constant k is calculated from the number of transitions N for microMg and nanoMg and eq S3. The errors for N and k are obtained from block averaging. ΔF* is the free energy difference between the top and the first minimum (Figure B). Values for parameters in TIP3P are taken from ref (26).

Transferability of Ion Parameters between Different Water Models

We tested the transferability of the recently developed Mg2+ parameters[26]microMg(TIP3P) and nanoMg(TIP3P), that were optimized in combination with the TIP3P water model. Figure and Table S3 show that the transferability of the parameters to other water models is limited and deviations from the experimental solvation free energy are observed. While for the 3-site water models (TIP3P-fb and SPC/E) the deviations are small (3–4 kJ/mol), significant deviations up to 80 kJ/mol are observed for the 4-site models (TIP4P/2005, TIP4P-Ew, and TIP4P-D). The influence of the different water models on the solvation free energy ΔGsolv, the radius R1 of the first hydration shell, and the coordination number n1 is also reflected in the free energy isolines which display differences between 3- and 4-site water models (Figure S3). Since the transferability of the Mg2+ parameters is limited, in agreement with a similar previous study on metal cations,[44] we have systematically optimized the parameters in combination with 5 different water models.

Figure 2

Difference in solvation free energy to the experimental value of −2532 kJ/mol[64] for parameters transferred to different water models. The microMg and nanoMg parameters[26] that have been optimized in TIP3P are used with one the following water models W (SPC/E,[46] TIP3P-fb,[47] TIP4P/2005,[48] TIP4P-Ew,[49] and TIP4P-D[50]), resulting in the converted parameter sets microMg(TIP3P|W) and nanoMg(TIP3P|W). Additional single-ion properties are shown in Table S7.

Optimization of Solvation Free Energy, Radius of Hydration Shell, and Coordination Number

For all 5 water models, we selected parameter combinations for σio and εio that accurately reproduce ΔGsolv, R1, and n1 simultaneously (Figure A and Table ). Hereby, the key to a successful parametrization was to take polarizability into account implicitly by considering a larger range of possible LJ parameters.[26] As shown exemplary for SPC/E water in Figure B (for the other waters see Figure S4), the interaction potentials with our optimized parameter sets microMg and nanoMg are more long-ranged compared to the standard parameter sets from the literature[15,20,21] and similar to 12-6-4 potentials, which explicitly consider ion–dipole interactions via additional terms in the interaction potential.[22] Thus, our optimized parameters take polarization effects into account implicitly, and the interaction between Mg2+ and water is more attractive and long-ranged compared to the standard force fields from the literature. This in turn allows us to reproduce the experimental results for ΔGsolv, R1, and n1 simultaneously and to achieve agreement comparable to 12-6-4 interaction potentials[22] (Figure A) without introducing additional parameters in the functional form of the interaction potential (eq ).

Figure 3

Comparison of the optimized parameter sets microMg (blue) and nanoMg (green) of various water models (SPC/E,[46] TIP3P-fb,[47] TIP4P/2005,[48] TIP4P-Ew,[49] and TIP4P-D[50]) with force fields from the literature[15,20−22,32] and experimental data. (A) Solvation free energy ΔGsolv for neutral MgCl2 pairs in correlation with the inverse of the Mg2+–oxygen distance of the first hydration shell 1/R1. The gray area indicates the experimental results from refs (64 and 87). The different marker shapes indicate Mg2+ models of different water models. The colors indicate the different Mg2+ models. The values for Åquist were taken from ref (89). Solvation free energy values for parameters by Li–Merz (both 12-6 and 12-6-4 based) and Åquist were combined with Cl– values by Marcus.[87] (B) Lennard-Jones interaction potential VLJ as a function of the Mg2+–oxygen distance rMgOx for different Mg2+ force fields and exemplary SPC/E water. Similar plots for the microMg and nanoMg parameters of the other water models can be found in Figure S4.

Optimization of the Water Exchange Rate

In the next step of our optimization, we selected two sets of parameters from the data that reproduced ΔGsolv, R1, and n1 in the previous step based on the calculated water exchange rate. In order to provide an accurate estimate of the rate constant, the rate was calculated from 1 μs long straightforward simulations since transition state theory, albeit computationally less demanding, can overestimate the true rate by more than 2 orders of magnitude.[28] The rates for the two sets of parameters for all water models are shown in Figure A and Table . The first parameter sets, microMg, yield water exchange on the microsecond time scale in agreement with experimental results.[73,74] The second parameter sets, nanoMg, yield accelerated water exchange with 106 to 108 exchanges per second dependent on the water model used (Table ). In all cases, the parameters for nanoMg were chosen such that they yield the highest rate while still reproducing all other experimental properties. Yet, the different properties of the water models,[36] in particular the dielectric constant of the models, lead to distinct differences in the parameter range that reproduces ΔGsolv, R1, and n1 (Figure S3). Consequently, the maximum exchange rate that could be achieved differs for the different water models. NanoMg(SPC/E) yields the highest exchange rate with 108 exchanges per second, similar to the acceleration achieved for TIP3P.[26]NanoMg(TIP4P-Ew) yields the next highest rate with 107 s–1 followed by nanoMg(TIP3P-fb), nanoMg(TIP4P/2005), and nanoMg(TIP4P-D) with 106 exchanges per second (Table and Figure A).

Figure 4

Water exchange in the first hydration shell of Mg2+. (A) Kinetic rate coefficients k (eq S3) of microMg (blue) and nanoMg (green) of different water models (SPC/E,[46] TIP3P-fb,[47] TIP4P/2005,[48] TIP4P-Ew,[49] and TIP4P-D[50]) and for TIP3P from the literature.[26] The gray vertical line indicates the experimental values[73,74] (Table ). (B) One-dimensional free energy profiles as a function of the distance between Mg2+ and the leaving water molecule RMg–Ox for different force fields in combination with SPC/E water. (C) The snapshots show representative conformations in the two stable states: (i) Before exchange: Leaving water (shown in green) is part of the first and entering water (shown in blue) is part of the second hydration shell. (ii) After exchange: Leaving water is in the second hydration shell and the entering water molecule filled the void in the first hydration shell. The snapshots were taken using microMg(SPC/E).

The improvement of the exchange rate compared to the standard sets from the literature is shown exemplary for SPC/E in Figure B (for the other waters, see Figure S5). The barrier height between the two stable states (i) and (ii), which corresponds to replacing one water molecule in the first hydration shell with a water molecule from the second shell (green and blue water molecules shown in Figure C), differs by more than 10 kBT for different force fields (Figure B). The different barrier heights provide a qualitative explanation of the different water exchange rates obtained with different force fields from the literature. Force fields with slow exchange kinetics (Åquist,[15] Mamatkulov–Netz,[20] and Li–Merz[21]) have high free energy barriers while force fields with a more attractive and long-ranged interaction potential (microMg and Li–Merz[12-6-4][22]) have lower energetic barriers and yield exchange rates close to the experimental result. Finally, nanoMg has the lowest barrier, resulting in the fastest exchange kinetics that is particularly useful for enhanced sampling of Mg2+ binding, as will be discussed further below.

The microscopic mechanism of water exchange is, however, more complex than the one-dimensional free energy profiles suggest. The exchange involves the concerted motion of two water molecules[28,45] and is captured more realistically by two-dimensional free energy profiles (Figures S6 and S7).

In addition, the self-diffusion coefficients D0 were calculated for each water model. Unlike the results in TIP3P water, where the diffusion coefficient matched the experimental value without optimization, D0 is underestimated for all water models investigated in the present study (Table ). Likely, the better agreement in TIP3P is due to the higher self-diffusion constant of TIP3P.[36]

Optimization of the Activity Derivative

Subsequently, we balance ion–water and ion–ion interactions by matching with experimental activity derivatives.[34,35,53,75,76] Similarly to previous work,[20,23,25,26,35] the standard combination rules (eq ) (with λσ,εCl– = 1.0) overestimate the anion–cation interaction and consequently underestimate the activity derivative over the full parameter range. Since the standard combination rules are valid only in idealized cases, polarization and charge transfer can lead to deviations. Introducing scaling factors in the combination rule allows us to take these effects into account and to provide closer agreement with experimental results without changing ΔGsolv, R1, n1, D0, or k. In all cases, increasing the effective size, σMgCl, and decreasing the cation–anion interaction energy, εMgCl, reproduces the experimental activity coefficient derivative of MgCl2 solutions over a broad concentration range (Figure ). In particular, we found one set of universal scaling factors (λσCl = 1.59 and λεCl = 0.1) that leads to reasonable agreement with the experimental values for all five water models (Table , Figure ). Note, however, that, at high concentrations, individually adjusted combination rules for water models such as TIP4P-D or SPC/E might lead to better agreement.

Figure 5

Optimization of ion–ion interactions. Activity derivative acc as a function of the MgCl2 salt concentration with the optimized scaling factors for microMg and nanoMg of the various water models (SPC/E,[46] TIP3P-fb,[47] TIP4P/2005,[48] TIP4P-Ew,[49] and TIP4P-D[50]) (Table ) and from experiments.[88] The inset shows a representative snapshot of Mg2+ in interaction with Cl–, including their first hydration shell. This snapshot was taken using microMg(SPC/E) and Cl– parameters from Smith–Dang.[54]

Optimization of the Binding Affinity to Nucleic Acids

So far, the optimization was based only on bulk ion properties. However, the transferability of such bulk-optimized ion parameters to describe the interactions of Mg2+ with biomolecules turned out to be limited.[25,26,30] Therefore, in a last step, we optimized the interaction between Mg2+ and nucleic acids by matching with the experimental binding affinity and distance to one of the nonbridging phosphate oxygens.[77,78] Similar to previous work,[19,26,30,57,58] a dimethylphosphate (DMP) is used to mimic the phosphate backbone, the most important inner-sphere Mg2+ binding sites on larger RNA molecules.[5,77,79−81] Here, the standard combination rules (eq with λσ,εRNA = 1.0) significantly overestimate the binding affinity toward the phosphate oxygen reflecting that the Mg2+–RNA interactions are too attractive (data not shown). As before, this problem can be solved by increasing the effective diameter via λσRNA and simultaneously reducing the interaction energy via λεRNA (Table ). Note that, in this case, the scaling parameters for the different water models are similar but not identical since, unlike for chloride, the force field parameters of the RNA were not adjusted to the different water models.

Figure shows the binding affinity ΔGb0 and binding distance Rb for different water models and force field parameters from the literature.[19,20,22,30] Clearly, the bulk-optimized parameters (Mamatkulov–Netz(SPC/E), Allnér–Villa(TIP3P)) significantly overestimate ΔGb0 and underestimate Rb. The Panteva–York[m12-6-4] parameters[30] provide improvement. However, in the optimization the 4-fold access of the phosphate oxygen binding site on the backbone compared to the nucleobase binding sites[77] was taken into account, which is only applicable in the context of the modular model used in the experimental study. The affinity is thus overrated. Also, the binding distance Rb was not explicitly considered. With the optimized ion-RNA scaling factors for microMg and nanoMg, the results for ΔGb0 and Rb agree within errors with the experimental results for all 5 water models (Figure , Table , and Tables S9 and S10).

Figure 6

Optimization of ion–RNA interactions. (A) Binding affinity ΔGb0 in correlation with the inverse of the Mg2+–phosphate oxygen distance 1/Rb for microMg and nanoMg in the water models SPC/E,[46] TIP3P-fb,[47] TIP4P/2005,[48] TIP4P-Ew,[49] or TIP4P-D.[50] The inset is a zoomed in view of the experimental area. The experimental values (gray) are taken from refs (77 and 78). (B, C) One dimensional free energy profiles along the distance between one of the nonbridging phosphate oxygens of the dimethylphosphate (DMP) and the Mg2+. The free energy profiles from Allnér–Villa and Panteva–York were obtained from refs (19 and 30) with permission. (D) Representative snapshot of the DMP in interaction with Mg2+, including the first hydration shell of Mg2+, in (i) inner-sphere and in (ii) outer-sphere contact. Snapshots were taken using microMg(SPC/E).

The Mg2+ force field parameters and the water model have a significant influence on the ligand exchange kinetics.[45] The effect is illustrated by the one-dimensional free-energy profiles as a function of the distance between Mg2+ and the phosphate oxygen for different force fields and water models (Figure B,C). The free energy profiles show two stable states (i) and (ii), corresponding to the inner-sphere conformation (direct contact between Mg2+ and phosphate oxygen) and the water-mediated outer-sphere conformation (Figure D). The bulk-optimized parameters (Mamatkulov–Netz(SPC/E), Allnér–Villa(TIP3P)) significantly overestimate ΔGb0 and underestimate Rb as reflected by the shift and depth of the first minimum (Figure B). For microMg and nanoMg, the inner-sphere and outer-sphere minima are identical, reflecting that the parameters have identical thermodynamic and structural properties. In addition, clear differences are observed for the barrier height, which corresponds in large part to the free energy cost of exchanging one water molecule from the first hydration shell with a water from the second shell. Force fields that underestimate the rate of water exchange (Mamatkulov–Netz(SPC/E), Allnér–Villa(TIP3P)) have high energetic barriers, Mg2+ association/dissociation too slow, binding affinity too high, and binding distance too small (Figure B). MicroMg and nanoMg provide significant improvement (Table ). At the same time, the free energy barrier for nanoMg is up to 4 kBT lower compared to microMg depending on the water models used (Figure B, Table ).

Reproducing ΔGb0 and Rb is crucial to correctly describe the structure and thermodynamics of specifically bound Mg2+ ions. Since microMg and nanoMg reproduce the properties of specifically bound ion as well as important bulk properties, both parameter sets are equally suited to calculate the distribution of Mg2+ around biomolecules. In simulations that target the binding kinetics of the metal cations, the microMg parameter set should be used. To enhance the sampling of binding events, the nanoMg set can be used since the exchange kinetics is accelerated, thereby improving the sampling. However, the acceleration of the binding kinetics strongly depends on the water model used. For the largest speed-up, we recommend using nanoMg in combination with the TIP3P or SPC/E water model (Figure A).

Conclusions

The distinct role of magnesium in biological systems has driven the development of force fields for molecular dynamics simulations. Despite considerable efforts, Mg2+ force fields based on the 12–6 Lennard-Jones (LJ) potentials showed significant shortcomings in reproducing a broad range of structural, thermodynamic, and kinetic properties. Since Mg2+ strongly polarizes its environment, the lack of polarizability and charge transfer effects in classical simulations likely causes such deviations. Recently, progress was made by an optimization procedure that implicitly accounts for polarizability. Considering an enlarged range of possible LJ parameters for Mg2+–water interactions and optimizing the Lorentz–Berthelot combination rule for Mg2+–Cl– and Mg2+–RNA interactions yielded simple and accurate force fields in TIP3P water.[26]

However, a large variety of water models exist, and models with improved properties are increasingly used in biomolecular simulations. Since simulations rely on available water models and ion force fields that provide physically meaningful results when combined,[82] the transferability of the Mg2+ parameters to other water models needs to be assessed and the parameters need to be optimized if necessary. Our results reveal that the Mg2+ parameters developed in combination with TIP3P show limited transferability to SPC/E, TIP3P-fb, TIP4P/2005, TIP4P-Ew, and TIP4P-D. While for the 3-site water models (TIP3P-fb and SPC/E) the deviations from the experimental solvation free energy are small, significant deviations are observed for the 4-site models (TIP4P/2005, TIP4P-Ew, and TIP4P-D).

To provide improvement, we have systematically developed improved Mg2+ parameters. The optimized parameters (Table ) reproduce the solvation free energy, the distance to oxygens in the first hydration shell, the hydration number, the activity coefficient derivative in MgCl2 solutions, and the binding affinity and distance to the phosphate oxygens on nucleic acids (Tables , 5). In order to provide consistent and robust results for the solvation free energy and the activity derivative, Cl– was chosen as a reference ion and its parameters were optimized in a preceding step for each water model (Table ). Here, an alternative and promising approach for future work is the simultaneous optimization of anion and cations parameters, which eliminates the necessity to select a reference ion.[83] In addition, including experimental binding affinities toward specific ion binding sites on biomolecules has turned out to be essential to capture the structure of specifically bound ions[26] as well as to reproduce the structural properties of large nucleic acids.[42]

Two parameter sets are presented for each water model: The first sets, microMg, yield water exchange on the microsecond time scale in agreement with experimental results.[73,74] For the second sets, nanoMg, the parameters were chosen to yield the highest exchange rate possible while still reproducing all other thermodynamic and structural properties. Since the water models have different properties,[36] including different dielectric constants and diffusion coefficients, the maximum achievable exchange rate differs and ranges from 106 to 108 exchanges per second (Table ).

In summary, the Mg2+ parameters presented here provide simple, computationally efficient, and highly accurate models for the simulation of Mg2+ ions in aqueous solutions of SPC/E, TIP3P-fb, TIP4P/2005, TIP4P-Ew, and TIP4P-D water. For simulations targeting the kinetics of ion pairing or ion binding, the microMg is recommended. For simulations targeting the distribution of Mg2+ around nucleic acids, proteins, or lipids the nanoMg is recommended to enhance the sampling of binding events. Here, the parameters for TIP3P[26] or SPC/E yield the largest enhancement.