Next Article in Journal
Effect of Heat Treatment on the Microstructure and Corrosion Resistance of Al0.75CoCr1.25FeNi High-Entropy Alloys
Previous Article in Journal
Microstructure Evolution and Tensile Properties of Medium Manganese Steel Heat Treated by Two-Step Annealing
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Design of Point Charge Models for Divalent Metal Cations Targeting Quantum Mechanical Ion–Water Dimer Interactions

by
Yongguang Zhang
,
Binghan Wu
,
Chenyi Lu
and
Haiyang Zhang
*
Department of Biological Science and Engineering, School of Chemistry and Biological Engineering, University of Science and Technology Beijing, Beijing 100083, China
*
Author to whom correspondence should be addressed.
Metals 2024, 14(9), 1009; https://doi.org/10.3390/met14091009
Submission received: 24 July 2024 / Revised: 30 August 2024 / Accepted: 2 September 2024 / Published: 3 September 2024

Abstract

:
Divalent metal cations are of vital importance in biochemistry and materials science, and their structural and thermodynamic properties in aqueous solution have often been used as targets for the development of ion models. This study presented a strategy for designing nonbonded point charge models of divalent metal cations (Mg2+ and Ca2+) and Cl by targeting quantum mechanics (QM)-based ion–water dimer interactions. The designed models offered an accurate representation of ion–water interactions in the gas phase and showed reasonable performance for non-targeted properties in aqueous solutions, such as the ion–water oxygen distance (IOD), coordination number (CN), and density and viscosity of MgCl2 and CaCl2 solutions at low concentrations. Our metal cation models yielded considerable overestimates of the hydration free energies (HFEs) of the ions, whereas the Cl model displayed good performance. Together with the overestimated density and viscosity of the salt solutions, these results indicated the necessity of re-optimizing ion–ion interactions and/or including polarization effects in the design of ion models. The designed Mg2+ model was capable of maintaining the crystal metal-binding networks during MD simulation of a metalloprotein, indicating great potential for biomolecular simulations. This work highlighted the potential of QM-based ion models to advance the study of metal ion interactions in biological and material systems.

1. Introduction

Metal cations of magnesium (Mg2+) and calcium (Ca2+) are widely distributed in nature and play crucial roles in biology, chemistry, and materials science. For instance, Mg2+ often serves as a cofactor for hundreds of enzymes, promoting their ability to achieve corresponding biological functions [1,2]. Similarly, Ca2+ ion is indispensable in cellular signaling pathways [3,4,5], muscle contraction [6], and bone formation [7,8]. The unique coordination properties of these ions have also been used to develop advanced materials, such as biocompatible implants [9,10] and responsive materials [11,12,13]. Given their diverse roles, a deep understanding of their behaviors in aqueous solution is imperative for interpreting the underlying interaction mechanism. Computational approaches offer a powerful tool to study ions, providing detailed insights into their interactions at the atomic level. The predictive capability is highly dependent on the accuracy of potential functions (known as force fields) used to model molecular interactions in computational simulations [14]. However, the development of accurate force field models for divalent metal ions remains a challenge due to complex electronic structures and diverse coordination environments.
In recent years, several strategies have been in use for ion modeling in molecular dynamics (MD) simulations, including the bonded model [15,16], nonbonded dummy model [17,18,19,20,21], and nonbonded point charge model [22,23,24,25]. The bonded model, by constraining the chemical bonds between metal ions and ligands, is capable of accurately reproducing the geometric structure of the metal center. However, this model does not permit ligand exchange and needs to be reparametrized when it comes to a new metal center. For the nonbonded dummy model, the dummy atoms are covalently bound to the metal center in a predefined geometry; this model avoids the formation of covalent bonds with ligands, thereby permitting both ligand exchange and coordination number changes. The nonbonded point charge model is particularly favored in MD simulations due to its simplicity and computational efficiency. It treats the metal ion as a point charge that interacts with other atoms via electrostatic and van der Waals (vdW) interactions, avoiding the complexity of explicitly modeling chemical bonds. The vdW interactions are often modeled by the 12–6 Lennard–Jones (LJ) potential (Equation (1)).
U L J r i j = ε i j ( R m i n , i j r i j 12 R m i n , i j r i j 6 )
where r i j is the distance between atoms i and j , while ε i j and R m i n , i j are the potential well depth and the distance at which the potential energy between atoms i and j reaches the minimum, respectively. The parameters of R m i n , i j and ε i j commonly satisfy the Lorentz–Berthelot (LB) mixing rule ( R i j = ( R i i + R j j ) / 2 ; ε i j = ε i i ε j j ). Only the LJ parameters of R and ε need to be optimized for this model. Therefore, the simplicity of this approach allows efficient simulation of complex systems, making it popular in force-field-based predictions.
Many efforts have been dedicated to developing the nonbonded point charge models for divalent metal ions, with the aim of accurately capturing the essential ion properties and interactions with the environment. Selecting appropriate target properties for guiding the development of metal ion models is of vital importance. Structural characteristics such as the ion–water oxygen distance (IOD) and coordination number (CN), as well as dynamics and thermodynamic properties such as the ion–water exchange rate, diffusion coefficient, interaction energy (IE), and hydration free energy (HFE), are often used as the benchmarks. These target properties ensure that the developed models are able to reproduce the corresponding experimental data accurately. Åqvist [26] was the first to develop 12–6 LJ parameters for alkaline earth metal ions (such as Mg2+ and Ca2+), which are capable of accurately reproducing the experimentally observed HFE. Babu and Lim optimized the vdW parameters for 24 divalent metal cations with known experimental HFE [27]. Villa et al. parameterized the LJ parameters for Mg2+ by fitting the activation energy in the ion–water system, and the designed set accurately reproduced the experimental exchange rate and improved the kinetic description of Mg2+ and phosphate ions [28]. Merz and coworkers reported three sets of 12–6 LJ models for 24 divalent metal ions, comprising one set targeting HFE, one targeting IOD, and a compromise (CM) set targeting a balance of both properties [22]. Our recent work has developed novel 12–6 LJ models by increasing the well depth of the LJ potential; these models were able to accurately reproduce both HFE and IOD in the first hydration shell for 16 divalent metal cations simultaneously [25]. Other methods were also reported with a double exponential function [29] or 12–6–4 LJ [23,24,30,31,32] for modeling vdW interactions as well as the use of scaled ionic charges in conjunction with the 12–6 LJ potential [33,34,35,36].
Current research mainly relies on the choice of existing experimental values as targets to optimize ion models. However, some of the experimental observations are still debated at present, such as the real single-ion hydration free energies (HFEs) [37]. Quantum mechanics (QM)-based calculation, such as density functional theory (DFT) [38,39], provides a rigorous description of the electronic and energetic properties of molecular systems, enabling the computation of accurate interaction energies and geometric parameters. There are few reports in the literature on the use of QM-derived data as a target for the design of force field models. Herein, we developed force field models for divalent metal cations compatible with the widely used TIP3P water model [40], specifically focusing on Mg2+ and Ca2+. A systematic comparison with the corresponding models in the Amber force field [41] was presented in terms of the HFE, IOD, CN, and ion–water interactions. We also developed vdW parameters for the Cl ion to evaluate the performance of the cation models in simulating the density and viscosity of salt solutions. The metalloprotein form of the phosphohydrolase OxsA (PDB ID: 5TKA) [42], was simulated to demonstrate the capability of our models in biomolecular simulations. In this work, we aim to show the feasibility of using gas-phase QM data as a target to develop ion models for aqueous simulations. We show this by using the gas-phase DFT calculations to tune the LJ parameters and applying the designed models to MD simulations in the liquid phase. It would be valuable for the further design of ion models and computational simulation of metal-containing systems.

2. Computational Methods

2.1. QM Calculations of Ion–Water Dimers

The geometries of ion–water dimers, namely, Mg2+–H2O, Ca2+–H2O, and Cl–H2O, were optimized using the popular B3LYP−D3(BJ) exchange-correlation functional [43,44] with the def2-TZVP basis set [45]. The C2v symmetry was considered for all three geometries of ion–water dimers, and for Cl–H2O, the Cs symmetry was also examined, as tested in the work by Joung and Cheatham [46]. A cluster model of ion–water dimers was used in gas phase (no periodic boundary condition) during the optimization.
Interaction energies (IEs) between ions and water were computed using the ωB97X-V density functional (ω = 0.3) [47] in combination with the large def2-QZVPPD basis set [48]. This level of theory showed good accuracy in IE predictions for these ion–water dimers compared to the highly accurate CCSD(T)/CBS references [49]. The IE values and the distances between ions and water oxygen (IODs) computed using the optimized geometries of ion–water dimers were used as targets to design the force field ion models as stated in the following.
We also performed a scan of ion–water interactions as a function of the IOD at the ωB97X-V/def2-QZVPPD level of theory. The optimized geometries of ion–water dimers in the C2v symmetry were rotated to make the vector connecting the ion and water oxygen parallel to the x-axis, and the water molecule was then moved toward or away from the ion along the x-axis in steps of 0.1 Å. A step size of 0.01 Å was used for the IOD values within 1 Å of the QM-optimized optimal distance. For the Cl–H2O dimer in the Cs symmetry, similarly, the vector connecting the ion and the neighboring hydrogen of water was chosen for the rotation. The geometry operations on the ion–water dimers were performed using the Multiwfn software (version: 3.8) [50] for the scan running in a batch mode. All of the QM calculations were carried out with the Beijing Density Functional (BDF) package (version: 2024A) [51,52,53,54].

2.2. LJ Parameter Scanning

The nonbonded point charge model was used to model Mg2+, Ca2+, and Cl ions, and the 12–6 LJ potential (Equation (1)) was adopted to describe vdW interactions. Thus, two LJ parameters (R and ε) could be adjusted to reproduce specific properties (i.e., target values). For Mg2+ and Ca2+ ions, we selected R values in a range from 0.1 to 2 Å with a step size of 0.1 Å (20 values in total). The ε values ranged from 0.01 to 3168.2 kcal/mol; for convenience, a form of –lg(ε) was used to represent ε. The corresponding –lgε values ranged from –3.5 to 2 kcal/mol, namely, –3.5, –3.25, –3, –2.5, –2, –1.5, –1, –0.5, 0, 0.5, 1, and 2 (12 values in total). A total of 20 × 12 = 240 LJ combinations were utilized for parameter scanning of the cations. For the chloride ion (Cl), the R values varied from 1.5 to 3.0 Å in increments of 0.1 Å, and the –lgε ranged from –2 to 4 kcal/mol in steps of 0.5, thereby resulting in 15 × 13 = 195 different Rε combinations.

2.3. Parameter Generation

Similar to the interaction energy (IE) scan mentioned in Section 2.1, we computed ion–water IE profiles as a function of IOD using a given Rε combination. The IE was a sum of the van der Waals (Equation (1)) and Coulomb interactions ( 1 4 π ϵ 0 Q i Q j r i j ). The IOD value was changed from 0.2 to 10 Å with a step size of 0.01 Å. For each Rε, a pair of IE and IOD corresponding to the minimum of the computed IE profile was recorded. The discrete points of IE and IOD obtained for different Rε combinations were used to construct the surface of IE and IOD as a function of R and ε. Then, the contour lines corresponding to the targeted IE and IOD values (i.e., QM predictions in gas phase in Section 2.1) were plotted on a 2D graph with −lg(ε) as the x-axis and R as the y-axis. The intersection point (if any) of the two contour lines was the designed LJ parameters, which would be able to simultaneously reproduce both target properties. An in-house Python script was used for performing 2D interpolation and generating contour plots [19,25]. We also tried the use of gas-phase IE and liquid-phase IOD (IODliq) for the development of ion models. The IOD values in aqueous solution were taken from our previous work [25]. These designed LJ parameters were validated by the following molecular dynamics (MD) simulations of salt solution and protein systems.

2.4. General MD Simulation Setup

All MD simulations were performed with GROMACS software (version 2024.1) [55]. The particle–mesh Ewald (PME) method [56,57] was used to handle electrostatic interactions with a cutoff distance of 10 Å, and the vdW interactions were truncated with the same cutoff. The temperature was maintained at 298.15 K through the velocity-rescaling algorithm [58] with a coupling time constant of 1 ps. For pressure coupling, the Parrinello–Rahman method [59,60] was adopted, utilizing a time constant of 5 ps. To mimic infinite solutions, periodic boundary conditions (PBCs) were applied in all directions during MD simulations. The SETTLE algorithm [61] was employed to maintain the structural rigidity of the TIP3P water model [40]. The chemical bonds related to hydrogen atoms were constrained using the LINCS algorithm [62], allowing a time step of 2 fs. As a preliminary step, an energy minimization was conducted for all MD simulations to avoid bad contacts. Subsequently, two simulations were performed: an initial 100 ps NVT equilibration, followed by a 400 ps NPT equilibration. Finally, production simulations were implemented under NPT conditions for data acquisition and further analysis.

2.5. Hydration Free Energy Calculation

Following our previous protocol [21,25,63,64], thermodynamic integration (TI) was used to compute the hydration free energy (HFE) of the ions tested. This approach utilized an alchemical reaction coordinate, λ, ranging from 0 to 1 in 25 discrete steps (λ = 0, 0.02, 0.04, 0.07, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5, 0.55, 0.6, 0.65, 0.7, 0.75, 0.8, 0.85, 0.9, 0.93, 0.96, 0.98, and 1) for decoupling electrostatic interactions [65]. Subsequently, a set of 16 λ values (λ = 0. 0.05, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.65, 0.7, 0.75, 0.8, 0.85, 0.9, 0.95, and 1) was employed to decouple LJ interactions [66], ensuring a more accurate treatment of nonbonded interactions. Each λ value was subjected to 500 ps MD simulations at a temperature of 298.15 K and a pressure of 1 bar. Each simulation box contained an ion and 600 water molecules. The initial 100 ps was discarded to ensure equilibration of the system, and the rest was used for data analysis.
The HFEs of ions were calculated from the TI statistics with the Bennett acceptance ratio (BAR) method [67]. The computed HFE derived from typical PME/PBC simulations only considered the bulk interaction between ions and their local media, thereby constituting an intrinsic free energy component ( G i n t r ). To obtain the physically real free energy of ionic hydration (Equation (2)), it was necessary to add a surface contribution ( G s u r f ) that accounted for the ion traversing the vacuum/bulk interface [66,68,69,70].
G , r e a l = G i n t r + G s u r f = G i n t r + q G
Here, q represented the net charge of ions (+2 e for Mg2+ and Ca2+, –1 e for Cl), and G denoted the solvent Galvani potential (–12.3 kcal/mol/e for TIP3P water model) [71].

2.6. IOD and CN Calculations

Production simulations of a single ion immersed in 600 water molecules were run for 4 ns, and the trajectories were recorded every 50 simulation steps (0.1 ps). As a result, 40,000 frames in total were utilized for calculating the radial distribution functions (RDFs) of the oxygen atoms (OW) of water molecules in the vicinity of the ions. This analysis enabled the derivation of the ion–water oxygen distance (IOD) and coordination number (CN) in the first hydration shell of the ion. The built-in “gmx rdf” tool in the GROMACS software was employed to carry out this task with a bin width of 0.01 Å [55].

2.7. Model Evaluation in Salt Solution

The developed cation and anion models were employed to investigate the properties of salt solutions (MgCl2 and CaCl2). The densities and viscosities of salt solutions were simulated in different systems comprising 555 water molecules and a certain number of ions at desired salt concentrations. For instance, to achieve a concentration of 1 molarity (m) (defined as the number of moles of salt dissolved per kilogram of water) [36], 10 MgCl2 (i.e., 10 Mg2+ and 20 Cl) were added to the simulated system with a simulation box of 26 × 26 × 26 Å3. A series of simulations were conducted for salt solutions with a variety of concentrations from 0.1 to 5 m (namely, 0.1, 0.5, 1, 2, 3, 4, and 5 m). To calculate the densities of the salt solutions, the production simulations were extended to 50 ns under NPT conditions, and the trajectories in the last 40 ns were utilized for data collection and analysis. For viscosity analysis, a 5 ns simulation was carried out at NPT to equilibrate the system, followed by a non-equilibrium NVT simulation of 10 ns using an acceleration amplitude value of 0.05 nm·ps−2 [72]. Each system was repeated three times using different initial velocities.

2.8. Model Evaluation in a Metalloprotein

The metalloprotein of the HD-domain phosphohydrolase OxsA (PDB ID: 5TKA) [42] was selected as an example to evaluate the usability of the designed ion models. This protein has a Mg2+ binding site where one Mg2+ ion is coordinated with four amino acid residues (His31, His66, Asp67, and Asp132) and two water molecules, as presented in Figure 1. The figures of 3D structures in this work were made in PyMOL (https://www.pymol.org; version 2.5.4; accessed on 24 July 2024). Each residue contributes a single atom to the coordination sphere of metal center (that is, in a monodentate binding mode). The default Mg2+ model in the Amber force fields (adapted from Åqvist’s model [26]) was also used for comparison. The Amber ff14SB force field [73] and TIP3P water model [40] were used to model protein and water molecules, respectively. Crystal water molecules located within 5 Å of Mg2+ ion were preserved. The simulation box was ca. 70 × 70 × 70 Å3 and contained one protein, 3 Mg2+, and 9345 water molecules. During the equilibration stages for the NVT and NPT ensembles, position restraints were imposed on the protein backbone atoms using a harmonic potential with a force constant of 1000 kJ/(mol nm2). Production simulations were run for 50 ns under NPT conditions without any position restraints, enabling a comprehensive exploration of Mg2+ dynamics and interactions with the protein residues.

3. Results and Discussion

3.1. Targeted Properties from Gas-Phase QM Calculations

Optimized geometries of ion–water dimers at the B3LYP-D3(BJ)/def2-TZVP level of theory in the C2v and Cs symmetries are given in Figure 2. During optimization, the water structure was constrained to its experimental geometry in the gas phase with an O–H bond length of 0.9572 Å and a H–O–H angle of 104.52 degrees [74], as depicted in Figure 2a; this geometry was the same as the TIP3P water model [40]. The ion–water oxygen distances (IODs) for Mg2+, Ca2+, and Cl in the C2v symmetry amounted to 1.92 (Figure 2a), 2.23 (Figure 2b), and 3.16 Å (Figure 2c), in good agreement with the previously reported values of 1.91, 2.22 [49], and 3.20 Å [75], respectively. The IOD for the Cl–H2O dimer in the Cs symmetry was 3.12 Å, and the distance between Cl and one hydrogen atom of water was 2.18 Å (Figure 2d). Both distances agreed well with the highly accurate CCSD(T)/CBS reference values of 3.12 and 2.15 Å [49,76].
The interaction energies (IEs) computed at the ωB97X-V/def2-QZVPPD level of theory for Mg2+, Ca2+, and Cl in the C2v symmetry were −82.14, −56.67, and −11.76 kcal/mol, respectively. For the Cs symmetry, the IE between Cl and a single water amounted to −13.76 kcal/mol, representing stronger binding than that in the C2v symmetry. This indicated a more thermodynamically stable configuration of the Cl–H2O dimer in the Cs symmetry than in the C2v case. Our computed IE values were in line with previous work [49,76] with errors ranging from 0.3 (for Mg2+ in the C2v symmetry) to 1.1 (for Cl in the Cs symmetry) kcal/mol. The differences in the IOD and IE were mainly ascribable to the water geometry used and/or the level of theory. Table 1 summarizes the IE and IOD values obtained from our QM-optimized geometries of ion–water dimers in the gas phase, and these predictions were then used as target properties for guiding the development of ion models.
We also attempted to design ion models targeting the experimental IOD in aqueous solution (liquid phase, IODliq), as given in Table 1. Experimental HFE and CN values in Table 1 were listed as well and were used to evaluate the model performance in predicting the ion properties that were not targeted in the model design. Note that G o , c o n v represented the conventional HFE (experimentally measurable) taken from Ref. [77], and it was possible to use Equation (3) to obtain the real HFEs of ions.
G * , r e a l ( i o n ) = G o , c o n v ( i o n ) + q G o H + G o *
where q was the net charge of ions, G o H + denoted the real HFE of the proton H+, and G o * represented a standard-state correction of 1.9 kcal/mol. The reported experimental HFEs have been highly debated due to the choice of different G o H + , and we chose −259.8 kcal/mol in this work, as recommended by Hünenberger and Reif [78].

3.2. Designed Model Parameters

As stated in Section 2.3, the intersection (if existent) of the contour lines of IE and IOD corresponding to the target values was the optimal LJ parameter for a given ion. Two sets of models were attempted: one (model I) targeted the gas-phase IE and the gas-phase IOD (IODgas), and the other (model II) targeted the gas-phase IE and the liquid-phase IOD (IODliq). Contour plots depicting the two models for Mg2+, Ca2+, and Cl are presented in Figure 3 (panels a–c for model I and panels d−f for model II). Intersection points existed for Mg2+ and Ca2+ ions (Figure 3a,d). For the Cl ion in the C2v symmetry, however, there was no intersection between the IE and IOD contours, and the two contours seemed to be parallel for an increasing R and a decreasing ε, as indicated by Figure 3b,e. When considering the Cs symmetry, the intersections existed, and we were able to obtain the optimal Rε parameters for the Cl anion (Figure 3c,f).
The final LJ Rε parameters of the tested ions for use with the TIP3P water model are listed in Table 2 along with the calculated IE, IODgas, and HFE and the IODliq and CN values in the first hydration shell of the ions. The force field parameters derived from the target properties of gas-phase IE and IODgas are denoted as model I, while model II represents the ion models targeting gas-phase IE and IODliq. For comparison, Mg2+ and Ca2+ models (adapted from Åqvist’s work [26]) in the Amber force fields were tested as well; the Cl model by Joung and Cheatham was chosen, denoted as Cl(JC), and this model was tuned for a balanced performance of crystal and solution properties [46].
For the metal cations, our models appeared to have a smaller R and a larger ε than Åqvist’s models (Table 2). The designed models accurately reproduced the corresponding target properties of IE, IODgas (model I), and IODliq (model II), whereas they gave considerable overestimates (i.e., more negative values) of the ionic HFE by 50−100 kcal/mol (Table 2). Mg(II) and Ca(II) yielded IODliq values of 2.09 and 2.46 Å (Table 2 and Figure 4), respectively, in good agreement with the experimental observations (Table 1); however, both models gave a larger coordination number (CN) of water molecules in the first hydration shell of the ions than the experiments, probably due to the overestimated HFEs. When targeting the gas-phase IOD, model I for Mg2+ and Ca2+ tended to yield a relatively small IOD in aqueous solution (i.e., IODliq), despite its good performance in CN prediction (Figure 4). As expected, Mg(Åqvist) targeted the HFE and gave a value close to the experiment; however, Ca(Åqvist) yielded an HFE value of −331.98 kcal/mol, underestimating the experimental value by 46 kcal/mol. It should be noted that in the original publication [26], the Mg2+ and Ca2+ models showed a reasonable predictions for both HFE and IOD. After adaptation for compatibility with the Amber force fields, the Ca2+ model appeared not to work well.
For the Cl anion, R and ε values were on the same order of magnitude as the JC model (Table 2). Surprisingly, although model Cl(I) did not target the liquid-phase IOD and HFE, it reproduced both properties accurately (Table 2). Cl(JC) appeared to give a relatively small IODliq of 3.12 Å compared to the experimental observation (3.187 Å), as given in Table 2 and Figure 4. Comparisons of the predicted properties with gas-phase QM calculations and liquid-phase MD simulations are given in Tables S1 and S2, respectively, in the Supplementary Materials.

3.3. Gas-Phase IE Comparison

The interaction energies (IEs) of ion−water dimers and the IOD in the gas phase (IODgas) were well reproduced when both properties were targeted in the design of ion models (i.e., model I in Table 2), as indicated by the minimum of the IE profiles (red dashed lines) in Figure 5. In the force field models, the IE between the ion and a single TIP3P water for different ion models was composed of the Coulomb and vdW contributions; the corresponding QM predictions were performed at the ωB97X-V/def2-QZVPPD level of theory (Section 2.1). Model II targeted the liquid-phase IOD (IODliq) and tended to give a large IODgas for the metal cations (blue dashed lines in Figure 5a,b), whereas Cl(II) displayed good performance in predicting the distance between Cl and a single water with the C2v and Cs symmetries in the gas phase (Table 2 and green lines in Figure 5c,d). Note that we only used the Cs configuration of the Cl−H2O dimer for the model development (Figure 5d). Although the resulting Cl model showed a good performance for IODgas in the C2v symmetry, it greatly overestimated the interaction energies compared to the QM predictions (Table 2 and Figure 5c).
Targeting IE and IOD in the gas phase, model I appeared to give much narrower IE profiles for the metal cations than the QM predictions, as indicated by the red dashed lines and solid lines in Figure 5a,b, respectively. For the Cl ion, however, the IE profiles from force field and QM predictions were almost identical (Figure 5d). This finding indicated that the use of a simple 12–6 LJ potential appeared accurate enough to model Cl−water interactions, whereas a more complex form might be necessary for accurate modeling of divalent metal cations.

3.4. Density and Viscosity of Salt Solutions

Salt solutions of MgCl2 and CaCl2 in different concentrations (0−5 m) were used to evaluate the performance of the designed models and to compare with the existing models of M2+ by Åqvist [26] and Cl by Joung and Cheatham [46]. The densities of MgCl2 and CaCl2 solutions at a temperature of 298.15 K and a pressure of 0.1 MPa are presented in Figure 6. At low concentrations (<1 m), the densities of both salt solutions were well reproduced. For high concentrations of >1 m, our designed Mg2+ models overestimated the solution densities by less than 3%, whereas the Mg(Åqvist) model gave slight underestimates (Figure 6a). The Mg(I) model targeting the gas-phase IOD showed the best performance among the tested three models. The Ca(I) and Ca(II) models appeared to yield considerable overestimates of the densities of salt solutions at high concentrations, while the Ca(Åqvist) showed good performance (Figure 6b).
Similar to the performances for the densities of salt solutions, the three models tested were able to reproduce the solution viscosities with good accuracy at low concentrations of up to 3 m (Figure 7). For high concentrations, the viscosities of MgCl2 and CaCl2 salt solutions were significantly overestimated, which might be due to the overestimated densities of these solutions for our designed models (I and II). Although Åqvist’s model yielded a lower density than the experiment (Figure 6), all three models still greatly overestimated the viscosities at high concentrations (Figure 7). These finding indicated a problem in the modeling of ion−ion interactions, in particular for the solutions at high salt concentrations.

3.5. Metalloprotein Simulation

To evaluate the performance of the designed models in complex systems such as proteins, the metalloprotein of phosphohydrolase OxsA (PDB ID: 5TKA) [42] was simulated for 50 ns. Three Mg2+ models were utilized to model the metal center, namely, Mg(I), Mg(II), and Mg(Åqvist). Figure 8 presents the root-mean-square deviations (RMSDs) of the protein backbone and metal-binding residues from the crystal structure as a function of simulation time. The protein structures tended to reach equilibrated states after 30 ns simulations, and the trajectories from the last 20 ns were thus used for data analysis. The RMSDs of the protein backbone amounted to 4.1, 3.1, and 4.4 Å for the Mg(I), Mg(II), and Mg(Åqvist) models, respectively. The corresponding RMSDs of the Mg2+-binding residues were 0.5, 0.9, and 0.7 Å, respectively.
Notably, a significant increase in the RMSD of metal-binding residues was observed for the Mg(II) model during the early stage of MD simulations, reflecting a large change in the structure of the Mg2+ coordination environment. In the crystal structure, only one oxygen atom of each Asp side chain (OD2 in Asp67; OD1 in Asp132) was coordinated with the Mg2+ ion (Figure 1), forming a monodentate binding mode. During MD simulations, the Mg(I) and Mg(Åqvist) models supported the monodentate mode, and they also coordinated two histidine residues (His31 and His66) and two water models, giving rise to a coordination number (CN) of six, as indicated by Figure 9a,c, respectively. However, both oxygen atoms (OD1 and OD2) of the Asp side chains participated in the metal binding in the Mg(II) model, forming a bidentate binding mode (Figure 9b). As a result, only one water molecule was able to coordinate in the Mg(II) model, and the resulting CN of the Mg2+ ion was changed to seven (Figure 9b). This issue might be related to the identical charges of OD1 and OD2 atoms in the Asp residue of the Amber ff14SB force field and could be solved by modifying the atomic charges of the oxygen atoms in Asp [17], using the nonbonded dummy atom model of metal ions [21], or choosing appropriate ion models such as the point charge models presented in Åqvist’s work [26], our previous work [83], and this work.
The interatomic distances between Mg2+ and metal-binding residues were computed using the last 20 ns of the MD simulations and listed in Table 3 along with the crystal distances. In the crystal structure, Mg2+ bound to two histidine residues (His31 and His66), two aspartic residues (Asp67 and Asp132), and two water molecules, resulting in a coordination number (CN) of six (Figure 1). The Mg(I) and Mg(Åqvist) models were able to maintain the crystal metal-binding networks, while they tended to overestimate the distances between Mg2+ and amino acids (His and Asp) by up to 0.2 Å and to underestimate the Mg2+−water distances by ca. 0.2−0.3 Å (Table 3). These two models showed similar performances in predicting the distances between Mg2+ and oxygen atoms in Asp and water, and Mg(Åqvist) displayed better performance for the distance between Mg2+ and nitrogen atoms in His than the Mg(I) model, as indicated by Table 3.

4. Concluding Remarks

The development of ion models has depended on the availability of reliable experimental observations. Thermodynamic properties have often been used as key target properties (for instance, hydration free energy, HFE). However, the ionic HFE is inherently tied to the HFE of the proton G H + , yet there is no consensus on this quantity in the literature. Experimental observations for ionic HFEs originated mostly from three sources: Marcus [77], Noyes [84], and Rosseinsky [85]. The former one adopted a G H + value of −252.3 kcal/mol, whereas the latter two were based on a value of −260.5 kcal/mol. This resulted in discrepancies in the referenced HFEs [18,22,23,27,29,66], thereby influencing the calibration and validation of ion models and introducing inconsistencies between different sets of ion models.
In this work, we showed that the gas-phase quantum mechanics (QM) calculations were able to guide the design of classical force field models for ions and that the designed models can be used for MD simulations in liquid phase. We attempted to target QM-based ion−water dimer interactions in the gas phase; more precisely, the interaction energies (IEs) between the ion and a single water molecule and the distance between the ion and water oxygen (IODgas) served as target properties. The resulting model was denoted as model I. This strategy considered the energetics of ion binding and had the advantage of avoiding the use of debated HFE references. Surprisingly, the designed Cl model, Cl(I), showed good performance for some non-targeted properties in aqueous solutions, such as HFE, IOD, and CN. Our Mg2+ model, Mg(I), yielded reasonable IOD and CN in the liquid phase as well, and it appeared to work well for simulating metalloproteins in water. However, Mg(I) considerably overestimated the HFE, which might be ascribed to the fact that the polarization effects in bulk media [86,87,88,89] were not fully captured in the model design. The designed models were able to effectively reproduce the density and viscosity of salt solutions at low concentrations, but a large degree of overestimation was observed for high concentrations. This indicated the necessity of re-optimization of ion−ion interactions, which might be related to the polarization effects as well. Addressing this issue will require incorporating complicated polarization models [90,91,92,93] or using scaled ionic charges [33,34,35,36,94,95].
Targeting gas-phase IE and liquid-phase IOD was also attempted. Despite good performance for the IOD (IODliq) in the aqueous phase, the designed model (model II) showed a similar or even worse performance for the tested properties in solutions and they were not able to maintain the crystal structure of the metal-binding center (Mg2+) during MD simulations of the metalloprotein tested. Therefore, this strategy (i.e., model II) is not recommended for the development of ion models on the basis of our tests. This work has implications for the further design of ion models and will be valuable for MD simulations of ion-containing systems.

Supplementary Materials

The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/met14091009/s1, Table S1: Comparison of each calculated interaction energy (IE, kcal/mol) and ion–water oxygen distance (IOD, Å) in the gas phase with the corresponding DFT reference value (ref) for the designed ion models; Table S2: Comparison of each calculated interaction energy (IE, kcal/mol) and ion–water oxygen distance (IOD, Å) in the gas phase with the corresponding experimental observation (exp.) for the designed ion models.

Author Contributions

Conceptualization, Y.Z. and H.Z.; Data Curation, Y.Z., B.W. and C.L.; Formal Analysis, Y.Z., B.W. and C.L.; Investigation, Y.Z.; Methodology, Y.Z. and H.Z.; Software, H.Z.; Supervision, H.Z.; Validation, Y.Z., B.W. and C.L.; Visualization, Y.Z.; Writing—Original Draft, Y.Z.; Writing—Review and Editing, H.Z. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the National Natural Science Foundation of China (grant: 21977013 and 21606016).

Data Availability Statement

The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author.

Acknowledgments

We thank Yang Gao, Rongding Lei, and Shaoqin Zhang from HZWTECH for technical assistance regarding QM calculations with the BDF software.

Conflicts of Interest

The authors declare no conflicts of interest.

References

  1. Al Alawi, A.M.; Majoni, S.W.; Falhammar, H. Magnesium and Human Health: Perspectives and Research Directions. Int. J. Endocrinol. 2018, 2018, 9041694. [Google Scholar] [CrossRef]
  2. Berta, D.; Buigues, P.J.; Badaoui, M.; Rosta, E. Cations in Motion: QM/MM Studies of the Dynamic and Electrostatic Roles of H+ and Mg2+ Ions in Enzyme Reactions. Curr. Opin. Struct. Biol. 2020, 61, 198–206. [Google Scholar] [CrossRef] [PubMed]
  3. Eshra, A.; Schmidt, H.; Eilers, J.; Hallermann, S. Calcium Dependence of Neurotransmitter Release at a High Fidelity Synapse. eLife 2021, 10, e70408. [Google Scholar] [CrossRef]
  4. Song, Z.; Wang, Y.; Zhang, F.; Yao, F.; Sun, C. Calcium Signaling Pathways: Key Pathways in the Regulation of Obesity. Int. J. Mol. Sci. 2019, 20, 2768. [Google Scholar] [CrossRef] [PubMed]
  5. Oh, B.-C. Phosphoinositides and Intracellular Calcium Signaling: Novel Insights into Phosphoinositides and Calcium Coupling as Negative Regulators of Cellular Signaling. Exp. Mol. Med. 2023, 55, 1702–1712. [Google Scholar] [CrossRef] [PubMed]
  6. Rall, J.A. Discovery of the Regulatory Role of Calcium Ion in Muscle Contraction and Relaxation: Setsuro Ebashi and the International Emergence of Japanese Muscle Research. Adv. Physiol. Educ. 2022, 46, 481–490. [Google Scholar] [CrossRef] [PubMed]
  7. Morrell, A.E.; Brown, G.N.; Robinson, S.T.; Sattler, R.L.; Baik, A.D.; Zhen, G.; Cao, X.; Bonewald, L.F.; Jin, W.; Kam, L.C.; et al. Mechanically Induced Ca2+ Oscillations in Osteocytes Release Extracellular Vesicles and Enhance Bone Formation. Bone Res. 2018, 6, 6. [Google Scholar] [CrossRef]
  8. Qiao, W.; Pan, D.; Zheng, Y.; Wu, S.; Liu, X.; Chen, Z.; Wan, M.; Feng, S.; Cheung, K.M.C.; Yeung, K.W.K.; et al. Divalent Metal Cations Stimulate Skeleton Interoception for New Bone Formation in Mouse Injury Models. Nat. Commun. 2022, 13, 535. [Google Scholar] [CrossRef]
  9. Meenachi, P.; Subashini, R.; Lakshminarayanan, A.K.; Gupta, M. Comparative Study of the Biocompatibility and Corrosion Behaviour of Pure Mg,Mg Ni/Ti, and Mg 0.4Ce/ZnO2 Nanocomposites for Orthopaedic Implant Applications. Mater. Res. Express 2023, 10, 056503. [Google Scholar]
  10. Shan, Z.; Xie, X.; Wu, X.; Zhuang, S.; Zhang, C. Development of Degradable Magnesium-Based Metal Implants and Their Function in Promoting Bone Metabolism (a Review). J. Orthop. Transl. 2022, 36, 184–193. [Google Scholar] [CrossRef]
  11. Zhao, C.; Tan, J.; Li, W.; Tong, K.; Xu, J.; Sun, D. Ca2+ Ion Responsive Pickering Emulsions Stabilized by Pssma Nanoaggregates. Langmuir 2013, 29, 14421–14428. [Google Scholar] [CrossRef] [PubMed]
  12. Li, N.; Cao, Y.; Liu, J.; Zou, W.; Chen, M.; Cao, H.; Deng, S.; Liang, J.; Yuan, T.; Wang, Q.; et al. Microenvironment-Responsive Release of Mg2+ from Tannic Acid Decorated and Multilevel Crosslinked Hydrogels Accelerates Infected Wound Healing. J. Mater. Chem. B 2024, 12, 6856–6873. [Google Scholar] [CrossRef]
  13. Bril, M.; Fredrich, S.; Kurniawan, N.A. Stimuli-Responsive Materials: A Smart Way to Study Dynamic Cell Responses. Smart. Mater. Med. 2022, 3, 257–273. [Google Scholar] [CrossRef]
  14. Dauber-Osguthorpe, P.; Hagler, A.T. Biomolecular Force Fields: Where Have We Been, Where Are We Now, Where Do We Need to Go and How Do We Get There? J. Comput. Aided Mol. Des. 2019, 33, 133–203. [Google Scholar] [CrossRef] [PubMed]
  15. Hu, L.H.; Ryde, U. Comparison of Methods to Obtain Force-Field Parameters for Metal Sites. J. Chem. Theory Comput. 2011, 7, 2452–2463. [Google Scholar] [CrossRef]
  16. Li, P.; Merz, K.M., Jr. Mcpb.Py: A Python Based Metal Center Parameter Builder. J. Chem. Inf. Model. 2016, 56, 599–604. [Google Scholar] [CrossRef] [PubMed]
  17. Duarte, F.; Bauer, P.; Barrozo, A.; Amrein, B.A.; Purg, M.; Åqvist, J.; Kamerlin, S.C.L. Force Field Independent Metal Parameters Using a Nonbonded Dummy Model. J. Phys. Chem. B 2014, 118, 4351–4362. [Google Scholar] [CrossRef]
  18. Jiang, Y.; Zhang, H.; Feng, W.; Tan, T. Refined Dummy Atom Model of Mg2+ by Simple Parameter Screening Strategy with Revised Experimental Solvation Free Energy. J. Chem. Inf. Model. 2015, 55, 2575–2586. [Google Scholar] [CrossRef]
  19. Jiang, Y.; Zhang, H.; Tan, T. Rational Design of Methodology-Independent Metal Parameters Using a Nonbonded Dummy Model. J. Chem. Theory Comput. 2016, 12, 3250–3260. [Google Scholar] [CrossRef]
  20. Liao, Q.; Pabis, A.; Strodel, B.; Kamerlin, S.C.L. Extending the Nonbonded Cationic Dummy Model to Account for Ion-Induced Dipole Interactions. J. Phys. Chem. Lett. 2017, 8, 5408–5414. [Google Scholar] [CrossRef]
  21. Peng, J.; Zhang, Y.; Jiang, Y.; Zhang, H. Developing and Assessing Nonbonded Dummy Models of Magnesium Ion with Different Hydration Free Energy References. J. Chem. Inf. Model. 2021, 61, 2981–2997. [Google Scholar] [CrossRef]
  22. Li, P.; Roberts, B.P.; Chakravorty, D.K.; Merz, K.M., Jr. Rational Design of Particle Mesh Ewald Compatible Lennard-Jones Parameters for +2 Metal Cations in Explicit Solvent. J. Chem. Theory Comput. 2013, 9, 2733–2748. [Google Scholar] [CrossRef]
  23. Li, P.; Merz, K.M., Jr. Taking into Account the Ion-Induced Dipole Interaction in the Nonbonded Model of Ions. J. Chem. Theory Comput. 2014, 10, 289–297. [Google Scholar] [CrossRef] [PubMed]
  24. Li, Z.; Song, L.F.; Li, P.; Merz, K.M. Systematic Parametrization of Divalent Metal Ions for the OPC3, OPC, TIP3P-FB, and TIP4P-FB Water Models. J. Chem. Theory Comput. 2020, 16, 4429–4442. [Google Scholar] [CrossRef]
  25. Zhang, Y.; Jiang, Y.; Peng, J.; Zhang, H. Rational Design of Nonbonded Point Charge Models for Divalent Metal Cations with Lennard-Jones 12-6 Potential. J. Chem. Inf. Model. 2021, 61, 4031–4044. [Google Scholar] [CrossRef] [PubMed]
  26. Åqvist, J. Ion-Water Interaction Potentials Derived from Free Energy Perturbation Simulations. J. Phys. Chem. 1990, 94, 8021–8024. [Google Scholar] [CrossRef]
  27. Babu, C.S.; Lim, C. Empirical Force Fields for Biologically Active Divalent Metal Cations in Water. J. Phys. Chem. A 2006, 110, 691–699. [Google Scholar] [CrossRef]
  28. Allnér, O.; Nilsson, L.; Villa, A. Magnesium Ion-Water Coordination and Exchange in Biomolecular Simulations. J. Chem. Theory Comput. 2012, 8, 1493–1502. [Google Scholar] [CrossRef] [PubMed]
  29. Man, V.H.; Wu, X.; He, X.; Xie, X.Q.; Brooks, B.R.; Wang, J. Determination of Van Der Waals Parameters Using a Double Exponential Potential for Nonbonded Divalent Metal Cations in TIP3P Solvent. J. Chem. Theory Comput. 2021, 17, 1086–1097. [Google Scholar] [CrossRef]
  30. Jafari, M.; Li, Z.; Song, L.F.; Sagresti, L.; Brancato, G.; Merz, K.M., Jr. Thermodynamics of Metal–Acetate Interactions. J. Phys. Chem. B 2024, 128, 684–697. [Google Scholar] [CrossRef]
  31. Li, Z.; Song, L.F.; Sharma, G.; Koca Fındık, B.; Merz, K.M., Jr. Accurate Metal–Imidazole Interactions. J. Chem. Theory Comput. 2023, 19, 619–625. [Google Scholar] [CrossRef] [PubMed]
  32. Koca Fındık, B.; Jafari, M.; Song, L.F.; Li, Z.; Aviyente, V.; Merz, K.M., Jr. Binding of Phosphate Species to Ca2+ and Mg2+ in Aqueous Solution. J. Chem. Theory Comput. 2024, 20, 4298–4307. [Google Scholar] [CrossRef]
  33. Nikitin, A.; Del Frate, G. Development of Nonbonded Models for Metal Cations Using the Electronic Continuum Correction. J. Comput. Chem. 2019, 40, 2464–2472. [Google Scholar] [CrossRef]
  34. Zeron, I.M.; Abascal, J.L.F.; Vega, C. A Force Field of Li+, Na+, K+, Mg2+, Ca2+, Cl, and SO42− in Aqueous Solution Based on the TIP4P/2005 Water Model and Scaled Charges for the Ions. J. Chem. Phys. 2019, 151, 134504. [Google Scholar] [CrossRef] [PubMed]
  35. Blazquez, S.; Conde, M.M.; Abascal, J.L.F.; Vega, C. The Madrid-2019 Force Field for Electrolytes in Water Using TIP4P/2005 and Scaled Charges: Extension to the Ions F, Br, I, Rb+, and Cs+. J. Chem. Phys. 2022, 156, 044505. [Google Scholar] [CrossRef] [PubMed]
  36. Blazquez, S.; Conde, M.M.; Vega, C. Scaled Charges for Ions: An Improvement but Not the Final Word for Modeling Electrolytes in Water. J. Chem. Phys. 2023, 158, 054505. [Google Scholar] [CrossRef]
  37. Duignan, T.T.; Baer, M.D.; Schenter, G.K.; Mundy, C.J. Real Single Ion Solvation Free Energies with Quantum Mechanical Simulation. Chem. Sci. 2017, 8, 6131–6140. [Google Scholar] [CrossRef]
  38. Jones, R.O. Density Functional Theory: Its Origins, Rise to Prominence, and Future. Rev. Mod. Phys. 2015, 87, 897–923. [Google Scholar] [CrossRef]
  39. Verma, P.; Truhlar, D.G. Status and Challenges of Density Functional Theory. Trends Chem. 2020, 2, 302–318. [Google Scholar] [CrossRef]
  40. Jorgensen, W.L.; Chandrasekhar, J.; Madura, J.D.; Impey, R.W.; Klein, M.L. Comparison of Simple Potential Functions for Simulating Liquid Water. J. Chem. Phys. 1983, 79, 926–935. [Google Scholar] [CrossRef]
  41. Cornell, W.D.; Cieplak, P.; Bayly, C.I.; Gould, I.R.; Merz, K.M., Jr.; Ferguson, D.M.; Spellmeyer, D.C.; Fox, T.; Caldwell, J.W.; Kollman, P.A. A Second Generation Force Field for the Simulation of Proteins, Nucleic Acids, and Organic Molecules. J. Am. Chem. Soc. 1995, 117, 5179–5197. [Google Scholar] [CrossRef]
  42. Bridwell-Rabb, J.; Kang, G.; Zhong, A.; Liu, H.-W.; Drennan, C.L. An HD Domain Phosphohydrolase Active Site Tailored for Oxetanocin-a Biosynthesis. Proc. Natl. Acad. Sci. USA 2016, 113, 13750–13755. [Google Scholar] [CrossRef] [PubMed]
  43. Stephens, P.J.; Devlin, F.J.; Chabalowski, C.F.; Frisch, M.J. Ab Initio Calculation of Vibrational Absorption and Circular Dichroism Spectra Using Density Functional Force Fields. J. Phys. Chem. 1994, 98, 11623–11627. [Google Scholar] [CrossRef]
  44. Grimme, S.; Ehrlich, S.; Goerigk, L. Effect of the Damping Function in Dispersion Corrected Density Functional Theory. J. Comput. Chem. 2011, 32, 1456–1465. [Google Scholar] [CrossRef]
  45. Weigend, F.; Ahlrichs, R. Balanced Basis Sets of Split Valence, Triple Zeta Valence and Quadruple Zeta Valence Quality for H to Rn: Design and Assessment of Accuracy. Phys. Chem. Chem. Phys. 2005, 7, 3297–3305. [Google Scholar] [CrossRef]
  46. Joung, I.S.; Cheatham, T.E., III. Determination of Alkali and Halide Monovalent Ion Parameters for Use in Explicitly Solvated Biomolecular Simulations. J. Phys. Chem. B 2008, 112, 9020–9041. [Google Scholar] [CrossRef]
  47. Mardirossian, N.; Head-Gordon, M. ωb97x-V: A 10-Parameter, Range-Separated Hybrid, Generalized Gradient Approximation Density Functional with Nonlocal Correlation, Designed by a Survival-of-the-Fittest Strategy. Phys. Chem. Chem. Phys. 2014, 16, 9904–9924. [Google Scholar] [CrossRef]
  48. Rappoport, D.; Furche, F. Property-Optimized Gaussian Basis Sets for Molecular Response Calculations. J. Chem. Phys. 2010, 133, 134105. [Google Scholar] [CrossRef]
  49. Mao, Y.; Demerdash, O.; Head-Gordon, M.; Head-Gordon, T. Assessing Ion–Water Interactions in the Amoeba Force Field Using Energy Decomposition Analysis of Electronic Structure Calculations. J. Chem. Theory Comput. 2016, 12, 5422–5437. [Google Scholar] [CrossRef]
  50. Lu, T.; Chen, F. Multiwfn: A Multifunctional Wavefunction Analyzer. J. Comput. Chem. 2012, 33, 580–592. [Google Scholar] [CrossRef]
  51. Liu, W.; Hong, G.; Dai, D.; Li, L.; Dolg, M. The Beijing Four-Component Density Functional Program Package (BDF) and Its Application to Euo, Eus, Ybo and Ybs. Theor. Chem. Acc. 1997, 96, 75–83. [Google Scholar] [CrossRef]
  52. Zhang, Y.; Suo, B.; Wang, Z.; Zhang, N.; Li, Z.; Lei, Y.; Zou, W.; Gao, J.; Peng, D.; Pu, Z.; et al. BDF: A Relativistic Electronic Structure Program Package. J. Chem. Phys. 2020, 152, 064113. [Google Scholar] [CrossRef] [PubMed]
  53. Liu, W.; Wang, F.; Li, L. The Beijing Density Functional (BDF) Program Package: Methodologies and Applications. J. Theor. Comput. Chem. 2003, 02, 257–272. [Google Scholar] [CrossRef]
  54. Liu, W.; Wang, F.; Li, L. Relativistic Density Functional Theory: The BDF Program Package. In Recent Advances in Relativistic Molecular Theory; World Scientific Publishing: Singapore, 2004; pp. 257–282. [Google Scholar]
  55. Abraham, M.J.; Murtola, T.; Schulz, R.; Páll, S.; Smith, J.C.; Hess, B.; Lindahl, E. GROMACS: High Performance Molecular Simulations through Multi-Level Parallelism from Laptops to Supercomputers. SoftwareX 2015, 1, 19–25. [Google Scholar] [CrossRef]
  56. Darden, T.; York, D.; Pedersen, L. Particle Mesh Ewald: An N.Log(N) Method for Ewald Sums in Large Systems. J. Chem. Phys. 1993, 98, 10089–10092. [Google Scholar] [CrossRef]
  57. Essmann, U.; Perera, L.; Berkowitz, M.L.; Darden, T.; Lee, H.; Pedersen, L.G. A Smooth Particle Mesh Ewald Method. J. Chem. Phys. 1995, 103, 8577–8593. [Google Scholar] [CrossRef]
  58. Bussi, G.; Donadio, D.; Parrinello, M. Canonical Sampling through Velocity Rescaling. J. Chem. Phys. 2007, 126, 014101. [Google Scholar] [CrossRef]
  59. Parrinello, M.; Rahman, A. Polymorphic Transitions in Single Crystals: A New Molecular Dynamics Method. J. Appl. Phys. 1981, 52, 7182–7190. [Google Scholar] [CrossRef]
  60. Nose, S.; Klein, M.L. Constant Pressure Molecular Dynamics for Molecular Systems. Mol. Phys. 1983, 50, 1055–1076. [Google Scholar] [CrossRef]
  61. Miyamoto, S.; Kollman, P.A. Settle—An Analytical Version of the Shake and Rattle Algorithm for Rigid Water Models. J. Comput. Chem. 1992, 13, 952–962. [Google Scholar] [CrossRef]
  62. Hess, B.; Bekker, H.; Berendsen, H.J.C.; Fraaije, J.G.E.M. Lincs: A Linear Constraint Solver for Molecular Simulations. J. Comput. Chem. 1997, 18, 1463–1472. [Google Scholar] [CrossRef]
  63. Zhang, Y.; Jiang, Y.; Qiu, Y.; Zhang, H. Rational Design of Nonbonded Point Charge Models for Highly Charged Metal Cations with Lennard-Jones 12-6 Potential. J. Chem. Inf. Model. 2021, 61, 4613–4629. [Google Scholar] [CrossRef] [PubMed]
  64. Qiu, Y.; Jiang, Y.; Zhang, Y.; Zhang, H. Rational Design of Nonbonded Point Charge Models for Monovalent Ions with Lennard-Jones 12-6 Potential. J. Phys. Chem. B 2021, 125, 13502–13518. [Google Scholar] [CrossRef] [PubMed]
  65. Zhang, J.; Tuguldur, B.; van der Spoel, D. Force Field Benchmark of Organic Liquids. 2. Gibbs Energy of Solvation. J. Chem. Inf. Model. 2015, 55, 1192–1201. [Google Scholar] [CrossRef] [PubMed]
  66. Misin, M.; Fedorov, M.V.; Palmer, D.S. Hydration Free Energies of Molecular Ions from Theory and Simulation. J. Phys. Chem. B 2016, 120, 975–983. [Google Scholar] [CrossRef] [PubMed]
  67. Bennett, C.H. Efficient Estimation of Free Energy Differences from Monte Carlo Data. J. Comput. Phys. 1976, 22, 245–268. [Google Scholar] [CrossRef]
  68. Lamoureux, G.; Roux, B. Absolute Hydration Free Energy Scale for Alkali and Halide Ions Established from Simulations with a Polarizable Force Field. J. Phys. Chem. B 2006, 110, 3308–3322. [Google Scholar] [CrossRef]
  69. Zhang, H.; Jiang, Y.; Yan, H.; Yin, C.; Tan, T.; van der Spoel, D. Free-Energy Calculations of Ionic Hydration Consistent with the Experimental Hydration Free Energy of the Proton. J. Phys. Chem. Lett. 2017, 8, 2705–2712. [Google Scholar] [CrossRef]
  70. Zhang, H.; Jiang, Y.; Yan, H.; Cui, Z.; Yin, C. Comparative Assessment of Computational Methods for Free Energy Calculations of Ionic Hydration. J. Chem. Inf. Model. 2017, 57, 2763–2775. [Google Scholar] [CrossRef]
  71. Zhang, H.; Yin, C.; Jiang, Y.; van der Spoel, D. Force Field Benchmark of Amino Acids: I. Hydration and Diffusion in Different Water Models. J. Chem. Inf. Model. 2018, 58, 1037–1052. [Google Scholar] [CrossRef] [PubMed]
  72. Song, Y.; Dai, L.L. The Shear Viscosities of Common Water Models by Non-Equilibrium Molecular Dynamics Simulations. Mol. Simulat. 2010, 36, 560–567. [Google Scholar] [CrossRef]
  73. Maier, J.A.; Martinez, C.; Kasavajhala, K.; Wickstrom, L.; Hauser, K.E.; Simmerling, C. Ff14sb: Improving the Accuracy of Protein Side Chain and Backbone Parameters from Ff99sb. J. Chem. Theory Comput. 2015, 11, 3696–3713. [Google Scholar] [CrossRef] [PubMed]
  74. Benedict, W.S.; Gailar, N.; Plyler, E.K. Rotation-Vibration Spectra of Deuterated Water Vapor. J. Chem. Phys. 1956, 24, 1139–1165. [Google Scholar] [CrossRef]
  75. Kim, J.; Lee, H.M.; Suh, S.B.; Majumdar, D.; Kim, K.S. Comparative Ab Initio Study of the Structures, Energetics and Spectra of X·(H2O)n=1–4 [X=F, Cl, Br, I] Clusters. J. Chem. Phys. 2000, 113, 5259–5272. [Google Scholar] [CrossRef]
  76. Bajaj, P.; Götz, A.W.; Paesani, F. Toward Chemical Accuracy in the Description of Ion–Water Interactions through Many-Body Representations. I. Halide–Water Dimer Potential Energy Surfaces. J. Chem. Theory Comput. 2016, 12, 2698–2705. [Google Scholar]
  77. Marcus, Y. Thermodynamics of Solvation of Ions. Part 5.—Gibbs Free Energy of Hydration at 298.15 K. J. Chem. Soc. Faraday Trans. 1991, 87, 2995–2999. [Google Scholar] [CrossRef]
  78. Hünenberger, P.; Reif, M. Single-Ion Solvation: Experimental and Theoretical Approaches to Elusive Thermodynamic Quantities; Royal Society of Chemistry: Cambridge, UK, 2011. [Google Scholar]
  79. Gates, J.A.; Wood, R.H. Densities of Aqueous Solutions of Sodium Chloride, Magnesium Chloride, Potassium Chloride, Sodium Bromide, Lithium Chloride, and Calcium Chloride from 0.05 to 5.0 mol kg-1 and 0.1013 to 40 MPa at 298.15 K. J. Chem. Eng. Data 1985, 30, 44–49. [Google Scholar] [CrossRef]
  80. Laliberté, M.; Cooper, W.E. Model for Calculating the Density of Aqueous Electrolyte Solutions. J. Chem. Eng. Data 2004, 49, 1141–1151. [Google Scholar] [CrossRef]
  81. Laliberté, M. A Model for Calculating the Heat Capacity of Aqueous Solutions, with Updated Density and Viscosity Data. J. Chem. Eng. Data 2009, 54, 1725–1760. [Google Scholar] [CrossRef]
  82. Laliberté, M. Model for Calculating the Viscosity of Aqueous Solutions. J. Chem. Eng. Data 2007, 52, 321–335. [Google Scholar] [CrossRef]
  83. Fan, K.; Zhang, Y.; Qiu, Y.; Zhang, H. Impacts of Targeting Different Hydration Free Energy References on the Development of Ion Potentials. Phys. Chem. Chem. Phys. 2022, 24, 16244–16262. [Google Scholar] [CrossRef] [PubMed]
  84. Noyes, R.M. Thermodynamics of Ion Hydration as a Measure of Effective Dielectric Properties of Water. J. Am. Chem. Soc. 1962, 84, 513–522. [Google Scholar] [CrossRef]
  85. Rosseinsky, D.R. Electrode Potentials and Hydration Energies. Theories and Correlations. Chem. Rev. 1965, 65, 467–490. [Google Scholar]
  86. Fong, K.D.; Sumić, B.; O’Neill, N.; Schran, C.; Grey, C.P.; Michaelides, A. The Interplay of Solvation and Polarization Effects on Ion Pairing in Nanoconfined Electrolytes. Nano Lett. 2024, 24, 5024–5030. [Google Scholar] [CrossRef] [PubMed]
  87. Park, S.; Park, H.; Park, C.B.; Sung, B.J. The Effects of Polarization on the Rotational Diffusion of Ions in Organic Ionic Plastic Crystals. J. Chem. Phys. 2022, 157, 144501. [Google Scholar] [CrossRef]
  88. Phan, L.X.; Lynch, C.I.; Crain, J.; Sansom, M.S.P.; Tucker, S.J. Influence of Effective Polarization on Ion and Water Interactions within a Biomimetic Nanopore. Biophys. J. 2022, 121, 2014–2026. [Google Scholar] [CrossRef] [PubMed]
  89. Zhang, Z.; Zofchak, E.; Krajniak, J.; Ganesan, V. Influence of Polarizability on the Structure, Dynamic Characteristics, and Ion-Transport Mechanisms in Polymeric Ionic Liquids. J. Phys. Chem. B 2022, 126, 2583–2592. [Google Scholar] [CrossRef]
  90. Lemkul, J.A.; Huang, J.; Roux, B.; MacKerell, A.D., Jr. An Empirical Polarizable Force Field Based on the Classical Drude Oscillator Model: Development History and Recent Applications. Chem. Rev. 2016, 116, 4983–5013. [Google Scholar] [CrossRef]
  91. Bedrov, D.; Piquemal, J.-P.; Borodin, O.; MacKerell, A.D., Jr.; Roux, B.; Schröder, C. Molecular Dynamics Simulations of Ionic Liquids and Electrolytes Using Polarizable Force Fields. Chem. Rev. 2019, 119, 7940–7995. [Google Scholar] [CrossRef]
  92. Kognole, A.A.; Aytenfisu, A.H.; MacKerell, A.D. Balanced Polarizable Drude Force Field Parameters for Molecular Anions: Phosphates, Sulfates, Sulfamates, and Oxides. J. Mol. Model. 2020, 26, 152. [Google Scholar] [CrossRef]
  93. Villa, F.; MacKerell, A.D., Jr.; Roux, B.; Simonson, T. Classical Drude Polarizable Force Field Model for Methyl Phosphate and Its Interactions with Mg2+. J. Phys. Chem. A 2018, 122, 6147–6155. [Google Scholar] [CrossRef] [PubMed]
  94. Kostal, V.; Jungwirth, P.; Martinez-Seara, H. Nonaqueous Ion Pairing Exemplifies the Case for Including Electronic Polarization in Molecular Dynamics Simulations. J. Phys. Chem. Lett. 2023, 14, 8691–8696. [Google Scholar] [CrossRef] [PubMed]
  95. Duboué-Dijon, E.; Javanainen, M.; Delcroix, P.; Jungwirth, P.; Martinez-Seara, H. A Practical Guide to Biologically Relevant Molecular Simulations with Charge Scaling for Electronic Polarization. J. Chem. Phys. 2020, 153, 050901. [Google Scholar] [CrossRef] [PubMed]
Figure 1. Crystal structure of the HD-domain phosphohydrolase OxsA (PDB code: 5TKA) and an enlarged view of the metal-binding center. Mg2+ coordinated with four amino acid residues (His31, His66, Asp67, and Asp132) and two water molecules.
Figure 1. Crystal structure of the HD-domain phosphohydrolase OxsA (PDB code: 5TKA) and an enlarged view of the metal-binding center. Mg2+ coordinated with four amino acid residues (His31, His66, Asp67, and Asp132) and two water molecules.
Metals 14 01009 g001
Figure 2. Optimized structures of (a) Mg2+, (b) Ca2+, and (c,d) Cl ions interacting with one single water in the C2v and Cs symmetries at the B3LYP-D3(BJ)/def2-TZVP level as simulated using the BDF software. Mg2+, Ca2+, and Cl ions are shown as green, yellow, and cyan balls, respectively. Water is depicted in red (oxygen) and white (hydrogen), and its experimental geometry is presented in panel (a).
Figure 2. Optimized structures of (a) Mg2+, (b) Ca2+, and (c,d) Cl ions interacting with one single water in the C2v and Cs symmetries at the B3LYP-D3(BJ)/def2-TZVP level as simulated using the BDF software. Mg2+, Ca2+, and Cl ions are shown as green, yellow, and cyan balls, respectively. Water is depicted in red (oxygen) and white (hydrogen), and its experimental geometry is presented in panel (a).
Metals 14 01009 g002
Figure 3. Contour plots of gas-phase IE and gas-phase IOD (ac) and of gas-phase IE and liquid-phase IOD (df) for Mg2+, Ca2+, and Cl ions in the C2v and Cs symmetries. The intersection points are the designed LJ parameters for the corresponding ion models.
Figure 3. Contour plots of gas-phase IE and gas-phase IOD (ac) and of gas-phase IE and liquid-phase IOD (df) for Mg2+, Ca2+, and Cl ions in the C2v and Cs symmetries. The intersection points are the designed LJ parameters for the corresponding ion models.
Metals 14 01009 g003
Figure 4. Radial distribution functions (RDFs) (left axis) and coordination numbers (CNs) of water molecules around the ions of (a) Mg2+, (b) Ca2+, and (c) Cl. The gray dotted line indicate the experimental IOD values (Table 1).
Figure 4. Radial distribution functions (RDFs) (left axis) and coordination numbers (CNs) of water molecules around the ions of (a) Mg2+, (b) Ca2+, and (c) Cl. The gray dotted line indicate the experimental IOD values (Table 1).
Metals 14 01009 g004
Figure 5. Interaction energy (IE) profiles between the ion and a single TIP3P water in the C2v (ac) and Cs (d) symmetries as a function of ion–water oxygen (rion−OW; IOD) distance. The black lines represent the QM predictions at the B97X-V/def2-QZVPPD level as calculated using BDF software.
Figure 5. Interaction energy (IE) profiles between the ion and a single TIP3P water in the C2v (ac) and Cs (d) symmetries as a function of ion–water oxygen (rion−OW; IOD) distance. The black lines represent the QM predictions at the B97X-V/def2-QZVPPD level as calculated using BDF software.
Metals 14 01009 g005
Figure 6. Calculated densities as a function of concentration for (a) MgCl2 and (b) CaCl2 salt solutions. The solid lines (in black) indicate the experimental measurements taken from Refs. [79,80,81].
Figure 6. Calculated densities as a function of concentration for (a) MgCl2 and (b) CaCl2 salt solutions. The solid lines (in black) indicate the experimental measurements taken from Refs. [79,80,81].
Metals 14 01009 g006
Figure 7. Calculated viscosities of (a) MgCl2 and (b) CaCl2 solutions with different ionic concentrations. The solid lines (in black) indicate the experimental observations taken from Refs. [81,82].
Figure 7. Calculated viscosities of (a) MgCl2 and (b) CaCl2 solutions with different ionic concentrations. The solid lines (in black) indicate the experimental observations taken from Refs. [81,82].
Metals 14 01009 g007
Figure 8. Root-mean-square deviations (RMSDs) of the protein backbone and Mg2+-binding residues from the crystal structure of the phosphohydrolase OxsA (PDB ID: 5TKA) as a function of simulation time for the three Mg2+ models of (a) Mg(I), (b) Mg(II), and (c) Mg(Åqvist). The RMSDs were calculated without consideration of hydrogen atoms.
Figure 8. Root-mean-square deviations (RMSDs) of the protein backbone and Mg2+-binding residues from the crystal structure of the phosphohydrolase OxsA (PDB ID: 5TKA) as a function of simulation time for the three Mg2+ models of (a) Mg(I), (b) Mg(II), and (c) Mg(Åqvist). The RMSDs were calculated without consideration of hydrogen atoms.
Metals 14 01009 g008
Figure 9. Interatomic distances between Mg2+ and carboxyl oxygen atoms (OD1 in black and OD2 in red) of Asp67 and Asp132 as a function of simulation time and the representative metal-binding centers for the (a) Mg(I), (b) Mg(II), and (c) Mg(Åqvist) models, respectively.
Figure 9. Interatomic distances between Mg2+ and carboxyl oxygen atoms (OD1 in black and OD2 in red) of Asp67 and Asp132 as a function of simulation time and the representative metal-binding centers for the (a) Mg(I), (b) Mg(II), and (c) Mg(Åqvist) models, respectively.
Metals 14 01009 g009
Table 1. QM predictions of interaction energies (IEs, kcal/mol) and ion–water oxygen distances (IODs, Å) for the tested ion–water dimers in the gas phase as well as the conventional HFE (ΔGo,conv), real HFE ( G * , r e a l , kcal/mol), IOD, and coordination number (CN) values of the ions in aqueous solution.
Table 1. QM predictions of interaction energies (IEs, kcal/mol) and ion–water oxygen distances (IODs, Å) for the tested ion–water dimers in the gas phase as well as the conventional HFE (ΔGo,conv), real HFE ( G * , r e a l , kcal/mol), IOD, and coordination number (CN) values of the ions in aqueous solution.
IonsGas Phase Liquid Phase
SymmetryIEIOD G o , c o n v G * , r e a l IODCN
Mg2+C2v−82.141.92 65.5−456.02.0906
Ca2+C2v−56.672.23 143.1−378.42.4608
ClC2v−11.763.16 −335.6−77.73.1876–8.5
ClCs−13.763.12 −335.6−77.73.1876–8.5
ΔGo,conv, IOD, and CN are taken from Refs. [25,64,77]; ΔG*,real was computed using Equation (3).
Table 2. Designed point charge models of the ions and the simulated IE (kcal/mol), IOD (Å), HFE (kcal/mol), and CN values for Mg2+, Ca2+, and Cl with TIP3P water.
Table 2. Designed point charge models of the ions and the simulated IE (kcal/mol), IOD (Å), HFE (kcal/mol), and CN values for Mg2+, Ca2+, and Cl with TIP3P water.
IonsModel aTarget Properties bR (Å)ε (kcal/mol)IE (kcal/mol)IODgasIODliqHFECN
C2vCsC2vCs
Mg2+Mg(I)IE, IODgas0.3293697.6676−82.141.931.98−539.60 ± 0.536.0
Mg(II)IE, IODliq0.39271611.5925−81.862.032.09−560.57 ± 0.146.7
Mg(Åqvist) cHFE0.79260.8947−69.561.891.98−457.75 ± 0.176.0
Ca2+Ca(I)IE, IODgas0.782569.7699−56.672.232.34−432.02 ± 0.098.0
Ca(II)IE, IODliq0.8081445.7715−56.952.382.46−467.90 ± 0.089.0
Ca(Åqvist) cHFE1.71310.4598−38.162.612.70−331.98 ± 0.128.9
ClCl(I)IE, IODgas2.34420.1365−13.59−13.763.193.123.18−77.42 ± 0.077.7
Cl(II)IE, IODliq2.30210.1892−13.61−13.763.203.123.18−77.64 ± 0.027.6
Cl(JC) dIE, HFE, LE, LC2.51300.0356−13.67−13.953.153.073.12−76.87 ± 0.077.3
a Models I and II were designed in this work. The existing models of Mg(Åqvist), Ca(Åqvist), and Cl(JC) in the Amber force fields were examined for comparison. b Targeted properties in the model development, including interaction energy (IE), ion−water oxygen distance (IOD) in gas or liquid phases, and hydration free energy (HFE) as well as the crystal properties of lattice energy (LE) and lattice constant (LC). c Adapted from Ref. [26]. d Taken from Ref. [46].
Table 3. Interatomic distance (Å) between Mg2+ and binding residues from the last 20 ns of the MD simulations for the Mg2+-binding center of the phosphohydrolase OxsA (PDB code: 5TKA).
Table 3. Interatomic distance (Å) between Mg2+ and binding residues from the last 20 ns of the MD simulations for the Mg2+-binding center of the phosphohydrolase OxsA (PDB code: 5TKA).
ModelH31 NE2H66 NE2D67 OD1D67 OD2D132 OD1D132 OD2Water1Water2
Mg(I)2.092.093.301.861.863.682.012.01
Mg(II)2.242.152.001.982.011.992.10none
Mg(Åqvist)2.172.193.581.861.863.742.012.02
crystal2.242.213.302.092.123.931.941.92
A distance smaller than 2.5 Å was considered to indicate a coordinated ligand; the standard deviations of the calculated distances for the coordinated and non-coordinated ligands amounted to 0.03−0.09 and 0.12−0.18 Å, respectively, and are not given here; the term “none” indicates that a second water coordinated with Mg2+ did not exist.
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Zhang, Y.; Wu, B.; Lu, C.; Zhang, H. Design of Point Charge Models for Divalent Metal Cations Targeting Quantum Mechanical Ion–Water Dimer Interactions. Metals 2024, 14, 1009. https://doi.org/10.3390/met14091009

AMA Style

Zhang Y, Wu B, Lu C, Zhang H. Design of Point Charge Models for Divalent Metal Cations Targeting Quantum Mechanical Ion–Water Dimer Interactions. Metals. 2024; 14(9):1009. https://doi.org/10.3390/met14091009

Chicago/Turabian Style

Zhang, Yongguang, Binghan Wu, Chenyi Lu, and Haiyang Zhang. 2024. "Design of Point Charge Models for Divalent Metal Cations Targeting Quantum Mechanical Ion–Water Dimer Interactions" Metals 14, no. 9: 1009. https://doi.org/10.3390/met14091009

APA Style

Zhang, Y., Wu, B., Lu, C., & Zhang, H. (2024). Design of Point Charge Models for Divalent Metal Cations Targeting Quantum Mechanical Ion–Water Dimer Interactions. Metals, 14(9), 1009. https://doi.org/10.3390/met14091009

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop