Affordable virtual screening with Blaze: Benchmarks


We released BlazeGPU a couple of years ago, allowing the full power of the Blaze virtual screening system to be used on a few consumer graphics cards rather than a full-scale Linux cluster. Since then, graphics cards and CPUs have only got faster, so we decided that it was time to update our benchmarks and see how well all of the new hardware performs.

For these benchmarks we took a random subset of 4,000 molecules from our in-house Blaze data set and searched with a medium-sized query molecule. The molecules in the data set average 80 conformers each. We’ve run with three different search conditions: the full slow-but-accurate simplex algorithm, the standard clique algorithm and the new fastclique algorithm. All of these were run with 50% fields and 50% shape.

CPU performance

Firstly, the CPU benchmarks. All of these are single-core performance, but with all cores loaded so that we’re not benefitting from Intel Turbo Boost. In most cases Blaze will be saturating all cores, so this is representative of real-world performance. Note that the vertical axis is on a log scale.

CPU benchmarks

As can be seen, there’s a significant performance difference between the older CPUs at Cresset (such as the Q6600) and the newer Ivy Bridge i7-3770K chips, but not nearly as much as you would expect given that the Q6600s are around 7-8 years old at this point. The significant speed improvements of the fastclique algorithm are clearly visible with the throughput being more than 4x greater than the original clique algorithm. The last set of columns on the graph are from an Amazon c4.xlarge instance and show that the performance of each core on those systems is roughly the same as the Sandy Bridge i3-2120.

GPU Performance

Moving on to the GPUs, we’ve tested the throughput on a variety of different systems. Firstly, we’ve tested a variety of GTX580s on different motherboards and processors. As you would expect, for the most part the performance is governed by the GPU, but the exception is the fifth test system which is noticeably slower than the others. That card is sitting in a much older chassis with an older motherboard and hence is probably suffering from lack of backplane bandwidth to the GPU.

GPU benchmarks

The newer GTX960s perform extremely well on the Blaze calculations. We weren’t sure if they would, after the disappointment of the GTX680 which was noticeably slower than the 580 (data not shown). The difference is noticeable in the clique stages, but really stands out in the simplex calculations where a GTX960 is 50% faster than the GTX580s. By contrast, the high-end Tesla hardware is not a great performer on the Blaze OpenCL kernels. By all accounts the Tesla hardware is significantly faster than the consumer hardware on double precision workloads, but the Blaze code is all single precision and in that realm the cheap consumer hardware has an unbeatable price/performance advantage.
Finally, the GRID K520 is the hardware found on the Amazon g2.2xlarge and g2.8xlarge instances. As can be seen, it’s not a brilliant performer on the Blaze workload, being around the same speed as the Tesla on the fastclique algorithm but noticeably slower than all of the other cards tested on the simplex workload. However, it provides a nice test of GPU scaling: when running on a 4 times larger data set on all 4 GPUs of a g2.8xlarge instance, we observed substantially the same throughput as running the original data set on a single K520 GPU, showing that we can parallelise across multi-GPU systems with no loss of performance.

Cost efficiency on Amazon

Converting the throughput shown above, we can look at the cost of screening on the Amazon cluster with Blaze. The raw cost to screen a million molecules is shown in the table. Note that the actual costs will be somewhat higher, due to job overheads and data transfer costs.

Cost efficiency on Amazon

The Amazon GPU solutions are noticeably cheaper for fastclique jobs, roughly cost-competitive for the clique runs, but the poor performance of the K520 on the simplex task means that it is significantly more expensive there. As a result, at the moment there’s no real impetus to use the Amazon GPU resources unless you can get them significantly more discounted than the CPU instances on the spot market.


New hardware is significantly faster at running Blaze than old stock as would be expected. However, the speed increases are much lower than they have been in the past, with CPUs that are well past their best still performing adequately. On the GPU side, Blaze performs particularly well on commodity graphics cards leaving few reasons for us to invest in dedicated GPU co-processing cards.

The cost of running a million molecule virtual screen on the Amazon cloud has never been cheaper. If tiered processing is used as is the default for Blaze then these screens can be performed for a very low cost indeed – less than $15 per million molecules for the processing costs.

Contact us for a free evaluation to try Blaze on your own cluster, or Blaze Cloud.

Bioisosterism – harder than you might think

A common concept in medicinal chemistry is the bioisostere. The exact definition of a bioisostere is rather fuzzy, but the Wikipedia entry is as good as any: ‘bioisosteres are chemical substituents or groups with similar physical or chemical properties which produce broadly similar biological properties to another chemical compound.’ Bioisosteric replacements come in two general classes: core replacements, where the centre of a molecule is changed, creating a new chemical series; and leaf replacements where a replacement is sought for part of the molecule on the periphery, keeping it in the same series. Core replacements can get you out of a difficult ADMET or IP situation or can be a useful stepping stone to developing a backup series. Leaf replacements are more common as part of the day-to-day work of lead optimization – having discovered a substituent that works, the obvious thing to do is to try more things like it!

There are many software products available to search for bioisosteric replacements, both commercial and freely-available. The proliferation of methods is largely because searching for bioisosteres looks so easy. Look at the piece you want to replace, search a database for something largely the same size and geometry, and present the results. Simple, right? However, if you delve into it, the whole process is much more complex than it appears …

The first issue for a bioisostere replacement method involves deciding whether a fragment is physically the right size – are the connection points in the right place and at the right angle? Again this seems a trivially easy question to address until we consider that this is a 3D question. You need to consider the 3D geometry of the fragment – what conformation is it in? We chose to use a mixed approach to this: we provide both fragments derived from the CSD where the conformation comes directly from the small molecule crystal structure and large databases of fragments where we provide a small distribution of low energy conformations. As a result, you can just search experimental small molecule crystal conformations if you wish, but you can augment that search with a much larger and richer conformation database if you so desire.

The second criterion is pharmacophoric: a bioisosteric replacement must have ‘broadly similar biological properties’ and hence its interactions with the protein must be similar to the original molecule. Again, there’s a subtlety here. If I’m looking for a replacement for a triazolothiazole, then when I’m presented with a candidate fragment the obvious thing to do is ask “how similar is it to triazolothiazole?” This approach is flawed because the properties of the candidate fragment (and the initial triazolothiazole!) depend on their environment i.e. the rest of the molecule. This is especially true if you are interested in molecular electrostatics, as the electrostatic potential of a molecule is a global property that cannot generally be piecewise decomposed. That is why Spark performs all scoring in product space, not in fragment space (Figure 1). However, even if you are using a more crude method such as shape or pharmacophore similarity, treating the fragment in isolation can be very misleading.

Figure 1: Merging a new fragment can (subtly) change the electrostatics of the rest of the molecule. Red=positive, Blue=negative regions.

Working in product space brings advantages in steric considerations as well as electronic. Even though a particular fragment might have a beautiful geometric fit, and match the original core very well in terms of shape/electrostatics/color/pharmacophores/whatever, you still need to assess whether the fragment is compatible with the conformation of the original molecule. That assessment cannot be done in fragment space! In Spark we handle this in two ways. Firstly, all product molecules are minimized prior to scoring. If, due to steric or electronic effects, the candidate fragment just is not compatible with the desired conformation of the original molecule, then this minimization will distort the invariant parts and lead to a very low similarity score. Fragments which cause steric clashes or lone pair-lone pair repulsions are thus automatically filtered out (Figure 2). Secondly, we follow this up with an explicit assessment of strain energy around the newly-formed bonds, so the user can immediately see whether there are any causes for concern.

Toluyl is similar to phenyl except when it causes a steric clash
Figure 2: Toluyl is similar to phenyl, except when it causes a steric clash. Whether it does or not depends on what R is!

Another thing that can be drastically affected by the environment is the hybridization and charge state of the fragment (and the rest of the molecule). The main culprit here is nitrogen. Is an amine basic or not? Is the nitrogen pyramidal or trigonal planar? Both these questions depend strongly on what the nitrogen is attached to and hence you cannot assess the suitability of a fragment without reference to the chemical environment that you are going to place it in. These questions turn out to be surprisingly difficult to handle in a completely robust way: there is a large amount of code in Spark which assesses the local chemical environment around the newly-formed bonds and determines whether any hybridization or formal charge changes are required. It turns out you cannot just recharge the product molecule, as in some cases the user may have assigned particular charge states to parts of the molecule (based on experimental knowledge of the charge state when bound to the protein, for example) and you don’t want to undo that.

The final complication to performing the scoring properly (i.e. in product space) is accounting for the available degrees of freedom. If you are doing a core replacement, then the positon of the core is generally completely determined by the attachment points. However, if you are doing a leaf replacement, then the question arises as to what rotamer to choose around the attachment bond. The only solution is to perform a limited conformation search around that bond and score each conformer, which has the potential to get very slow although there are some computational tricks you can do to reduce the search space. It’s not just leaf groups that have to be searched: replacing the centre of a molecule might require a rotation scan if there are two attachment points which happen to be collinear.

So, our nice and simple bioisosterism calculation has become rather more complicated. Rather than just see if a fragment fits into the hole in the original molecule, we need to merge it in properly, assess any hybridization or formal charge state changes that are required, minimize the resulting molecule, perform a rotation scan around any underspecified degrees of freedom, recompute the electrostatic potentials, align to the original molecule and calculate the steric and electrostatic similarities. Only then can we decide whether the fragment is any good or not! Luckily, Spark does all of that for you, so as a user you can concentrate on the results, rather than on the complicated calculations required to produce those results.

Try it for yourself – download a free evaluation of Spark.

Electrostatics of the Topliss tree

Since its publication on J. Med. Chem in 1972, J.G. Topliss’ operational scheme for aromatic substitution,1 commonly known as ‘the Topliss tree’, has been widely used by medicinal chemists to explore the substitution pattern on a benzene ring.

This decision tree approach was designed for maximum simplicity of application by the medicinal chemists, avoiding use of computers and statistical procedures. It is based on a stepwise selection of substituents designed to progress as rapidly as possible to the most potent compounds in a series.

The Topliss tree is based on the concept, pioneered by Hansch,2,3 that the observed changes in activity following the introduction of a substituent on a benzene ring depend on its lipophilic, electronic and steric properties, which in Topliss’ original publication were represented in a quantitative manner by the physico-chemical parameters µ, σ and Es.

To complement this traditional drug design approach, we have used Torch,4 Cresset’s powerful molecular design tool for medicinal and synthetic chemists, to visualize and give new insight to the deep changes in electrostatic and hydrophobic properties introduced by the substitution patterns recommended by the Topliss tree:

Electrostatics_Topliss tree

Download ‘Electrostatics of the Topliss tree’

1. J. Med Chem. 1972, 15, 1006
2. J. Amer. Chem. Soc. 1964, 86, 1616
3. ‘Drug Design’, Vol I 1971, 271

Is this compound worth making?


In deciding if a compound is worth making medicinal chemists need to consider a number of factors including:

  • does it fit with the known structure activity relationship (SAR)?
  • does it explore new interactions or regions of space?
  • does it have good physico-chemical properties?

Summarizing the SAR content of large compound series in a single informative picture is challenging. Traditional 3D-QSAR techniques have been used for this purpose, but are known to perform poorly where the structure-activity landscape is not smooth. 3D Activity cliff analysis has the potential to explore such rugged SAR landscapes: pairs of compounds with a high similarity and a large difference in activity carry important information relating to the factors required for activity. However, looking at pairs of compounds in isolation makes it hard to pinpoint the changes which consistently lead to increase potency across a series.

Similarly, summarizing the regions and interactions that have been explored in a single picture is not trivial, especially if the electrostatic character of compounds is to be considered (e.g., replacement of electron rich with electron poor rings).

Here we present Activity Atlas, a new technique to analyse a series of compounds and derive a global view of the structure-activity relationship data.



The Regions explored summary gives a comprehensive 3D picture of the regions explored in electrostatic and shape space. A novelty score is assigned to each compound, enabling to predict whether newly designed candidates are likely to contribute additional SAR knowledge, and are thus worth making.

Average of actives summarizes electrostatic and shape properties which consistently lead to potent ligands across the series.

The Activity cliff summary gives a global view at the activity cliff landscape highlighting the regions where either more positive/negative electrostatic potential, or larger steric/hydrophobic potentials increase activity.

A traditional 3D-QSAR model2 was built on the same data set (q2 = 0.7). While 3D-QSAR seems better at extracting information where SAR is continuous, Activity Atlas gives more definition in regions where SAR requirements are critical.

Activity Atlas gives more definition in regions where SAR requirements are critical

Accounting for uncertainties in the alignment

Activity Atlas takes into account multiple alignments, weighted by their probability to be correct:

Accounting for uncertainties in the alignment1
Based on the distance from the top result

Accounting for uncertainties in the alignment2
Based on the absolute similarity value (CCDC-AZ dataset)

Analyzing multiple alignments per molecule enables the technique to account for uncertainties in the orientation of flexible side chains not represented in the reference(s).

Application to selectivity

Application to selectivity

Application to a set of adenosine receptor antagonists with activities against A1, A2a and A3 receptors3 demonstrates the utility of the technique in locating critical regions for selectivity. Examination of the steric and electrostatic maps for the three subtypes clearly shows which regions should be targeted in order to enhance subtype selectivity.

In the example above, the right hand side of the molecules can be used to discriminate between A3 and the other two subtypes, while A1 and A2a can be separated by increasing steric bulk and positive charge around the top of the molecules.


The Activity Atlas technique is a powerful way of summarizing SAR data in 3D. By combining information across multiple pairs, it enables a global view of the critical points in the activity landscape. The Average of actives summary captures in one picture the 3D requirements for potency, while the Regions explored summary enables prioritizing compounds which add crucial SAR information over trivial analogues.


1. J. Chem. Inf. Model. 2006, 46, 665-676
3. J. Chem. Inf. Model. 2011, 51, 258-266

Rapid technique for new scaffold generation II: What is the best source of inspiration?


Scaffold hopping and R-group replacement remain central tasks in medicinal chemistry for generating and protecting intellectual property. Spark is a bioisostere replacement tool (available as a desktop software application) for rapidly generating reasonable yet novel scaffold and R-group replacements using Cresset’s molecular field points.

Cresset’s field technology condenses the molecular fields down to a set of points around the molecule, termed ‘field points’. Field points are the local extrema of the electrostatic, van der Waals and hydrophobic potentials of the molecule.

field points

Spark workflow

The Spark approach uses a database of molecule fragments, or available reagents, to suggest replacements that maintain the shape and electrostatic character of a known active molecule. The user identifies the region of a known active
molecule that they wish to replace, and this piece is removed.

region of known active molecule to replace
The number of bonds broken is recorded together with the distance and angle between any pair of broken bonds. This information is used to search a database of fragment conformations for replacement moieties.

search for replacement moieties
The product molecule is energy minimized and then scored as a replacement. Scoring is performed using an average of field and shape similarity on the product molecule. Scoring the product (rather than the fragment) allows the electronic changes induced in the rest of the molecule to be taken into account.

scored as a replacement
By default, the scoring reflects the change relative to the original molecule, but the user can choose to add other molecules that can be used in the scoring. In this way compounds with sub-optimal interactions can be improved by mimicking other known actives.

Fragment sources in Spark

Spark generates bioisosteres from databases of fragments derived from:

  • Commercially available, real compounds and reagents (ZINC)
  • Theoretical aromatic rings (VEHICLe)
  • Literature reports of bioactive compounds (ChEMBL)
  • Fragments from the Cambridge Structural Database (CSD) of small molecule crystal structures

In this case study we investigate which of the fragment sources available in Spark is the best source of inspiration.

best source of inspiration
If you have access to significant proprietary chemistry, to specialized reagents, or want to consider fragments from reagents that you have in stock, then the creation of custom databases with the Spark Database Generator will enable you to exploit your own proprietary chemistry to generate and protect intellectual property.

R-group replacement to D3 antagonists

The ChEMBL ‘common’, ‘rare’ and ‘very rare’, ZINC ‘very common’, ‘common’, ‘less common’, ‘rare’, ‘very rare’, ‘singleton’ and the VEHICLe fragment databases were searched using ‘Accurate But Slow’ calculation settings. Compounds with piperazine scaffolds were filtered out as these are very well known in the literature.

known scaffolds
Known D3 scaffolds were found in ChEMBL or ZINC (commercially available compounds) databases. Novel solutions were found in the ChEMBL database.

An analysis of the chemical diversity of the known D3 scaffolds retrieved from each database clearly shows that the less common fragments derived from the literature database are a precious source of potentially
useful chemical diversity. Note that these less common fragments may be associated with more complex and less documented synthetic routes.

analysis of chemical dversity

Scaffold hopping application to Sildenafil

The ChEMBL and VEHICLe fragment databases were searched using ‘Accurate But Slow’ calculation settings. The protein structure for 1UDT was used as an excluded volume, constraining the field points associated with the interaction with glutamine (Gln817) in the 1UDT protein.

known actives were found
Known actives were found in ChEMBL and VEHICLe databases. Novel but highly plausible solutions were found in the VEHICLe database.


Spark provides both known active scaffolds and novel solutions that represent opportunities for scaffold hopping and R-group replacement.

The nature of the experiment appears to dictate the best source of fragments. It is therefore important to have a wide range of fragment sources to choose from for each experiment, to provide a balance between novelty and synthetic accessibility.

The creation of fragment databases from proprietary collections of compounds can be a powerful way of increasing the chemical diversity available to Spark.


Examining the diversity of large collections of building blocks in 3D


2D fingerprint-based Tanimoto distances are widely used for clustering due the overall good balance between speed and effectiveness. However, there are significant limitations in the ability of a 2D fingerprint-based method to capture the biological similarity between molecules, especially when conformationally flexible structures are involved. Structures which appear to largely differ in functional group decoration may give rise to quite similar
steric/electrostatic properties, which are what actually determine their recognition by biological macromolecules.

In BioBlocks’ Comprehensive Fragment Library (CFL) program, we were confronted with clustering a very large collection of scaffolds generated from first principles. Due to the largely unprecedented structures in the set and our design aim to populate the 3D ‘world’, using the best 3D metrics was critical. The structural diversity of the starting collection of about 800K heterocyclic scaffolds with variable functional group decoration was not adequately captured by 2D ECFP4 fingerprint Tanimoto distances, as shown by the rather flat distribution of 2D similarity values across the set, and by their lack of correlation with the 3D similarity metrics.

The initial step of any clustering procedure is the computation of an upper triangular matrix holding similarity values between all pairs of compounds. This step becomes computationally demanding when using 3D methods, since an optimum alignment between the molecules needs be found taking into account multiple conformers.

The presentation covers the methodological and technical solutions adopted to enable 3D clustering of such a large set of compounds. Selected examples will be presented to compare the quality and the informative content of 3D vs 2D clusters.


See presentation ‘Examining the diversity of large collections of building blocks in 3D‘ given at 250th ACS national meeting.

Is it worth making? Assessing the information content of new structures


We have recently presented a method of summarizing the information obtained from 3D activity cliff analysis: examination of all pairs of molecules can distinguish between apparent cliffs that are outliers, or due to measurement error, and those which consistently point to particular electrostatic and steric features having a large impact on activity. To do this it has proved essential to allow for alignment noise: no 3D alignment technique is perfect, so we apply a Bayesian analysis to correct for potential misalignments and for the case where a molecule is aligned correctly except for a flexible substituent whose conformation is under-constrained. We use the recent AZ/CCDC alignment validation data set to determine valid estimates for the Bayesian priors.

As an extension of this technique, it is possible to mine the data for a simple picture of explored pharmacophoric space, corrected for the conformational and alignment flexibility of each molecule. This provides an invaluable picture to the chemist of which parts of property space around a molecule have been adequately explored. When considering a new molecule for synthesis, it is possible to compute the amount that this would increase the explored pharmacophoric space and hence present an ‘information content’ score for the new molecule: if we made and tested this new molecule, how much would it actually increase the structure activity relationship (SAR) information content of the data set?

The combination of this, with the activity cliff summary data, allows a simple qualitative evaluation of the SAR of a data set in 3D, alongside guidance on which parts of pharmacophoric space have been mined out and which remain underexplored. We present the application of these techniques to several literature data sets.


See presentation ‘Is it worth making? Assessing the information content of new structures‘ given at the 250th ACS National Meeting.

Scaffold hopping into new DPP-IV protease inhibitors


Cresset’s powerful scaffold hopping and fragment replacement software can be used to mimic an existing ‘active’ molecule’s electrostatic patterns to quickly, and efficiently, generate new molecular designs. The outcome, based on 3D molecular electrostatic similarity, is more biologically relevant than that from other similarity metrics. Here, we show how Spark can be used to rapidly generate potentially valuable chemical ideas for new DPP-IV inhibitors.


Two peptide hormones, GLP-1 and GIP, mediate lowering of blood glucose level through stimulation of insulin release and inhibition of Glucagon release. DPP-IV cleaves a dipeptide from the N-terminus of both GLP-1 and GIP hormones to give their inactive forms, thus abolishing their glucose lowering action. DPP-IV inhibitors have been shown to be important agents useful for treating type II diabetes.

Scaffold hopping methodology

The method involves:

  • An initial bioactive conformation: Preferably a potent target molecule of interest. This can either be extracted from an X-ray or modeled (e.g. from a pharmacophore or alignment / binding hypothesis).
  • Field pattern generation: Cresset’s proprietary XED force field used to generate electrostatic properties.
  • Database of molecule fragments: Spark uses an internal database of fragments derived from ChEMBL and ZINC. Custom sets can be generated from proprietary chemistry.
  • Automated reconstruction of the new 3D molecules: Can specify linking chemistry and replacement sites (e.g. scaffold or decoration).
  • Aligns and scores output as full minimized molecules: Using ’field similarity’ and shape similarity relative to the starting template(s).
  • Protein target can be used as excluded volume: Can take account of protein pocket by penalising examples with steric clashes.
  • Filter output on physiochemical properties.

DPP-IV X-ray ligand electrostatic and shape similarity analysis

A plethora of molecules are already known which inhibit DPP-IV and have useful anti-diabetic properties.

Hierarchical clustering of DPP-IV x-ray ligands by electrostatic and shape similarity
Figure 1. Hierarchical clustering of DPP-IV x-ray ligands by ‘electrostatic and shape similarity’.

Superimposition of over 20 published x-ray crystal structures, followed by hierarchical clustering using all-by-all field similarity on the ligands, reveals two main clusters and a distinct mode of binding for a fluoro-olefin example (Figure 1). This unique clustering, performed in Cresset’s ligand-focused workbench, Forge, allowed the selection of three distinct inhibitors for further scaffold hopping work.

Experiment and results

Alogliptin (1), Omarigliptin (2) and the fluoro-olefin (3, PDB: 3C45) represent some of the most ligand efficient examples from these clusters. Two experiments were performed in Spark using (1) and (2) in a simple scaffold hopping exercise. A final experiment was a chemotype merging experiment: a truncated (2) was used, with (3) as a second template, to find molecules bridging the two series. Workflow as shown in Figure 2.

Example results (Table 1) shows the diverse range of output suggestions provided for new chemistry and validates the method by providing examples which already have precedents in patents and the literature.

Scaffold hopping and merging_input structures fields field points and output examples
Figure 2. Scaffold hopping and merging – input structures, fields, field points and output examples.

Diverse range of output suggestions provided for new chemistry
Table 1. Scaffold hopping and merging examples of output 2D structures.


The Spark searches output a wide range of diverse chemotypes that included active or very close architectures to known active frameworks. A number of these had been discovered through HTS rather than through rational design. Tight control over the chemistry ensures that feasible chemistry is provided from known fragments whilst maintaining the features necessary for activity.


Spark is a powerful molecular modeling tool for the rapid virtual elaboration of scaffold ideas; either in scaffold hopping, merging, fragment growing or linking experiments. Applying Spark to chemotypes bound to the active site of DPP-IV provided a range of interesting and synthetically-feasible suggestions.


Green et al, Diabetes and Vascular Disease Research 2006, 3 No3, p159-65.

Protein pictures were rendered using open source Pymol from Delano Scientific.

[This content was presented as a poster at Proteinase 2015.]

Applying the XED molecular mechanics force field to the binding mechanism of GPCRs (part 3 of 3)

Agonist and antagonist differences emerge for the β2 Adrenergic GPCR

Part 3 reports on work carried out up to January 2015 and continues from Part 1 (posted March 2014) and Part 2 (posted April 2014). Since January 2015, further work has revealed more information that is being written up for an academic paper.

The story so far

In Part 2 of this series, published online in April 2014, I described the ‘DRY’ lock mechanism and proposed that it was the best tool for monitoring the conversion of the β2AR from the active to inactive form. I also undertook to report how I would investigate the mechanism of the G-protein entry and exit from the active receptor.

However, “it is a truth universally acknowledged, that a scientific report distributed whilst relevant work is still in progress will be soon outdated!” My research is ongoing and my views have changed as a consequence. In reflection of this I have decided to continue reporting as I proceed.

Since posting Parts 1 and 2, I have done more work that showed that my initial results were neither consistent nor reproducible. I realize that this is because more preparatory care was needed on the input structures. I have carried out extensive work to address this problem and this work has revealed results that begin to suggest that we can distinguish agonists from antagonists and explain the way neutral antagonists work.

The continuing work also suggests that the ‘DRY’ lock may not be as important as I first thought, at least for the β2AR. Indeed, the definition of Asp130 and Arg131 (both on the receptor) as the accepted ‘DR’ part of the DRY lock is open to some confusion in the literature. I concluded that the Asp130/Arg131 ion pair is not appropriate as an indicator of active/passive receptor states and that a DR ion pair between Arg131 and Glu392 (on the G protein) may be more useful.

In this blog I report on the tentative finding that agonists stabilize the full complex and fight to keep the Arg131/Glu392 closed whilst antagonists destabilize the complex and encourage the ion pair to open and release the G-protein.

Making sure the input data set is sound

The majority of the effort involved in protein computation is to do with making sure the quality of the input structure is chemically sound. Most input data errors stem from incomplete and poorly refined X-ray coordinates.

Molecular mechanics can be very forgiving about yielding answers from poor structures but can fail if atom types are badly defined or bits of chemistry are missing. For example, a molecule of water and a piece of protein can be so close in the X-ray that they are inside a known physical barrier. Furthermore, along with other abstruse errors, many X-ray structures:

  • Are missing entire residues
  • Give a limited idea of which way protonic hydrogens may point
  • Include the alpha carbons of named residues without the residue coordinates
  • Have poorly defined terminations
  • Cannot distinguish the orientation of terminal amide groups (ASN and GLN).

The most serious problems stem from the inability of most protein X-ray electron density maps to ‘see’ hydrogen atoms. Among other things, this means that the positions of hydrogen bonds, which are a major part of protein structure and action, have to be guessed, along with tautomeric states and formally chargeable ions. Group pKa values can vary drastically when imbedded in ‘protein space’ and distorting effects from neighbouring molecules in a crystal are rarely taken into account.

Where possible, correction routines have been developed to resolve these errors or at least highlight their presence. Plainly, entire missing residues cannot sensibly be replaced. General corrections to formal charge states can be applied with respect to the extent of solvation but pKa variations from known aqueous values are not yet tractable.

The general strategy for this work was to ‘clean’ the X-ray structure as far as possible and fully optimize every β2AR/G-protein (R*) complex using XED molecular mechanics before any analysis was attempted. Essentially, this strategy is a means of ‘normalizing’ all structures.

Molecular mechanics issues

Unlike quantum mechanics which is ab initio, molecular mechanics provides estimates of the strain in a structure based on empirical data from experimental observations. It returns a ‘strain energy’ (referred to as just the ‘energy’ in the bulk of this report). This ‘energy’ reflects how far the structure is from optimum but has no real physical counterpart.

Because entropic contributions are not included, the ‘strain energy’ should only be used for relative comparisons, e.g., the energies of two conformers of the same structure would be expected to have the same entropic contribution which, we assume, cancels on taking the energy difference. The same approximations apply to optimising structures by mathematically minimising their strain energy and comparing the difference in their final optimized energy.

The value of this ‘strain energy’ is unique to each molecular mechanics implementation. Since its usefulness rises with the quality of the molecular mechanics force field, I am keen to know if the XED force field can return valid conclusions from proteins and other large structures given its ability to take account of electron polarizability and improved handling of conjugated and aromatic systems.

The chosen ligands

A modified version of the table of chosen ligands from Part 2 is presented in Supplementary Material Appendix 1. Two main types of G-protein couple with the β2AR:

  • Gi, which results in adenylate cyclase (AC) inhibition
  • Gs, which stimulates it.

Current pharmacology strongly suggests the existence of at least three active β2AR* conformations that can either couple with both the Gi and Gs proteins or couple with one or the other.

The ‘action’ column in the table in Appendix 1 shows that propranolol inhibits Gs coupling without seemingly affecting Gi. Clenbuterol seems to agonise through Gi control only. ICI-118551 inhibits Gs and stimulates Gi. Since Gi stimulation acts as an AC controlling feedback mechanism, the stimulation of both Gi and Gs by some agonists and antagonists makes the agonist/antagonist distinction a matter of degree rather than definitive.

The main target for β2AR is the adenylate cyclase pathway. However, this versatile receptor stimulates other paths (MAPK, certain kinases and β-arrestin) but probably independently of the G-proteins and in some cases only after phosphorylation.

Improving the starting ligand conformations

As mentioned above, the early work could not be reproduced. This was apparent whenever small changes were applied to the conformation of the GR* complex or its ligand. The problem arises when minimising a large system with very many possible energy minima. Because minimisation merely moves energy downwards, the resulting optimized structure reflects the starting conformation. Therefore, if the starting conformation is wrong, the final structure will be wrong.

Simulated annealing techniques may overcome this problem somewhat, and dynamics should do even better, but the principle applies across the board, especially if the quality of the molecular mechanics is poor. It seemed reasonable to assume that more effort would be needed to get the starting conformations right.

Between 45 and 80 randomly different conformations for each ligand were placed in the active site of the β2AR* crystal structure (truncated 3SN6.pdb coordinates) and the complete complex was minimized to a 0.02 rms exit limit. Supplementary Material Appendix 2 shows the plots of the binding energy of each ligand versus the total complex energy.

In principle, the conformation with the lowest total complex energy was deemed to be the desired starting conformation for each ligand. Using this criterion, new starting conformations for carazolol, propranolol, ICI-118551, and salbutamol were isolated from their conformational packs and placed into a ‘starting’ set.

However, two discrete conformational regions were identified for POG, adrenaline and clenbuterol and the lowest complex energy conformation was kept from each of the clusters. Two clusters were clearly discernible from the neutral antagonist ‘Chem 7929193’ (7929193-Chembridge) conformational collection but appeared much more loosely ordered than other ligands. Two relevant conformations were added to the starting set.

Analysing all minimized conformations

I now had a set of conformations of the fully optimized total complex (β2AR*+ligand) for each of the chosen ligands. The next question to ask was, are the conformers of this starting set related by structure?

Figure 7 summarizes the comparison of the Cα backbone of each conformer against all others. It may be inferred from this comparison that the inverse agonists share one unique conformation and the agonists share another.

Figure 7 The Cα RMSD of each conformer across all others. The cut-off to discriminate ‘same’ conformers from ‘different’ conformers was taken as 0.7rms.
Figure 8 summarizes the analysis of these clusters with respect to their relative complex energies, the experimental pharmacology and the way each ligand affects the distance between the ‘DR’ (Arg131/Glu392) ion pair lock. This ‘DR’ lock (reported in Part 2) forms an ion pair between the G-protein and the receptor and is obviously going to be closed when the G-protein is bound. There may be variation in the other ion pairs that the two residues form after the G-protein leaves the receptor.

Summary of analysis of collected confs
Figure 8. Summary of the analysis of Figure 7. The ligands in red are inverse agonists, those in green are agonists and ‘MT’ is the apo-complex (no ligand). The blue ligand is a neutral antagonist having no observable effect on the complex other than blocking its active site. Etot is the fully minimized complex energy for each ligand conformer started from the 3SN6.pdb truncated to include the G-protein surrogate and the receptor chains only. The estimated error attached to the full complex energy values is between ±10 and ±20 kcal/mole for approximately 6000 atoms and 6000 ‘xeds’ (polarisable extended electron constructs unique to XED molecular mechanics). The distance (Å) in column 3 is not the usual Arg131/Asp130, which shows little variation, but the alternative DR ‘lock’ across Arg131on the receptor and Glu392 on the G-protein.
The results shown in Figure 8 suggest that:

  1. The inverse agonists Carazolol, Propranolol and ICI-118551 destabilize the complex by ~20 kcal/mole over the empty complex and show the largest potential lock monitor distance
  2. The agonists POG, Adrenaline, Salbutamol and Clenbuterol decrease the energy over the empty complex by ~50 kcals/mole and show a tendency to decrease the potential lock monitor distance
  3. The neutral antagonist Chem 7929193 can stabilize two conformers with little change in binding energy, one of which decreases the stability over the empty complex by >100kcals/mole in line with inverse agonist action whilst the other stabilizes the complex by ~40kcal/mole in line with the agonists.

Chem 7929193 has been shown to be a neutral antagonist. Received wisdom, based on the ‘two state’ receptor theory, concludes that this category of antagonist rarely occurs as it necessitates agonism and inverse agonism to exactly cancel out. Our investigations of both the low and high energy forms of Chem 7929193 reveal conformers that seem to agonise and antagonise at the same time! Surely this must confuse the receptor and maybe bring it to a standstill whilst stopping any other ligand from entering the receptor.

In summary

If we assume that the conformations of the eight β2AR active complexes are correct, this work shows that there seems to be a relationship between the ligand type and its affect on G-protein/receptor complex stability. Why agonists appear to stabilize the complex more than inverse agonists is not yet clear. The ability of agonists to hold the Arg131/Asp130 DR lock more open than inverse agonists is intuitively attractive but inconsistently maintained across other GPCRs and unconfirmed in the present study.

Further work will concentrate on expanding the ligand base to see if the principle holds up, namely that G-protein/receptor stability is enhanced by agonists and diminished by inverse agonists.

Despite the lack of surrounding structure (membranes, water, sugars, etc.), the good correlation of the pharmacology of eight diverse β2AR* ligands with energies calculated using molecular mechanics is remarkable. Justification must reside in the fact that we are using a single system which is probably allowing cancellation of entropy, enthalpy, solvation and other neglected effects.

Supplementary Material

Appendix 1 – The β2AR ligands chosen for this study


Appendix 2 – Binding and Total Ligand Energies

Between 45 and 80 randomly different conformations for each ligand were placed in the active site of the β2AR crystal structure (truncated 3SN6.pdb coordinates) and the complete complex minimized to a 0.02 rms exit limit.


Summarizing activity and selectivity in one picture

Summarizing activity and selectivity in one picture was presented by Dr Mark Mackey, Cresset at the Cresset European User Group Meeting 2015.


3D-QSAR based on molecular interaction potentials can provide a wealth of information about the exact molecular characteristics required for activity. However, current techniques have a number of issues such as alignment noise, sampling errors and descriptor choice which can make it difficult to reliably produce effective models. We have presented in the past techniques for solving the sampling problem and shown that using accurate electrostatics combined with simple shape descriptors often gives meaningful models. However, there are still times when it is not possible to obtain a statistically valid linear regression model.

One useful qualitative data analysis method that is being increasingly used is activity cliffs analysis. In this technique, pairs of compounds are located that are similar (in some sense), but have different activities. Traditionally activity cliff analysis has used a 2D definition of similarity, but extension to 3D similarity metrics gives additional information that is very useful to locate the source of and reason for the activity differences.

An extension of 3D activity cliff analysis is to mine the entire data set for corresponding cliffs and use this to build a model for activity. Analysis of the data set to locate activity cliffs locates the pairs of molecules with the highest information content. However, this needs to be tempered with an analysis of how likely it is that the molecules are aligned correctly, as only properly-aligned molecules contain any information. We apply Bayesian corrections to the activity cliff data to obtain a map of the electrostatic and shape characteristics that seem to locally correlate with improving activity. The resulting model is semi-quantitative in that it attempts to describe the entire data set without building a linear regression model. This technique provides a valuable fall back to the computational chemist for information extraction from ligands in 3D.