Chelator-Based Parameterization of the 12-6-4 Lennard-Jones Molecular Mechanics Potential for More Realistic Metal Ion-Protein Interactions.

Metal ions are associated with a variety of proteins and play critical roles in a wide range of biochemical processes. There are multiple ways to study and quantify protein-metal ion interactions, including molecular dynamics simulations. Recently, the AMBER molecular mechanics forcefield was modified to include a 12-6-4 Lennard-Jones potential, which allows for a better description of nonbonded terms through the additional pairwise Cij coefficients. Here, we demonstrate a method of generating Cij parameters that allows parametrization of specific metal ion-ligating groups in order to tune binding energies computed by thermodynamic integration. The new Cij coefficients were tested on a series of chelators: ethylenediaminetetraacetic acid, nitrilotriacetic acid, egtazic acid, and the EF1 loop peptides from the proteins lanmodulin and calmodulin. The new parameters show significant improvements in computed binding energies relative to existing force fields and produce coordination numbers and ion-oxygen distances that are in good agreement with experimental values. This parametrization method should be extensible to a range of other systems and could be readily adapted to tune properties other than binding energies.

Introduction

Metal ions are thought to be associated with around 50% of all proteins.[1] They play important roles in a multitude of biologically significant processes, such as enzyme catalysis and signal transduction.[2−4] Inclusion of metal ions into molecular mechanics (MM) force fields used for molecular dynamics (MD) simulations poses specific challenges because of the high polarizability of metal ions, which can take part in charge transfers, coordination number (CN) changes, and ligand swapping.[5−7] Several approaches have been developed for simulating metal ions, but these typically lack transferability between systems or simulate one property accurately at the expense of others. For example, explicitly bonded models treat metal–ligand coordination as immutable bonds, which does not allow for ligand and CN changes.[8] Such methods therefore cannot provide insight into coordination dynamics, and multiple simulations with different bonding parameters are required to determine optimal CN in a given system.[7] Another commonly used approach is to rely on nonbonded terms that treat metal–ligand coordination as nonbonded interactions facilitated by Lennard-Jones (LJ) and Coulombic terms. The 12-6LJ nonbonded model is the default in most MM force fields and defines the standard LJ potential using two pairwise parameters for each type of atom–atom interaction (eq ). This approach allows for the metal ions to switch both CN and ligands but lacks charge transfer effects and polarizability. Nonetheless, this approach is commonly used for MD simulations as it is computationally efficient and does not require assumptions about the ligands coordinated to the metal ion, which usually require previous knowledge about the simulated system.

Recently, an array of new 12-6LJ parameters were developed to describe divalent metal ions. The “12-6LJ HFE” parameters accurately reproduced hydration free energies (HFEs), the “12-6LJ IOD” parameters reproduced ion-oxygen distances (IODs), while the “12-6LJ CM” parameters were designed as a compromise between both.[6] As the 12-6LJ potential does not accurately reproduce both HFE and IOD experimental values with a single set of parameters, this led to the development of the new 12-6-4LJ potential.[6,9] This contains an additional pairwise parameter, which accounts for the charge induced dipole and dipole induced dipoles that are neglected in the 12-6LJ potential (eq , below). This addition to the force field made it possible to use a single set of parameters to adequately describe HFE, IOD, and CN values in three commonly used water models, TIP3P, SPC/E, and TIP4Pew.[9−11] However, the C coefficients of the new C/r4 term were parameterized only for the interactions with the water oxygen atoms, whereas parameters for other types of atoms were adopted from molecular polarizability tensor calculations.[12] In a subsequent study of metal ion–nucleotide interactions, it was shown that the adopted C coefficients may not produce accurate binding energies as they overestimated the binding energies to adenine, guanine, and phosphate by up to 20 kJ mol–1.[13] These authors were able to generate a new set of C coefficients that increased the accuracy of estimated binding energies (±0.5 kJ mol–1 relative to experimental values) to nucleotides and phosphate without affecting the metal–water interactions, as the metal ion parameters contain separate C coefficients for each atom type. This suggests that further parameterization of these C coefficients with amino acids offers a reasonable approach to improve the accuracy of MD simulations of metal ions.

In this study we have developed a new set of pairwise C coefficients for Ca2+, Mg2+, Y3+, and La3+ ions using the 12-6-4LJ nonbonded potential. These metals were chosen to investigate metal binding to a recently discovered lanthanide selective protein, lanmodulin (LanM), which is homologous to the calcium binding protein calmodulin (CaM).[14,15] Ethylenediaminetetraacetic acid (EDTA) and nitrilotriacetic acid (NTA) were used to determine the parameters, as these contain metal ion ligating carboxylate groups similar to those found in glutamate and aspartate (Figure and Figure S1) and accurate binding affinities are available (Figure , Table S1). A linear equation system based on the ratios of ligating oxygens and tertiary nitrogen atoms was developed to deconvolute the contribution from each ligating atom type. The new C coefficients were tested against another chelator, egtazic acid (EGTA), and the EF1 loop peptides of LanM and CaM. They give binding energies that are significantly closer to experimental values than the 12-6LJ CM/IOD parameter sets or the 12-6-4LJ potential with default C coefficients and reproduce the observed metal ion selectivity in LanM.

Figure 1

Metal ion coordination in the EF1 loops of CaM (left) and LanM (right). Ca2+ is bound to CaM (PDB 1CLL) and Y3+ is bound to LanM (NMR conformer 1 of PDB 6MI5).[15,16]

Figure 2

Experimental binding energies for metal binding to NTA (blue) and EGTA (red) plotted against those for EDTA. Black line shows the diagonal (x = y). NTA data are fitted to a linear function with a slope of 0.68 and a y intercept of +3 kJ mol–1 and the EGTA data are fitted to a second-order polynomial to guide the eye. Binding affinities are given in Table S1.

Computational Methods

MD Simulations and Thermodynamic Integration

The nonbonded potentials used by the AMBER ff14SB force field are described by eqs and 2, which contain 12-6 and 12-6-4LJ terms, respectively:[17,18]Essentially, A/r12 is a repulsive term that prevents the attraction from becoming too strong at short distances, B/r6 is an attractive term derived from London dispersion forces, C/r4 accounts for ion-induced dipoles and the e2Q/r Coulombic term accounts for the electrostatic interactions between the atoms.[9,19]

Thermodynamic integration (TI) was used to calculate Ca2+, Mg2+, Y3+, and La3+ ion HFEs and binding energies, according to the standard thermodynamic cycle shown in Figure S2. All simulations were carried out using AMBER18 with the ff14SB force field and TIP3P solvation extended at least 15 Å from the solute.[17,18,20] Chelator parameters were determined using the Antechamber package from AmberTools19[17] with charges calculated using the AM1-BCC charge model.[21] Charges are given in Table S2. Particle Mesh Ewald was used for long-range nonbonded electrostatic interactions with a cut off distance of 10 Å.[22] The SHAKE algorithm was used to constraint all covalent bonds that involve hydrogen atoms and a 2 fs timestep was used.[23] The temperature was maintained at 300 K using the Langevin thermostat with a 2.0 ps–1 collision frequency and the pressure at 1 atm using the Berendsen barostat.[17,24] After initial TI testing it was observed that fewer intermediate λ states were required to reach a converged ΔGVdW value than for the ΔGele+pol term, so 9 λ steps were used to calculate ΔGVdW and 12 for ΔGele+pol with Gaussian integration.[17] The simulations to determine ΔGele+pol were run in triplicate to ensure sufficient sampling and to allow error estimation. The lengths of EM, NVT, NPT and MD sampling per intermediate λ state are outlined in Table S3. More details on the simulation parameters and protocols are available in the Supporting Information.

We tested three existing parameter sets: 12-6LJ HFE, 12-6LJ CM/IOD, and the 12-6-4LJ standard parameters. The 12-6LJ HFE parameters should provide accurate ion HFEs but are known to produce IODs that are shorter than experimental values by an average of 0.27 and 0.29 Å for divalent and trivalent ions, respectively.[10] The 12-6LJ CM parameters should provide a compromise between IOD and HFE, but these are not available for trivalent ions. The 12-6-4LJ parameters were reported to reproduce accurate IODs and HFEs.[9] Even though 12-6LJ parameters are expected to produce poorer results, they were included as a benchmark.

In order to calculate average IODs and CNs, data were taken from the TI simulations with λ = 0.00922, which reproduced the known IOD and CN of fully charged Ca2+ and Mg2+ ions in water (Table S4). For full-length LanM, a separate 10 ns MD simulation was used. IODs and CNs were calculated from radial distribution functions (RDFs), which were calculated to a resolution of 0.01 Å using VMD.[25,26] The IODs were taken as the peak of a quadratic fit applied to ±0.1 Å of the first peak of the RDF, and CN (taken as the number of atoms within the first coordination sphere around the metal ion) was calculated by integrating the RDF from 0 to the first minimum.

Experimental Benchmarking

The 12-6LJ and 12-6-4LJ parameters were obtained using experimental HFE values, so we used the same values for our comparisons.[6,9,10,27] EDTA and NTA metal ion stability constants (K1) were taken from NIST reports.[28,29] For EGTA, no NIST values were found, so values from the Dojindo metal chelate affinity report were taken.[30] Dissociation constants (Kd) for LanM and CaM were taken from refs (14,31−33). Note that these will give average binding energies for the four EF loops in each protein. All experimental values were converted to binding energies in kJ mol–1 at 298 K (Table S1) assuming:

Results and Discussion

Initially, EDTA binding energies for Ca2+, Mg2+, Y3+ and La3+ were computed by TI using three established parameter sets:[6,9,10] 12-6LJ HFE and 12-6-4LJ for all metals and 12-6LJ CM for Ca2+, Mg2+, or 12-6LJ IOD for Y3+ and La3+ (12-6LJ CM parameters not available for trivalent ions). The difference between the computed EDTA-metal ion binding energies (ΔGbsim) and the experimental values in Figure (ΔGbexp) is shown in Figure . HFEs are given in Table S5 and absolute binding energies are given in Table S6. In most cases, all three parameter sets significantly overestimate the binding energy (producing more negative values), although they produced relatively accurate binding energies for Ca2+. While 12-6-4LJ produced the best results, there is significant room available for improvement as this showed errors ranging from +4.8 to −52.8 kJ mol–1.

Figure 3

Left, the structures of EDTA, NTA, and EGTA and below, snapshots taken from MD simulations of the Ca2+ bound chelators performed with 12-6-4LJ default C coefficients. Right, the difference between experimental and computed EDTA-metal ion binding energies obtained using four sets of parameters. 12-6LJ CM was used for Ca2+ and Mg2+ while 12-6LJ IOD was used for Y3+ and La3+. “12-6-4LJ Ch-BE” denotes the parameters generated during this study.

To increase the accuracy of the TI calculations for the binding energy between chelators and metal ions, we chose to re-parameterize specific C coefficients of the 12-6-4LJ parameter set (eq ) as this allows the modification of the pairwise interactions between the metal ion and specific atom types (i.e., ligating atoms) without affecting other interatomic interactions and the HFE.[13] EDTA, NTA, and EGTA can coordinate metal ions using both carboxylate groups and tertiary nitrogen atoms (Figure ). As the ligating groups are in similar chemical environments in the three chelators, we reasoned that a single set of ligating oxygen and nitrogen parameters can be shared between these molecules. EGTA also contains ether groups, which we did not specifically parameterize in this study. C coefficients for the carboxylate oxygen and the tertiary nitrogen were parameterized using TI simulations of EDTA and NTA using the 12-6-4LJ standard parameters. EGTA was then used for benchmarking.

Figure (top panel) shows ΔGbsim values for Ca2+ coordination by EDTA and NTA, which were computed at different ligating oxygen and nitrogen C values while holding all other C parameters to their default 12–6-4-LJ values, that is, C(O) was set to default when varying Cij(N) and vice versa. Data for other metal ions are shown in Figures S4–S6. These data all show linear dependences of ΔGbsim on C(O), so they were fitted to a linear function, with the gradients, m given in Table . While there is not a strong dependence of ΔGbsim on C(N) for some metals, these data were also fitted to a linear function for consistency. Binding free energies were also computed with C values for both ligating oxygen and nitrogen atoms set to 0 and these ΔGb(C = 0,0) values are also given in Table . In principle, the following relationship should then describe the experimental binding energy:

Figure 4

Top, Computed binding energies for the chelation of Ca2+ by EDTA (blue) and NTA (red) as a function of the ligating oxygen and nitrogen 12-6-4LJ C values. Filled symbols are for oxygen and open symbols for nitrogen. Solid lines are linear fits to the data and the horizontal dashed lines are the experimental ΔGbexp values from Table . Default C values are labeled. Bottom, 2D plot of ΔGbsim for the chelation of Ca2+ by EDTA versus the ligating oxygen and nitrogen 12-6-4LJ C values. Contour corresponding to the ΔGbexp value is shown in black and the new Ch-BE values from Table are shown as red dashed lines.

Table 1

Metal Binding Energies (kJ mol–1) for EDTA and NTA and Gradient (m) Values (kJ mol–1C–1) for the Data in Figures , S4–S6

	ΔGbexp[28,29]		ΔGbsim(Cij default), ΔGbsim(Cij = 0,0)a		m(O)		m(N)
metal	EDTA	NTA	EDTA	NTA	EDTA	NTA	EDTA	NTA
Ca2+	–60.8	–37.5	–56.0, −41.3 (−24.0)	–32.7, −17.6	–0.593	–0.442	–0.177	–0.081
Mg2+	–50.1	–30.6	–88.1, −37.6 (−21.8)	–84.9, −52.5	–1.028	–0.842	–0.124	+0.077
Y3+	–103.2	–65.5	–156.0, −81.2 (−82.6)	–148.9, −107.1	–0.752	–0.691	–0.058	+0.008
La3+	–87.6	–59.1	–123.3, −85.0 (−79.0)	–115.4, −86.7	–0.632	–0.491	–0.056	–0.028

Computed using standard parameters and with C for ligating oxygen and nitrogen groups both set to 0. The EDTA values in (parenthesis) are the ΔGbsim′(C = 0,0) values determined using eq .

Table 2

Default and Newly Derived “Ch-BE” 12-6-4LJ Pairwise Oxygen and Nitrogen C Coefficients (kcal mol–1 Å–4) for Ca2+, Mg2+, Y3+, and La3+ Ligation by EDTA and NTA

metal	default		Ch-BE
	O	N	O	N
Ca2+	34.4	65.9	29.2	110.6
Mg2+	52.4	100.3	12.3	126.1
Y3+	85.1	163.0	14.9	161.9
La3+	59.9	114.7	9.5	46.3

A 2D plot of ΔGbsim versus C(O) and Cij(N) for the chelation of Ca2+ by EDTA is shown in the bottom panel of Figure . The ΔGbexp value corresponds to a contour in this plot, so to find a unique solution for C(O) and C(N),eq must be solved simultaneously for two or more different chelators with different binding energies and m(O) and m(N) values. However, when calculating the binding energies for EDTA and NTA with a range of C values, it became apparent that ΔGbsim(C = 0,0) values for NTA with Mg2+, Y3+ and La3+ were significantly larger than the experimental values (Table ). Since larger C coefficients increase the binding energy, this meant that in order to satisfy eq at least one of the C coefficients would have to become negative, which is not physically realistic. Further, some computed NTA binding energies are larger than the EDTA values for the same metal ion, in opposition to the experimental values (Table , Figure ). This suggests that ΔGbsim(C = 0,0) are unreliable as computed. Instead, an equivalent binding free energy ΔGbsim′(C = 0,0) was determined by back-extrapolating from ΔGbsim(C) determined using default 12-6-4-LJ values:

The EDTA ΔGb′(C = 0,0) values were determined using eq and are given in Table . Equivalent values can be determined for NTA by scaling the EDTA values by the ratio of the experimental binding constants:

This approach ensures that the ΔGbsim′(C = 0,0) values for different chelators follow the experimental binding energy trends.

As seen in Table , for Mg2+ and Y3+, the NTA m(N) values are not physically realistic (non-negative). Visual inspection of the relevant MD trajectories suggested that these metal ions did not interact with the tertiary nitrogen atoms in the same way as EDTA (Figure S7), which likely leads to underestimation of these m(N) values.[34] From the EDTA and NTA Ca2+ simulations, we noted that the ratios of the m(O) and m(N) values are approximately 3/4 and 1/2 for oxygen and nitrogen, respectively, thus mirroring the maximum CNs available (Figure ).[35] Based on this observation, we can modify eq to scale the EDTA m(O) and m(N) values by 3/4 and 1/2, respectively, in order to describe the NTA binding energy. Making use of this approximation, and by substituting eqs and 6 into eq , we can determine unique C values using only the EDTA TI simulation data and the experimental binding energies for EDTA and NTA in eq . The resulting ligating oxygen and nitrogen C coefficients are given in Table and the Ca2+ values are plotted on the contour plot in Figure to show they intersect at ΔGbexp. If different CNs are observed[36] or suspected, eq can be readily modified to use ratios other than 3/4 and 1/2. Examples of the analysis performed with alternative ratios is given in the Supporting Information.

Validation of New Cij Coefficients

Our new C coefficients for the 12-6-4LJ nonbonded potential, which we denote “12-6-4LJ Ch-BE,” resulted in significant improvements to the TI-computed binding energies compared to the 12-6LJ CM/IOD, 12-6LJ HFE, and default 12-6-4LJ parameter sets. These are shown in Figure and the absolute binding energy values are given in Table S6. Our new binding energies are all within 8.6 kJ mol–1 (∼2 kcal mol–1) of the experimental values, while the previously best results using the 12-6-4LJ default set differed by up to ∼50 kJ mol–1 from the experimental data. The average absolute error was reduced to 5.6 kJ mol–1 compared to 32.8 kJ mol–1 for the default 12-6-4LJ set.

To ensure that the new C coefficients were not overfitted for EDTA, they were also tested on NTA (the NTA TI values were not used for parametrization; see eq ) and a similar chelator EGTA (Figure ). The difference in computed and experimental binding energies are shown in Figure and the absolute binding energy values are given in Table S7. In each case, a clear improvement was observed. Overall, the “12-6-4LJ Ch-BE” parameter set reduced the absolute average binding energy errors from 49.7 to 27.1 kJ mol–1 for NTA and from 30.9 to 17.2 kJ mol–1 for EGTA, relative to the default 12-6-4LJ parameters. This represents ∼1.8-fold increase in accuracy for both NTA and EGTA, although it should be noted that there was significant variability in the EGTA values with Y3+ using both forcefields. This arose due to EGTA adopting another conformation in some simulations, which introduces additional hydrogen bonding between EGTA and solvent molecules. More detailed analysis of these results is discussed in the Supporting Information.

Figure 5

Difference between computed NTA (A) and EGTA (B) (ΔGbsim) and experimental (ΔGbexp) binding energies obtained using TI with default 12-6-4LJ (green) and 12-6-4-LJ Ch-BE (yellow) O and N C coefficients for the ligating oxygen and nitrogen groups. Energies are given below or above the bar and the red and blue crosses in (B) indicate the approximate energies for EGTA Y3+ in two different chelator conformations (Figure S8).

Next, we compared the 12-6-4LJ Ch-BE parameters to the 12-6-4LJ default and 12-6LJ CM/IOD parameter using TI simulations of metal ion binding to the CaM and LanM 12-residue EF1 loop peptides (Figures , S1). The 12-6LJ HFE set was not selected as it was the worst performing set of parameters with EDTA (Figure and Table S6), and it is known to produce inaccurate IODs.[6,10] Although 12-6LJ CM/IOD produced poorer results in EDTA-metal ion simulations than 12-6-4LJ (Figure ) it was still included here as 12-6LJ potentials are more widely used. As the EF-hand peptides do not possess tertiary amines, only the new metal ion-oxygen pairwise C coefficients were used. These were used for the carboxylate oxygens, while the coordinating oxygen atom from the backbone of Thr 7 retained the default C coefficient as it is part of an amide bond which has not yet been re-parameterized. The length of the constant-pressure equilibration (NPT) and MD trajectories were optimized using LanM EF1 peptide with Ca2+ (Figure S10). Two nanoseconds of NPT and 5 ns of MD sampling were used as they give reasonably converged binding energies. The van der Waals contribution to the chelator binding energies is relatively small and does not vary significantly; for Ca2+, ΔGVdW = 9.5 kJ mol–1 in EDTA, 9.2 kJ mol–1 in NTA, 9.1 kJ mol–1 in EGTA and 9.1 kJ mol–1 in LanM EF1. Consequently, we chose to use a fixed −ΔGVdW contribution of 9 kJ mol–1 for all metal ions in both the LanM EF1 and CaM EF1 systems to reduce computational cost. The ΔGele+pol portion of the TI calculations was run in triplicate for each system to improve sampling and allow error estimation. The van der Waals contribution (−ΔGVdW) to ΔGbsim was not included in the error analysis. The average binding energies for LanM EF1 and CaM EF1 are shown in Table .

Table 3

Mean Binding Energies (kJ mol–1) for Metal Binding to the LanM and CaM EF1 Loop Peptides

metal	exp[14,31−33]	12-6LJ CM/IODa	12-6-4LJ	12-6-4LJ Ch-BE
LanM
Ca2+	–18.0	–67.4	–52.3	–55.5 ± 2.7
Mg2+	NAb	–119.2	–84.3	–39.9 ± 8.6
Y3+	–61.4	–205.6	–151.7	–73.2 ± 14.1
La3+	–64.3	–159.6	–127.9	–95.8 ± 7.8
CaM
Ca2+	–33.9	–53.6	–71.5	–54.1 ± 1.9
Mg2+	–22.8	–121.1	–100.3	–57.8 ± 3.2
Y3+	NAb	–189.5	–148.3	–74.7 ± 7.5
La3+	–45.3	–152.3	–126.7	–89.1 ± 3.2

12-6LJ CM was used for divalent and 12-6LJ IOD for trivalent ions.

Not available.

As seen in Table , all simulations significantly overestimate the binding energies. As the experimental values are the average affinity for the 3–4 EF loop binding sites in LanM and CaM, some of this error may reflect differences in affinity between the different EF loops in each protein. The 12-6LJ CM/IOD parameter sets performed the worst in each case, except for Ca2+ binding to CaM EF1. The 12-6-4LJ Ch-BE parameters performed significantly better than the default 12-6-4LJ parameters, but still significantly overestimate the binding energy for all cases. The average errors in binding energy are 84.6, 63.1 and 29.0 kJ mol–1 for 12-6LJ CM/IOD, 12-6-4LJ, and 12-6-4LJ Ch-BE, respectively, and there is a 2.9-fold and 2.2-fold increase in accuracy for the 12-6-4LJ Ch-BE parameter set compared to 12-6LJ CM/IOD and 12-6-4LJ, respectively. The observed standard deviations in ΔGbsim values (between triplicate ΔGele+pol simulations) are relatively low, rarely exceeding 9 kJ mol–1, indicating that these simulations give consistent results when starting from the same input geometry.

The 12-6-4LJ Ch-BE parameters were the only parameters that reproduced the experimentally observed order of binding affinities of LanM EF1 with three metal ions: La3+ > Y3+ > Ca2+.[14] Although there are no reported LanM EF-Mg2+ binding energies, the data in Table predict that the affinity is weaker than that for Ca2+. This is usually the case for CaM EF-hands.[31,32] For CaM EF1, all three parameter sets failed to predict the correct order of affinities for the three metal ions: La3+ > Ca2+ > Mg2+ (Table ). Instead, they predict the same binding affinity order of: La3+ > Mg2+ > Ca2+. Nevertheless, the 12-6-4LJ Ch-BE parameters performed significantly better than the other parameter sets. There is no experimental data for Y3+ binding to CaM, but the predicted binding energy is 14.4 kJ mol–1 weaker than that for La3+, suggesting CaM EF1 has a binding preference for La3+.

To gain better insight into EF-hand metal ion coordination, the IODs and CNs were calculated using the same three parameter sets and these values are given in Table . Additionally, to compare the metal ion coordination of isolated EF-hands to the EF-hand motifs in a full-length protein and to further test the 12-6-4LJ Ch-BE parameters, 10 ns MD simulations on full-length LanM were performed to determine the IODs and CNs for each of the EF1, EF2 and EF3 binding sites with bound Ca2+, Mg2+, Y3+, and La3+ (Tables and S8). In general, the new 12-6-4LJ Ch-BE parameters produced very similar results to the existing 12-6-4LJ parameters.

Table 4

Ion-Oxygen Distances (IODs in Å) and Coordination Numbers (CNs) of Metal Ions in Full-Length LanM and the EF1 Loop Peptides of LanM and CaM

	12-6LJ CM/IODa		12-6-4LJ		12–6-4LJ Ch-BEb
metal	CN	IOD	CN	IOD	CN	IOD
LanM EF1 peptide
Ca2+	8.3–8.7c	2.43	8.1	2.39	8.0–8.5c	2.40
Mg2+	6.0	1.94	7.0	2.08	6.0–6.7c	2.07
Y3+	9.0	2.28	9.0	2.30	9.0–9.3c	2.31
La3+	10.0	2.45	10.0	2.50	10.0	2.49
LanM EF1 in full-length LanM
Ca2+	7.7–7.8c	2.41	8.0–8.1c	2.39	7.8–8.0c	2.38
Mg2+	6.0	1.94	6.5–6.8c	2.07	6.1–6.3c	2.06
Y3+	9.0	2.29	9.0	2.29	9.0	2.29
La3+	9.9	2.44	10.0	2.46	10.0	2.48
CaM EF1 peptide
Ca2+	8.2	2.42	8.0	2.40	7.7–8.0c	2.40
Mg2+	6.0	1.94	7.0	2.09	6.0	2.05
Y3+	9.0	2.29	9.1	2.31	9.0	2.30
La3+	10.0	2.45	10.0	2.48	10.0	2.49

12-6LJ CM for divalent and 12-6LJ IOD for trivalent ions.

IODs were calculated as the average from triplicate simulations, whereas for CNs, if different values were obtained between simulations, the full observed range is shown.

No clear minima were observed so the range is given.

As seen in Table , Ca2+ and Mg2+ ion coordination in LanM and CaM EF-hands is not represented by a single binding mode as the CN values are noninteger. Y3+ and La3+ ions displayed consistent CN numbers in both LanM EF1 and CaM EF1 being ∼9 and ∼10, respectively. Inspection of the simulations revealed that some coordinating residues were fluctuating between monodentate and bidentate ligating modes (Figure S11). This is likely to be physically realistic (e.g., as is observed in CaM, PDB 1CFF). However, if a fixed integer CN is required while retaining the ability to switch ligands this could be achieved by applying the cationic dummy atom model (CDAM) in combination with re-parameterization of the ligating atom C terms. The CDAM partitions the metal ion charge between itself and dummy atoms surrounding the central metal ion in a predefined geometry, which interacts with the surrounding atoms and can freely exchange ligands.[37] Recently, CDAM was combined with 12-6-4LJ potential for Mg2+, Fe3+, Al3+, and Cr3+ ions in a predefined octahedral geometry and was shown to reproduce HFE, IOD, and CN values in a water solution.[38] Alternatively, it has recently been suggested that modified angle and torsion parameters on ligands combined with the double decoupling method[39] might improve the TI simulation performance for transition metal binding to coordination complexes and proteins.[40] This approach should be applicable to the study of metal binding to CaM and LanM and in future work it would be interesting to see if this improves the computed binding energies and/or CN values.

The computed IODs from the Ca2+, Mg2+, Y3+, and La3+ simulations with the EF1 peptides and full-length LaM do not significantly differ between the 12-6-4LJ and 12-6-4LJ Ch-BE parameters. These IODs varied within ±0.02 Å in cases where CN remained the same (Table ). The 12-6LJ CM/IOD parameter set usually gave rise to lower IODs for Mg2+, Y3+, and La3+, while the Ca2+ IODs were relatively similar to those computed with the 12-6-4LJ forcefield. All three parameter sets produced relatively accurate average IODs for Y3+ in LanM EF1 (∼2.3 Å) compared to experimental structures in PDB 6MI5 (2.2–2.4 Å depending on the ligating group).[15] The Ca2+ IODs were also consistent with experimental CaM structures. Experimental IODs of 2.42 Å for monodentate and 2.41 Å for bidentate ligands in EF1 of PDB 1CLL are in good agreement with the 12-6-4LJ Ch-BE average IODs of 2.40 Å for two bidentate and four monodentate ligands.[16] For Mg2+, the 12-6-4LJ Ch-BE average IOD of 2.05 Å is comparable with the experimental monodentate ligation distance of 2.11 Å in EF1 of PDB 3UCW.[41] No crystal structures of LanM or CaM with La3+ are available for comparison.

Conclusions

We have demonstrated a method of generating C parameters for 12-6-4LJ MM forcefields that allows parametrization of specific ligating groups in order to tune binding energies computed by TI. The new C coefficients were tested on a series of chelators: EDTA, NTA, EGTA, and EF1 loop peptides from LanM and CaM proteins and showed significant improvements in computed binding energies relative to existing forcefields. The new parameters also produce CN and IOD values that are in good agreement with experimental values. The parametrization method should be extensible to a range of other systems and could be readily adapted to tune properties other than binding energies.