Computer-Aided Drug Design 2003-2005

BY BERNARD COUPEZ, HENRIK MOBITZ AND RICHARD A. LEWIS

Novartis Institutes for Biomedical Research, Basel CH-4002, Switzerland

1 Introduction

The themes for this review again have been driven strongly by the need of the Pharmaceutical industry to make the discovery process quicker and more reliable. Virtual screening in all its forms is at the heart of most research, from bioavailability niters through to rigorous estimations of the free energy of binding. Two areas of relative heat have been docking/scoring, and ADME/ Tox. On the other hand, 3D-QSAR and pharmacophores have become quiet. Part of the reason for this may arise from the successes in high-throughput crystallography, delivering more targets and complexes, the relative failure of HTS, and the increase in the amount of high quality data coming from late-phase research/early-phase development concerning the fate of clinical candidates. These trends look set to continue in the future, and the next two years should yield many new breakthroughs.

2 ADME/Tox and Druggability

There has been a fresh impetus to the modelling of ADME, Toxicity and druggability phenomena, partly driven by a desire to understand why such complex phenomena can, apparently, be described so simply, and partly to see if better models can be built, to improve the attrition rate in medicinal chemistry still further.

2.1 Druggability and Bioavailability. – In the continuing debate over what physicochemical properties are required for bioavailability, Vieth et al. have surveyed 1729 marketed drugs with respect to their route of administration, h-bonding capability, lipophilicity and flexibility. One conclusion they draw is that these properties have not varied substantially over time, implying that oral bioavailability is independent of target or molecular complexity. Compounds with lower molecular weight, balanced lipophilicity and less flexibility tend to be favoured. Leeson and Davis claim that molecular weight, flexibility, the number of O and N atoms and hydrogen-bond acceptors have risen, by up to 29%. This may be partly due to the choice of 1983 as the reference year, or the advent of more complex targets with greater selectivity needs (e.g. kinases). In the same vein, a study re-examined the correlation of flexibility and polar surface area (PSA) with bioavailability proposed by Veber et al. One conclusion is that there are significant differences in the ways of denning flexibility and PSA, and the correlations depend markedly on the method used (this is not surprising, as neither quantity is precisely definable). A second conclusion was that the limits denned (Number of rotatable bond< 10, PSA< 140 Å2 excluded a significant number of compounds with acceptable rat bioavailability. In the authors' words, "This observation underscores the potential danger of attempting to generalise a very complicated endpoint and of using that generalisation in a prospective selection application". Despite this, another bioavail-ability score has been devised, to predict the probability that a compound has > 10% bioavailability in the rat. Compounds are grouped by ionisation class (anions, cations, neutral). It was found that the standard rule-of-5 does well for cations and neutrals (88% of the compounds predicted to have low bioavailability are observed as such). Anionic compounds were better described by PSA limits. Some simple rules are given to compute the bioavailability score. In Abbott laboratories, this score is now routinely computed for all compounds and is used for hit-list triaging. It will be interesting to see if the results can be repeated on other data sets; the paper has certainly sparked much interest in the modelling community. Wegner provides support for the idea that human intestinal absorption correlates with PSA, by generating a classification model. The justification is that the error in the experimental data is 25%, and 80% of the observations occur in the top and bottom quartiles, that is, the data is more binary than evenly spread. In addition to PSA, other descriptors that reflect the electronic character of atoms and their environment also came to the fore.

2.2 Metabolism, Inhibitors and Substrates. – The field of cytochrome modelling is becoming more mature as we begin to understand the limitations of the experimental data and the subtleties of the mechanisms (the whole field of cytochrome P450 modelling, including homology, pharmacophore and 3D-QSAR models has been reviewed in detail recently). Empirical models are still preferred, especially for rapid evaluation of large libraries. In one case, use of a jury system improved prediction accuracy to over 90%. Chohan et al. have developed 4 models for Cytochrome P450 (Cyp) 1A2 inhibition, and identified the expected descriptors as being important to the QSAR (lipophilicity, aromaticity, HOMO/LUMO energies). Perhaps a more interesting result in this paper was the use of the k index to assess predictive powers of the models using test data.

k = observed agreement-chance agreement/total observed-chance agreement

This index should prove useful for data sets that are diverse and noisy. The validity of QSAR model predictions has also been studied by Guha and Jurs. The protocol is quite straightforward. The initial QSAR models were built, and the residuals of the compounds in the training set were used to classify the trains set predictions into good and bad. The threshold for the classification is arbitrary. Test compounds were predicted, and the predictions were grouped by substructural similarity to the nearest neighbour in the training set. It was seen that test compounds that had neighbours with low/good residuals were themselves well-predicted, with the reverse being the case for neighbours with high residuals. The success rate for classifying the strength of the prediction was 73% to 94%. The Merck group performed a retrospective study of in-house data sets, and concluded that the distance to the nearest neighbour, and the number of nearest neighbours (local density) were the two most useful measures for predicting prediction quality. They also concluded that distance does not have to be measured in the same descriptor space as was used to build the QSAR model. Topological descriptors combined with a Dice coefficient worked equally well.

A number of groups have been active in the prediction of the most likely sites of metabolism of molecules that are substrates for cytochromes. Singh et al. developed a semi-quantitative method based on the energy barrier to the creation of hydrogen radicals as calculated by AMI. Using a set of 50 substrates for Cyp 3A4, they were able to show that only hydrogens with a solvent-accessible surface area over 8 Å2 are susceptible to attack. The expensive quantum mechanic calculations could be approximated by local neighbourhood descriptors which could be well correlated to the energies (R 2 = 0.98), offering a fast and practical method for screening large libraries. An extension of this concept is embodied in the MetaSite program, which uses propensity to react, accessibility and GRID molecular interaction fields as descriptors. The methodology is more general, and can be applied to any cytochrome structure: in validation experiments, an accuracy of 80% is claimed. It is also important to be able to predict which compounds will be inhibitors as well as substrates, to avoid drug-drug interactions. A classifier based on a support vector machine (SVM) has been created that correctly predicts compounds into high, medium and low affinity at 70% accuracy, even with simple 2D descriptors. The improved accuracy was obtained through a systematic variation and optimisation of the SVM parameters.

Considering the success of surprisingly simple, semiempirical methods in ADME modelling, it is interesting to see whether more advanced methods could bring further improvements. A recent paper of Beck provides a link to the rich literature of DFT studies of hemes and cytochromes. The author uses Fukui functions to gauge the site of highest nucleophilicity of a number of known drugs. The predictions give mixed results and demonstrate that the implicit assumption of Fukui functions, i.e. an isotropic electrophilic attack, is flawed, not to mention that their MO-like shape does not allow a ranking of single atoms. In conclusion, the study suggests that it is more important to have an accurate description of the cytochrome-ligand complex than to invest in a high-level description of the chemical reactivity. De Visser et al. have used DFT on 10 C–H barriers with reference to bacterial cytochromes, and claim an excellent correlation between bond energy and observed activation energy barriers, so there is still some mileage in this approach.

2.3 Toxicity. – Unacceptable toxicity is still a key source of compound failure in clinical trials. Several groups have developed tools and programs for predicting toxicity for use in early phase, but the question arises about the accuracy of these models, and the levels of false positives and negatives that are acceptable. In research, an overly strict model with no false negatives may cause the discarding of a perfectly reasonable lead series. In development, missing a toxic alert which shows up in a later phase is unacceptable. Similarly, any program that is used by regulatory authorities to screen compounds must be very unforgiving of any flaw. In a recent study by the FDA on maximum human therapeutic dose, rules-based programs managed 64% accuracy, not much over random, giving an indication of the pitfalls in this field; Helma has given an overview of this area. Clearly, the domain of the models is critical and this has been addressed explicitly for QSARs that make toxicity predictions. Another route to predicting ADME properties is to use screening results, as exemplified by the Bioprint approach. 1198 drugs have been assayed against 130 screens, to give an activity fingerprint. QSAR models are then derived using pharmacophore descriptors. New compounds can be run through the models to predict binding affinity in all the screens, compared to the nearest neighbours in the database and finally fingerprinted themselves for confirmation. Using the affinity fingerprint alone, one can again identify similar molecules (sometimes with surprising results) and extrapolate to the potential side-effect profiles. This is very useful when selecting one from several lead series for optimisation.

3 Docking and Scoring

3.1 Ligand Database Preparation. – The ligand database is the basis for virtual screening (VS). Special care must be taken at this stage; accurate and physically relevant tautomeric and protonation states need to be assigned. Often compounds are registered in a database as a tautomer that is not necessarily the most probable state of the molecule and it is difficult to assign the correct state, so all relevant states should be generated. Similarly, as the stereochemistry of chiral centres is often not known, one must generate all stereoisomers. A recent article reveals the impact of pre-processing a database containing both known actives and inactives, where multiple protonated, tautomeric, stereochemical, and conformational states have been enumerated. The authors show that the interplay between 2D representations, stereochemical information, protonation states, and ligand conformation ensembles has a profound effect on the success rates of VS and conclude that the enrichment is highly dependent on the initial treatments used in database construction. In a paper that is bound to become a citation classic for the service that it has provided to the academic modelling community, Irwin and Shoichet describe the creation of the ZINC database of commercially available compounds, available via the web. The resource can be used in virtual screening studies, as the authors have taken care to provide compounds in multiple protonation and tautomer states, even multiple conformations. The paper provides a useful recipe for creating such a database for general use.

3.2 Target Preparation. – Thanks to high throughput crystallography and structural genomics, we have the X-ray structures of many targets of therapeutic interest, with the obvious exception of membrane proteins. When no experimental structure is available, it is possible to generate a 3D structure based on a template protein of similar sequence and a known structure, for example the model of CDK10, based on the CDK2 crystal structure, that was successfully used for a docking study.

If several structures of the target are available, which structures should be used: the apo form, a holo complex, or a homology model? This issue was examined by McGovern et al. They docked a large number of small molecules against 10 targets using the apo-, holo-, and modelled forms of the binding site. Using enrichment rates, they found that the holo form gave the best results (70% enrichment) followed by the apo (20%) and then the modelled form (10%). However, the holo form can be over influenced by the ligand in holo complex, if the active site has "collapsed" around the ligand. Then one would get a lower retrieval rate of similar but larger ligands, due to the increased steric constraints; the apo form of an active site can be markedly different from holo form. The conclusion is that VS using any form of the target will do better than random, but the holo form will give a best enrichment. This was also confirmed by Erickson et al. which show that the docking accuracy decreases dramatically if one uses an average or apo structure. Another approach is to use softened repulsive terms in the Lennard-Jones potential, to allow a closer approach of ligand and protein atoms that could be later resolved by minimisation. The T4 lysozyme system was used, with the ACD database as the source of ligands. The soft function was worse than the hard function, if multiple protein conformations were used, and vice versa for a single model. It was concluded that soft potential favour the decoys as much as the true ligands, so needs to be used with care.

Like the ligand preparation, the preparation of the target also requires great care. Incorrect protonation states or tautomers of histidines can lead to serious docking errors. For example Polgar et al. demonstrate the importance of protonation states in virtual screening for β-secretase (BACE1) inhibitors. They observed improvement of enrichment rates when they assigned different protonation states to catalytic Asp32 and Asp228 residues. Some docking methods require the addition of hydrogens. It is recommended that after the addition of hydrogen atoms to the protein, the positions of the hydrogens are relaxed by energy minimization to avoid any steric clashes. The positioning of hydrogen atoms on hydroxyl groups in the active site should also be checked and changed if necessary. In some instances, hydrogen bonds to crystallographic waters might need to be maintained for the docking.

However increasing the degrees of flexibility also increases the computational complexity and cost. Different methods have been described in the literature to tackle this critical issue (for a review see ref. 29). Often these methods model the flexibility in the binding site exclusively, by sampling the protein conformational space using molecular dynamics or Monte Carlo calculations or rotamer libraries. Another way of treating protein flexibility is to use an ensemble of protein conformations, rather than a single one. In a recent paper, Barril and Morley use all the X-ray structures of cyclin dependant kinase 2 (CDK2) and heatshock protein 90 (HSP90) to assess the performance of flexible receptor docking. They observe that flexible receptor docking performs much better in binding-mode prediction than rigid receptor docking. However, they also noticed that for library screening, ensembles of cavities often result in worse hit rates than rigid docking. This trend can be reversed by selecting those ligands that bind consistently well to many cavities in the ensemble.

3.3 Water Molecules. – Another challenge in protein-ligand docking is the modelling of the water molecules in protein ligand recognition. Water can form hydrogen bonds between the protein and the ligand or can be displaced by the ligand. Recently a new approach that allows this was implemented by Verdonk et al. in GOLD. The method allows water molecules to switch on and off and to spin. The explicit inclusion of water molecules in a docking program improves the binding mode when a ligand interacts with a water molecule. A distinction can also be made between the compounds that can displace a water molecules and the compounds that cannot. They claim that their algorithm correctly predicts water mediation/displacement in 93% of their tests and they observe some slight improvements in binding mode quality for water-mediated complexes. Similar results were reported by De Graf et al. for cytochrome P450s. The waters were either removed, or the crystallographic waters retained, or waters in GRID minima were used. Surprisingly, the last scenario gave the best results by up to 20% in the number of correct poses that scored highest.