Force field independent metal parameters using a nonbonded dummy model.

The cationic dummy atom approach provides a powerful nonbonded description for a range of alkaline-earth and transition-metal centers, capturing both structural and electrostatic effects. In this work we refine existing literature parameters for octahedrally coordinated Mn(2+), Zn(2+), Mg(2+), and Ca(2+), as well as providing new parameters for Ni(2+), Co(2+), and Fe(2+). In all the cases, we are able to reproduce both M(2+)-O distances and experimental solvation free energies, which has not been achieved to date for transition metals using any other model. The parameters have also been tested using two different water models and show consistent performance. Therefore, our parameters are easily transferable to any force field that describes nonbonded interactions using Coulomb and Lennard-Jones potentials. Finally, we demonstrate the stability of our parameters in both the human and Escherichia coli variants of the enzyme glyoxalase I as showcase systems, as both enzymes are active with a range of transition metals. The parameters presented in this work provide a valuable resource for the molecular simulation community, as they extend the range of metal ions that can be studied using classical approaches, while also providing a starting point for subsequent parametrization of new metal centers.

Introduction

From organometallic catalysis to biology, metal ions are ubiquitous in both natural and synthetic processes. Catalytic metal centers are found in all different classes of enzymes, accounting for 44% of oxidoreductases, 40% of transferases, 39% of hydrolases, 36% of lyases, 36% of isomerases, and 59% of ligases.[1] These metal centers can play a variety of roles during the catalytic step, which include activating nucleophiles,[2] acting as redox centers,[3] and facilitating optimal positioning of the reacting system.[4] Metal ions also appear to play an important role in determining enzyme specificity, as demonstrated by a number of studies showing that metal substitutions can not only change the observed catalytic activity but even generate completely new activities in an enzyme.[5,6] Therefore, understanding the factors that determine metal binding specificity and their function within the context of enzyme reactivity is a complex problem.

Over the past decades, the availability of more accurate experimental and computational techniques[7−10] has helped to elucidate various structural and electronic aspects of the role of metals in biomolecular systems. However, while providing valuable information, each of these techniques has its own limitations. From an experimental point of view, correct assignment of the identity of the metal ion and its coordination pattern still remain a significant challenge.[11] On the other hand, while computational approaches such as molecular dynamics and Monte Carlo simulations have made it possible to study metal ions not only in solution,[12−15] but also in biomolecular systems,[16−19] there are often still difficulties with obtaining stable and physically meaningful descriptions of metal centers in biomolecular simulations.[19] Ideally, one would want to move to a full quantum mechanical (QM) description of the metal ion, as is being increasingly done in a range of QM/MM studies of biochemical reactions.[9,16,20] Nevertheless, this is not a trivial issue, as, despite significant advances in this area, there are still substantial errors in reproduction of physicochemical properties associated with current quantum mechanical treatments of metal centers.[21,22] Additionally, the very high computational cost makes a QM description untenable with increasing system size, particularly if one wants to perform substantial configurational sampling in free energy calculations.

A number of strategies have been adopted in order to incorporate metal ions into classical simulations for various force fields. These can be broadly classified into three distinct strategies: the nonbonded soft-sphere model,[14,15,23,24] the covalent bond approach,[25−28] and the dummy-model approach.[17,29−32] The simplest of these is the nonbonded soft-sphere model, in which the metal–ligand interactions are simply described through electrostatic and van der Waals potentials. This approach has been successfully used to describe alkali and alkaline-earth ions,[14,24,32] and also modified to include the effects of polarization and metal-to-ligand charge transfer on Zn2+ systems.[33] Nevertheless, it appears to be inadequate when it comes to more complex situations such as systems containing multinuclear metal centers with closely located metal ions,[31] or for the correct treatment of transition metals.[32] In the latter case, the challenge becomes to obtain a parameter set that can simultaneously reproduce both solvation free energies and metal–water distances.[32,34]

On the other hand, covalent or bonded approaches, which have been widely used to model, for instance, Zn-metalloenzymes,[27,35] suffer from the fact that they include predefined covalent bonds between the metal and ligands, thus not allowing for ligand exchange and/or interconversion between different coordination geometries.[27] Additionally, this approach includes challenges such as the large number of parameters to be optimized (which also makes it highly system specific), dependence on the choice of internal coordinates, and double counting of the electrostatic and Lennard-Jones interactions.[36]

The final strategy, namely the dummy-model approach, provides a promising solution to these problems.[17,32] In this approach, the metal center is described by a set of cationic dummy atoms connected around a central atom in the specific coordination geometry to be attained (Figure 1). The dummy model was specifically introduced to solve the problems with modeling transition-metal ions, where manifestations of crystal field effects lead to a more complicated pattern of solvation energies for the transition-metal ions compared to alkali and alkaline-earth ions.[37] This makes it challenging to obtain both the correct solvation free energy and metal–oxygen distances using a standard soft-sphere model.[34] The fact that this is a nonbonded model means that while the positions of the dummies are normally adjusted to the metal coordination in the relevant binding site, this is a transferable model that can be used in sites with coordination geometries other than that which was originally intended, and where changes in the ligand coordination occur, without the need for any further parametrization (see, e.g., ref (17) and the Ca2+ model presented in this work).

Figure 1

(A) Schematic representation of the dummy model used in this work. (B) Representative active site of human GlxI where the active site metal has been replaced by an octahedral dummy model to represent Zn2+. The central atom and the dummy atoms are shown in gray and white, respectively. This figure has been adapted from ref (38). Reproduced by permission of the PCCP Owner Societies.

In the present work, we extend the scope of the original model[31,32] to three biologically relevant transition metals, namely, Ni2+, Co2+, and Fe2+, as well as refining the existing parameters in the literature for Mn2+,[32] Zn2+,[17] Mg2+,[31] and Ca2+.[39] In all cases, our parameters are able not only to reproduce the relevant radial distribution functions (RDF) in aqueous solution (i.e., the metal–oxygen distances) but also to consistently reproduce experimental solvation free energies. We have tested our model using two different popular water models (SPC[40] and TIP3P[41]) and show that our parameters are easily transferable to any force field that describes nonpolar interactions using a Lennard-Jones potential. Finally, we demonstrate our parameters “in action” by modeling metal substitution in the active site of the human and E. coli variants of glyoxalase I. We believe that the parameters presented in this work provide a valuable resource for the molecular simulation community, while also providing a foundation for the subsequent parametrization of other metal centers.

Theoretical Background

Octahedral Dummy Model

As our starting point for the present work, we have used the octahedral dummy atom model originally presented by Åqvist and Warshel[17,32] in 1990. The original octahedral dummy model consists of six particles with fractional positive charges, henceforth referred to as “dummy atoms”, placed around a central particle in an octahedral geometry (see Figure 1). Each of the dummy particles possesses a charge of +δ, while the center of the system possesses a charge n–6δ. Therefore, the total complex retains the net charge of n+ possessed by the metal center of interest. Such charge delocalization away from the center is particularly advantageous in the case of systems with multinuclear metal centers,[31] as re-distributing the charge prevents excessive repulsion between the metal centers. For this reason, this model has been shown to also be useful for maintaining crystallographic structures in the case of alkaline-earth metals such as Ca2+[39] and Mg2+,[31] and the parameters for these metals have been further refined in the present work.

The geometry of the dummy complex itself is kept rigid by the imposition of large force constants on the metal–dummy bonds (see Table 1). However, as there are no bonds between the dummy complex and the surrounding ligands, overall rotation of the six-center frame about the nucleus is allowed, and no internal forces are associated with such rotation.[32] Here, it should be noted that as the dummy complex is allowed to freely rotate around the metal center, the coordination geometry is not constrained to the geometry of the dummy model used, but rather, the system is free to exchange ligands on the relevant time scale, provided the simulations are run for a sufficiently long time. An example of such flexible ligand coordination was demonstrated, for example, in the case of carbonic anhydrase,[17] where an octahedral dummy model was used to describe a zinc center with tetrahedral coordination.

Table 1

Force Field Parameters for the Dummy Models Used in This Work and Addtional Parameters for Each Atom

Force Field Parameters
bond typea	Kb	r0
M–D	800.0	0.900
Di–Dj≠i	800.0	1.273

U = K(b – b0)2, where Kb is in kcal mol–1 Å–2 and r0 is in Å.

Uθ = (1/2)kθ(θ – θ0)2, where Kθ is in kcal mol–1 rad–2 and θ 0 is in degrees.

Where Lennard-Jones parameters are given in units of [kcal1/2 mol–1/2 Å–6] for Ai and [kcal1/2 mol–1/2 Å–3 ] for Bi.

Overview of the Parametrization Procedure

Geometric parameters that successfully maintain the octahedral conformation of the dummy model have been previously presented in the literature,[17,31,32,39] and they have been used, with some modifications, in the present work. Additionally, as in previous work,[31] negligible van der Waals parameters were used on the dummy atoms (Table 1); therefore, the remaining parameters that required adjustment are the van der Waals parameters on the metal center, as well as the distribution of charges between the metal center and the peripheral dummy atoms.

The derivation of the relevant van der Waals parameters for the metal centers under study has been done using the corresponding metal–aquo complexes, [M(H2O)], where x indicates the number of water molecules in the first coordination sphere and n+ indicates the total charge of the system. Metal–aquo complexes have been extensively studied using both experimental[37,42−44] and theoretical[45−48] approaches. Additionally, once the metal center has been rigorously parametrized in aqueous solution, it is then possible to transfer the parameters to a different environment such as an enzyme active site.

In atomistic force fields, the simplest description of the nonbonded interactions between atoms involves an electrostatic term, expressed by a Coulomb potential, as well as a van der Waals term, expressed by a Lennard-Jones potential[15] (note, however, that more complex functions to describe these interactions are being increasingly used in the literature; see, for instance, refs (49−51)). The potential function, U, used in this work has the following form:where A and B are the geometric Lennard-Jones parameters for the interaction between atoms i and j. The Lennard-Jones parameters are defined per atom type as A and B (Table 1) and are combined using the geometric rule: A = AA and B = BB, where A = A1/2 and B = B1/2.

In the present work, the van der Waals coefficients A and B were systematically refined in order to fit simulated properties to known experimental values. The total charge of the metal was distributed on both the central and dummy atoms following Warshel’s model for Zn2+ and Mg2+[17,31] (although other charge distributions have also been used in the literature[32,39]). Note that, with rigorous parametrization, either approach should work, provided that the relevant nonbonded parameters are fitted carefully and the system is tested for both stability and its ability to reproduce relevant observables.

In order to validate our parameter sets, we aimed to reproduce experimentally observed solvation free energies and M2+–O distances, as well as testing the stability of our parameters in selected enzyme active sites. These properties were chosen as they allow us not only to validate structural parameters for the specific model but also to validate our electrostatic treatments. Calculations of solvation free energies provide one of the most direct benchmarks for this, as they directly quantify electrostatic (ion–dipole) interactions of the metal with its surroundings.[52,53] Clearly, because the dummy model still uses a limited parameter set, there exist more than one combination of parameters that can reproduce individual observables; however, by iteratively refining our parameters to consistently reproduce a range of relevant observables, we obtain a more robust and reliable parameter set.

Evaluating Solvation Free Energies

The solvation free energy is calculated using the following expression:The first term, ΔGMFEP, refers to the free energy of charging the ion in water. The second term in eq 2, ΔGBorn, is a correction term that accounts for the error introduced by use of a finite interaction cutoff radius for the electrostatic interactions. This correction can be estimated from the Born formula (in kcal/mol), as was done in ref (13):where Q is the net charge of the solute, and RBorn is the radius of the cavity in the macroscopic medium with dielectric constant ε(T) in which the charge is embedded.[24] ΔGcav refers to the “cavitation energy”, which corresponds to the cost of creating a cavity in the solvent for the solute. This value is expected to be in the region of +2.5 kcal/mol, as discussed in refs (13 and 24). Finally, ΔGcorrstd.state corresponds to the (−1.85 kcal/mol) correction to the experimentally cited solvation free energy associated with moving the metal ion from the gas phase (standard state of 1 atm) to solution (standard state of 1 mol/L). These last two terms cancel out to within ∼1 kcal/mol, and their effect on the total calculated energetics is therefore expected to be negligible.

ΔGMFEP is obtained from a standard free energy perturbation (FEP) simulation, mapping from Q = 0 to n+ in n discrete steps which represent intermediates between the initial and final states:[15]Here, ⟨···⟩ represents an ensemble average on the mixed mapping potential (U):where U0 and U are the initial and final states, respectively; and λ is a mapping parameter that changes from 0 to 1 in fixed increments during the simulation.

Now while accurately reproducing experimental solvation free energies is a very important benchmark for electrostatics, the challenge in this case is that there can be significant variation in the results obtained from different experiments;[54−56] this issue has also been commented on by other authors.[57−60] For example, even in the simple case of Mg2+, Marcus,[55] Rosseinsky,[56] and Noyes[54] provide solvation free energies of −437.4, −455.5, and −452.2 kcal/mol, respectively. As can be seen from Table 2, Noyes and Rosseinsky both estimate similar solvation free energies, whereas Marcus cites substantially different values due to the way the solvation free energy of the proton is treated. Considering these potentially large variations, we have decided, for consistency, to remain with using the values that are extensively tabulated by Noyes[54] as our frame of reference, although it should be pointed out that despite the large absolute values of these deviations, they are only a smaller percentage of the total solvation free energies (which are all in the range of several hundred kilocalories per mole).

Table 2

Comparison of Experimental Literature Values for the Experimental Solvation Free Energies of the Different Metal Ions Studied in This Worka

	ΔGhyd, experimental
M2+	Noyes[54]	Marcus[55]	Rosseinsky[56]
Mg	–454.2	–437.4	–455.5
Mn	–436.4	–420.7	–437.8
Fe	–451.8	–439.8	–456.4
Co	–481.0	–457.7	–479.5
Ni	–492.8	–473.2	–494.2
Zn	–483.3	–494.7	–484.6
Ca	–379.5	–359.7	–380.8

Data are based on values provided by Noyes,[54] Marcus,[55] and Rosseinsky.[56] All values are in kcal/mol.

Simulation Setup

Parametrization of the Dummy Model for Different Metals

As outlined in section II, our starting point is an octahedral dummy model, using existing geometric parameters available in the literature.[31] These parameters were slightly modified in order to obtain a tighter dummy model (Table 1). The resulting metal dummy model was then complexed with six water molecules and immersed in an 18 Å water sphere, with both SPC[40] and TIP3P[41] water models being used for comparison purposes. To reproduce the experimental density and polarization of water molecules close to the sphere boundary, radial and polarization restraints were used according to the SCAAS algorithm as implemented in the program Q. The SHAKE[61] procedure was applied to all solvent molecules, and the nonbondend pair list was updated every 30 steps. A weak harmonic restraint of 0.5 kcal/(mol·Å2) was applied to the dummy model to keep it close to the center of the sphere. The time step of the simulation was set to 1 fs, a nonbonded cutoff of 10 Å was used, and electrostatic interactions beyond this cutoff were treated with the local reaction field (LRF)[62] method. No cutoffs were applied to the dummy model during the FEP run. The temperature was controlled using the Berendsen algorithm.[63] All simulations were performed with the Q simulation package[64] (version 5.0).

Once solvated, the system was first subjected to a short initial molecular dynamics (MD) equilibration of 5 ps at 20 K, 5 ps at 100 K, and 50 ps at 300 K, in order to relax the water sphere around the dummy. Subsequently, FEP calculations were performed in 101 λ steps. At each step, the system was simulated for 50 ps and potential energies were saved every 0.005 ps. In the calculation of the energy ensemble, the first 5% of each trajectory was discarded as equilibration and the energy was estimated as an average of a FEP calculation in both forward and backward directions. In order to obtain statistically meaningful data and assess the convergence of the calculations, in all cases, five replicas were obtained using different random number generator seeds. Finally, an independent 1 ns MD simulation was carried out at 300 K in order to obtain data to calculate the corresponding radial distribution function and metal coordination number, during which data were collected every 0.5 ps for analysis, discarding the first 10% of the simulation as equilibration.

Molecular Dynamics Simulations of the Glyoxalase I Variants

Crystal structures of both the human[65] and E. coli(66) variants of glyoxalase I (GlxI) were obtained from the Protein Data Bank[67] (PDB IDs 1QIP(65) and 1F9Z,[66] respectively). For the simulation of the human variant, the crystal structure used was obtained in complex with the product analogue S-p-nitrobenzyloxycarbonylglutathione (NBC-GSH) at 1.72 Å resolution. In this system Gln33, Glu99, Glu172, His126, and two water molecules coordinate the native zinc ion. For the MD simulations, the product analogue was removed from the active site.

For the simulations of the E. coli variant of GlxI, the native E. coli Ni2+–GlxI structure bound protein (1.50 Å resolution) was used as a starting point for simulations with all relevant metal ions. Note that crystal structures of this enzyme in complex with Co2+, Cd2+, and Zn2+ have also been reported from the same authors; however, these structures are at lower resolution than that of native enzyme.[66] In the Ni2+-complexed structure four protein residues (His5 and Glu56 from one monomer and His74 and Glu122 from the other) and two water molecules are coordinated to the metal. Note that, for both the human and E. coli variants, the first seven residues from each chain were not included, as they were not accurately traced from the electron density maps. Additionally, one of the two water molecules coordinating the metal center was replaced by a hydroxide ion, due to the expected low pKas of transition metal coordinating water molecules.

The MD simulations of the protein systems were performed using the Q[64] simulation package and the OPLS-AA[68] force field. The system was solvated using a spherical droplet of TIP3P[41] water molecules with a 24 Å radius centered between the two subunits for both human and E. coli variants of glyoxalase (allowing us to capture both active sites in our simulation sphere). In order to prevent Glu residues from artificially doubly coordinating to the catalytic metal center due to the identical charge on both carboxylate oxygens used in most force fields, the charges of the oxygen atoms of coordinating Glu side chains were modified to −0.918/–0.750. This captures some of the true polarization of these residues in the “real” system, and the degree of polarization of each residue was based on examining charge distributions obtained from calculations using model systems involving acetate bound to a metal center coordinated by five extra water molecules, averaged over all metals considered in this work. The atomic charges were fit to the electrostatic potential at points selected according to the Merz–Singh–Kollman scheme[69,70] using the B3LYP density functional[71−73] and 6-31G* basis set, as well as implicit solvation using the polarizable continuum model (PCM).[74] These calculations were performed using Gaussian 09.[75] Once the enzyme had been prepared for simulation, the relevant dummy model was placed into each active site such that the central atom overlaid with the metal center in the original crystal structure. After this, the system was heated from 1 to 300 K in a stepwise manner with initial random velocities taken from a Maxwell–Boltzmann distribution (increasing step size from 0.1 to 1.0 fs). The systems were then equilibrated for 20 ns using a 1 fs time step and a weak harmonic restraint of 0.1 kcal/(mol·Å2) on the heavy atoms. Net charges were assigned to ionizible residues located within 18 Å of the center of the sphere. The first 10% of each simulation was considered the equilibration time, and thus removed from the final analysis of the system. Unless stated otherwise, the other MD parameters were used as described in the previous section.

Results and Discussion

As our starting point in our parametrization protocol, we tested the performance of the currently available parameters for Mn2+,[32] Zn2+,[17] Mg2+,[31] and Ca2+.[39] The data from these simulations are shown in Table 3. While the original parameters perform well, reproducing M2+–O distances for Ca2+, Zn2+, and Mg2+, they showed larger differences for Mn2+. In terms of the free energy of solvation, the parameters for Zn2+ provide slightly lower values compared to the experiment (by 4 kcal/mol), while for Mn2+, Ca2+, and Mg2+ values of 4, 24, and 30 kcal/mol more negative are obtained, respectively. Finally, we also examined the performance of a recent dummy model of Ca2+,[39] as shown in Table 3. Note as an aside that comparison to standard soft-sphere models is complicated for the transition metals presented in this work, as it is not possible to simultaneously reproduce the complicated solvation patterns of transition metals and the physical M2+–O distance with the same parameter set (see the discussions in refs (32 and 34)). However, for comparison, we calculated solvation free energies of the standard soft-sphere models for Ca2+,[15] Mg2+,[15] and Zn2+[23] using the TIP3P water model (see Supporting Information Table S1 and Figure S1).

Table 3

Comparison of Calculated and Observed Solvation Free Energies (ΔGhyd, kcal/mol), and Ion–Water Oxygen Distances (M2+–O, Å) Employing Parameters Developed by Åqvist and Warshel[17,32] for Mn2+ and Zn2+ Dummy Models, by Oelschlaeger et al.[31] for Mg2+ Dummy Model, and by Saxena and Sept[39] for Ca2+ Dummy Model. For Comparison, the Results Using Soft-sphere Models are also Presented in Table S1.a

								TIP3P		experimental
	e	AiM	BiM	mM	AiD	BiD	mD	ΔGhyd	M2+–Oa	ΔGhyd[54]	M2+–O[52]
Mnb	–0.1	145	25.0	48.94	0.00	0.00	1.0	–440.9 ± 0.1	2.14 ± 0.05	–436.4	2.20
Znb	–1.0	136	41.0	59.38	0.00	0.00	1.0	–479.5 ± 0.1	2.07 ± 0.04	–483.3	2.08
Mgc	–1.0	70.0	41.0	6.3	0.05	0.00	3.0	–484.9 ± 0.1	2.08 ± 0.03	–454.2	2.10
Cad	0.0	233.2	35.5	33.1	0.05	0.00	1.0	–408.6 ± 0.2	2.32 ± 0.02	–379.5	2.39–2.46[76,77]

All values are averages and standard deviations over five trajectories, as outlined in the main text. M–O distances for all of the water molecules bound to the metal were monitored along the simulation. Only for calcium, which shows a rapid water exchange, the M–O distance was directly taken from the peak of the RDF (see Supporting Information Figure S1).

Kb = 1600 (kcal mol−1Å−2) and Kθ = 250 (kcal mol−1rad−2) and no bond between dummies.

Kb = 640 (kcal mol−1Å−2) and Kθ = 55 (kcal mol−1rad−2).

Kb = 540 (mol−1Å−2) and Kθ =55 (kcal mol−1rad−2) and no bond between dummies.

The parameters used for the calculations in Table 3 were subsequently refined, in order to reproduce both M2+–O distances and solvation free energies, and used as starting points for the parametrization of the other metals. As mentioned before, the existing geometric parameters available in the literature for Mg2+[31] were slightly modified for the new dummy models in order to obtain a tighter dummy model with better structural agreement with experiment while maintaining reasonable thermodynamic properties (Table 1).

The calculated solvation free energy values as well as the M2+–O distances for the new dummy models are shown in Table 4. The corresponding M2+–O radial distribution functions (g(r)), as well as the coordination numbers for the first coordination shell obtained from integration over the first M2+–O peak (N[g(r)]), are shown in Figure 2.

Table 4

Comparison of Calculated and Observed Solvation Free Energies (ΔGhyd, kcal/mol), and Ion–Water Oxygen Distances (M2+–O, Å) Using Our Parameter Set from Table 1 for Different Metalsa,b

	TIP3P		SPC		experimental
	ΔGhyd	M2+–O	ΔGhyd	M2+–O	ΔGhyd[54]	M2+–O[52]
Fe	–451.9 ± 0.2	2.13 ± 0.04	–450.0 ± 0.1	2.14 ± 0.04	–451.8	2.12
Ni	–492.7 ± 0.1	2.06 ± 0.03	–490.8 ± 0.2	2.07 ± 0.03	–492.8	2.06
Co	–480.5 ± 0.1	2.08 ± 0.03	–478.7 ± 0.2	2.09 ± 0.03	–481.0	2.08
Zn	–483.4 ± 0.1	2.08 ± 0.03	–481.4 ± 0.2	2.09 ± 0.03	–483.3	2.08
Mn	–436.9 ± 0.2	2.19 ± 0.03	–433.1 ± 0.1	2.20 ± 0.04	–436.4	2.20
Mg	–454.4 ± 0.1	2.12 ± 0.04	–452.3 ± 0.1	2.13 ± 0.04	–454.2	2.10
Ca	–379.9 ± 0.2	2.38 ± 0.02	–377.9 ± 0.2	2.39 ± 0.01	–379.5	2.39–2.46[76,77]

Calculations were performed using both SPC[40] and TIP3P[41] water models.

All values are averages and standard deviations over five trajectories, as outlined in the main text. M–O distances for all of the water molecules bound to the metal were monitored along the simulation. Only for calcium, which shows a rapid water exchange, the M–O distance was directly taken from the peak of the RDF (see Supporting Information Figure S1).

Figure 2

Radial distribution functions (g(r), y1-axis) and coordination number (N[g(r)], y2-axis) corresponding to the first hydration shells of (A) Ni2+, (B) Co2+, (C) Mg2+, (D) Mn2+, (E) Zn2+, and (F) Fe2+, obtained as outlined in the main text using the dummy-model parameters presented in Table 1. In all cases, the M2+–O g(r) are represented by solid lines and the N[g(r)] by dashed lines.

We would like to remind the reader that, as outlined in section II, it is difficult to select the appropriate metal solvation free energy against which to parametrize our systems due to the large deviations in the experimental values reported by different workers. For consistency between our different systems, we keep here to the data presented by Noyes in ref (54), which includes thermodynamic parameters for a wide range of metal centers, thus capturing the relative effect of the different metals. As can be seen from Table 4, in all cases, we obtain extremely good agreement with experimental values using the TIP3P water model, which is the model for which our systems were originally parametrized. Upon changing the water model from TIP3P to SPC we still are within the range of experimentally measured values; however, a systematic deviation of ∼2 kcal/mol is obtained in the calculated solvation free energy for all metals. This error corresponds to less than 1% of the total value. On the other hand, this shift to lower values is also in line with the change in charge distribution when moving from TIP3 to SPC water model, with the SPC water molecules being slightly less polarized than their TIP3P counterparts (oxygen partial charges for SPC and TIP3P are 0.820 and 0.834, respectively). Therefore, our parameters provide a good transferrable starting point, which could be fine-tuned if necessary for a different water model.

Despite the slight deviation in solvation free energy for the SPC model, the calculated M2+–O distances lie within the range of observed values, with an estimated standard error of 0.04 Å. Additionally, for all of the systems examined in this work, except calcium, integration of the RDF gives a coordination number of 6, in agreement with the existence of an octahedral arrangement of water molecules around the dummy model. In the case of Ca2+, the integration gives a coordination number of 7, and, as discussed below, experimental values fall in the range of 6–8 for this metal,[45,77] further highlighting the geometric flexibility of our dummy model. Here, it is important to emphasize that this octahedral arrangement was not user-defined, but rather, the nonbonded parameters provided in Table 1 are sufficient to attain this configuration both in aqueous solution and enzyme active sites after just a few picoseconds of relaxation from any arbitrary starting position of the dummy model. Therefore, the water molecules surrounding the metal automatically reorganizes in order to create the correct coordination sphere around the metal.

We would also like to note that because the model presented here is completely nonbonded, this should allow for ligand exchange around the metals which is important if one wants to study the long time scale dynamical behavior of metal ions in biological systems. However, due to the long time scales of ligand exchange for the systems presented here, with the exception of calcium (which will be described later) no ligand exchange was observed around the first coordination shell of the metal, and the metal maintained stable coordination for the duration of the simulations, in agreement with experimental results which show that this process takes place on average in the microsecond time scale.[78]

For calcium a different behavior is observed. The coordination geometry for this metal is quite flexible and strongly influenced by the second coordination sphere, with coordination number ranging from six to eight[45,77] and Ca2+–O distances in the range of 2.39–2.46 Å depending on the coordination number.[76,77] For this model, the coordination number was found to be seven as is often the case in biological systems (see Figure S2) and water exchange, taking place via a substitution mechanism, was observed in the picosecond time scale, in agreement with experimental reports.[78,79] As mentioned above, recently a seven-coordinated Ca2+ dummy model was also presented by Saxena and Sept, which was successfully tested in the calbindin system.[39] In the present work, in addition to structural properties, we also provide better thermodynamic properties, extending the scope of this model to systems where the metal ion is directly involved in chemical reactivity and therefore where correct solvation of the metal ion is crucial.

Exploring the Effect of Metal Substitution in Different Gloxalase I Variants

In order to test the parameters presented in this work in an actual biological system, we have performed molecular dynamics simulations of our dummy models for different metal centers in the active site of both human and E. coli variants of glyoxalase I (GlxI, EC 4.4.1.5) as a representative system. GlxI is a member of the metalloglutathione transferase superfamily, which catalyzes the first of two reaction steps in the detoxification of cytotoxic methylglyoxal (MG) via the conversion of nonenzymatically produced HG-GSH hemithioacetals to S-d-lactoylglutathione (Figure 3), thus playing a critical detoxification role in cells.[80] This enzyme has an absolute requirement for metal ions, showing a completely different response to the binding of different metals, depending on the species from which the enzyme originates.[66] Despite the low sequence identity (36%), both enzymes are structurally related, with three of four metal ligands conserved (the fourth ligand in E. coli has been assigned to His5, replacing Gln33 in Homo sapiens (H. sapiens); Figure 4). GlxI variants from H. sapiens,[81]Saccharomyces cerevisiae (S. cerevisiae)[81] and Pseudomonas putida (P. putida)[82] are zinc-dependent enzymes. In these systems, replacement of the native zinc ion with other divalent metal ions (such as Mn2+ and Co2+) yields an enzyme that, while active, has generally slightly impaired catalytic activity. Here, interestingly, only Mg2+ can fully restore the enzyme’s catalytic activity.[83,84] In contrast, the GlxI variant from E. coli is completely inactive in the presence of Zn2+, but fully active in the presence of Ni2+ and partially active with Co2+, Cd2+, and Mn2+.[85] The fact that both variants from S. cerevisiae(81) and P. putida(82) are active with Zn2+ and also feature the replacement of Gln by His suggests that swapping this residue is not the only culprit for the critical difference in activity.[86,87] Computational studies on the human GlxI variant[88] have shown that the metal ion plays a key role in the stabilization of the enolate intermediate (Figure 3). Other experimental studies have also suggested that the role for the metal ion is activating bound water molecules for nucleophilic attack.[66,89]

Figure 3

Proposed reaction mechanism for glyoxalase I. The mechanism involves a base extracting a proton from the C1 atom of the hemithioacetal of glutathione followed by reprotonation at C2. This proposed mechanism is based on that from ref (90).

Figure 4

Superposition of the E. coli Ni2+–GlxI structure from (blue) on the H. sapiens GlxI Zn2+–GlxI (yellow). Two residues from each domain form the active site, which is situated in a barrel formed only on dimerization. The metal and its coordinating residues are shown in a ball and stick representation with the zinc colored yellow and nickel blue (top right). This figure was created from the atomic coordinates deposited as PDB entries 1QIP and 1F9Z and is partially adapted from ref (66).

In order to test our parameters, we have studied both the human and E. coli variants of GlxI by means of molecular dynamics simulations in order to explore not just the structural stability of the system but also potential structural rearrangements around the different metal centers. GlxI provides an ideal test system for our parameters, as it is an enzyme that shows some level of activity with almost all of the metal centers presented in this work. Therefore, we have performed 20 ns of molecular dynamics simulations with each metal center as outlined in section III.2, to show that the system is stable without the need for additional constraints or artificial bonds (the RMSD of the backbone atoms is shown in Supporting Information Figures S3 and S4 to illustrate this fact). Here, we have taken both the native forms of the human (Zn2+) and E. coli (Ni2+) enzymes, as well as replacing the catalytic metal centers with Mn2+, Co2+, and Mg2+ in the human form and Mn2+, Co2+ and Zn2+ in the E. coli variant. Because all metals ions under consideration can substantially lower the pKa of the water molecules bound to them (pKa: Ni2+ < Co2+ < Zn2+ < Mn2+ < Mg2+), this greatly increases the probability of having one water molecule in its deprotonated form. Therefore, we have modeled one of the two water molecules as a hydroxide ion, as outlined in the Simulation Setup (note that the cost of deprotonating a second water molecule in the presence of the additional negative charge of the hydroxide ion would be expected to be substantially higher).

E. coli Glyoxalase I Variant

In the Glx I variant from E. coli, the coordination sphere of the active site metal ion is composed of four protein residues (His5 and Glu56 from one monomer and His74 and Glu122 from the other) and two water molecules. Table 5 compares the calculated backbone RMSD between the time-averaged dynamics structure and the crystal for the different metals in the two active sites (distinguished by the subscripts A and B). The results for the different metal ions are generally very similar, with the two active sites showing almost identical results. As can be seen from the backbone RMSD values for the different metals, the protein structure is well preserved, and the octahedral conformation remains stable during the simulation time. Additionally, no ligand exchange is observed on these time scales.

Table 5

Time Averages of the Root Mean Square Deviation (Å) of the Protein Backbone Atoms (RMSDbackbone) and of the Metal-Binding Residues (RMSDmetal) of E. coli Glyoxalase for Different Metalsa

systemb	RMSDbackbone	RMSDmetal
MnA	0.44 ± 0.03	0.59 ± 0.07
MnB		0.55 ± 0.05
NiA	0.45 ± 0.03	0.57 ± 0.06
NiB		0.53 ± 0.05
CoA	0.46 ± 0.03	0.57 ± 0.06
CoB		0.53 ± 0.05
ZnA	0.45 ± 0.03	0.56 ± 0.06
ZnB		0.53 ± 0.05

RMSDs for the protein backbone have been calculated by taking into account only the atoms within 20 Å of the system center, i.e., those that are inside the solvent sphere and are not subject to any restraint.

Subscripts A and B refer to the metal centers in active sites A and B from the different monomers respectively.

Substitution of the different metal ions in the active site of GlxI (Table 6) does not change the overall structure of the protein, in agreement with structural analysis of different crystal structures of this enzyme in complex with other metal ions.[66] As expected, however, there are a number of subtle changes in the active site architecture upon metal substitution. Specifically, the ligand coordination distances from the Mn2+ center are larger than for Co2+ and Ni2+, which is also consistent with the trend in increasing metal radii as one moves across the series of metals studied (Mn2+ > Co2+ > Ni2+ ≈ Zn2+; see ref (71)). Additionally, in the case of Mn2+, longer average distances of about 2.3–2.4 Å are observed between the metal center and the two coordinating histidines (His74 and His131, located on the x–y plane and on the z-axis respectively; cf. Figure 5), while shorter metal–ligand distances are observed between the metal and the charged carboxylate groups of Glu122 and Glu182. Following from this, and as would be expected, the observed average M2+–O distance to the hydroxide ion is slightly shorter than that to the water molecule.

Table 6

Interatomic Distances in E. coli GlxI Calculated from Our Molecular Dynamics Simulationsa

system	ODx Glu122	ODx Glu182	His74	His131	OH	HOH
MnA	2.10 ± 0.03	2.10 ± 0.03	2.39 ± 0.08	2.31 ± 0.06	2.18 ± 0.04	2.23 ± 0.05
MnB	2.10 ± 0.03	2.10 ± 0.03	2.39 ± 0.08	2.32 ± 0.06	2.18 ± 0.03	2.23 ± 0.05
CoA	2.01 ± 0.03	2.01 ± 0.03	2.27 ± 0.09	2.19 ± 0.05	2.08 ± 0.03	2.11 ± 0.04
CoB	2.01 ± 0.03	2.01 ± 0.03	2.27 ± 0.09	2.19 ± 0.05	2.08 ± 0.03	2.11 ± 0.04
NiA	1.99 ± 0.03	1.99 ± 0.03	2.17 ± 0.05	2.13 ± 0.04	2.08 ± 0.03	2.08 ± 0.04
NiB	1.99 ± 0.03	1.99 ± 0.03	2.17 ± 0.05	2.13 ± 0.04	2.08 ± 0.03	2.08 ± 0.04
ZnA	2.01 ± 0.03	2.01 ± 0.03	2.25 ± 0.08	2.18 ± 0.05	2.08 ± 0.03	2.11 ± 0.04
ZnB	2.00 ± 0.03	2.01 ± 0.03	2.25 ± 0.08	2.18 ± 0.05	2.08 ± 0.03	2.11 ± 0.04
CoA,exp	2.2	2.1	2.4	2.3	2.4	2.3
CoB,exp	2.4	2.1	2.3	2.2	2.4	2.2
NiA,exp	2.1	2.1	2.3	2.2	2.2	2.1
NiB,exp	2.1	2.1	2.3	2.1	2.2	2.1

In the cases of Co2+ and Zn2+, the experimental distances were obtained from the corresponding M–GlxI crystal structure (PDB IDs, 1FA6 and 1FA5, respectively,[66] last four rows). Note that because these dimeric enzymes have two metal-binding sites, subscripts A and B refer to the two different sites, respectively.

Figure 5

Coordination sphere of the catalytic metal centers in the active site of the (A) E. coli and (B) H. sapiens GlxI variants, where the active site metal has been replaced by the octahedral dummy model. Shown here are models for the native Ni2+ and Zn2+ ions, respectively. The central atom and the dummy atoms are shown in blue/yellow and silver, respectively, and the surrounding ligands have been highlighted to show the stability of the metal coordination sphere after 20 ns of MD.

A similar trend is seen for all other metal ions examined, with overall distances contracted in line with the changing metal radii described above. In the two cases for which experimental distances were known from crystal structures (Co2+ and Ni2+, from PDB IDs 1FA6 and 1FA5, respectively,[66]), our calculated average distances are all within 0.1–0.2 Å of the experimental value. Interestingly, in the case of Co where there are differences in the M–O distances in the two active sites from the crystal structure, we were also able to reproduce this difference computationally. Finally, no substantial distortion is observed in the case of Zn2+, suggesting that the origin of the deactivation of this enzyme by zinc is simply not structural but rather more complex.

Human Glyoxalase I Variant

In the human GlxI variant, the metal coordination is virtually identical, with the exception of the fact that one of the histidines of the E. coli variant has been replaced by a glutamine (Gln26). As indicated in Tables 7 and 8, the initial distances between the Zn2+ ion and the atoms bound to it are around 2.0 Å, with the exception of one of the water molecules, for which the O–Zn2+ distance is 2.8 Å and for which a much lower density is found when analyzing the density map.

Table 7

Time Averages of the Root Mean Square Deviation (Å) of the Protein Backbone Atoms (RMSDbackbone) and of the Metal-Binding Residues (RMSDmetal) of Human GlxI for the Different Metalsa

systemb	RMSDbackbone	RMSDmetal
MnA	0.37 ± 0.02	0.74 ± 0.10
MnB		0.68 ± 0.06
CoA	0.38 ± 0.02	0.76 ± 0.12
CoB		0.68 ± 0.06
MgA	0.37 ± 0.02	0.71 ± 0.10
MgB		0.67 ± 0.06
ZnA	0.37 ± 0.02	0.63 ± 0.08
ZnB		0.66 ± 0.06

RMSD for the protein backbone have been calculated by taking into account only the atoms within 20 Å of the system center, i.e., those that are inside the solvent sphere and are not subject to any restraint.

Subscripts A and B refer to the metal centers in active sites A and B from the different monomers respectively.

Table 8

Interatomic Distances in Human GlxI Calculated from Molecular Dynamics Simulations (20 ns)a

systemb	O2 Gln33	O1 Glu99	His126	Glu172	HO(Oc)	HOH(OHc)
MnA	2.20 ± 0.05	2.11 ± 0.03	2.34 ± 0.07	2.10 ± 0.03	2.17 ± 0.03	2.22 ± 0.05
MnB	2.20 ± 0.05	2.10 ± 0.03	2.36 ± 0.07	2.11 ± 0.03	2.17 ± 0.03	2.22 ± 0.05
CoA	2.20 ± 0.05	2.11 ± 0.05	2.35 ± 0.07	2.11 ± 0.03	2.17 ± 0.03	2.22 ± 0.05
CoB	2.20 ± 0.05	2.09 ± 0.05	2.37 ± 0.08	2.09 ± 0.03	2.17 ± 0.03	2.22 ± 0.09
MgA	2.15 ± 0.05	2.05 ± 0.03	2.32 ± 0.09	2.04 ± 0.03	2.09 ± 0.03	2.15 ± 0.05
MgB	2.15 ± 0.05	2.03 ± 0.03	2.36 ± 0.11	2.04 ± 0.03	2.09 ± 0.03	2.15 ± 0.05
ZnA	2.07 ± 0.04	2.00 ± 0.03	2.21 ± 0.06	1.99 ± 0.03	2.07 ± 0.03	2.10 ± 0.04
ZnB	2.08 ± 0.04	1.99 ± 0.03	2.24 ± 0.08	2.01 ± 0.03	2.07 ± 0.03	2.10 ± 0.04
ZnA,exp[65]	2.1 (2.0)	2.0 (1.9)	2.0 (2.1)	2.0 (3.3)	2.6(2.1)	2.3 (2.1)
Zne,exp[65]	1.9 (2.1)	2.0 (1.9)	2.1 (2.1)	2.0 (3.3)	2.8(2.1)	2.2 (2.0)

Experimental distances were obtained from the corresponding Zn2+–GlxI crystal structure (PDB ID, 1QIP). The zinc ligand distances for the transition state analogue bound-protein structure (PDB ID, 1QIN) are shown in parentheses.b

Subscripts A and B after each metal refer to the relevant active site.

In the HIPC–GSH complex the oxygen atoms of the hydroxycarbamoyl group of the inhibitor replace the two oxygen atoms.

Similar to the E. coli variant, the human GlxI enzyme is also homodimeric, with both monomers being very close in structure to each other. This is particularly evident in the active site, where very similar interatomic distances are obtained from the two active sites. In this case the presence of the Mn2+ dummy model also leads to larger M–O distances throughout, compared to the other metals examined. The largest distance in all cases is observed for the interaction with His126 followed by Gln33, while the two Glu residues show a much tighter interaction, following the trend observed for the E. coli variant. These distances become shorter with Mg2+ and Zn2+. In this last case all distances are very similar, about 2.0 Å, and Gln33 shows a shorter distance compared to the analogous histidine in the E. coli variant.

As in the previous case, no major structural changes were observed upon metal substitution, and all of the metals studied for these variants showed a stable octahedral coordination during the simulation time. The main difference among the metals studied is found in the M2+–O/–N distances. Considering the effect of the metal center has been suggested to be the stabilization of the enolate intermediate,[88] these differences could already give some insight about the stabilization of this substrate when bound to the metal.

Conclusion

Metal ions are crucial to biology, fulfilling a range of chemical and structural functions that make life possible, but accurately describing metal centers in molecular simulations remains a significant bottleneck in computational chemistry and biology. Ideally, one would like to have a full quantum mechanical description of the metal center, particularly as this would allow one to correctly describe ligand-to-metal charge transfer, and changes of spin state in the case of transition metals. However, while complex metal centers are by now routinely modeled in computational organometallic chemistry[10] and in QM/MM studies of biological systems,[91] the high computational cost associated with using a density functional theory based description of the metal center makes use of such an approach untenable in long time scale biomolecular simulations. Additionally, while there have been substantial advances in classical force field based descriptions of metal ions, one still often has to resort to using either artificial restraints or bonds in order to maintain stable coordination geometry. While this is fully reasonable in shorter simulations where the metal center only plays a structural role, it does not allow for ligand exchange in long time scale simulations and also becomes problematic in studies of chemical reactivity, where a flexible coordination sphere around the metal can become essential.[92]

Åqvist and Warshel proposed a solution to this problem in 1990,[32] presenting a nonbonded octahedral dummy model to describe transition metals, which allowed one to capture fundamental properties of the metal center such as M2+–O distances and solvation free energies, without the need for any artificial constraints or bonds. This approach has been adopted by other workers since then[29,31,93,94] and extended to a range of metals with both octahedral and tetrahedral coordination geometries. However, the parameters presented have still been adjusted for specific force fields, and in the case of tetrahedral zinc[29,30,93] also specific coordination geometries without the consideration of solvation effects (although this is not possible in the case of a tetrahedral zinc model, which, while common in enzyme active sites,[95] has no counterpart in aqueous solution). We demonstrate that while a classical model will, per definition, lack a proper description of electronic properties, we can nevertheless reproduce important structural and physical properties such as metal–ligand coordination distances and solvation free energies. As a test of the reliability of our parameters to maintain correct coordination in enzyme simulations, we have also performed simulations of different dummy models in the active sites of the human and E. coli variants of GlxI, demonstrating that our models not only maintain stable coordination geometries, but also capture the subtle geometric changes that would be expected upon metal substitution without the need for any artificial bonds or constraints. We believe that the parameters described herein are an important and useful aid to the broader scientific community in the study of the role of metal ions in complex biological systems.