Assessing the known and unknown unknowns: WYSI(N)WYG

Ken Dill and David Mobley from UCSF have a really nice review in Structure on computational modeling of protein-drug interactions and the problems inherent in the process. I would strongly recommend anyone interested in the challenges of calculating protein-drug binding to read the review, if not for anything else for the copious references provided. The holy grail of most such modeling is to accurately calculate the free energy of binding. For doing this we frequently start with a known structure of a protein-ligand complex. The main point that the authors emphasize is that when we are looking at a single protein-ligand complex, deduced either through crystallography or NMR, we are missing a lot of important things.

Perhaps the most important factor is entropy which is not at all obvious in a single structure. Typically both the protein and the ligand will populate several different conformations in solution. Both will have to pay complex entropic penalties to bind one another. The ligand strain energy (usually estimated at 2-3 kcal/mol for most ligands) also plays an important role. The desolvation cost for the ligand also can prominently figure. In addition both protein and ligand will have some residual entropy even in the bound state. As if this were not enough of a problem, much of the binding energy can come from the entropic gain that the release of water molecules from active sites engenders. Calculating all these entropies for protein, ligand and solvent is important for accurately calculating the free energy of protein-ligand binding. But there are few methods that can accomplish this complex task.

Among the methods reviewed in the article are most of the important methods used currently. Usually the tradeoff for each method is between cost and accuracy. Methods like docking are fast but inaccurate although they can work well on relatively rigid and well-parameterized systems. Docking also typically does not take protein motion and induced-fit effects into account. Slightly better methods are MM-PBSA or MM-GBSA which as the names indicate, combine docking poses with an implicit solvent model (PBSA or GBSA). Entropy and especially protein entropy is largely ignored, but since we are usually comparing similar ligands, such errors are expected to cancel. Going to more advanced techniques, relative free-energy calculations use molecular dynamics (MD) to try to map the detailed potential energy surfaces for both protein and ligand. Absolute free-energy perturbation calculations are perhaps the gold standard in calculating free energies but are hideously expensive. They work best for ligands that are simple.

There is clearly a long way to go before calculation of ∆Gs becomes a practical endeavor in the pharmaceutical industry. There are essentially two factors that contribute to the recalcitrance of the problem. The first factor as indicated is the sheer complexity of the problem; assessing the thermodynamic features of protein, ligand and solvent in multiple configurational and conformational states. The second problem is a problem inherent in nature; the sensitivity of the binding constant to the free energy. As iterated before, the all-holy relation ∆G = -RT ln K ensures that an error of even 1 kcal/mol in calculation will translate to a large error in the binding constant. The myriad complex factors noted above ensure that errors of 2-3 kcal/mol already constitute the limit of what the best methods can give us. Recall that an error of 3 kcal/mol means that you are dead and buried.

But we push on. One equal temper of heroic hearts. Made weak by time and fate, but strong in will. To strive, to seek, to find, and not to yield. At some point we will reach 1 kcal/mol. And then we will sail.

Reference:
Mobley, D., & Dill, K. (2009). Binding of Small-Molecule Ligands to Proteins: “What You See” Is Not Always “What You Get” Structure, 17 (4), 489-498 DOI: 10.1016/j.str.2009.02.010

Strain Energies in Ligand Binding: Round Two- Fight!

Or why to be wary of ligands in the PDB, force field energies, and anybody who tells you not to be wary of these two



One of the longstanding questions in protein-ligand binding has been; what is the energy penalty that a protein has to pay in order to bind a ligand? Another question is; what is the strain energy that a protein pays in order to bind the ligand? Contrary to what one might initially think, the two questions are not the same. Strain energy is the price paid to twist the conformation of the ligand into the binding conformation. Free energy of binding is the energy that the protein has to pay in addition to the strain energy in order to bind the ligand.

A few years ago, this question shot into the limelight because of a publication in J. Med. Chem. by Perola et al. from Vertex. The authors did a meticulous study of hundreds of ligands in their protein-bound complexes, some from the PDB and others proprietary. They used force fields to estimate the difference between the energy of the bound conformation of the ligands and the nearest local energy minimum conformation- the strain energy penalty. For most ligands, they obtained strain energies ranging from 2-5 kcal/mol. But what raised eyebrows was that for a rather significant minority of ligands, the strain energies seemed to be more than 10 kcal/mol, and for some they seemed to be up to 20 kcal/mol.

These are extremely high numbers. To understand why this is so, consider a fact that I have frequently emphasized on this blog; the concentration of a particular conformation in solution is virtually negligible if the free energy difference between it and a stable conformation is about only 3 kcal/mol. For a conformation to pay that much of an energy penalty in order to transform itself into the bound conformation would already be a stretch, considering its low concentration. For a conformation to pay an energy penalty of 20 kcal/mol does not make sense at all in this light, since such a conformation should be non-existent. Plus, think about the fact that hydrogen bonds usually contribute about 5 kcal/mol and that energy at room temperature is itself about 20 kcal/mol- significantly greater than the rotational barriers in most molecules- and this number for the strain energy penalty starts looking humungous. Where exactly would it come from?

Perola's paper generated a lot of buzz- a good thing. It was discussed by speakers at a conference in March last year that I attended. Now, a paper in J. Comp. Chem. seems to clear up the air a little. In a nutshell, the authors conclude that the strain energies they have measured seldom, if ever, surpass 2 kcal/mol. Needless to say, this is a huge difference compared to the earlier studies.

Why such a startling difference? It seems that as always, the answer strongly depends on the method and the data.

First of all, the PDB is not as flawless as people assume it is. Most people who are crystallizing protein-ligand complexes are first and foremost interested in the structure of the protein. They often do a poor job of fitting ligands to the electron density; Gerard Kleywegt of the University of Uppsala has done some marvelous work on detecting errors in PDB ligands, and his review on this should be a must-read for all scientists even marginally connected with crystallography. Because of poor fits, conformations of ligands in the electron densities in the PDB can be completely unrealistic and at the very least, brutally strained. Amides can be cis or non-planar, and more rarely planar aromatic rings can be deformed. There can be severe steric clashes which are not easily apparent. Quite naturally, such conformations when refined would lead to huge drops in energy. Therein lies the first source of the unrealistically large strain energy differences.

The second factor has to do with the vagaries and inadequacies of force fields, often unknown to crystallographers but known to experienced computational chemists. Force fields are quite poor at determining energies and their results are especially skewed by an overemphasis on electrostatic interactions which the force fields are ill-equipped to damp. Now consider what happens when a ligand in a PDB that has a positively and negatively charged group in it is optimized. If you relax it to the nearest local energy minimum, these two groups would instantly snap together and form a very strong ionic bond. This would lead to a huge overstabilization of the conformation, thus again giving the illusion of a large strain energy difference between the PDB conformation and the local minimum.

Finally, the devil is in the details. In doing the initial refinement of the conformation, the earlier study used a constraint called the flat-bottom potential in optimizing the PDB ligands in their bound state. However the flat-bottom potential, which extracts no penalties for atomic movement within a certain short distance and suddenly ramps up the penalty, is not physically realistic. A better method might be to use a harmonic potential which continuously and smoothy extracts a penalty proportional to atomic displacement.

The present study takes all these factors into account and also substitutes the force field results with some well-established quantum chemical energy determinations at the B3LYP/6-31G* level. They use this method to calculate the energies of bound and local energy minimum conformations. Secondly, they use a well-established continuum solvation model (PCM) as incorporated in the latest version of the Gaussian program to incorporate damping effects due to solvation. Thirdly as indicated above, they use the harmonic potential for optimization. Fourthly and most importantly, for the cases where the strain energy seems unusually high (and even there they set the bar quite high- anything greater than 2 kcal/mol), the authors closely investigate the relevant PDB entries and find that indeed, the ligands were not fit well into the electron density and had unrealistically strained conformations.

Once they tackled these problems, the strain energies all fell down to between 0.5 and 2 kcal/mol, which seems to be a realistic penalty that a conformation with a respectable concentration in solution could pay. There is now a second question; what is the maximum strain energy penalty that a ligand can pay to be transformed into the bound conformation? The authors are working on this question, and we will await their answer.

But this study reiterates two important lessons that should be remembered by anyone dealing with structure at all times:
1. Don't trust the PDB
2. Don't trust force field energies

Better still, as old Fox Mulder said, trust no one and nothing.

References:
1. Keith T. Butler, F. Javier Luque, Xavier Barril (2009). Toward accurate relative energy predictions of the bioactive conformation of drugs Journal of Computational Chemistry, 30 (4), 601-610 DOI: 10.1002/jcc.21087

2. Emanuele Perola, Paul S. Charifson (2004). Conformational Analysis of Drug-Like Molecules Bound to Proteins: An Extensive Study of Ligand Reorganization upon Binding Journal of Medicinal Chemistry, 47 (10), 2499-2510 DOI: 10.1021/jm030563w

3. A Davis, S Stgallay, G Kleywegt (2008). Limitations and lessons in the use of X-ray structural information in drug design Drug Discovery Today, 13 (19-20), 831-841 DOI: 10.1016/j.drudis.2008.06.006