Hi!!

If you do not find in the literature will have to use some software ab initio to calculate these parameters.

I'll just mention something to clarify an issue - it's not "in GROMACS" that the parameters are defined; it's in a force field. The key to properly parametrizing the ion is to do so in a manner that is compatible with the force field you have chosen. QM is the way to go, but realize that transition metal ions dramatically impact their local environment and have complex coordination geometry and charge transfer effects that MM force fields generally do not reproduce very well. Several force fields have parameters for prosthetic groups containing iron, like heme and Fe-S clusters. You may want to look those up as a starting point.

Recent parameters for +2 metal cations:

Li P, Roberts BP, Chakravorty DK, Merz KM. 2013. Rational design

of particle mesh ewald compatible Lennard-Jones parameters for +2

metal cations in explicit solvent. J Chem Theory Comput 9(6):2733–

2748.

Please, remember about Gromacs versus above paper units! I have used this parameters for Fe2+ binding helical motif and TIP3P and got very good agreement of binding affinities with experiment conducted in parallel. We have recently published:

Zhao, Y., Marzinek, J. K., Bond, P. J., Chen, L., Li, Q., Mantalaris, A., Pistikopoulos, E. N., Noro, M. G., Han, L. and Lian, G. (2014), A Study on Fe2+ – α-Helical-Rich Keratin Complex Formation Using Isothermal Titration Calorimetry and Molecular Dynamics Simulation. J. Pharm. Sci.. doi: 10.1002/jps.23895

Or below:

Article A Study on Fe 2+ - α-Helical-Rich Keratin Complex Formation ...

Thank you Jan, this is very interesting.

Thanx everone, especially Jan for the very interesting paper.....

Congratulation for the paper........

To parameterize your force field of choice do some digging and find out how it was generated. For example, the OPLS force field (available in GROMACS) uses parameters from ab inito calculations, calculated at the B3LYP/6-31G(d) level of theory.

At a minimum you'll want to start there and generate your DFT model, but you may need to augment the method for the metal centre by applying an effective core potential basis set, such as LANL2DZ, in the event that you don't get the right coordination geometries out.

Your DFT model should look like the binding site you want to simulate. In general as well you're better off trying to simulate a static structure with MM methods (as opposed to exchangeable ligand binding at the metal or something similar). For example, your DFT model may be of a Fe2+ heme with one or two additional donor ligands in the axial positions, or it could be something simpler like two CH3S- ligands to represent deprotonated Cys residues and two 4-mmethylimidazole ligands representing coordinating His residues.

You'll want to optimize your geometry first, as this will provide you with bond lengths, bond angles and dihedrals (assuming they're not available from another technique). Be warned that protein crystallography is useful for providing useful starting coordinates for the amino acids but metal binding environments can suffer from non-physical distortions in the final structure. If EXAFS has been done on your protein, and it has been done correctly, then that should be the most accurate source of metal-to-ligand bond lengths to parameterize your force field.

With your optimized geometry you'll also want to perform a harmonic frequency calculation for two reasons 1) for metal sites there may be a number of low-lying local minimal or your geometry optimization may have found a saddle point, the frequency calculation will tell you if you're in a local minimum or not because the calculated frequencies will be positive (i.e. no imaginary/negative frequencies), 2) the software you use should also spit out force constants for each mode. You'll want a program like JMol (free) to visualize the frequency output and identify the mode that best represents the Fe-Ligand stretching mode you're after. You'll also likely have to scale the force constant to put it into units meaningful for the force field in GROMACS.

You will also want to calculate the partial atomic charges (for OPLS this requires the CHELPG method). For Fe you will likely also need to find a suitable radius to use as it may not be included in your program by default or it may not be appropriate to the oxidation state you want. For eg, Gaussian requires the use of pop=(chelpg,readatradii) for metal centres.

GROMACS has Lennard-Jones parameters for Fe2+, but they may or may not be what you want, but they'll do in a pinch because your metal binding site is likely going to be surrounded by bound ligands anyway and the Lennard-Jones potential will only be a slight perturbation to those local interactions.

You will have to put all of this data into several places. Depending on your model, if you need to make non-standard amino acids or new ligands you'll need to make a new entry in the residue topology file, as well as specify the H-bonding in the hdb file, and perhaps also list new atoms/atomtypes in the bonded and non-bonded parameters files (these are located in a sub directory within the GROMACS directory, but it is recommended you copy them out into a local working directory where you can modify them). Gromacs looks if your current directory first for parameters before defaulting back to its usual place.

Finally you'll need to enter all of the bond, pair, angle, dihedral, etc. interactions. I do this in the .top file that is generated by GROMACS, and then I usually do stability testing by running the simulation for a while and closely monitoring the energy, structure, etc. Just because the simulation doesn't blow up doesn't mean you don't have a problem with how you parameterized your model.

It's been about a year or so since I built anything, so I may have missed something, but this is most of what you want to do. Have fun!!

I should hasten to add that if you have amino acids coordinating your metal you should also modify their parameters in your .top file that is generated for you, as their charge distribution will be very different when they are bound to a metal. Generally for things like His or deprotonated Cys residues I modify their charges all the way back to the beta carbon and its hydrogens, but I leave the backbone alone unless it is also directly involved in metal coordination.

You could use a semi-empirical or ab initio calculation.

The MOPAC software (semi-empirical) has included a keyword to work with large systems like proteins, or you can use the Quantum Mechanics/Molecular Mechanic (QM/MM) mix. Also the Fragment Molecular Orbital (FMO) is capable to do calculations for large systems using DFT, HF, etc. (it uses a "divide and conquer" technique).

Any of them will work much, much better than docking (even for ligand/protein).

Camps

@M. Pushie

Could you give a reference on CHELPG in OPLS?