Conquering the curse of morphine, one docked structure at a time

Morphine, one of the most celebrated molecules in history, also happens to be one of the deadliest. It belongs to the same family of opioid compounds that includes oxycodone and heroin. So does fentanyl, a weaponized form of which is rumored to have caused havoc and killed dozens in a terrorist standoff in Russia a few years ago. Potent painkillers which mitigate pain that is impossible to cure using traditional painkillers like ibuprofen, morphine and its cousins are also devils in disguise.

Constipation is the least concerning side effect of opioid use. Even in moderate doses they have a high risk of causing death from respiratory failure; in the jargon of drug science, their therapeutic index - the difference between a beneficial and a fatal dose - is very low. To get an idea of how potent these compounds are, it's sobering to know that in 2014, more than 28,000 deaths were caused by opioid overdoses, and half of these were because of prescription opioids; the addiction potential of these substances is frighteningly high. The continuing oxycodone epidemic in places like New Hampshire bears testimony to the hideous trail of death and misery that opioids can leave in their wake.

Quite unsurprisingly, finding a chemical agent that has the same unique painkiller profile as morphine without the insidious side effects has been a longstanding goal in drug discovery, almost a holy grail. It would be hard to overestimate the amount of futile resources that pharmaceutical companies and academic labs alike have spent in this endeavor. It was long thought that it might be impossible to decouple morphine's good effects from its bad ones, but studies in the last few decades have provided important insights into the molecular mechanism of morphine's action that may potentially help us accomplish this goal. Morphine binds to the so-called µ-opioid receptor (MOR), one of a family of several opioid receptors involved in pain, addiction and reward pathways. The MOR is a GPCR, a family of complex seven-helix protein bundle which are ubiquitously involved in pretty much every important physiological process, from vision to neuronal regulation.

Almost 30% of all drugs work by modulating the activity of GPCRs, but the problem until now has been the difficulty of obtaining atomic-level crystal structures of these proteins that would allow researchers to try out a rational drug design approach against them. All that changed with a series of breakthroughs that has allowed us to crystallize dozens of GPCRs (that particular tour de force bagged its inventors - Brian Kobilka and Robert Lefkowitz - the 2012 Nobel Prize, and I had the pleasure of interviewing Kobilka a few years ago). More recently they have crystallized the various opioid receptors, and this definitely represents a promising opportunity to find drugs against these proteins using rational drug design. When activated or inhibited by small molecules, GPCRs bind to several different proteins and engage several different biochemical pathways. Among other insights, detailed studies of GPCRs have allowed us to determine that the side effects of opioid drugs emerge mainly from engaging one particular protein called ß-arrestin.

A paper published this week in Nature documents a promising effort to this end. Using computational docking, Brian Shoichet (UCSF), Bryan Roth (UNC) and their colleagues have found a molecule called PZM21 that seems to activate the MOR without also activating the molecular mechanisms which cause dangerous side effects like respiratory depression. The team started with the crystal structure of the MOR and docked about 3 million commercial lead-like molecules against it. The top 2500 entries were visually inspected, and using known data on GPCR-ligand interactions (especially knowledge of the interaction between the ligands and an aspartate - Asp 147) and sound physicochemical principles (discarding strained ligands), they narrowed down the set to about 25 ligands which gave binding affinities ranging from 2 µM - 14 µM; interestingly, along with the Asp 147 interaction, the model located another unique hydrogen bond with Asp 147 that has never been seen before in GPCR-ligand binding (compound 7 in figure above). Using traditional structure-based design that allowed the group to access novel interactions, the affinity was then improved all the way to single digit nanomolar, which is an excellent starting point for future pharmacokinetic studies. The resulting compound hit all the right endpoints in the right assays; it had greatly reduced recruitment of ß-arrestin and subsequently reduced respiratory arrest and constipation, it mitigated pain efficiently in mice and showed much lesser addiction potential. It also showed selectivity against other opioid receptor subtypes, which may make it a more targeted drug.

This paper is noteworthy for several aspects, and not just for the promising result it achieved. Foremost among them is the fact that a relatively straightforward docking approach was used to find a novel ligand with novel interactions for one of the most coveted and intractable protein targets in pharmacology. But to me, equally noteworthy is the fact that the group examined 2500 docked ligands to pick the right ones and discard bad ones based on chemical intuition and existing knowledge; the final 25 ligands were not just selected from the top 50 or 100. That exercise shows that no modeling can compensate for the virtues of what chemist Garland Marshall called "a functioning pair of eyes connected to a good brain", and as I have found out myself, there's really no substitute for patiently slogging through hundreds of ligands and using your chemical intuition to pick the right ones.

There is little doubt that approaches such as this one will continue to find value in discovering new ligand chemistry and biology in valuable areas of pharmacology and medicine. As is becoming clear in the age of computing, the best results emerge not from computers or human beings alone, but when the two collaborate. This here is as good an example of that collaboration as I have recently seen.

Hit picking parties, eviscerating weaknesses and other vistas from the world of molecular docking

Here’s a good review of the pitfalls and promises of molecular docking by John Irwin and Brian Shoichet (UCSF) which is worth your time, especially if you are a non-specialist who wants a summary of what’s happening in the field. As the review notes, there are many first principles-based reasons why docking should not work – poor calculation of ligand conformations, poor treatment of protein and ligand electrostatics and desolvation, non-existent consideration of protein movement, wishful treatment of water molecules, sloppy representation of x-ray structures...the list just goes on. When docking 10⁷ or so ligands involving 10¹³ total configurations, any one of these “maddening details” can doom your study.

And yet as the authors note both through general considerations and a few case studies, incremental but steady improvements in docking methodology have now made the technique respectable in most structure-based drug design campaigns (which interestingly have provided more drug candidates than HTS). Several reasons have contributed to this respectability. The first is the sheer throughput; no experimental technique can possibly screen 10 million ligands in a few days or weeks, so even with its flaws docking rises up at least as a potential complement to experiment. As the review notes, even a 10% success rate in finding new binders to a good target would be an improvement, and with well-defined binding sites the success rate can surpass high throughput screening (thus, the correct question to ask is not whether your method gives you false positives and negatives but whether this error rate is enough to overwhelm the experimental discovery rate for true positives). The second factor is the existence of massive databases like ZINC and ChEMBL containing millions of annotated ligands which could serve as starting points for ligand discovery.

Thirdly, while the holy grail of docking would be to correctly predict the absolute affinities of your ligands or at least to rank them, the more modest goal is to try to separate binders from non-binders and to discover novel chemotypes. Docking has been reasonably successful in meeting this goal, and the review presents several case studies that discovered interesting ligands which were dissimilar to known ones. In turn however, what really makes the discovery of novel chemotypes interesting is that it could lead to novel biology. It is this ability to potentially “break out of medicinal chemistry boxes” that makes docking attractive. For instance you could potentially find agonists by docking when you only found antagonists before, or – in what is one of the more interesting examples illustrated in the article – you could ‘deorphanize’ an enzyme by docking potential substrates to it and predicting its reaction profile. I still find this evidence anecdotal, but it's at least a good starting point for trying out things.

The rest of the review is also useful, not in the least because - given that it's from the Shoichet lab - there's also an instructive checklist of caveats to keep in mind while experimentally screening (PAINS, aggregators etc.) ligands. There is also a discussion of using homology models for docking. Using homology models is tricky even for lead optimization, so I would be wary in applying them too widely for high-throughput docking. While the review does illustrate an interesting case involving a GPCR, I want to note that I did blog about this case when it was published and described how – new ligands notwithstanding – the calculation seemed to use enough computing power and models to light up a startup, along with copious expert input.

It’s this last point in the review which is really the crux of the matter. When the servers have cooled down and the electrons have stopped flowing, the most important equipment that one can bring to bear on a docking study is a good pair of eyes connected to an experienced brain. Even with small error rates "the scum can rise to the top" ("The scum is out there" could be a good tagline for a X-Files episode about HTS). There is little substitute for careful inspection of top docked hits and looking at things like strain, abnormal charged interactions and wrong tautomeric states; otherwise staring down a computer screen would present the same risks as staring down a gun barrel.

The Shoichet lab has occasional "hit-picking parties" where teams of medicinal and computational chemists examine docked structures, and I suspect these parties are more common in other places than you think (although probably not as common as they should be). It’s only when the high throughput-low accuracy domain of docking meets the low throughput-high accuracy domain of the human mind that docking will continue to be successful. Given what we have seen so far I think there are grounds for hope.

Drugable.com "ranks billions of drug interactions"? Hold your horses.

Now here's a study that should make most seasoned molecular modelers cringe. Nature News reports on an effort by website Drugable.com that docked 600,000 compounds to 7,000 protein targets and predicted which ones would show activity against these targets based on docking scores:

Predicting how untested compounds will interact with proteins in the body, as Drugable attempts to do, is more challenging. In setting up the website, Cardozo’s group selected about 600,000 molecules from PubChem and the European Bioinformatics Institute’s ChEMBL, which together catalogue millions of publicly available compounds. The group evaluated how strongly these molecules would bind to 7,000 structural ‘pockets’ on human proteins also described in the databases. Computing giant Google awarded the researchers the equivalent of more than 100 million hours of processor time on its supercomputers for the mammoth effort.

But mammoth computing resources do not translate to carefully constructed protocols or correct predictions. In its current incarnation, docking is best for finding the binding pose, that is, the orientation of a drug bound into a protein pocket. Ranking compounds is far more difficult, and predicting absolute binding affinities is a very distant, currently unachievable third goal.

Anyone who has tried to run a hit to lead or lead optimization project based on docking scores would know how riddled with problems and qualifications any prediction based on these highly subjective numbers is. For starters, every modeling program gives you its own docking scores. Absolute values of these numbers (which ideally should reflect the free energy of binding but which seldom do) are almost always useless. If you are dealing with a congeneric series of molecules and are fairly confident about the binding orientation (usually confirmed by x-ray crystallography or some other technique) then maybe you could get some help from the scores in ranking the compounds, but even then mostly in terms of trends rather than quantitative differences.

Unfortunately the news piece says nothing about what method was used to generate the poses, whether there was any clustering or whether only the top pose was considered, what the false positive rate was, and most importantly, whether there was any experimental verification whatever of the ranking. The website is also not helpful in this regard. It also does not tell us if the protein structures used for docking were well-resolved or refined or whether they were homology models. In the absence of all this information the ranking of the compounds is tenuous at best and useless at worst and as it stands the study sounds little better than throwing darts in the dark and hoping some of them will stick. Ranking often fails even for similar compounds, so how well (or badly) it would work for 600,000 diverse compounds bound to 7,000 diverse protein targets is anyone's guess.

The report also compares the study to a similar activity prediction study by Brian Shoichet in which drug similarity was used to predict activity against unexpected targets. But that was a very different kettle of fish; it was a compound similarity - not docking - study so it did not have to deal with the complexities of error-ridden protein crystal structures or homology models, it verified a lot of the predictions using carefully constructed assays, and even then it gave a hit rate which did not exceed about 50%.

Either the Drugable.com study itself has failed to validate its predictions or the news report is woefully incomplete. Maybe I am wrong and in fact the study has laid the careful groundwork and validation that is necessary for trusting docking. As it stands however, the purpose of the report mainly seems to be to highlight the fact that Google generously donated 100 million hours of its computing power to the docking. This heightened, throw-technology-at-it sense of wonder and optimism is exactly what the field does not need. I would be the first one to welcome reliable predictions of drug-protein affinity based on wholesale docking of compounds to targets, but I don't think this work achieves that goal at all.

Protein-ligand crystal structures: WYSI(N)WYG

Crystal structures of proteins bound to small molecules have become ubiquitous in drug discovery. These structures are routinely used for docking, homology modeling and lead optimization. Yet as several authors have shown over the years, many of these structures can hide flaws that don't become apparent until they are actually revealed through analysis. Worse still, there may be flaws that never become apparent because nobody takes the trouble to look at the original data.

A recent paper from the group at OpenEye has a pretty useful analysis of flaws in crystal structures. They carry on a tradition most prominently exemplified by Gerard Kleywegt at Uppsala. The authors describe common metrics used for picking crystal structures from the PDB and demonstrate their limitations. They propose new metrics to remedy the problems with these parameters. And they use these metrics to analyze about 700 crystal structures used in structure-based drug design and software validation and find out that only about 120 structures or so pass the rigorous criteria used for distinguishing good structures from bad ones.

The most important general message in the article is about the difference between accuracy and completeness of data, and the caveat that any structure on a computer screen is a model and not reality. Even very accurate looking data may be incomplete and this flaw is often neglected by modelers and medicinal chemists when picking crystal structures. For instance, the resolution of a protein structure is often used as a criterion for selecting one among many structures of the same protein from the PDB. Yet the resolution only tells you how far apart atoms can be distinguished from each other. It does not tell you if the data is complete to begin with. So for instance, a crystallographer can acquire only 80% of the theoretically maximum possible data and present it with a resolution of 1.5 A, in which case the structure is clearly incomplete and possibly flawed for use in structure-based drug design. Another important metric is the R-free factor which is obtained by omitting certain parts of the data and refitting the rest to the model. A difference between R-free and the R-factor (a factor denoting the original difference between the full set of data and the model) of more than 0.45 for a structure with resolution 3.5 A or more is a red flag.

The OpenEye authors instead talk about a variety of measures that provide much better information about the fidelity of the data than resolution. The true measure of the data is of course the actual electron density. Any structure that is seen on the computer screen results from the fitting of a model to this electron density. While proteins are often fit fairly well to the density, the placement of ligands is often more ambiguous, partly because protein crystallographers are not always interested in the small molecules. The paper documents several examples of electron density that was either not fit or incorrectly fit by ligand atoms. In some cases sparse or even non-existent density was fit by guesswork, as in the illustration above. All these badly fit atoms show up in the final structure but only an expert would know this. The only way to overcome these problems is to take a look at the original electron density, which is unfortunately a task that most medicinal chemists and modelers are ill equipped for.

The worst problem with these crystal structures is that because they look so accurate and nice on a computer screen, they may fall into Donald Rumsfeld's category of "unknown unknowns". But even with "known unknowns" the picture is a disturbing one. When the authors apply all their filters and metrics to a set of 728 protein-ligand structures frequently used in benchmarking docking programs and in actual structure-based drug design projects, they find that only 17% or so (121) of the structures make it through. They collect these structures together and call the resulting database Iridium which can be downloaded from their website. But the failure of the 580 or so common structures used by leading software vendors and pharmaceutical companies to pass important filters leaves you wondering how many resources we may have wasted by using them and how much uncertainty is out there. 

Something to think about, especially when considering the 50,000 or so protein structures in the PDB. What you see may not be what you get.

Constancy of the discodermolide hairpin motif

Our paper on the conformational analysis of discodermolide is now up on the ACS website. The following is a brief description of the work.

Discodermolide (DDM) is a well-known highly flexible polyketide that is the most potent microtubule polymerization agent known. In this capacity it functions very similar to taxol and the epothilones. However the binding mode of DDM will intimately depend on its conformations in solution.

To this end we have performed multiple force field conformational searches on DDM and the first surprising thing we noticed was that all four force fields located the same global minimum for the molecule in terms of geometry. This is surprising because, given the dissimilar parameterization criteria used in different force fields, minima obtained for flexible organic molecules are usually different for different force fields. Not only that, but all the minima closely superimposed on the x-ray structure of DDM which we call the "hairpin" motif. This is also surprising since the solid state structure of such a highly flexible molecule should not generally bear resemblance to a theoretically calculated global minimum.

Next, we used our NAMFIS methodology that combines parameters from conformational searches to coupling constants and interproton distances obtained from NMR data to determine DDM conformations in two solvents, water and DMSO. We were again surprised to see the x-ray/force field global minimum structure existing as a major component of the complex solution conformational ensemble. In many earlier studies, the x-ray structure has been located as a minor component so this too was unexpected.

However, this same structure has also been remarkably implicated as the bioactive conformation bound to tubulin by a series of elegant NMR experiments. To our knowledge, this is the first tubulin binder which has a single dominant preferred conformation in the solid-state, as a theoretical global minimum in multiple force field conformational searches, in solution as well as in the binding pocket of tubulin. In fact I personally don't know of any other molecule of this flexibility which exists as one dominant conformation in such extremely diverse environments; if this happened to every or even most molecules, drug discovery would suddenly become easier by an order of magnitude since all we would have to do to predict the binding mode of a drug would be to crystallize it or to look at its theoretical energy minima. To rationalize this very pronounced conformational preference of DDM, we analyze the energetics of three distributed synthons (methyl-hydroxy-methyl triads) in the molecule using molecular mechanics and quantum chemical methods; it seems that these three synthons modulate the conformational preferences of the molecule and essentially override other interactions with solvent, adjacent crystal entities, and amino acid elements in the protein.

Finally, we supplement this conformational analysis with a set of docking experiments which lead to a binding mode that is different from the earlier one postulated by NMR (as of now there is no x-ray structure of DDM bound to tubulin). We rationalize this binding mode in the light of SAR data for the molecule and describe why we prefer it to the previous one.

In summary then, DDM emerges as a unique molecule which seems to exist in one dominant conformation in highly dissimilar environments. The study also indicates the use of reinforcing synthons as modular elements to control conformation.

Jogalekar, A., Kriel, F., Shi, Q., Cornett, B., Cicero, D., & Snyder, J. (2009). The Discodermolide Hairpin Structure Flows from Conformationally Stable Modular Motifs Journal of Medicinal Chemistry DOI: 10.1021/jm9015284

New ligands for everyone's favorite protein

A landmark event in structural biology and pharmacology occurred in 2007 when the structure of the ß2-adrenergic receptor was solved using xray crystallography by Brian Kobilka's and Raymond Stevens's groups at Stanford and Scripps respectively. The structure was co-crystallized with the inverse agonist carazolol. Until then the only GPCR structure available was that of rhodopsin and all homology models of GPCR were based on this structure. The availability of this new high resolution structure opened new avenues for structure-based GPCR ligand discovery.

The ß2 binding pocket is especially suited for drug design since it is tight, narrow and lined with mostly hydrophobic residues with polar residues well-separated. Two crucial residues, an Asp and a Ser bind to the ubiquitous charged amino nitrogen present in most catecholamines and the aromatic section of the molecule docks deep into the hydrophobic pocket. These particular features also make computational docking more facile; a mix of polar and non-polar features with bridging waters can make docking and scoring more challenging.

Since the ß2 structure has been published, attempts are being made to use it as a template to build homology models of other GPCRs. A couple of months back I described an interesting proof-of-principle paper by Stefano Costanzi that sought to investigate how well a homology model based on the ß2 would perform. In that study carazolol itself was used as a ligand for docking into the homology model. Comparison with the original crystal structure revealed that while the ligand docked more or less satisfactorily, an important deviation in its orientation could be explained by a counterintuitive orientation of a Phe residue in the binding site. The study indicated that the devil is in the details when one is considering homology models.

However, finding ligands for the ß2 itself is also an important and interesting endeavor. Virtual screening could help in such studies. To this end Brian Shoichet, Brian Kobilka and their group have used the DOCK program to virtually screen one million lead-like ligands from their ZINC database against the ß2. Out of the 1 million ranked poses, they chose and clustered the top 500 compounds (0.05% of the database) into 25 unique chemotypes, a choice also guided by visual inspection of the protein-ligand interactions and commercial availability. They then tested these 25 compounds against the ß2 and found 6 compounds with IC50s better than 4 µM. One of these compounds with an IC50 of 9 nM is perhaps the most potent inverse agonist of the ß2 known. The binding poses revealed substantial overlap of similar functional groups with the carazolol structure. Two compounds turned out to have novel chemotypes and bore very little similarity with known ß2 ligands. A negative test was also run where a known predicted binder was chemical modified so that it would not bind.

Interestingly all the compounds found were inverse agonists. The ZINC library is somewhat biased against aminergic ligands as is most of chemical space. The catecholamine scaffold is one of the favourite scaffolds in medicinal chemistry. However, subtle difference in protein structure can sometimes turn an inverse agonist into an agonist. In this case, small changes in the orientation of the crucial Ser residue near the mouth of the binding pocket. In a past study for instance, slightly changing the rotameric features of the Ser residue thus resulting in a different orientation of the hydroxyl was sufficient to retrieve agonists.

The study thus shows the value of virtual screening in the discovery of new ß2 ligands and indicates the effect of library bias and protein structure on such ligand discovery. Many factors can contribute to the success or failure of such a search; nature is a multi-armed demon.

Reference:
Kolb, P., Rosenbaum, D., Irwin, J., Fung, J., Kobilka, B., & Shoichet, B. (2009). Structure-based discovery of ß2-adrenergic receptor ligands Proceedings of the National Academy of Sciences DOI: 10.1073/pnas.0812657106

Post-docking as a post-doc, and some fragment docking

I am now ready to post-doc. I am also now ready to post-dock, that is, engage in activities beyond docking. Sorry, I could not resist cracking that terrible joke. It's been a long journey and I have enjoyed every most moments of it. Thanks to everyone in the chemistry blogworld who regaled, informed, provoked and entertained on this blog. I am now ready to move on to the freakingly chilly Northeast. Location not disclosed for now, but maybe later.



Speaking of docking, here is a nice paper from the Shoichet group in which they use fragment docking to divine hits from a large library for a beta-lactamase. Fragment docking can often be tricky compared to "normal" docking since fragments being small usually demonstrate promiscuity, low-affinity and non-selectivity in binding. Fragment docking thus is not yet a completely validated technique.

In their study, the present authors screen their ZINC library for fragments binding to the ß lactamase CTX-M by docking using the program DOCK. They also screen a lead-like library for larger molecules. The top hits from the fragment docking results were assayed and showed micromolar inhibition against the lactamase. These included tetrazole scaffolds not seen before. Importantly, five of these hits could be crystallized and the high-res crystal structures validated the docking modes.

What was interesting was that the same tetrazole scaffolds in the larger lead-like library were ranked very low (>900) and would not have ever been selected had their tetrazole fragments not showed up at the top in the fragment docking results. These compounds, when assayed showed sub-milimolar to micromolar activity against the lactamase. Thus, the protocol essentially demonstrated that fragment docking can reveal hits that can be missed by docking larger lead-like molecules. One of the reasons DOCK succeeds in this capacity is because of its use of a physics-based scoring function that has no bias against fragments. It also helps that the active site of CTM-X is relatively rigid with little protein motion.

The fragments were also assayed against another lactamase for Amp C. Usually, hits for CTM-X and Amp C are mutually exclusive. What was seen was that the higher the potency of the fragments for CTX-M, the higher the specificity for CTX-M, not surprising considering that increased potency translates to a much better complementary fit of the fragments for CTX-M.

Fragment docking can be messy since fragments can bind non-selectively and haphazardly to many different parts of many different proteins. But this study indicates that fragment docking is not an uninteresting strategy to possibly find hits from other lead-like libraries that may be otherwise concealed.

The potencies of the compounds found may look pretty weak, but because there are extremely few molecules inhibiting these medicinally important lactamases, such advances are welcome. Lactamases are of course an important target for overcoming resistance in antibiotic treatment.

Reference:
Chen, Y., & Shoichet, B. (2009). Molecular docking and ligand specificity in fragment-based inhibitor discovery Nature Chemical Biology DOI: 10.1038/nchembio.155

Water-Inclusive Docking with Remarkable Approximations

The role of water in mediating protein-ligand interactions has now been well-recognized by both experimentalists and modelers. However it's been relatively recently that modelers have actually started taking the unique roles that water plays into account. While the role of water in bridging ligand and protein atoms is obvious, a more subtle but crucial role of water is to fill up hydrophobic pockets in proteins. Such waters can be very unhappy in such pockets because of both unfavourable entropy (not much movement) and enthalpy (inability to form a full complement of 4 hydrogen bonds). If one can design a ligand that will displace such waters, significant gains in affinity would be obtained. One docking approach that does take such properties of waters into consideration is Schrodinger's Glide, with a recent paper attesting to the importance of such a method for Factor Xa inhibitors.

Clearly the exclusion of water molecules during docking and virtual screening (VS) will hamper enrichment factors, namely how well you can rank actives above inactives. Now a series of experiments from Brian Shoichet's group illustrates the benefits of including waters in active sites when doing virtual screening. These experiments seem to work in spite of two approximations that should have posed significant problems, but surprisingly did not.

To initiate the experiments, the authors chose a set of 24 targets and their corresponding ligands from their well-known DUD ligand set. This is a VS data set in which ligands are distinguished by topology but not by physical properties such as size and lipophilicity. This feature makes sure that ligands aren't trivially distinguished by VS methods on the basis of such properties alone. Importantly, the complexes were chosen so that the waters in them are bridging waters with at least two hydrogen bonds to the protein, and not waters which simply occupy hydrophobic pockets. Note that this would exclude a lot of important cases where affinity comes from displacement of such waters.

Now for the approximations. Firstly, the authors treated each water molecule separately in multiple configurations. They then scored the docked ligands against each such configuration as well as the rest of the protein. The waters were treated as either "on" or "off", that is, either displaced or not displaced. Whether to keep a water or not depended on whether the score improved or not when it was displaced by a ligand. The best scored ligands were then selected and figured high on the enrichment curve. This is a significant approximation because the assumption here is that every water contributes to ligand binding affinity independently of the other waters. While this would be true in certain cases, there is no reason to assume that it would generally hold.

The second approximation was even more important and startling. All the waters were regarded as energetically equivalent. From our knowledge of protein-ligand interactions, we know that the reason why evaluating waters in protein active sites is such a tricky business is precisely because each water has a different energetic profile. In fact the Factor Xa study cited above takes this profile into consideration. Without such an analysis it would be difficult to tell the medicinal chemist which part of the molecule to modify to get the best binding affinity from water displacement.

The most important benefit of this approximate approach was a linear increase in computational time instead of an exponential one. This was clearly because of the separate-water configuration approximation. The calculation of individual water free energies would also have added to this time.

In spite of these crucial approximations, the results indicate that the ability to distinguish actives from inactives was considerably improved for 12 out of 24 targets. This is not saying much, but even 50% sounds like a lot in the face of such approximations. Clearly an examination of the protein active site will also help to evaluate which cases will benefit, but it will also naturally depend on the structure of the ligand.

For now, this is an encouraging result and indicates that this approach could be implemented in virtual screening. There are probably very few cases where docking accuracy decreases when waters are included. With the sparse increases in computational time, this would be a quick and dirty but viable approach for virtual screening.

Reference:
Niu Huang, Brian K. Shoichet (2008). Exploiting Ordered Waters in Molecular Docking Journal of Medicinal Chemistry, 51 (16), 4862-4865 DOI: 10.1021/jm8006239

Multiconformational MMGBSA Rescoring; Advancing On Mount Free Energy

Blogging has been a little slow lately mainly because there have been exciting new developments with one of the projects I have been involved in and I was in meetings related to this. One of the topics that was discussed at the conference I was at last week was the accurate prediction of free energies of binding, one of the holy grails of drug discovery. Free-energy perturbation (FEP) still remains the gold standard to get relative free energies of binding, but the procedure is very computer intensive and therefore can be carried out only with small changes in congeneric series of inhibitors. The goal remains elusive and extremely challenging.

A poor man's way of quickly obtaining such ∆Gs is MMGBSA (Molecular Mechanics Generalized Born Surface Area). The GBSA model is well-established as a continuum solvation model for taking solvation into account. What MMGBSA does is take a docked ligand structure and then calculate the free energy of binding as the difference between the bound and unbound states using a force field, including implicit solvation.

Therefore, it calculates

∆G (binding) = ∆G (protein-ligand complex) - ∆G (protein) - ∆G (ligand)

Clearly it has to calculate the energies of the free ligand and free protein. Much of the challenge lies in these two terms. For starters, one has to calculate the strain energy penalty that the protein has to pay in order to bind the ligand. The binding energy that we see experimentally emerges after the protein has paid this strain penalty. How much this strain energy can be has been a controversial topic recently and I will get into it in another post. Suffice it to say that it's a challenging calculation that is not always handled well by MMGBSA. This is because in calculating the ligand free energy, MMGBSA essentially uses a force field to relax the ligand from the bound conformation to the nearest local energy minimum. However, a complex ligand exists in several local energy minima in solution and this force field local minimum may not correspond to any of them. Thus, one has to consider the global strain penalty that the protein has to pay. For this the method also has to consider the multiple conformations that a ligand adopts in solution. Sadly there are very few techniques that will deconvolute the Boltzmann population of a ligand's real conformations in solution and give us the global minimum. This problem in calculating strain energies remains an important drawback of the method.

Calculating ∆G (protein) is also not a trivial matter. We need to consider the entropy of the protein. One can get this from time-consuming MD simulations but it's not certain if the force field is parametrized well and if conformational space has been sampled comprehensively. Another uncertain factor is the induced fit effects involved in binding. A lot of these effects can be subtle and may extend to second shell amino acid residues.

Given these drawbacks, MMGBSA has nonetheless been quite successful in improving agreement with experiment. One of the reasons it works so well is that when you are dealing with congeneric series of ligands for a given target, many of the terms like conformational entropy and protein reorganization energy are the same or very similar and cancel, although there can be surprises. It seems now that at least one of the problems in MMGBSA- not considering the multiple conformations of the ligand in solution- can be tackled. A simple way to get multiple conformations of a ligand in solution is to do a conformational search. Assuming that the search is "complete", one can then calculate the conformational entropy penalty that the ligand has to pay in order to sacrifice all conformations except one in which it binds to the protein. There has been an implicit way to take this into account- many docking programs include a penalty of 0.65 kcal/mol per frozen rotatable bond. But clearly this penalty may be quite less if there are hundreds of conformations in solution that would lead to a large conformational penalty.

Now a group from Amgen has done such multi-conformational MMGBSA rescoring for four important targets and their ligands- CDK2, Thrombin, Factor Xa and HIV-RT. They compare scores obtained with Schrodinger's GlideXP routine with experimental binding affinities. Then they compare scores obtained with MMGBSA rescoring either with a single ligand conformer representation or with a multiple conformer representation that takes ligand conformational entropy into account. The comparison between single and multiple conformers gives somewhat mixed results and sometimes the single conformer representation also does fairly well; however, one thing is strikingly clear, that MMGBSA rescoring can radically improve correlation with experimental affinities compared to simple GlideXP scoring. In some cases the correlation coefficient jumps from essentially 0.00 to a whopping (by current standards) 0.75. There is a lot of interesting methodology described in the paper worth taking a look at. But it's quite clear how including some of the explicit physical effects involved in protein-ligand binding can substantially improve correlation with experiment. In this case the extra effort expended is a fraction of the cost involved in FEP calculations and the methods can also tackle more diverse ligands.

Even if we are not close to conquering the free energy fort, at least we seem to be getting concrete footholds on it.

Guimaraes, C.R., Cardozo, M. (2008). MM-GB/SA Rescoring of Docking Poses in Structure-Based Lead Optimization. Journal of Chemical Information and Modeling, 48(5), 958-970. DOI: 10.1021/ci800004w

Computational modeling of GPCRs: not too bad

GPCRs constitute one of the most important family of proteins in our body, both for their innate importance in signal transduction and neurotransmission, and as important targets for drugs. Many of the important drugs on the market today target GPCRs. And yet there is an unusual gap between knowledge and application when it comes to this important family. That's because only two crystal structures of GPCRs are known. And one of them was derived last year, so there's been a real dearth of structural information about GPCRs for a long time.

We do know something about many GPCRs, however. We know that they are 7-TM receptor-spanning proteins. And the two structures we do know about shed valuable insight on GPCR function. One is rhodopsin which has been around for a while. Then there was big news last year about the second important GPCR whose structure was determined- the ß-2 adrenergic receptor.

Given the paucity of structural information and the availability of two structures, a logical question is whether computational modeling can teach us something new about GPCRs whose structure is unknown. To this end, Stefano Costanzi at the NIH did a nice set of experiments which he published in J. Med. Chem. He attempted to build a homology model of the adrenergic receptor based on the sequence and structure of rhodopsin. Since we now have a crystal structure of the adrenergic GPCR, we have something concrete to compare modeled structures and ligand orientations to.

Costanzi was particularly interested in knowing how a small molecule-carazolol- binds to the modeled GPCR. This is important both from a structural and functional drug-discovery point of view. His results indicate that we can do pretty well. In essence, he built two models of the receptor, one of them de novo. While the models were similar to rhodopsin in the conserved regions, the important differences were with respect to a loop that flaps on top of the protein. In one model the loop was buried inside the binding pocket, and in the other one it was open. Docking of carazolol into the buriled-loop model using the Glide program from Schrodinger gave a binding pose in which the ligand was, not surprisingly, buried deeper into the cavity compared to the crystal structure. This was naturally the effect of the loop blocking part of the pocket. The other model in which the loop was not buried gave much better results. Curiously, the ligand was buried a little deep in the pocket even in this model, even though it was much less buried compared to the previous one. It still misaligned considerably with the experimental pose. Inspection revealed that there was a Phe in the pocket which was anti in the model but +gauche in the crystal structure. Since the corresponding residue in rhodopsin was Ala, there was no way this unusual conformation could have been predicted ab initio. Fixing the conformation of this residue to +gauche suddenly gave excellent alignment with the ligand orientation in the crystal structure.

An instructive piece of work that shows that homology modeling and docking of ligands into GPCRs of unknown structure can be fruitful. However, it also indicates caveats like the Phe conformation which are hard to account for de novo. However, since structures of members in this important family of proteins are unavailable anyway, even some predictive ability might be welcome in this area.

Costanzi, S. (2008). On the Applicability of GPCR Homology Models to Computer-Aided Drug Discovery: A Comparison between In Silico and Crystal Structures of the ß2-Adrenergic Receptor. Journal of Medicinal Chemistry DOI: 10.1021/jm800044k

A relatively rare example of docking-based virtual screening

Many studies published in the last few years have demonstrated that in general, ligand-based methods are better for virtual screening compared to structure-based docking methods. For example, a 2007 Merck study showed that 2-D similarity searching methods are quite good for finding similar leads, while 3-D methods can do some scaffold hopping and find new families of structures. Both methods are generally superior to docking. One of the reasons for this is that docking is not really designed for virtual screening; docking is much more valuable for prediction of crystallographic conformations and most importantly, predicting binding affinity, which is the holy grail of the industry. The latter task is still extremely challenging, although dents have been made in tackling it.

In any case, so this group from Vertex tackled a kinase inhibitor search problem for Pim-1 kinase using docking, and this seems to be one of those cases where docking with Schrodinger's Glide program helped complement and indeed improve upon HTS. The group screened a large database enriched in kinase inhibtors by HTS and got only a 0.3% hit rate. They decided to find out if VS could do better. They used Glide to screen a corporate collection that was less enriched in kinase inhibitors, to avoid bias. They used Glide not in the VS mode but the regular docking mode which takes more time but is more accurate. They used some astute filters to avoid getting false hits from large molecules that fit better in the site. They also used an C-H aromatic hydrogen bond constrain in the docking.

After screening out compounds that were too large and hydrophobic, they got 4 compounds (a 4% hit rate) with activities ranging from 90 nM to 550nM. Two of these could be crystallised and it was confirmed that the experimental conformation was very close to the predicted binding conformation. Glide also picked up the "weak" C-H aromatic hydrogen bond. The authors conjecture that the reason why Glide chose this H-bond is because the traditional hinge region of Pim-1 kinase is more hydrophobic than that in other kinases because of a proline residue. The study demonstrates how VS can serve as a valuable complement to HTS.

Pierce, A.C., Jacobs, M., Stuver-Moody, C. (2008). Docking Study Yields Four Novel Inhibitors of the Protooncogene Pim-1 Kinase. Journal of Medicinal Chemistry DOI: 10.1021/jm701248t

The BI guys should have used Glide XP Dock

In the last post, I was wondering why the BI guys did not use any docking to conjecture the binding mode of their new p38 inhibs. While docking is not a foolproof predictor, it can shed light on possible anomalous modes. Particularly interesting was this observation about the benzothiophene sulfur reversing its trans preference to become cis to the adjacent amide oxygen. This counterintuitive observation was later explained by alluding to the glutamate that would have had an unfavourable electrostatic interaction with the sulfur had it been trans to the oxygen. The observation was revealed by the crystal structure. However, before the crystal structure was obtained, they went through the design and synthesis of two compounds based on the reasonable hypothesis that the sulfur would be trans. It was only after their puzzlement with the failure of these designed compounds that the situation became clear through crystallography.

Well, the observation is also revealed by Glide XP docking, and I think this could have saved them some time. I first duplicated the X-ray pose of the earlier BIRB compound in p38 (1KV2 PDB id), and then docked the new compound in. The result; the binding score was still good, if not as good as for BIRB. But the important fact was that all the top 5 docking poses from Glide indicated the sulfur to be cis to the oxygen and away from the Glu, just as in the crystal. In fact, even the fragment docked in the same position as the crystal structure indicated, with the sulfur cis. Not totally surprising if electrostatic interactions are well parametrized and recognised by the scoring function. A further advancement where polarization effects are taken into account during docking would help a lot (This is in the works)

This was one case in which docking could have said basically the same thing that the crystal structure did. In this case, docking could have saved them the design and synthesis of two extra compounds based on a misleading hypothesis, and perhaps additional head-scratching validated by crystallography. On a related note, Glide is pretty well-parametrized on a couple of kinases, including Lck, CDK2, and p38.

How do you choose a good crystal structure for docking?

The first step in much of SBDD, including docking, is the selection of a good crystal structure if it exists. The crystal structure is used as the starting point for seeking new leads and optimizing them. Consider any docking method evaluation paper in J. Med. Chem. and one will come across a benchmarking set of protein structures that are used as starting models for testing the docking protocols.

Now crystal structures are frequently as close as you can get to "reality", but even they are models and should be treated with some skepticism. But the more obvious question for such a study when multiple crystal structures of a protein are available is, which crystal structure among those should you use?

The short answer to this question is, choose one with good resolution (preferably 2.0 A or less), which does not have missing portions, and which is preferably also unencumbered by the presence of a whole lot of counterions, stabilizing molecules, and other ligands.

But is that really all? Maybe not. Recently, I was playing around with docking some molecules into kinase crystal structures. I was trying to see if docking scores can correlate with the selectivity for one related protein over the other. Usually they don't, but I was going to look at similar proteins and similar structures, so I though it may be worth a shot. I was particularly looking at cyclin-dependent kinases (CDKs) which share a lot of homology especially in their ATP binding pocket. CDK2 is probably the most well-characterised CDK among the CDKs, and there are at least four to five different high-resolution CDK2 structures in the PDB. Also, I was more keen on using CDK2, because it was one of the proteins used for benchmarking the docking program.

So I decided upon two structures, both of high resolution. One had ATP docked into it, the other one had Staurosporine. I took an inhibitor which was known to be selective for another CDK over CDK2. First I docked it into that other CDK, and into the CDK2 structure that had ATP bound to it (without the ATP of course). I noted that the score for the other CDK was higher (which actually means more negative, since it is supposed to reflect the free energy of binding). That was consistent with the experimental data, which showed that the inhibitor was in fact more selective for the other kinase. But then, I docked it into the other CDK2 structure, and now the score was much better than for the other kinase. So the two docking runs gave two opposite results for the same protein. One predicted that the inhibitor would be less selective for CDK2, and the other one predicted that it would be more selective.

Now one of the things this says is that you cannot trust docking scores much. But this still was weird, because the question persists; which CDK2 structure should I use if I am going to do some SBDD and selectivity studies? I don't know the answer to this question, but I took a look at the two structures to try to figure out. In the one with the ATP, the adenine region of ATP nicely made two hydrogen bonds with the hinge region of the kinase, and so did my inhibitor which was supposed to be an ATP mimic. In the other one however, the backbone carbonyl that was supposed to form the hydrogen bond to the inhibitor was rotated by almost 90 degrees upwards. It did not form a bond with stauroporine, and it did not have to, because staurosporine does not "look" like ATP. And needless to say, it could not hydrogen bond with my inhibitor too. That's why the docking score was much worse.

What's the solution for circumventing such a problem? One quick answer that comes to my mind is; if you are docking a ligand that is "similar" to ATP, use the protein structure that has ATP bound to it. However, "similarity" can be a tricky concept, and should be considered carefully. Also, it may be slightly easy for kinase inhibitors, because there are literally hundreds of very typical planar, heterocyclic amino-pyrimidine based kinase inhibitors that share some very obvious similarity to ATP (or not...)

But probably the best message to take home from this from a computational standpoint is that rigid protein docking not surprisingly can get you into some bad trouble. Not allowing the protein to move means that you are going to preconstrain the protein based on its preconstrained conformation in the crystal. To test this thought, I did an induced-fit docking run on both structures with the inhibitor. Gratifingly, both the runs converged on the same protein-ligand structure.

Choosing a PDB x-ray structure may not be as easy as we think, and may have to be done critically. And more importantly as usual, what we put in is what we get out. Rigid docking is ok if there's only one crystal structure, and then only because there's no other choice. But in other circumstances, always allow the protein to move. That's closer to nature.