Discussion U. Ryde Paper: Do Quantum Mechanical Energies Calculated for Small Models of Protein-Active Sites Converge?

We discussed the paper by Hu et al.

Bottom line:

• The paper discusses how the energy difference between an HIP state and an HID state is affected by the residues included in the respective model
• Let $E_B$ be the energy of the product state HIP for system $s$
• Let $E_A$ be the energy of the reactant state HID for system $s$
• Let $E_B'$ be the energy of the product state HIP for a system ($s$ + residue $i$), called $s'$
• Let $E_A'$ be the energy of the reactant state HID for a system ($s$ + residue $i$), called $s'$
• Reported is the $\Delta \Delta E := (E_B' - E_A') - (E_B - E_A)$
• Vacuum calculations are essentially worthless because even if the size of the system is set such that the difference in energy appears converged (above 15 residues, Fig. 2, 4, 5), still the energy differences obtained from such a system can change dramatically, once charged residues are added on top of this system (Fig. 7)

• This is less dramatic but still present when using continuum solvent models ($\epsilon =4$), Fig. 7
• It is now a question how much this effect could be further reduced with a higher dielectric constant
• It is however also possible that once the system reaches a certain size, the specific value of the dielectric constant does not affect the energy difference any longer (p.11796)
• One conclusion is that only including every position avoids convergence difficulties due to system size
• This however then most likely requires QM/MM approaches or linear scaling techniques
• If the study is to be carried out using only fragments of the protein back bone, charged residues should be included only together with their counter charge pair if available
• Charged residues in the interior of the protein are usually in counter-ion pairs
• Quantum only studies as the one presented however remain useful in mechanistic studies where one mechanism can be ruled out in favor for another mechanism because their respective activation energies are very different
• For such studies, it is recommended to include only a minimum model consisting of the groups immediately involved in the reaction

That's Classic!

I've very recently had the joy of setting up classical simulations for use in my research developing the next generation of force-fields. Since I've dealt with quantum mechanics and classical mechanics derived from quantum mechanics, my need for classical simulations on proteins and ligands has been void. Until now.

I only know of two tools to actually calculate Molecular Dynamics simulations on entire proteins: The Tinker package and the GROMACS package. I've tinkered(!) around with the first package and kept a good distance to the latter, but times have changed since I needed to include a ligand. GROMACS (and Tinker for that matter) are tuned to either proteins (in solution) or small molecules (in solution)*. If you want another molecule to use for the MD, you must first obtain the force-field parameters somehow. Since the most of us have no time for that, or that the next generation force-field is on its way but not quite there yet, what do you do?

To start from nothing and finish two days later, I did what any (in)sane scientist would do if no-one in the department has apparent experience with a particular piece of software: I Googled and I browsed the GROMACS "guides". The first is good if you know what to look for, and the latter is good to get inspiration, but it does not quite get you there since information always seem to be outdated just a bit.

Here is how I did it:

1. Find the structure of interest to you on the Protein Data Bank. Build your ligand into that system using your favorite build tool. You should now have a PDB file with a protein and a ligand.
2. Use SwissParam (more info below) to make a topology file (actually it is a force-field) of your ligand.
3. Follow the SwissParam Gromacs guide until step 6 (you must include this step!)
4. Follow the Protein-Ligand Complex guide by Justin Lemkul from the solvation step (it is good to actually read the entire guide). This will get you to make a production run MD of 1 ns of your own system.
5. Analyze your MD! This is out of scope for this post, but man was I happy when the MD actually successfully finished.

Note on SwissParam
The SwissParam tool takes a small molecule of interest and generates what is know in GROMACS-speak as a topology file, but in regular office-chatter, this is known as a force-field. Here is an exert from the main page of SwissParam

The data are derived from the Merck Molecular ForceField (MMFF). Dihedral angle terms are taken as is, while only the harmonic part of the bond, angle and improper terms are retained. Charges are taken from MMFF. Van der Waals parameters are taken from the closest atom type in CHARMM22.

The paper of SwissParam(doi:10.1002/jcc.21816) has been published so give it a read if you're really feeling geeky.

* I should note that if one chooses the MMFF 94 force field in Tinker, one should potentially be able to use it for everything, but it is a bit messy to actually set up and run not to mention that I found absolutely no way to automate the atom type-setting.