Purity in the atomic force-fields of molecular dynamics simulations

In the protein modeling community, there is a split between those who pursue a more "physics-based" approach, and those pursue a more "empirical-based" approach. Nothing highlights this better than the acrimonious debate between atomic force-fields and statistical potentials. These force-fields are used in everything from protein folding studies, ligand-binding studies and protein-protein interaction studies.

From a brutally pragmatic point of view, the successes in protein folding and homology modeling show that empirical (or statistical) force-fields have galloped far ahead. Still, the purists leer down on empirical force-fields, as they pray themselves to sleep at night that the test of time will expose empirical force-fields for the dirty statistic aberrations that they are.

The purists do have a point. We've arrived at a place where molecular dynamics (MD) has become such an easy-to-use tool that any biologist with a souped-up laptop can carry out MD simulations. The whole point of MD is to extend simulations into the unknown, and so one should use chemically pure and reliable force-fields.

Except they aren't. Chemically 'derived' force-fields are neither pure nor reliable. For those of you who haven't had the pleasure of diving into the guts of your typical atomic force-field, you will not have seen the morass of compromises and shady half-steps in the guts of every force-field. And once you've opened up the sausage, you may not want to take another bite.

What does Molecular Dynamics (MD) do?

Molecular Dynamics (MD) is a way of shoe-horning the chemical properties of a molecule into a Newtonian framework. Now, if you didn't know any better, you might think MD is a straightforward first order approximation of quantum mechanics. Although you can derive Newton's equations from Schrödinger's equation, such derivations are more sophistry than useful methodology.

The fact is, Schrödinger's equation is only a one-electron equation. It has taken decades and several Nobel prizes to build solvable quantum mechanics (QM) systems for multi-electron molecules. Trying to extract a Newtonian approximation in terms of MD out of the QM systems has turned out to be a process that one might charitably describe as ad-hoc. The essential approximations of MD are:

1. Bonds can't be broken and are approximated as springs
2. Electrostatic interactions to electron orbitals are to be approximated with point charges.
3. Pauli exclusion principle to be replaced by the van der Waals terms.
4. The torsion degree of freedom are modeled as Fourier terms

But the most important loss in MD is that molecules in MD can't make or break chemical bonds – real chemistry such as catalysis is not allowed to happen in MD. Instead, MD ends up to be useful in studying how large molecules, such as proteins and DNA, move around to get in place for chemical reactions. It's like watching a grid-iron game where you can see all the players running around for the play but the actual throw has been edited out.

How atomic force-fields are derived

Every MD package depends on an atomic force field, a collection of parameters for a bunch of generic quasi-atomic parameters that could be applied to any conceivable organic molecule. The principal type of parameters are: masses, bond lengths, bond angles, spring terms, atom types, partial charges, vDW terms, and torsion parameters.

Some of these parameters are relatively straightforward and brook little argument. Masses of atoms are taken from standard physics measurements. Bond lengths and angles are taken from crystal structures of small molecules. Spring terms are derived from the vibration frequencies observed in molecular spectra of small molecules. There's not a lot you can do with these, and they are pretty much similar across different force-fields.

Everything else, however, is up for grabs.

Watching how atom types are defined is a little like watching your little brother collect Magic the Gathering cards. He starts his hand with a few useful cards, but then spends the rest of his life collecting an endless number of new cards, each boasting more and more combinations. Every new release of a force-field sees an expansion of basic atom types, not only to allow the MD package to model more types of molecules, but also to fix problems found in standard molecules. To transfer these parameters to the actual molecule you want to model, you have to follow arcane recipes. It is a process guaranteed to make your hair fall out.

And then there are the cases where some atom-types are defined purely to save some calculations. For years force-fields were split between those that model hydrogen atoms separately or those that fatten up the attached heavy atom, then dropping the hydrogen atom altogether. This is the united atom approach, and it still infects some modern-day force-fields, especially with models for explicit water.

For the vDW parameters of carbon and hydrogen atoms, monte-carlo liquid simulations of ethanols and alcohols are used to fit the vDW to calorimetric measurements. Now, we can also use highly detailed liquid measurements, such as density curves of liquids. Others optimize the vDW parameters to reproduce the geometry of small molecule crystals. Still, there is plenty of scope as to whether you want to fit it to gas-phase or liquid-phase experiments, or in the case of crystals, solid-phase experimental data. And these can be taken from all sorts of different temperatures.

There are many different options in fitting partial charges. Early on, the AMBER guys developed the RESP algorithm, which basically generates an electrostatic potential from QM calculations for small molecules, then fits the partial charges onto the nuclei of the atoms that best reproduces the electrostatic potential. The AMBER guys were pretty excited because it offers a pipeline that, on the surface, appeared to be completely automatic. The mystery of course is to decide what conformation of the molecule that you want to use, and in what QM package. And then later how to transfer the partial charges to your generic atom types. As well, you have to decide how the partial charges derived from a chemical fragment must be truncated at the ends to maintain overall electronegativity, and to allow the cobbling of polymers from your fragments. But the most egregious error (as we shall see later with the Ramachandran plot) is that these partial charges are taken from QM calculations done in vacuum.

Hydrogen bonds have also proven to be a mine-field of parameterization. Completely different types force-fields have been tried, including ones with exponential terms, others with a distance parameter. A popular parameterization has been to leave out the hydrogen vDW parameters, allowing only the electrostatic interaction to provide the interaction. This was justified by the fact that the simulations displayed no material difference to the simulation with more complex parameterizations. Still the usefulness of hydrogen bonds have been severely curtailed by their inability to identify bifurcated hydrogen bonds and hydrogen bonds with lone-pair electrons.

And perhaps, the most abject performance of atomic force-fields are in the torsion terms, which I'll discuss next.

The incommensurability between force-fields and protein crystallography

There is a tension between organic chemists and protein crystallographers. In the world of chemistry, it's the partial charges and vDW parameters that are important. After all, it's what affects calculations of calorimetry, and the calculations of the crystal structures of small organic molecules. In the world of protein crystallography, however, it's the torsion angles that matter. Unfortunately, it's the chemists that have pulled the hardest in keeping force-fields pure, and this has probably the biggest obstacle to long-time MD simulations.

When deriving a protein structure from an electron density, the backbone conformation is the first thing you work out. Everything else is subordinate to that. The backbone conformation is described by the phi/psi torsion angles, and their 2D experimental distribution, known as the Ramachandran Plot, is one of the most intensely studied constructs in protein crystallography (Lovell et al 2003 Proteins 50:437; Laskowski et al 1993 J Appl Cryst 26:283; Kleywegt and Jones 1996 Structure 4:1395). I know because I've written 3 papers on it. Major errors in protein structure determination have been eradicated by carefully following the structure reconstruction in accordance with the known distribution of Ramachandran plot.

If you look at force-field development, you can see that the force-field community have more-or-less ignored the massive work done by the crystallography community. In study after study, they have singularly failed at reproducing the Ramachandran plot of crystal structures in their simulations of the protein backbone. The reason is that force-field developers have attempted to derive torsional terms of the backbone from QM models of the peptide fragments. Unfortunately, these QM models have mostly been done in vacuum. Only recently, has their been studies where QM studies have got a grip on the Ramachandran plot of the protein crystallographer. Hu et al (2003 Proteins 50:451) simulated a dipeptide QM model in a bath of MD water. Tsai et al (2009 J Phys Chem B 113:309) made very good calculations with a dipeptide QM calculation using a sophisticated solvent model. But alas, these are not the studies used to parameterize the phi/psi terms in atomic force-fields.

Not only is the Ramachandran plot important for crystallography, but it is absolutely vital for protein folding. When the protein is folded in its native state (the starting point for a lot of MD simulations), there are many mutually reinforcing interactions that maintain the folded structure. In the early stages of protein folding, however, these reinforcing interactions don't exist, and the backbone is 'unstructured'. The main determinant of its spatial orientation is the backbone conformational biases, i.e. the Ramachandran plot. If the Ramachandran plot is wrong in your model, you will have biased your simulation in the wrong geometrical space at the beginning (and most crucial part) of your simulation. This problem applies to any kind of simulation involving unstructured backbone: loop motions, hinge motions of domains, and protein unfolding intermediate studies.

Given the failure of reproduction of the Ramachandran plot, is it then any surprise that initial efforts in protein folding have been so terrible?

The only test that matters is protein folding

Although MD simulations is often used to simulate huge protein/DNA complexes, surprisingly little optimization is made for such simulations in the developments of atomic force-fields. In early force-field papers, the testing involves running MD for a few nanoseconds on a few small folded proteins, and showing that the proteins did not fall apart. Even in 2005, my then lab published a paper merely to show that several nanoseconds of replica-exchange of several small proteins were stable.

Simulating folded structures is hardly a challenge, as there are so much redundant interactions in a folded structure. The hard task is to fold a protein, arguably the most spectacular chemical reaction that does not involve the breaking or making of chemical bonds – and thus theoretically possible in MD. It's taken a long time to really test this, as the fastest proteins fold in microseconds. It's only in the last few years that various groups have been able to run MD simulations for the time long enough to fold these small proteins.

Probably the most embarrassing simulation paper I've had the pleasure of reading is one of the first microsecond MD folding paper of the WW domain (Reddolino et al. 2008 Biophys J 10:L75) using the CHARMM22/CMAP force-field. Yes, the protein did fold in microseconds, but no, it folded to the wrong structure. Happily, a more recent paper was able to rectify this (Shaw et al 2010 Science 350:341) using a modified AMBER ff99SB force field. Apart from showing that getting the right answer is the difference between a Science article and a Biophys J article, it is interesting to note that the Shaw group had to introduce a whole bunch of ad-hoc atom types and parameters based on fitting rotamer directly to the PDB.

In the end, atomic force-fields are statistical potentials too

Here we are in the world of microsecond MD simulations, and we have the ability to seriously test the veracity of the atomic force-fields. Such tests have forced the chemists who created the atomic force-fields to deal pragmatically with the short-comings, in particular the modeling the backbone Ramachandran plot.

Indeed, arguments have raged in the literature over the intrinsic superiority of pure chemistry (MD with pure atomic force-fields) over mongrel protein statistics (such David Baker's Rosetta package with its salad of statistical potentials). Yet, for the purposes of pragmatic protein folding studies, Rosetta has triumphed over all physics-based protein folding systems. Moreover, Rosetta actually encodes better chemistry insight then any existing atomic force-field. Rosetta includes a hydrogen-bonding term that accounts for lone-pair electrons and bifurcated hydrogen-bonds (Kortemme et al 2003 JMB 4:1239).

Nevertheless, the purists can claim solace in the hope that once computers can provide sufficient oomph, MD simulations with pure atomic force-fields will turn the tables. Unfortunately, as computers have improved, MD folding simulations have exposed structural flaws in atomic force-fields. The irony is that, to compensate, force-field developers have been forced to hand-code statistics-based terms in the force-fields. The approach in CHARMM22 is to add the CMAP term, which is a grid-based potential of the protein backbone, ultimately based on PDB statistics (Buck et al 2005 Biophysics J 90:L36). Similarly, for the Shaw group, ad-hoc terms were added to the AMBER force-field that fixed the worst deviations from known PDB rotamer distributions. Another direct introduction of statistics is the Yasara Nova force-field that directly optimizes the AMBER force-field against the refinement of an extensive bunch of crystal structures (Krieger et al 2002 Proteins 47:393).

It turns out that building a useful atomic force-field is every bit as dirty as building a statistical force-field. By adding statistical terms to atomic force-fields, the statistical barbarians have breached the gates, and the two have started to converge. Purity, alas, has become a dirty word.