Magnesium Force Fields for OPC Water with Accurate Solvation, Ion-Binding, and Water-Exchange Properties: Successful Transfer from SPC/E

Magnesium plays a vital role in a large variety of biological processes. To model such processes by molecular dynamics simulations, researchers rely on accurate force field parameters for Mg2+ and water. OPC is one of the most promising water models yielding an improved description of biomolecules in water. The aim of this work is to provide force field parameters for Mg2+ that lead to accurate simulation results in combination with OPC water. Using twelve different Mg2+ parameter sets, that were previously optimized with different water models, we systematically assess the transferability to OPC based on a large variety of experimental properties. The results show that the Mg2+ parameters for SPC/E are transferable to OPC and closely reproduce the experimental solvation free energy, radius of the first hydration shell, coordination number, activity derivative, and binding affinity toward the phosphate oxygens on RNA. Two optimal parameter sets are presented: MicroMg yields water exchange in OPC on the microsecond timescale in agreement with experiments. NanoMg yields accelerated exchange on the nanosecond timescale and facilitates the direct observation of ion binding events for enhanced sampling purposes.


I. INTRODUCTION
Molecular dynamics simulations rely on accurate force field parameters for biomolecules, water molecules, and ions. It seems tempting to combine the force fields of the most promising water models with the most successful ion force fields in order to utilize the strengths of each parameter set. However, even for simple cations, the transferability of the ion parameters to different water models is limited. Different water models have a significant effect and can alter the thermodynamic and kinetic properties of the electrolyte solution considerably 1-3 . It is therefore crucial to assess whether the transfer of ion parameters to a different water model yields physically meaningful results. The aim of this work is to determine parameters for Mg 2+ in OPC water that leverage the strengths of both force fields, reproduce a broad range of experimental properties, and lead to accurate simulation results of biomolecular systems.
Water constitutes the major part in simulations of membranes, proteins, and nucleic acids. Due to its distinct role, a large variety of water models exists which differ in their complexity, accuracy, and computational efficiency 4,5 . One of the most promising recent water models is the 4-site OPC model 5,6 . OPC water was developed to accurately reproduce the electrostatic properties of water. It is quoted to improve simulations of intrinsically disordered proteins 7,8 . Balancing the interactions between the water molecules and amino acids is particularly important for disordered proteins since most models tend to favor overly compact and collapsed structures 7,9 . Tian et al. recommend to employ OPC water in combination with their recently developed protein force field ff19SB 10 . In addition, the description of hydration of small molecules 6 as well as the thermodynamics of ligand binding 11 improves with OPC water. Moreover, small RNA fragments 12 or central properties of DNA 13,14 show significant improve-ment when simulated in OPC water.
The apparent success of OPC water in simulating biological system raises the immediate question which ion force field should be used to obtain reliable results. In particular, accurate parameters for Mg 2+ are essential since these ions play a vital role in a large variety of physiological processes such as ATP hydrolysis 15 , cellular signaling 16,17 , or the catalytic activity of ribozymes [18][19][20] . Due to the prominent role of Mg 2+ , a variety of parameters for different water models exists today 1,21-34 . However, the development of accurate parameters for Mg 2+ turned out to be notoriously difficult leading to various shortcomings of the parameters: (i) Simultaneously capturing the solvation free energy and the structure of the first hydration shell failed unless polarization effects were included explicitly 28,35 . (ii) Without further optimization, the existing parameters led to unrealistically slow exchange in the first hydration shell of Mg 2+ rendering the simulation of ion binding events impossible 36 . (iii) The binding affinities to ion binding sites on biomolecules were typically overrated and needed further optimization 1,32,34,37 .
Recently, progress was made in tackling this problems by using a larger Lennard-Jones (LJ) parameter space and by modifying the standard combination rules 1, 34 . In most biomolecular simulations, the pair potential between atoms i and j is modeled as sum of the Coulomb term and the 12-6 Lennard-Jones (LJ) potential, where q i is the charge, r ij the distance, and ε 0 the dielectric constant of vacuum. The LJ term contains the interaction strength ε ij and the diameter σ ij . Increasing the LJ parameter space renders the interaction between water and Mg 2+ more attractive. This in turn allowed us to implicitly include polar-ization effects and to simultaneously reproduce the solvation free energy and the structure of the first hydration shell for 6 different water models 1,34 .
In addition, polarization effects provoke shortcomings of the standard combination rules for describing ion-ion interactions or the interaction of ions and biomolecules. By using scaling factors in the Lorentz-Berthelot combination rules [38][39][40][41] , the deviations from the standard combination rule can be taken into account. In particular, scaling parameters allowed us to reproduce activity derivatives over a broad range of MgCl 2 concentrations and to correct the excessive binding of the Mg 2+ ions to negatively charged groups on biomolecules 1,34 .
In the following, we systematically evaluate the transferability of those Mg 2+ force field parameters, that were optimized previously in combination with different water models, to OPC water. Our results show that two parameter sets, called microMg and nanoMg, that were initially optimized in combination with SPC/E water, perform best in reproducing a broad variety of experimental properties in OPC water.
In the following, we refer to these parameter sets as TIP3P-, TIP3P-fb-, TIP4P/2005-, TIP4P-Ew-, or TIP4P-D-optimized parameter sets. The corresponding Mg 2+ force field parameters are available from github https://github.com/bio-phys/ Magnesium-FFs and https://github.com/bio-phys/ optimizedMgFFs or from Table S2 in the SI. Note that for each water model two parameter sets exist, labeled microMg and nanoMg. The microMg parameter set reproduces the experimental water exchange rate on the microsecond time scale while the nanoMg parameter set yields accelerated exchanges on the nanosecond timescale. The parameters of the water models are listed in Table S1 of the SI.
The results upon transferring those force field parameters to OPC are compared to the Mg 2+ force field parameters by Li and coworkers 33 that were recently optimized in combination with OPC water. In the following, we refer to them as Li-Merz OPC-optimized parameters. In the work by Li et al. 33 , two different parameter sets are given. The first set of parameters is based on the commonly used 12-6 Lennard-Jones interaction potential. The second set of parameters is based on the so-called 12-6-4 potential which introduces an additional term to mimic polarization effects 28 .

B. Combination rules
To describe the Mg 2+ -water interactions, we use unmodified Lorentz-Berthelot combination rules (2) Note that the parameters ε Ow and σ Ow correspond to the ones for OPC water in this current work (see Table S1). The values for ε io and σ io therefore deviate from the values with the respective original water models 1,34 while the values for ε ii and σ ii remain unchanged (Table I).
To model the Mg 2+ -Cland Mg 2+ -RNA interactions, we introduce adjustable scaling parameters λ X σ and λ X ε in the Lorentz-Berthelot combination, similar to previous works 1, 32,34,[38][39][40][41] . With this, the Lorentz-Berthelot combination rules have the following form where X represents the Clanion or the atoms of an RNA. The scaling parameters allow us to take some of the effects of polarizablity into account which can cause deviations from the standard combination rule. At the same time, the Mg 2+ -water interaction remains unchanged such that the solvation free energy, the structural properties of the first hydration shell, and the rate of water exchange are not affected.

C. Simulation setup
All simulations were done with GROMACS 47 (version 2020) with the exception of the 12-6-4 interaction potentials since this interaction form is not readily available in GRO-MACS. Therefore, simulations with the 12-6-4 interaction potential were done with AMBER 48 (version 2018). We used different simulation setups to calculate the different physical properties. The details of the different setups are described in Section S3 of the SI.

D. Transferability of Mg 2+ parameters to OPC water
The simulations to evaluate the transferability of different force field parameters to OPC water were done in three consecutive steps.

Solvation free energy, radius, and coordination number
In the first step, we calculated all single-ion properties for the twelve transferred and for the two Li-Merz parameter sets. In order to provide neutral systems, we used Clas counter ion. In particular, we used the Clparameters by Mamatkulov-Schwierz 30 for TIP3P, Smith-Dang 49 for SPC/E, and Grotz-Schwierz 1 for TIP3P-fb, TIP4P/2005, TIP4P-Ew, and TIP4P-D.
We calculated the solvation free energy of neutral MgCl 2 ion pairs which includes the enthalpic and entropic contribution of the ions in water. In the computation, we took correction terms for finite size 50 and pressure 51 into account. Note that the term for interfacial crossing 51 cancels for neutral ion pairs. For the OPC-optimized 12-6 and 12-6-4 parameters, the experimental value for Cl -52 was added to the reported literature values 33 to obtain a neutral ion pair. Further details on the calculation of the solvation free energy can be found in Section S4 of the SI.
In addition, we considered the radius and coordination number as the most important structural parameters to characterize the first hydration shell. In the computation of radius and coordination number, the Sengupta-Merz 12-6 and 12-6-4 parameters 53 , that were optimized in combination with OPC, were considered as corresponding Clparameters for the Li-Merz Mg 2+ parameters.

Water exchange rate
In the second step, we selected the 8 parameter sets that performed best in the previous step and computed the water exchange rates in the first hydration shell. The exchange rates were obtained by counting transitions of water molecules between the first and second hydration shell from long trajectories in MgCl 2 . Note, however, that for some force field parameters water exchange becomes so rare that the rate could not be calculated from straight forward simulations 36 . While this limitation did not affect the transferred parameters, exchanges with the 12-6 and 12-6-4 Li-Merz parameters were rare. In both cases, water exchange was observed to be unrealistically slow compared to the experimental results and the rate could not be determined with good statistics  or not at all (12-6-4). More details on these computations can be found in Section S5 of the SI.

Activity derivative and binding affinity
In the last step, the activity derivatives in MgCl 2 and the binding affinity toward one of the non-bridging phosphate oxygens of RNA were calculated. Note that in both cases the respective modified Lorentz-Berthelot combination rules (eq 3) need to be used for accurate results. The corresponding scaling factors, λ Cl ε , λ Cl σ , λ RNA ε , and λ RNA σ are given in Tables I and S2. The activity derivatives were obtained from straight forward simulations in MgCl 2 solutions with different concentrations and Kirkwood-Buff theory (see Section S6 in SI for more details).
To calculate the binding affinity of Mg 2+ , we used the dimethylphosphate (DMP) as a simple model system containing the two non-bridging phosphate oxygen atoms. The GAFF 54 and AMBER RNA force field parameters (parmBSC0χ OL3 55-57 ) were used. Further details on the force field parameters for DMP can be found in ref. 34 . For the final two parameter sets the binding affinity was obtained from two different methods, namely by integrating potentials of mean force and via alchemical transformations. More details on the calculations of the binding affinities can be found in Section S7 of the SI.

E. Free energy profiles
To gain further insight into water exchange and Mg 2+ binding to DMP, we calculated the one-dimensional free energy profiles. For water exchange, we used umbrella sampling to obtain the free energy profiles as a function of the Mg 2+water oxygen distance. For Mg 2+ binding, we used umbrella sampling to obtain the free energy profiles as a function of the Mg 2+ -phosphate oxygen distance. Further details can be found in Section S5 and S7.1 of the SI. The barrier heights are listed in Table S6 for water exchange and in Table S7 for Mg 2+ binding to DMP.

F. Diffusion coefficient
For the final parameter sets, the self-diffusion coefficient was calculated from an additional 10 ns NVT simulation of the single Mg 2+ ion in OPC water. The first nanosecond was excluded from the analysis for equilibration. The selfdiffusion coefficient was calculated from the slope of the mean-square displacement. After a brief initial period, the mean-square displacement grows linearly and the diffusion coefficient is estimated from a straight line fit. The diffusion coefficient was corrected for system size effects 58 , where L is the box length, D pbc the computed self-diffusion coefficient, D the diffusion coefficient for the infinite nonperiodic system, k B the Boltzmann constant, T the absolute temperature, and ξ ew = 2.837297 the self-term for a cubic lattice. The empirical parameter α was set to 1.0. Since the solvent viscosity η for OPC water has not been reported in the literature, we used the experimental value η = 8.91 · 10 −4 kg m −1 s −14 .

III. RESULTS
The aim of this work is to evaluate the transferability of different Mg 2+ force fields to OPC water. In particular, we calculated all physical properties that were targeted in the initial optimization including the solvation free energy, the distance to water oxygens in the first hydration shell, the hydration number, the activity coefficient derivative in MgCl 2 solutions, and the binding affinity and distance to the non-bridging phosphate oxygens on nucleic acids. Finally, we selected two parameter sets (Table I), microMg and nanoMg, that reproduce the broad range of experimental properties (Table II) and lead to accurate simulation results in OPC water and are recommended for conventional or enhanced sampling purposes, respectively.  (Table II). σ ii , ε ii , σ io , ε io are the ion-ion and ion-water LJ parameters. λ X σ and λ X ε are the scaling factors for the Lorentz-Berthelot combination rules (eq 3) for the interaction with Clor the RNA atoms. Note that the scaling factors are only valid in combination with the Clparameters from Smith-Dang 49 for SPC/E water, and the parmBSC0χ OL3 RNA parameters 55 To evaluate whether the Mg 2+ parameters, that were optimized in combination with different 3-and 4-site water models, are transferable, we calculated the solvation free energy of MgCl 2 ion pairs in OPC water. Note that the usage of neutral ion pairs is more robust since it does not rely on conflicting experimental results for the hydration free energy of a proton 52,65,66 . Figure 1A shows the deviation of the calculated solvation free energy from the experimental results ∆∆G solv for the transferred parameters. For comparison, the values for the OPC-optimized 12-6 and 12-6-4 Li-Merz parameters 33 are shown. The simulations with the parameters transferred from SPC/E to OPC water yield the smallest deviation from experiments and perform even slightly better than those optimized directly within OPC.
Interestingly, ∆∆G solv is smaller for the parameters transferred from 3-site waters (SPC/E, followed by TIP3P-fb and TIP3P). 4-site water models (TIP4P-D, followed by TIP4P-Ew, and TIP4P/2005) lead to larger deviations despite the fact that OPC is also a 4-site model (Table S4). This result has not been expected based on previous work which showed reasonable transferability of ion parameters within water models of the same complexity 1,2,67 .
The size of the first hydration shell, measured by the distance between Mg 2+ and the water oxygens R 1 , agrees well with the experimental results for all parameter sets transferred to OPC water ( Figure 1B, Table S4). On the other hand, the OPC-optimized 12-6 model by Li-Merz significantly underestimates R 1 . This shortcoming led to the development of the 12-6-4 model 28 . Not surprisingly, the OPC-optimized 12-6-4 model 33 provides much better agreement ( Figure 1B).
For all parameter sets, the coordination number in OPC water is 6 and precisely matches the experimental result ( Table S4). Similarly, the simulations with Clparameter from different water models transferred to OPC closely reproduce R 1 (Figure 1C, Table S5). On the other hand, the OPC-optimized 12-6 Sengupta-Merz parameters deviate noticeably from the experimental value while the 12-6-4 version yields good agreement.
In summary, the six Mg 2+ parameter sets that were previously optimized for the 3-site water models (SPC/E, followed by TIP3P-fb and TIP3P) and their corresponding three Clparameter sets yield the best results for the solvation free energy and structure of the first hydration shell when transferred to OPC water and are used in the subsequent steps.

B. Water exchange rate
Including water exchange rates in the evaluation of the transferability is important to correctly capture ion binding and to avoid shortcomings due to unrealistically slow exchange dynamics as observed previously 36 . In particular, we evaluated the capability of microMg in reproducing the experimental exchange rate and the capability of nanoMg to speed up the association or dissociation processes for application in enhanced sampling simulations. Figure 2 shows the water exchange rate for the six parameter sets transferred to OPC water in comparison to the OPC-optimized Li-Merz parameters and the experimental results 61,62 . All transferred parameter sets yield a lower exchange rate in OPC water compared to the original water reflecting the influence of the water model on the exchange kinetics 3 . Still, the TIP3P-fb and SPC/Eoptimized microMg parameters yield the same order of magnitude as in experiments and are therefore considered reasonably accurate (Table S6). In addition, the transferred microMg parameters perform significantly better compared to the 12-6 and 12-6-4 Li-Merz parameters. Note, however, that there is a large uncertainty in the rate for the 12-6 Li-Merz parameters as the exchange dynamics is very slow with only 12 exchanges in 2 µs. With the 12-6-4 parameters, the exchange is even slower and not a single exchange event could be observed (Table S6). Even though it is surprising that the 12-6-4 parameters lead to a slower exchange compared to the 12-6 parameters despite including polarization effects via the r −4 term, the behavior is well-reflected in the one-dimensional free energy profiles ( Figure S1D).
For the nanoMg parameter sets, significant differences are observed: With the SPC/E-optimized parameters in OPC, exchanges happen on the order of 10 7 /s and are therefore almost as fast as in SPC/E (10 8 /s) 1 . With the TIP3P-fb-optimized and TIP3P-optimized nanoMg parameters transferred to OPC, the maximum rate is on the same order of magnitude (10 6 /s) and two orders of magnitude slower (10 8 /s to 10 6 /s), respectively. The SPC/E-optimized nanoMg parameters in OPC are therefore most beneficial for enhanced sampling of ion binding in OPC water.
In summary, based on the water exchange rate, the SPC/Eand TIP3P-fb-optimized parameters yield the best agreement with experimental results (microMg sets) while the SPC/E-TABLE II. Results for single-ion, ion-ion and ion-RNA properties in OPC for the optimal Mg 2+ parameters (Table I). Solvation free energy of neutral MgCl 2 ion pairs ∆G solv , Mg 2+ -oxygen distance in the first hydration shell R 1 , coordination number of the first hydration shell n 1 , water exchange rate from the first hydration shell k, binding affinity toward the phosphate oxygen of DMP ∆G 0 b , Mg 2+ -phosphate oxygen distance in inner-sphere coordination R b , and a cc the activity derivative of a MgCl 2 solution at 0.25 M concentration. The experimental value for ∆G 0 b is derived from the log stability constant (log K = 0.45) given in ref. 59   optimized nanoMg parameters yield the highest acceleration of exchanges.
C. Activity derivative and binding affinity to RNA Finally, we evaluate the transferability based on the activity derivative and the binding affinity and distance to the phosphate oxygens on nucleic acids. The activity derivative a cc is important in the evaluation of the transferability since it gives information on the balance between ion-ion and ion-water interactions 38 . Figure 3A shows that the TIP3P-, SPC/E-and TIP3P-fb-optimized parameters reproduce a cc over a broad concentration range in OPC water.
The binding affinity ∆G 0 b and distance R b provide insight into the accuracy of the Mg 2+ -RNA interactions and are shown in Figure 3B. The SPC/E-optimized parameters pre-cisely match ∆G 0 b and R b in OPC water, whereas parameters for TIP3P overestimate and those for TIP3P-fb underestimate ∆G 0 b .
In summary, the SPC/E-optimized microMg and nanoMg parameters yield the best agreement with the experimental results in OPC water. Figure 3C provides additional insight into the ion binding process from the free energy profiles along the Mg 2+phosphate oxygen distance for the SPC/E-optimized parameters in OPC. The two minima correspond to the inner-and outer-sphere coordination of Mg 2+ and are identical for mi-croMg and nanoMg, as expected. The energetic barrier, on the other hand, is significantly lower for nanoMg reflecting its enhanced ion binding kinetics.

IV. CONCLUSION
OPC water has proven to be one of the most promising 4site water models 6 due to its improved description of small molecules 6,11 , proteins 7,8,10 , and nucleic acids [12][13][14] . In this work, we evaluate the transferability of twelve different Mg 2+ force fields to OPC water based on solvation free energy, size of the first hydration shell, hydration number, water exchange rate, activity derivative, and binding affinity to the phosphate oxygens on nucleic acids. Our results show that the force field parameters, that were previously optimized with SPC/E 1 , are best suited to reproduce the broad range of solution properties in OPC water. In addition, the results show that the transferred parameters are compatible or better than OPC-optimized 12-6 and 12-6-4 parameters 33 in reproducing the experimental solvation free energy and water exchange rate.
Moreover, the SPC/E-optimized parameters closely reproduce the binding affinity of Mg 2+ toward the phosphate oxygen of RNA. Matching the binding affinity of metal cations to specific ion binding sites is particularly important since it improves the agreement between experiments and simulations for the structure and properties of larger nucleic acid systems 68 .
In summary, the SPC/E-optimized parameter sets microMg  in correlation with the inverse of the binding distance 1/R b toward the phosphate oxygen. The inset shows a simulation snapshot of the DMP molecule with one Mg 2+ ion in inner-sphere coordination. The gray area represents the experimental values from refs. 59,64 . (C) Free energy profiles F(R) along the Mg 2+ -phosphate oxygen distance for the optimal parameters sets in OPC water. The insets show simulation snapshots of Mg 2+ in the two stable states. and nanoMg provide an efficient and highly accurate model for the simulation of Mg 2+ in OPC water.
We recommend to use microMg to obtain accurate water exchange kinetics comparable to experimental results on the microsecond timescale. Further, we recommend to use the nanoMg parameters to yield accelerated exchange kinetics and ion-binding in enhanced sampling setups.

V. ACKNOWLEDGEMENTS
We acknowledge financial support from the DFG (Emmy Noether program, Grant No. 315221747 and the CRC902). GOETHE HLR is acknowledged for supercomputing access.

VI. SUPPORTING INFORMATION
Further details of the simulations and the analysis are provided in the supplementary material: Simulation setup, parameters of different water models and Mg 2+ force fields, details on the computation of the solvation free energy, water exchange rates, activity derivatives, radial distribution functions between Mg 2 and Cl 2 for the different models, and binding affinities. The Mg 2+ force field parameters are also available at https://github.com/bio-phys/MgFF_OPC.

VII. DATA AVAILABILITY
The data that support the finding of this study are available from the corresponding author upon reasonable request.