Bridging the 12-6-4 Model and the Fluctuating Charge Model

Metal ions play important roles in various biological systems. Molecular dynamics (MD) using classical force field has become a popular research tool to study biological systems at the atomic level. However, meaningful MD simulations require reliable models and parameters. Previously we showed that the 12-6 Lennard-Jones nonbonded model for ions could not reproduce the experimental hydration free energy (HFE) and ion-oxygen distance (IOD) values simultaneously when ion has a charge of +2 or higher. We discussed that this deficiency arises from the overlook of the ion-induced dipole interaction in the 12-6 model, and this term is proportional to 1/r 4 based on theory. Hence, we developed the 12-6-4 model and showed it could solve this deficiency in a physically meaningful way. However, our previous research also found that the 12-6-4 model overestimated the coordination numbers (CNs) for some highly charged metal ions. And we attributed this artifact to that the current 12-6-4 scheme lacks a correction for the interactions among the first solvation shell water molecules. In the present study, we considered the ion-included dipole interaction by using the 12-6 model with adjusting the atomic charges of the first solvation shell water molecules. This strategy not only considers the ion-induced dipole interaction between ion and the first solvation shell water molecules but also well accounts for the increased repulsion among these water molecules compared to the bulk water molecules. We showed this strategy could well reproduce the experimental HFE and IOD values for Mg2+, Zn2+, Al3+, Fe3+, and In3+ and solve the CN overestimation issue of the 12-6-4 model for Fe3+ and In3+. Moreover, our simulation results showed good agreement with previous ab initio MD simulations. In addition, we derived the physical relationship between the C 4 parameter and induced dipole moment, which agreed well with our simulation results. Finally, we discussed the implications of the present work for simulating metalloproteins. Due to the fluctuating charge model uses a similar concept to the 12-6 model with adjusting atomic charges, we believe the present study builds a bridge between the 12-6-4 model and the fluctuating charge model.


INTRODUCTION
Metal ion plays significant roles in various biological processes (Thomson and Gray, 1998;Andreini et al., 2004;Woodson, 2005;Waldron and Robinson, 2009;Kepp, 2012). Molecular dynamics (MD) simulations has become an important tool for studying biomolecules (Duan and Kollman, 1998;Lindorff-Larsen et al., 2011;De Vivo et al., 2016;Hollingsworth and Dror, 2018). It can provide information with atomic details and high time resolution. Meaningful MD simulations require reliable models and parameters, but it is challenging to develop accurate models for metal ions (Li and Merz, 2017). Different force field models have been developed for simulating metal ions in biological systems. For example, the bonded model, (Peters et al., 2010), the cationic dummy atom model, (Duarte et al., 2014), the nonbonded model, (Li et al., 2013), and the polarizable models (Sakharov and Lim, 2009;Zhang et al., 2012). Among these models, the nonbonded model is one of the most widely used models because of its simplicity and transferability.
In the past decade, Li, Merz, and co-workers found that the widely used 12-6 Lennard-Jones nonbonded model could not reproduce the experimental hydration free energy (HFE) and ionoxygen distance (IOD) simultaneously when metal ion has a charge of +2 or higher (Li et al., 2013;Li and Merz, 2014;Li et al.,2015a;Li et al., 2015b). They attributed this deficiency to that the 12-6 nonbonded model did not consider the ion-induced dipole interaction. Because the ion-induced dipole interaction is proportional to 1/r 4 , where r is the atomic distance, they proposed the 12-6-4 model and parameterized it for various atomic ions (Li and Merz, 2014;Li et al.,2015a;Li et al., 2015b). The 12-6-4 model could accurately reproduce the experimental HFE and IOD simultaneously after careful parameterization, which successfully solved the deficiency of the 12-6 model (Li and Merz, 2014;Li et al., 2015a;Li et al., 2015b). Moreover, Merz and co-workers have showed that the 12-6-4 model could well simulate the chelate effect (Sengupta et al., 2018) and the thermodynamics of ion binding in a metalloprotein system . However, it was found that for some highly charged metal ions such as Fe 3+ and In 3+ , the coordination number (CN) value was overestimated by the 12-6-4 model (Li et al., 2015a). For example, Fe 3+ has an experimental CN of 6, (Marcus, 1988), but the 12-6-4 model predicted its CN value as 6.8-6.9 when using the TIP3P, SPC/E, or TIP4P EW water model (Li et al., 2015a). Recently, Li, Merz, and co-workers have parameterized the 12-6-4 for various ions in conjunction with four new water models Li et al., 2021;Sengupta et al., 2021). But this issue still exists: the 12-6-4 model predicted a CN of 6.4-6.7 for Fe 3+ when using the OPC3, OPC, TIP3P-FB, or TIP4P-FB water model . It was proposed that this artifact was due to an inaccurate description of the water-water interactions inside the first solvation shell, which effect is not significant for monovalent and divalent metal ions but becomes more severe for the highly charged ions (Li et al., 2015a). To be more specific, a metal ion with a higher charge will cause larger induced dipoles of the first solvation water molecules, yielding a stronger repulsion between these water molecules. However, the current 12-6-4 model implementation does not take this effect into account, although it considers the ioninduced dipole interactions between the metal ion and surrounding water molecules.
In the present study, we tried to solve the CN overestimation issue by using a 12-6 model with increasing the dipole moments of the first solvation shell water molecules to account for the ioninduced dipole effect. Specifically, we studied the ion-aqueous systems containing Mg 2+ , Zn 2+ , Al 3+ , Fe 3+ , and In 3+ . We showed that this treatment could well reproduce the experimental HFEs and IODs for these ions, and solve the CN overestimation issue for Fe 3+ and In 3+ . The adjusted dipole moments of the first solvation shell water molecules for different ions are consistent with theory and ab initio molecular dynamics (AIMD) simulations. In addition, we derived the physical relationship between the C 4 parameter and induced dipole moment, which has excellent agreement with the simulation results we got. The concept and strategy of the current study can be further applied to biological systems such as metalloproteins. Considering the present study used a similar scheme to the fluctuating charge model, we believe it can bridge the 12-6-4 model and the fluctuating charge model.

COMPUTATIONAL METHODS
The potential function of the 12-6 Lennard-Jones nonbonded model in the AMBER force field (Cornell et al., 1995) is shown below: The potential function of the 12-6-4 nonbonded model (Li and Merz, 2014;Li et al., 2015a;Li et al., 2015b) is shown as follows: In these equations, U ij is the potential energy between atoms i and j. r ij is the distance between the two atoms. Q i and Q j are the particle charges of atoms i and j, respectively. e is the charge of the proton. C ij 12 /r 12 ij is the Pauli repulsive term, −C ij 6 /r 6 ij is the dispersion attractive term, while −C ij 4 /r 4 ij is the ion-induced dipole term. R min,ij is the distance between atoms i and j where the Lennard-Jones potential has its minimum, while ε ij is the well depth of the Lennard-Jones potential at this distance. The Berthelot-Lorentz (LB) combining rules are used in the AMBER force field: (Cornell et al., 1995): In the current study, the OPC3 water model was used to simulate the solvent as it represents an accuracy limit for the rigid three-point water models. (Izadi and Onufriev, 2016). The thermodynamic integration (TI) method (Kollman, 1993) was employed to calculate the ΔG values. Calculations of the ΔG values were performed under the NVT ensemble by assuming that the Gibbs free energy change is identical to the Helmholtz free energy change. The IOD and CN values were obtained based on the radial distribution function (RDF) between the metal ion and water oxygen atoms. Specific steps are described as follows.
For the 12-6 model with regular charges, we calculated the HFE, IOD, and CN values with the method described previously Sengupta et al., 2021). Among them the regular MD simulations were used to obtain the IOD and CN values.
For the 12-6 model with updated charges, we calculate the HFE, IOD, and CN values as described below. By treating the last snapshot of the regular MD simulations mentioned above as initial structure, we performed a set of TI calculations to obtain ΔG 1,ind for increasing the dipole moments of the six closest water molecules from the metal ion. Afterwards, we performed a set of TI calculations in the reverse manner to get ΔG 2,ind for decreasing the dipole moments of these water molecules to the original value. The final HFE was calculated as the sum of the HFE for the 12-6 model with regular charges and (ΔG 1,ind -ΔG 2,ind )/2. These TI calculations were performed with linear scaling: U(λ) (1 − λ)U 0 + λU 1 . With U 0 and U 1 represent the initial and final Hamiltonians, respectively. Seven λ windows (with λ 0.02544, 0.12923, 0.29707, 0.5, 0.70292, 0.87076, and 0.97455) were used for each set of TI calculations, with each window had 100 ps equilibration and 200 ps production. After finishing the TI calculations, the Gaussian quadrature was used to calculate ΔG 1,ind and ΔG 2,ind : ΔG i w i 〈zU/zλ〉 i . Herein, w i is the weight for window i, while 〈zU/zλ〉 i is the averaged zU/zλ value for window i. Considering the uncertainties of the calculated HFEs are ∼1 kcal/mol for monovalent and divalent metal ions and ∼2 kcal/mol for the highly charged metal ions when using the 12-6 model with regular charges, (Li et al., 2013;Li et al., 2015a;Li et al., 2015b), we estimated the uncertainties of the final HFE values are ∼3 kcal/mol or less based on the uncertainty propagation, by assuming the additional TI calculations provide an uncertainty of similar magnitude.
For Fe 3+ and In 3+ , a weak restraint was used for the additional TI calculations due to there was a CN switching: the CN value was eight or so when these water molecules have fixed charges of the OPC3 water model, while it changes to six for the equilibrated structure of the system when the six closest water molecules to the ion have enlarged dipole moments. To prevent these six water molecules leaving the first solvation shell (due to they have "induced dipoles" which are caused by their proximity to the ion), when calculating ΔG 1,ind and ΔG 2,ind , a one-side weak restraint of U k(r − r eq ) 2 with k equals 10 kcal/(mol·Å 2 ) was applied to each of these six ion-oxygen distances when r is bigger than 2.5 Å. The same weak restraint was applied to the initial and final states during the TI calculations. Because the CN switching was involved, we carried out longer simulations with each λ window had 200 ps equilibration and 800 ps production for the additional TI calculations of Fe 3+ or In 3+ .
We noticed that there were free energy changes by applying the restraints mentioned above. To calculate these free energy changes (that is the free energy change between the systems with and without restraints for each of the system with regular charges and the system with updated charges), the free energy perturbation (Zwanzig, 1954) equation was used (with assuming the ensembles A and B are the same): Where A represents the system without the restraints, and B represents the system with the restraints. Based on Eq. 5, we obtained the ΔG A → B value for the system with regular charges (herein termed as ΔG regularQ nores → res ), and the ΔG A → B value for the system with updated charges (termed as ΔG updateQ nores → res ). For each of the system with regular charges and the system with updated charges, by starting from the same final structure of the regular MD simulations mentioned above, 2 ns NPT equilibration plus 2 ns NVT production were performed, with the one-side restraints (mentioned in the last paragraph) applied to both the equilibration and production. The six ion-oxygen distances were collected every 10 fs, and they were used to calculate the E B − E A values afterwards, providing 200,000 data points for calculating each of ΔG regularQ nores → res and ΔG updateQ nores → res . Both these two ΔG values should be non-negative as it is always true that E B ≥ E A . Our calculations showed that the result of ΔG regularQ nores → res − ΔG updateQ nores → res is below 1 kcal/mol for Fe 3+ or In 3+ , which is below the uncertainty of the calculated HFEs, so we did not account these two ΔG values in the final reported HFE values.
Moreover, also by starting from the final structure of the regular MD simulations, we increased the dipole moments of the six closest water molecules of the metal ion by adjusting their atomic charges, and then we performed 2 ns NPT equilibration with reassigned initial velocities. Then we carried out 2 ns NVT production simulation with the snapshots saved every 500 fs. No restraints were applied in these MD simulations. Afterwards, the RDF between the metal ion and water oxygen atoms was calculated based on the production trajectory. Finally, the IOD and CN values were calculated based on this RDF as we obtained the IOD and CN values for the system with regular charges.
A quasi-cubic water box with a side length of 40 Å and periodic boundary conditions were used in our MD simulations (including the TI simulations). The particle mesh Ewald (PME) method (Darden et al., 1993;Essmann et al., 1995) was used to deal with the long-range electrostatic interactions. A cut-off of 10 Å was used for the nonbonded interactions in the real space. The size of the charge grid was set to 48 in each of the three dimensions in the reciprocal space. A time-step of 1 fs was used, and the three-point "SHAKE" algorithm (Miyamoto and Kollman, 1992) was employed to constrain the geometries of the water molecules with a tolerance of 10 −5 Å for each bond. The Langevin algorithm with a collision frequency of 2 ps −1 was used to control the temperature. The Brenden's barostat with a relaxation time of 1 ps was employed to control the pressure in the MD simulations in the NPT ensemble. All the simulations were performed using the pmemd. cuda (Salomon-Ferrer et al., 2013;Lee et al., 2020) program in the AMBER software package Frontiers in Chemistry | www.frontiersin.org July 2021 | Volume 9 | Article 721960 (Assisted Model Building with Energy Refinement (AMBER), RRID:SCR_014230, and AmberTools, RRID:SCR_018497) (Case et al., 2020).

Experimental Data
The experimental HFE, IOD, and CN values for the Mg 2+ , Zn 2+ , Al 3+ , Fe 3+ , and In 3+ ions are shown in Table 1. The Marcus HFE set (Marcus, 1991) was used is because it could be considered as a set of "intrinsic" HFEs, which only considers the free energy change caused by the ion-water interactions but not accounts the free energy change for ion crossing the gas-liquid interface. This is consistent with the free energy calculations performed in this study that used periodic boundary conditions, as discussed previously Sengupta et al., 2021).

CN Overestimation Issue of the 12-6-4 Model
Previously, we have shown that when in conjunction with the OPC3 water model, the 12-6 model could not simultaneously reproduce the experimental HFE and IOD values for each of Mg 2+ , Zn 2+ , Al 3+ , Fe 3+ , and In 3+ Li et al., 2021). Specifically, when the 12-6 model was able to reproduce the experimental HFE, it underestimated the IOD by 0.11, 0.22, 0.52, 0.48, and 0.31 Å for Mg 2+ , Zn 2+ , Al 3+ , Fe 3+ , and In 3+ , respectively Li et al.,2021). When the 12-6 model was able to reproduce the experimental IOD, it provided the HFE less negative than the experimental value by 33, 63, 152, 177, and 126 kcal/mol for Mg 2+ , Zn 2+ , Al 3+ , Fe 3+ , and In 3+ , respectively Li et al., 2021). Moreover, we have also parameterized the 12-6-4 model for these ions, which showed significant improvement than the 12-6 model: the 12-6-4 model could simultaneously reproduce the experimental HFE and IOD values for each of these ions Li et al., 2021). However, the 12-6-4 model overestimated the CN values for some highly charged ions (Li et al., 2015a;Li et al., 2021). In Table 2 we showed the 12-6-4 parameters for Mg 2+ , Zn 2+ , Al 3+ , Fe 3+ , and In 3+ in conjunction with the OPC3 water model, and the simulated HFE, IOD, and CN values from the previous studies Li et al., 2021). It can be seen that the 12-6-4 model overestimated the CN values of Fe 3+ and In 3+ as 6.7 and 8, which is significantly higher than the experimentally determined CN values of 6 (see Table 1). In addition, to better understand the influence of the C 4 term, we calculated the HFE, IOD, and CN values by ignoring the C 4 term in the potential (i.e. using the regular 12-6 model with the same Lennard-Jones parameters) in the current study. These results are shown in Table 2 as well. We can see that when ignoring the C 4 term, the simulated HFE changed to less negative by 46, 78, 201, 195, and

Solution of the CN Overestimation Issue
The CN overestimation issue of the 12-6-4 model not only exists for Fe 3+ and In 3+ but also for other highly charged metal ions such as Tl 3+ (Li et al., 2015a;Li et al., 2021). We have discussed that this artifact was caused by the fixed charge nature of the first solvation water molecules. As the ion charge increases, the first solvation shell water molecules should have larger induced dipole moments caused by the ion-induced dipole interaction, yielding a stronger repulsion between them. However, this behavior has not been considered by the current 12-6-4 scheme, which only takes into account the ion-induced dipole interaction between ion and water molecules but not considers its influence on the water-water interactions. Hence, to solve the CN overestimation issue, we increased the dipole moments of the first solvation shell water molecules explicitly to account for their enlarged dipole moments caused by the ion-induced dipole interaction. We understood that this strategy was not perfect as it did not consider the charge fluctuations along with each subtle configuration change inside the first solvation shell and did not change the water dipole moments in the second solvation shell or beyond. However, we believed that this strategy could well serve for a qualitative and demonstration purpose. For example, the ab initio MD simulations indicated that the water molecules in the second solvation shell of Zn 2+ , Al 3+ , and Fe 3+ had almost identical dipole moments as those of bulk water molecules Cauët et al., 2010). By adjusting the dipole moments of the water molecules in the first solvation shell (i.e., the six water molecules which are closest to the metal ion), we were able to reproduce the experimental HFEs within 1 kcal/mol by using the 12-6 model. Using the same water charge parameters, we calculated the IOD and CN values as well. The parameters and simulated results are shown in Table 3. The Lennard-Jones parameters of metal ions employed in these simulations are the same as described in Table 2. We can see that these parameters could well reproduce the experimental IOD and CN values with reasonable accuracy as well. Specifically, in order to reproduce the experimental HFEs, we need to increase the dipole moments of the first solvation shell water molecules by 0.34, 0.58, 0.81, 0.89, and 0.70 D for Mg 2+ , Zn 2+ , Al 3+ , Fe 3+ , and In 3+ respectively. This is consistent with the trend that ion with a higher charge will cause larger induced dipoles of the surrounding molecules, and Zn 2+ has a stronger such ability than Mg 2+ (by considering Zn 2+ has the same IOD but a more negative HFE than that of Mg 2+ , see Table 1).

Induced Dipole Moments Predicted by Theory and Quantum Calculations
To better understand the ion-induced dipole interaction, we performed the following calculations. The induced dipole of a particle caused by a point charge can be calculated based on the following equation: Where μ ind is the induced dipole, α is the polarizability of the particle, E is the electric field strength act on that particle, q is the charge magnitude of the point charge, r is the distance between the point charge and particle, ε 0 is the vacuum permittivity, ε r is the dielectric constant of the medium, α 0 is the polarizability volume of the particle which equals α/4πε 0 . The polarizability volume of a water molecule was determined to be 1.444 Å 3 experimentally (Eisenberg and Kauzmann, 1969). We can calculate its μ ind when there is a monovalent ion placed 2.05 Å away from it and ε r 1 (as used in the AMBER force field (Cornell et al., 1995)) based on Eq. 6: μ ind α 0 q ε r r 2 1.444 × (+1e) (2.05) 2 A · e 1.65 D The permanent dipole moment μ per of a water molecule was determined to be 1.85 D under gas phase (Dyke and Muenter, 1973). So this predicts that the total dipole moment μ tot (i.e., μ per + μ ind ) would be 3.5 D for the water molecule inside a monovalent ion-water complex with the ion-oxygen distance as 2.05 A, by assuming both μ ind and μ per point towards the monovalent ion. Because μ ind is proportional to q, so we can get the μ tot of the water molecule is 5.15 and 6.80 D in the analogous divalent ion-water complex and trivalent ion-water complex, respectively.
The above theoretical calculations assumed the water molecule was a single-point particle. Although it could provide some insights about the magnitude of μ ind , this description was not accurate. To further evaluate μ ind of the water molecule in an    . The 12-6-4 parameters and simulated HFE, IOD, and CN values for Al 3+ , Fe 3+ , and In 3+ are from Ref. .
Frontiers in Chemistry | www.frontiersin.org July 2021 | Volume 9 | Article 721960 ion-water complex, we performed quantum calculations for the dipole moment of a water molecule when placing a point charge on the central axis of the water molecule. The quantum calculations were performed using the B3LYP density functional (Becke, 1988;Lee et al., 1988;Becke, 1993;Stephens et al., 1994) with the 6-31G** basis set (Hehre et al., 1972;Hariharan and Pople, 1973). To keep consistent with the MD study, geometry of the OPC3 water model was used for the water molecule. These quantum chemistry calculations were performed using the Gaussian16 software (Gaussian, RRID:SCR_014897) (Frisch et al., 2019). Herein, we put the point charge 2.05 A away from the coordinating water oxygen atom. The calculated results are shown in Table 4. We can see that when the point charge has a magnitude of +1 e, the induced dipole of the water molecule is 0.78 D (with subtracting 1.96 from 2.74 D). This value is only about half of the value calculated by using Eq. 6. However, in terms of the correlation between μ ind and q, we can see a clear trend which is consistent with Eq. 6 that μ ind is proportional to q. Moreover, to test the dependence of the calculated values on density functionals, we also calculated these values by using the PBE0 density functional Perdew et al., 1997;Adamo and Barone, 1999) with the 6-31G** basis set. The results are shown in Table 4 as well. We can see these values are very close to what we obtained by using the B3LYP density functional.
Although having a quantum description of the water molecule, the above quantum calculations still use a pointcharge representation for the metal ion. Furthermore, we did a literature investigation about previous work of AIMD simulations for metal ions in aqueous phase. We found that the magnitudes of the dipole enlargements for the first solvation shell water molecules of Fe 3+ agreed well with our parameters shown in Table 3. Specifically, previous AIMD simulations of the Zn 2+ -aqueous system indicated the dipole moments of first solvation shell water molecules are 0.7 D larger than those of those of bulk water molecules (Cauët et al., 2010). Other studies indicated that this value is 0.2 D for the Mg 2+ ion (Lightstone et al., 2001) and 1.0 D for the Al 3+ ion . All of these values have excellent agreement with our parameters in Table 3. Previous AIMD simulations found that the first solvation shell water molecules of Fe 3+ had dipoles of 0.5 D higher than those of bulk water molecules, , which only qualitatively agreed with our value in Table 3. This may be due to the challenge of simulating the d-orbital effect by using the classical force field. In general, the strategy employed in the current study not only reproduced multiple physical properties simultaneously through a physically meaningful way but also showed well agreement with AIMD simulations.

Relationship Between the C 4 and μ ind Values
To elucidate the relationship between the C 4 parameter and μ ind , we performed the following derivation. The ion-induced dipole interaction potential can be described by the following equation: Herein, U(r) represents the energy needed to bring a neutral species from infinite away where the electric field strength equals zero to a certain position with an electric field strength of E.
Because the ion-induced dipole term equals −C 4 /r 4 in the 12-6-4 model, we can get: Moreover, we can obtain the following relationship based on Eq. 6: μ ind ε r r 2 α 0 q (10) By combining Eqs 9 and 10, we can obtain: We denote: Hence, we can obtain: k e is a constant of 8.988 × 10 9 J m/C 2 when ε r 1, so we can calculate μ ind based on C 4 , q and r for each of the metal ions investigated in the present study. The results are shown in Table 5. These values showed excellent agreement with the values in Table 3, validating the computational strategy employed in the present study.

Implications to Metalloproteins
Previously, we found that the 12-6 model could not reproduce the experimental HFE and IOD values simultaneously when metal ion has a charge of +2 or higher. We attributed this deficiency to the overlook of the ion-induced dipole interaction and developed the 12-6-4 model accordingly. Encouragingly, the 12-6-4 model could simultaneously reproduce the experimental HFE and IOD values for various metal ions (Li et al., 2013;Li and Merz, 2014;Li et al., 2015a;Li et al., 2015b). Moreover, we also found that the 12-6-4 model could provide increased stability than the 12-6 model when simulating a metalloprotein system . In addition, Sengupta et al. have shown that the regular 12-6 model could not well simulate the chelate effect between ethylenediamine and metal ions, while the 12-6-4 model was able to well model this effect after parameter optimization (Sengupta et al., 2018). Furthermore, Song et al. showed that the 12-6-4 model could well simulate metal ion binding to a protein active site after parameter refinement . These results indicated that the 12-6-4 model can overcome the deficiency of the 12-6 model and provide an improved representation of the metal ion containing systems.
Recently, Macchiagodena et al. developed a 12-6 model for fourcoordinated zinc sites in metalloproteins (Macchiagodena et al., 2019). In this model, they refitted the atomic charges of the coordinating groups in the ligating residues. In order to reproduce the mean distances between Zn 2+ and the ligating atoms, they adjusted their Lennard-Jones parameters. In their final parameters, the charges of the ligating sulfur atoms increased by 6-19% when comparing to the charge of the corresponding sulfur atom in the AMBER force field, while the charges of the ligating nitrogen atoms increased by 75-134% when comparing to charges of the corresponding nitrogen atoms in the AMBER force field. Their parameters successfully reproduced the experimentally determined mean distances between Zn 2+ and the ligating atoms and provided Zn 2+ binding affinities which correlated well with the experimentally determined dissociation constants. The strategy employed in their study is similar to that was used in the present work which considers the polarization effect by adjusting the atomic charges of the coordinating groups. And this strategy could improve the representation of the classical force field just as the 12-6-4 model did.

CONCLUSION
Our previous research found that the widely used 12-6 Lennard-Jones nonbonded model for ions could not reproduce the experimental HFE and IOD values simultaneously when ion has a charge of +2 or higher (Li et al., 2013;Li et al., 2015a;Li et al., 2015b). We discussed that this deficiency arises from the overlook of the ion-induced dipole interaction, which is proportional to 1/r 4 where r is the atomic distance (Li and Merz, 2014). Through an −C 4 /r 4 term which takes the ioninduced dipole interaction into account, we proposed and parameterized the 12-6-4 model for various ions (Li and Merz, 2014;Li et al., 2015a;Li et al., 2015b). Encouragingly, the 12-6-4 model could well reproduce the experimental HFE and IOD values at the same time (Li and Merz, 2014;Li et al., 2015a;Li et al., 2015b). However, we discovered a CN overestimation issue for some highly charged ions, and we discussed that this artifact arises from the fact that the current 12-6-4 model did not consider the increased repulsion among the first solvation water molecules of the metal ion (Li et al., 2015a;Li et al., 2021).
In the present study, we simulated multiple divalent and trivalent ions (Mg 2+ , Zn 2+ , Al 3+ , Fe 3+ , and In 3+ ) by using the 12-6 model with adjusting the atomic charges in the first solvation shell water molecules to consider the ion-induced dipole effect. We showed that this strategy could well reproduce the experimental HFE and IOD values simultaneously for these ions and solve the CN overestimation issue of the 12-6-4 model for Fe 3+ and In 3+ . In addition, we found that the enlarged magnitudes of the dipole moments of the first solvation shell water molecules in our simulations agreed well with previous AIMD simulations, further validating the current strategy. Besides, we derived the relationship between the C 4 parameter and the induced dipole moment according to the physical principles. Moreover, the induced dipole moments calculated based on the C 4 parameters showed good agreement with our simulation results. Finally, we discussed that both the 12-6-4 scheme and the 12-6 scheme with adjusting atomic charges could simulate the polarization effect in metalloprotein systems. Due to the fluctuating charge model also uses the 12-6 scheme with adjusting the atomic charges, we believe the present study can serve as a bridge between the 12-6-4 model and the fluctuating charge model.

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.

AUTHOR CONTRIBUTIONS
PL designed and performed the research; PL analyzed the data and wrote the paper.