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 . these metal centers 
 can play a variety of roles during the catalytic step , which include 
 activating nucleophiles , acting as redox 
 centers , and facilitating optimal positioning 
 of the reacting system . 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 . 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 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 . 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 , but also in biomolecular systems , there are often still difficulties 
 with obtaining stable and physically meaningful descriptions of metal 
 centers in biomolecular simulations . 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 . 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 . 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 , the covalent bond approach , and the dummy - model approach . 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 , and also modified to include 
 the effects of polarization and metal - to - ligand charge transfer on 
 zn systems . 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 , or for the correct 
 treatment of transition metals . in the 
 latter case , the challenge becomes to obtain a parameter set that 
 can simultaneously reproduce both solvation free energies and metal water 
 distances . on the other hand , covalent 
 or bonded approaches , which have been 
 widely used to model , for instance , zn - metalloenzymes , 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 . 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 . the final strategy , 
 namely the dummy - model approach , provides a 
 promising solution to these problems . 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 . this 
 makes it challenging to obtain both the correct solvation free energy 
 and metal  oxygen distances using a standard soft - sphere model . 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 ca model presented in this work ) . ( 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 
 zn . the central atom and the dummy atoms are shown in 
 gray and white , respectively . . reproduced by permission 
 of the pccp owner societies . in the present work , we extend the scope of the original 
 model to three biologically relevant transition 
 metals , namely , ni , co , and fe , as well as refining 
 the existing parameters in the literature for mn , zn , mg , and ca . 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 and tip3p ) 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 . as our starting 
 point for the present work , we have used the octahedral dummy atom 
 model originally presented by qvist and warshel 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 n6. 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 , 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 ca and mg , 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 . 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 , where an octahedral 
 dummy model was used to describe a zinc center with tetrahedral coordination . ub = kb(b  b0 ) , 
 where kb is in kcal mol  and r0 is 
 in  . 0 ) , where k is 
 in kcal mol rad and  0 is in degrees . where lennard - jones parameters are 
 given in units of [ kcal mol  ] for ai and [ kcal mol  ] 
 for bi . geometric parameters that successfully maintain the octahedral 
 conformation of the dummy model have been previously presented in 
 the literature , and they have been used , with some modifications , 
 in the present work . additionally , as in previous work , 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)x ] , 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 and theoretical 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 ( note , however , 
 that more complex functions to describe these interactions are being 
 increasingly used in the literature ; see , for instance , refs ( 4951 ) ) . the potential function , uij , used in this work has the following form:1where aij and bij are the 
 geometric lennard - jones parameters for the interaction between atoms i and j. the lennard - jones parameters are 
 defined per atom type as ai and bi ( table 1 ) and are combined using the geometric rule : aij = aiaj and bij = bibj , where ai = aii1/2 and bi = bii1/2 . in the present work , the van der waals coefficients ai and bi 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 zn and mg ( although other charge distributions have also been used in the 
 literature ) . 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 m  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 . 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 . the solvation free 
 energy is calculated using the following expression:2the first term , gmmfep , 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):3where qi 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 . 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 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 . gmmfep 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:4here , m represents an ensemble average on the mixed 
 mapping potential ( um):5where u0 and un are the initial and final states , respectively ; and m 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 ; this issue has also 
 been commented on by other authors . for example , even in the simple 
 case of mg , marcus , rosseinsky , and noyes 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 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 ) . as outlined in section ii , our 
 starting point is an octahedral dummy model , using existing geometric 
 parameters available in the literature . 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 and tip3p 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 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 ) 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 ) method . 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 . crystal structures of both the human and e. coli(66 ) variants of glyoxalase i ( glxi ) were obtained from the 
 protein data bank ( pdb ids 1qip(65 ) and 1f9z , 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 ni  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 co , cd , and zn have also been reported 
 from the same authors ; however , these structures are at lower resolution 
 than that of native enzyme . in the ni - 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 simulation package and the opls - aa force field . the system was solvated using a 
 spherical droplet of tip3p 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 kollman 
 scheme using the b3lyp density functional and 6 - 31 g * basis set , as well as implicit solvation using the polarizable 
 continuum model ( pcm ) . these calculations 
 were performed using gaussian 09 . 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 ) 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 . as our starting point in our parametrization 
 protocol , we tested 
 the performance of the currently available parameters for mn , zn , mg , and ca . the data from these simulations 
 are shown in table 3 . while the original parameters 
 perform well , reproducing m o distances for ca , zn , and mg , they showed larger 
 differences for mn . in terms of the free energy of solvation , the parameters for zn provide slightly lower values compared 
 to the experiment ( by 4 kcal / mol ) , while for mn , ca , and mg 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 ca , 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 m 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 ca , mg , and zn using the tip3p water model ( see supporting 
 information table s1 and figure s1 ) . all values are averages and standard 
 deviations over five trajectories , as outlined in the main text . 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 ) and k = 250 ( kcal molrad ) and no bond between dummies . kb = 
 640 ( kcal mol ) and k = 55 ( kcal molrad ) . kb = 
 540 ( mol ) and k = 55 ( kcal molrad ) and no bond between dummies . the parameters used for the calculations 
 in table 3 were subsequently refined , in order 
 to reproduce both m  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 mg 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 
 m  o distances for the new dummy models are shown 
 in 
 table 4 . o 
 radial distribution functions ( g(r ) ) , as well as the coordination numbers for the first coordination 
 shell obtained from integration over the first m o 
 peak ( n[g(r ) ] ) , 
 are shown in figure 2 . all values are averages and standard 
 deviations over five trajectories , as outlined in the main text . 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 ) . 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 m 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 ca , the integration gives a coordination number of 7 , and , as discussed 
 below , experimental values fall in the range of 68 for this 
 metal , 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 . radial distribution functions 
 ( g(r ) , y1-axis ) and coordination number ( n[g(r ) ] , y2-axis ) corresponding 
 to the first hydration shells of 
 ( a ) ni , ( b ) co , ( c ) mg , ( d ) 
 mn , ( e ) zn , and ( f ) fe , obtained 
 as outlined in the main text using the dummy - model parameters presented 
 in table 1 . in all cases , the m o g(r ) are represented 
 by solid lines and the n[g(r ) ] by dashed lines . for calcium a different behavior 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 and ca o 
 distances in the range of 2.392.46  depending on the 
 coordination number . 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 . as mentioned above , recently a seven - coordinated ca dummy model was also presented by saxena and sept , which was successfully 
 tested in the calbindin system . 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 . 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 . 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 . 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 , saccharomyces cerevisiae ( s. cerevisiae ) and pseudomonas putida ( p. putida ) are zinc - dependent enzymes . in these systems , 
 replacement of the native zinc ion with other divalent metal ions 
 ( such as mn and co ) yields an enzyme that , 
 while active , has generally slightly impaired catalytic activity . 
 here , interestingly , only mg can fully restore the enzyme s 
 catalytic activity . in contrast , the glxi variant 
 from e. coli is completely inactive in the presence 
 of zn , but fully active in the presence of ni and partially active with co , cd , and 
 mn . the fact that both variants 
 from s. cerevisiae(81 ) and p. putida(82 ) are active with zn 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 . computational studies 
 on the human glxi variant 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 . 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 ) . 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 
 ( zn ) and e. coli ( ni ) enzymes , 
 as well as replacing the catalytic metal centers with mn , co , and mg in the human form and mn , co and zn in the e. coli variant . because all metals ions under consideration can substantially 
 lower the pka of the water molecules bound 
 to them ( pka : ni < co < zn < mn < mg ) , 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 ) . superposition 
 of the e. coli ni  glxi 
 structure from ( blue ) on the h. sapiens glxi zn  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 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 . 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 . 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 mn center are larger than for co and ni , which is also consistent with the trend in 
 increasing metal radii as one moves across the series of metals studied 
 ( mn > co > ni  
 zn ; see ref ( 71 ) ) . additionally , in the case of mn , longer 
 average distances 
 of about 2.32.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 
 m  o distance to the hydroxide ion is slightly 
 shorter than that to the water molecule . in the cases of co and 
 zn glxi crystal 
 structure ( pdb ids , 1fa6 and 1fa5 , 
 respectively , 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 . 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 ( co and ni , from pdb ids 1fa6 and 1fa5 , 
 respectively , ) , our calculated average 
 distances are all within 0.10.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 zn , suggesting that the origin of the deactivation of this 
 enzyme by zinc is simply not structural but rather more complex . 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 . 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 . 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 zn ion and 
 the atoms bound to it are around 2.0  , with the exception of 
 one of the water molecules , for which the o zn distance is 2.8  and for which a much lower density is found 
 when analyzing the density map . 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 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 . all distances are very similar , 
 about 2.0  , and gln33 shows a shorter distance compared to the 
 analogous histidine in the e. coli variant . 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 . 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 m  o/n distances . considering the effect of 
 the metal center has been suggested to be the stabilization of the 
 enolate intermediate , these differences 
 could already give some insight about the stabilization of this substrate 
 when bound to the metal . experimental 
 distances were obtained 
 from the corresponding zn  glxi crystal structure 
 ( pdb i d , 1qip ) . the zinc ligand distances for the transition state analogue bound - protein 
 structure ( pdb i d , 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 . 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 and in qm / mm studies of biological systems , 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 . qvist and warshel proposed a solution to this problem 
 in 
 1990 , presenting a nonbonded octahedral 
 dummy model to describe transition metals , which allowed one to capture 
 fundamental properties of the metal center such as m 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 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 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 , 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 .
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 
 mn2 + , zn2 + , mg2 + , and ca2 + , as well as providing new parameters for ni2 + , co2 + , and fe2 + . 
 in all the cases , we are able to reproduce 
 both m2+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 .