Presentations and posters from 253rd American Chemical Society National Meeting

Download our presentation and posters from the 253rd ACS National Meeting. Also available is the presentation, and poster, by the University of Bristol which features Flare, our new structure-based design application.

Improving new molecule design using electrostatics

MEDI poster 135

Tim Cheeseright, Director of Products, Cresset

Electrostatics are critical to ligand binding and yet largely overlooked in new molecule design due to the difficulty in calculation and visualization of meaningful potentials. We have previously shown how electrostatics can be used effectively for scaffold hopping, virtual screening, ligand alignment and SAR interpretation. In this poster we will focus on ligand and protein electrostatics. We will show how considering the changes in ligand electrostatics improves the outcome for new molecule design. Going beyond traditional H-bonding based pharmacophore descriptors enables designers to map the effect of molecular changes on the full electrostatic potential of their molecules. Full exploitation of aromatic dipole moments, C-H hydrogen bonding, halogen bonds, and pi-pi interactions is only possible by understanding the electrostatic basis of these effects. Consideration of the protein electrostatics can inform ligand design to generate complementary patterns of electron rich and electron poor regions.

Putting electrostatics and water at the center of structure-based drug design

COMP poster 292

Tim Cheeseright, Director of Products, Cresset

The electrostatics of protein active sites can inform ligand design and SAR. However, no calculation of electrostatic potentials in the active site can be complete without also considering whether water molecules are tightly bound and contributing to the potentials and ligand binding. In this poster we will explore the effect of including water molecules shown to be energetically favorable on the electrostatic potential of individual proteins. We will present a new application to enable rapid and accurate calculation of water stability and protein interaction potentials. Combining these analyses with the popular Waterswap technique results in an in-depth knowledge of protein targets and ligand protein binding.

Combining protein interaction potentials with water analysis in structure-based design

COMP oral 497

Tim Cheeseright, Director of Products, Cresset

The XED molecular mechanics force field is unique in providing off-atom centered charges that enable studying complex charge distributions in ligands and proteins alike. We have previously described the use of the molecular interaction potentials derived from the XED force field for ligand similarity calculations which have, in turn, enabled virtual screening, scaffold hopping and SAR analysis to be performed using these effective and meaningful descriptors. In this talk we will describe in detail our application of the XED force field to proteins. We will focus on the calculation of protein interaction potentials showing how protein electrostatics provide deep insights into ligand binding that would otherwise be missed. In parallel to the calculation of the protein interaction potentials we describe using XED to minimize protein-ligand complexes and demonstrate the synergy between these methods.

Understanding protein-ligand binding at the molecular level: Using swap-based methods to visualise binding free energy components

COMP oral 103

Christopher Woods, University of Bristol

Binding free energy methods allow computational chemists to predict whether or not changes to a ligand increase or decrease its binding affinity to a medicinally important protein. While these methods can reveal whether or not a change to the ligand increases its binding affinity, they provide little information as to how this change affects individual molecular interactions. Such information would be extremely useful, as it could provide feedback to a drug designer that could inspire the next series of ligand modifications. For example, it would be valuable for the free energy calculation to reveal that addition of a hydroxyl group to a ligand strengthens its interaction with an active site aspartate residue, but simultaneously destabilises a water molecule that bridges between it and the protein. This feedback could inspire a drug designer to investigate ligand modifications that combine addition of the hydroxyl group with the addition of moeities that displace the bridging water molecule. We have developed a range of binding free energy methods, based on WaterSwap. These swap-based methods allow calculation of absolute and relative protein-ligand binding free energies using a single simulation over a single λ-coordinate. Use of a single coordinate allows components of the binding free energy (i.e. specific interactions between the ligand and individual active-site residues, or between the ligand and neighbouring water molecules) to be integrated across λ during the simulation. The resulting components are used to colour-code 3D views of the proteinligand system. This allows drug designers to easily visualise the affect of ligand modifications on the interactions between the ligand and individual residues in the protein, and individual water molecules in the binding site. These components are not true free energies. However, what they reveal is chemically intuitive, and provide a level of insight that allows drug designers to suggest modifications that lead to improved binding affinity. The components show cooperative affects, and also reveal how increasing the strength of interaction

Visualising the molecular drivers behind drug resistance

COMP Sci-Mix poster 347

Christopher Woods, University of Bristol

The huge expense of developing a new drug can be wasted if natural mutations of amino acid residues in the targetted protein lead to a loss of drug binding affinity. Rational drug design is a continual struggle, with evolution driving mutations that develop emerging drug resistance into widespread drug inefficacy. We have developed a new free energy method, based on WaterSwap, that can provide drug designers with the insight needed to understand how protein mutations affect drug binding. This method allows the change in binding free energy of the drug associated with the mutation of the protein to be predicted. In addition, as for other Swap-Based methods, this free energy change can be decomposed into components that can be used for visualisation. These components allow the change in binding free energy to be understood in terms of changes in specific ligand-protein interactions, or changes in the solvating water network in the active site. This method allows drug designers to pro-actively screen a new drug computationally against a range of likely protein mutants, thereby enabling the drug designer to get one step ahead of nature. Alternatively, it allows drug designers to investigate the molecular basis for reduced binding affinity in known drug-resistant mutants of a protein. This information would be useful for the development of the subsequent generations of the drug.

3D-RISM: Better water positions through improved electrostatics?

Professor Dr Paolo Tosco presented 3D-RISM: Better water positions through improved electrostatics? at the German Chemoinformatics Conference, Fulda, Germany. on 7th November.

Abstract

The Reference Interaction Site Model (RISM) is a modern approach to solvation based on the integral equation theory of liquids[1]. It has seen increasing use as a method to analyse the structure of water in and around protein active sites. The position and energetics of water molecules in and around the active site is known to be of crucial importance in understanding ligand binding, and in particular knowledge of which water molecules are tightly bound and which are energetically unfavourable can give valuable insights into structure-activity relationships.

The accuracy of 3D-RISM depends on the closure used and the functions used to compute the intermolecular potential between the protein and water molecules. In practise, the potential function consists of van der Waals interactions and electrostatic interactions. Increasing the accuracy of electrostatic calculations should therefore give more accurate 3D-RISM water placements and energies.

To date, 3D-RISM has been applied using traditional force fields such as AMBER[2] which utilise simple unpolarised atom-centred charges to represent electrostatic potentials. We have investigated whether a more accurate electrostatic model, XED[3], which provides both polarizability and electronic anisotropy, improves the effectiveness of the 3D-RISM method. As water energetics around a protein active site are not in general experimentally observable, we have in the first instance concentrated on computing water positions and energetics around a range of small molecule model systems, and validated these using DFT calculations with ONETEP[4] and interaction statistics from the CSD. Polarization and electronic anisotropy are found to be particularly important when assessing water structure around carbonyl groups, and the results from AMBER are found to be noticeably unrealistic in many cases. This has important implications for the accuracy of 3D-RISM calculations on protein active sites.

Read about progress in structure-based design at Cresset.

[1] Skyner, R.; McDonagh, J. L., Groom, C. R., van Mourik, T., Mitchell, J. B. O.; Groom, C. R.; Van Mourik, T.; Mitchell, J. B. O. A Review of Methods for the Calculation of Solution Free Energies and the Modelling of Systems in Solution. Phys. Chem. Chem. Phys 2015, 17 (9): 6174.

[2] J.A. Maier; C. Martinez; K. Kasavajhala; L. Wickstrom; K.E. Hauser; C. Simmerling. ff14SB: Improving the Accuracy of Protein Side Chain and Backbone Parameters from ff99SB. J. Chem. Theor. Comput. 2015, 11:3696-3713

[3] J.G. Vinter Extended Electron Distributions Applied to the Molecular Mechanics of some Intermolecular Interactions. J. Comput.-Aided Mol. Des. 1994, 8: 653-668

[4] C.-K. Skylaris; P. D. Haynes; A. A. Mostofi; M. C. Payne Introducing ONETEP: Linear-scaling density functional simulations on parallel computers. J. Chem. Phys. 2005, 122:084119


paolo-tosco-at-gcc-2016
Professor Dr Paolo Tosco, Cresset, Presenting at the German Chemoinformatics Conference
(Photo courtesy of Wendy Warr)

At the heart of medicinal chemistry design: The synergistic relationship between Forge and Spark

Abstract

We present an exercise in in silico medicinal chemistry: using Activity Atlas1 to assess and understand the features that drive activity for a collection of CDK2 inhibitors; using Spark2 to suggest new ideas to explore; then visually assessing those new ideas in Forge3 in the context of the qualitative 3D-SAR models generated by Activity Atlas.

Introduction

The efforts of medicinal chemists in their quest to design new, active ligands can be summarized by answering three deceptively simple questions:

  1. What to make next?
  2. Why should I make it?
  3. Has this chemical space been explored previously?

These questions are frequently addressed separately. However, there is much value in considering the questions simultaneously with the help of Forge and Spark.

In this case study, Cresset ligand-based design tools are used to address the questions above with respect to the well-studied target for oncology Cyclin-dependent kinase 2 (CDK2).

Method

Compound collection and alignment

A collection of 38 imidazopyridine containing CDK2 inhibitors was gathered from CHEMBL4 and curated to create a molecular dataset spanning a pIC50 range of 4-9.

Prior to assessing any relationship between activity and 3D structure, it is necessary to align the collection of molecules. In this case we chose to align the molecules to the ligand derived from the crystal structure 1OIT (CHEMBL73303: pIC50 =9.0). As many of the compounds in the dataset contained a solvent-exposed group that was not present in this reference an additional active ligand was identified (CHEMBL70808: pIC50 = 8.52) that maintained the core, but expanded beyond the terminal sulfonamide of the 1OIT ligand. The CHEMBL70808 ligand was aligned to the 1OIT ligand using manual and automated methods to ensure good correspondence of the imidazo[1,2-a]pyridine groups that are deeply embedded in the binding site. Further manual alignment was employed to ensure comparable orientation of the sulfonamide, despite being solvent-exposed. The pair of ligands (Figure 1) was then used as a combined, equally weighted, reference for the subsequent alignment of the 38 compounds gathered.

Figure 1_ Alignment of PDB1OIT ligand
Figure 1. Alignment of PDB:1OIT ligand (pink) and CHEMBL70808 (grey).

Compound collection and alignment

A collection of 38 imidazopyridine containing CDK2 inhibitors was gathered from CHEMBL4 and curated to create a molecular dataset spanning a pIC50 range of 4-9.

Prior to assessing any relationship between activity and 3D structure, it is necessary to align the collection of molecules. In this case we chose to align the molecules to the ligand derived from the crystal structure 1OIT (CHEMBL73303: pIC50 =9.0). As many of the compounds in the dataset contained a solvent-exposed group that was not present in this reference from an additional active ligand was identified (CHEMBL70808: pIC50 = 8.52) that maintained the core, but expanded beyond the terminal sulfonamide of the 1OIT ligand. The CHEMBL70808 ligand was aligned to the 1OIT ligand and subjected to additional manual alignment in order to better align the imidazo[1,2-

a]pyridine groups that are deeply embedded in the binding site. Further manual alignment was employed to ensure comparable orientation of the sulfonamide, despite being solvent-exposed. The pair of ligands (Figure 1) was then used as a combined, equally weighted, reference for the subsequent alignment of the 38 compounds gathered.

Within Forge, the collection of CDK2 inhibitors was aligned to the two equally weighted references, using 50% fields and 50% shape and a soft protein excluded volume to generate an alignment (Figure 2) and similarity score (range 0.59 through 0.85). As alignment is at the heart of every 3D-QSAR approach, the alignments were visually inspected and manually tweaked to ensure consistent orientation of all side chains.

Figure 2_ Alignment of 38 CDK2 inhibitors to PDB1OIT and CHEMBL70808
Figure 2. Alignment of 38 CDK2 inhibitors to PDB:1OIT and CHEMBL70808.

Modeling 3D SAR with Activity Atlas

To generate a qualitative assessment of the features driving activity, an Activity Atlas model calculation was performed.

Activity Atlas is available in Forge, and generates visually striking maps based on the aligned molecules. The models depict: the average electrostatics and shape for active compounds; the summary of activity cliffs in electrostatics, shape and hydrophobicity; and the regions explored. The models generated by Activity Atlas are intended to help answer the questions asked earlier and are shown in Figures 3-5.

In Figure 3, the average active molecule map is depicted. It indicates the common electrostatics and hydrophobic features of active molecules in the dataset – essentially molecules that do not have these features are unlikely to have any significant activity. The central H-bond donor (amine)/H-bond acceptor (pyrimidine) motif mapping the interaction with the hinge region of CDK2 is present in all actives. Negative regions near the H-bond acceptors of sulfonamide and imidazopyridine are also present on all actives.

Figure 3_ Activity Atlas map for average electrostatics and hydrophobics for active molecules
Figure 3. Activity Atlas map for average electrostatics and hydrophobics for active molecules. Blue = negative electrostatics; Red = positive electrostatics; Gold = hydrophobics.

Figure 4 shows the summary of activity cliffs in electrostatics and steric space – the fine detail on the average of actives map. Areas that are blue suggest that making this area more negative (or less positive) could enhance activity; areas that are red suggest that making the area more positive (or less negative) can drive activity; and finally, steric bulk is favorable in the green regions, whereas it is not tolerated in the pink regions.

Figure 4 (right) displays the activity cliff analysis in an orientation better suited to examining the imidazopyridine. Tilted on its side, the Activity Atlas model suggests that a more positive/less negative π-cloud should enhance activity, as well as noting that a stronger negative near the H-bond accepting ring nitrogen should also favor activity.

Figure 4_ Two views of the Activity cliff summary for electrostatics and sterics
Figure 4. Two views of the Activity cliff summary for electrostatics and sterics. Left: Electrostatics and sterics. Right: Electrostatics only in a view rotated around the x axis. Areas in blue suggest that negative electrostatics are favored; areas in red favor positive electrostatics. Areas in green are favorable sterically, while areas in pink are unfavorable for sterically bulky groups.

The assessment of the electrostatic regions explored suggests that at least 10 of 38 molecules have contributions as shown in Figure 5. This means that regions that are not explored may present opportunities for modification and further optimization.

The region explored analysis also performs a calculation of novelty for all the data set compounds as well and novel designs. The calculated novelty can be used as an indication for how much information would be gained from the molecule should it be made.

Figure 5_ Aligned ligands with a display of the regions of electrostatic space explored
Figure 5: Aligned ligands with a display of the regions of electrostatic space explored.

Generating ideas in Spark

With a SAR analysis in hand, Spark was used to find bioisosteric replacements for the highlighted portion of

the 1OIT ligand as shown in Figure 6. The fragments were sourced from the Spark reagent database of boronic acids derived from eMolecules.

Figure 6_ The Spark eMolecules boronic acids reagent database
Figure 6. The Spark eMolecules boronic acids reagent database was designated to identify bioisosteric replacements for the moiety shown in pink in the PDB:1OIT ligand.

Results and analysis

Four results were selected from the Spark experiment were selected for further analysis and are shown in Figure 7. Each of the selected Spark hits have been chosen to either enhance activity, and/or to challenge the model since it was built on a relatively small number of compounds.

Each of the proposed bioisosteric replacements are derived from commercially available reagents that are shipped from the supplier within 2 days to 4 weeks.

In (a), the pyrazolopyridine replacement would test the importance of the H-bond acceptor strength on the original imidazopyridine. Moving the heteroatom across the ring results in a weaker H-bond acceptor while maintaining a reasonable density above the ring. Note that this replacement results in a richer central pyrimidine ring.

The dihdropyrroloimidazole (b) has similar electrostatics to that starting ligand but tests the requirement for an electron deficient aromatic by saturating this ring, removing all pi-electrons from this region.

In (c), the benzofuran might be suitable for testing the model as it removes the strong H-bond acceptor altogether, replacing it with a weak feature. Whilst this might not be productive in activity the molecule would add knowledge to the model and the reagent is available immediately.

Finally, (d) suggests replacing the imidazopyridine with a quinoline. The change from 5,6 to 6,6 ring introduces a different geometry for the H-bond acceptor but maintains many of the electrostatic features of the starting ligand in a reagent that is available for immediate dispatch. If this substitution was tolerated then the spark results contain other, more exotic suggestions that include a H-bond acceptor in a 6 membered ring in this position.

Figure7
Figure 7. The starting ligand from 1OIT and selected Spark results (a-d). Top: in 2D. Center: in 3D showing positive and negative interaction potentials. Bottom: presented in the context of the Activity Atlas activity cliff summary for electrostatics.

Conclusion

This article presents our approach to answering the three questions at the heart of medicinal chemistry design:

1. What to make next?

Every chemist has a collection of tricks and substitutions that have been built up from their laboratory experience. Spark enhances the generation of ideas in an unbiased way, and in combination with a chemist’s intuition, can provide ideas about what compounds can be made that are in potentially open IP space, but still retain the characteristics of active molecules. The Spark experiment using eMolecules reagent databases also provides tier and ordering information for applicable reagents exploiting specific chemical reactions to improve laboratory efficiency.

2. Why should it be made?

Taking the Spark and chemist-designed ideas and examining them in the context of the Activity Atlas models is a way to ensure that new compounds meet the characteristics of active molecules, or conversely, are designed to test specific aspects of the models. The combination of average active molecule and activity cliff summary maps can provide the rationale needed to give confidence that taking a new design to the laboratory is worth the effort.

3. Has this chemical space been explored previously?

The region explored analysis combined with the novelty calculation within Activity Atlas can be used to assess whether the Spark and chemist-design ideas lie within the already explored chemical space, or whether they map regions of this space so far unexplored.

References and links

http://www.cresset-group.com/activity-atlas/

http://www.cresset-group.com/products/spark

http://www.cresset-group.com/products/forge

4 P. Bento, A. Gaulton, A. Hersey, L.J. Bellis, J. Chambers, M. Davies, F.A. Krüger, Y. Light, L. Mak, S. McGlinchey, M. Nowotka, G. Papadatos, R. Santos and J.P. Overington (2014) ‘The ChEMBL bioactivity database: an update.’ Nucleic Acids Res., 42 1083-1090.

Understanding torsions in Spark 10.4

Spark is the most advanced and complex bioisostere finding tool available today. While it excels at finding bioisosteric replacements that match the steric, pharmacophoric and electronic nature of the original molecule, we realize that the output of a Spark search is only the beginning. Spark suggests new molecules, so someone needs to make them. Given the time and effort that goes into even a simple synthesis, we need to make sure that Spark’s suggestions are as robust as they can be to maximize the changes of success for each Spark result.

Spark’s internal algorithms do a lot of complicated work to ensure that the result molecules are chemically sensible. When you are stitching a fragment into a new molecule there are many subtle effects that you have to take account of:

  • Does the hybridization state of any of the atoms change (especially important for nitrogens, which may shift from pyramidal to planar)?
  • Does the pKa of any ionizable groups change, and if so do we need to re-assign protonation states?
  • How do the electrostatic properties of the fragment change once it is stitched into the rest of the molecule (and vice versa)?
  • If any newly-created bonds are free to rotate, what is the rotamer that maximizes the similarity of the new product molecule to the starting point?
  • Is this conformation energetically sensible?

On this last point, Spark carries out a rotation around the newly-formed bond in order to estimate the amount of strain energy carried by that bond in the assigned torsion angle. For speed, this rotational scan is performed holding the parts of the molecule on each side of the bond rigid, and as a result the computed strain energy is only an estimate. However, even this estimated strain energy can be very useful to flag up cases where the new molecule is in an energetically unfavorable conformation and hence would need further investigation before it could be recommended for synthesis.

While this purely computational procedure works well, it would be nice to bring in some more empirical knowledge gleaned from the vast numbers of small molecule crystal structures that are available in the Cambridge Structural Database (CSD). Luckily, this analysis has already been done by Prof. Rarey’s group at the University of Hamburg, resulting in a pair of papers detailing a hierarchical set of rules to determine preferred values for torsion angles in small drug-like molecules.1, 2 This rule set, called the Torsion Library, has been incorporated into the most recent release of Spark (version 10.4).

Whenever a new product molecule is formed, Spark applies the Torsion Library rules to the newly-formed bonds, and highlights cases where the torsion angle is not one that is frequently observed in the CSD. This doesn’t automatically exclude that result from consideration (there may be a reason such as steric clashes for the uncommon torsion), but it does flag that result molecule as needing careful inspection. The Torsion Library results are automatically displayed for all results and can be used in Spark’s filters like any other result column.

As a quick example of the usefulness of the Torsion Library, we have re-examined the 2011 case study FieldStere V3.0 Example 2: Fragment Growing. In this study we were performing a fragment growing experiment in p38, using an existing ligand to guide the growth of a fragment towards the hinge binding region (Figure 1). If we had particular chemistries in mind, then instead of searching the general reagent databases, we could search specific reagent databases to find what commercially-available reagents we could use.

Fragment growing Spark experiment
Figure 1: Fragment growing Spark experiment.

In this case, we assume that we could grow the fragment either through addition of a thiol (leading to a sulfur linker) or via addition of an amine (leading to a nitrogen linker). Both searches lead to a variety of interesting-looking results which look potentially active. However, analysis of the Torsion Library result reveals significant differences. Looking at the amine linker results (Figure 2), 493 of the top 500 results have only a ‘Medium’ Torsion Library score, indicating that to get the amino substituents to reach to the hinge binding region you have to twist the amine into a moderately unfavorable conformation.

Amine linkers
Figure 2: Amine linkers.

However, when we look at the thioether result set (Figure 3), the majority (327/500) of the results are in the ‘High’ frequency category.

Thioether linkers
Figure 3: Thioether linkers.

Of course, it’s possible that good highly-active compounds could be obtained from either linking chemistry. However, the Torsion Library results clearly indicate that a thioether linkage is preferred here, purely because it orients the added fragments towards the hinge better. This is valuable knowledge if we were planning a small combinatorial library on this fragment expansion.

Adding the Torsion Library to Spark makes the results even more robust and useful, allowing you to see at a glance what the known experimental conformational preferences of small molecules say about the conformer quality in your Spark results. The new feature is available now as part of the Spark 10.4 releasetry a free evaluation today.

1 Shärfer  et al., J. Med. Chem. 2013, 56(5), 2016

2 Guba et al., JCIM 2016, 56, 1

What rings do medicinal chemists use, and why?

The vast majority of small molecule drugs contain at least one ring. The rigidity, synthetic accessibility and geometric preferences of rings mean that medicinal chemistry series are usually defined in terms of which ring or rings they have at their core. However, ring systems are more than just scaffolds waiting to be elaborated: the electrostatic and pharmacophoric properties of ring systems are usually crucial to the biological activity of the molecules that contain them.

We have conducted an investigation into the most common ring system and substitution patterns in the recent medicinal chemistry literature, as derived from the ChEMBL database. For each of these rings, the electrostatic potential has been calculated allowing the chemist to see at a glance the electronic properties of each system. In addition, applying the Spark bioisostere evaluation metric to the rings database reveals the best bioisosteric replacements for each ring system.

In the poster, What rings do medicinal chemists use, and why?, selected entries from the rings database are shown and discussed. The full data set is an invaluable aid to the medicinal chemist looking to understand the properties of their lead molecule and the opportunities for variation of its core.

Examining the diversity of large collections of building blocks in 3D as presented at Sheffield Chemoinformatics Conference 2016

Abstract

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.

Presentation

See presentation ‘Examining the diversity of large collections of building blocks in 3D‘ as presented at Sheffield Chemoinformatics Conference 2016.

User Group Meetings 2016 review

The Cresset 2016 North American and European User Group Meetings combined presentations and workshops in Cambridge, MA and UK respectively. A large and diverse group of scientists gathered to share ideas and experiences in the latest computational approaches to the design and optimization of small molecules.

As well as talks from invited speakers from the pharmaceutical and agrochemical industries, plus academia, Cresset staff also gave updates on new developments, including a preview of the new Cresset protein tool, due for alpha release this summer. You can see the presentations we have permission to publish below.

Workshops ran in parallel streams, focusing on increasing and deepening understanding of Cresset software. You can see the content covered in the workshops below.

Feedback from delegates on the scientific program and workshops was excellent, and I particularly enjoyed the chance to catch up with familiar faces and meet new people. These were our largest User Group Meetings to-date and we are already looking forward to celebrating 15 years of Cresset in 2017 – I hope many of you will join us.

Photos from User Group Meetings

Pros and cons of alignment-independent descriptors

Working with molecules in 3D is computationally expensive compared to most 2D methods. Most modern cheminformatics toolkits can do hundreds of thousands to millions of 2D fingerprint comparisons per second, with 3D similarity techniques being multiple orders of magnitude slower.

The computationally-expensive part of the 3D calculations usually involves aligning conformations to each other. The natural tendency therefore, is to see if we can skip this step and compute a set of properties that can tell us if two molecules are similar in 3D without actually having to align them. If this works we can get the best of both worlds: the speed of 2D comparisons combined with the accuracy and structural independence of 3D similarity functions.

Pharmacophoric descriptors

The earliest version of this idea is the simple pharmacophore. All you have to do is assign a few pharmacophoric points to each molecule (usually based on some sort of functional group pattern recognition), then generate a set of descriptors based on sets of these (usually 2 or 3). If two molecules share one or more pharmacophore descriptor, then they match.

Pharmacophore searches succeed on some counts: they are indeed very fast, and they do encode some 3D information. However, they involve a very crude binning of the wide array of possible intermolecular interactions into a few pharmacophore types, and they describe shape poorly giving them very limited predictive power.

If pharmacophores can’t describe shape well, are there other techniques that can? A number of different methods have been presented, such as those based on multipole moments or spherical harmonic coefficients (e.g. ParaSurf/ParaFit), as well as methods based on statistical moments such as Ultrafast Shape Recognition (USR). None of these has achieved widespread use: harmonic coefficients are not rotation-invariant, while the USR technique correlates poorly to more accurate measures of shape similarity. 1

Field descriptor distances

It would be nice if there was a way of providing alignment-independent descriptors which described both electrostatics and shape/pharmacophoric properties with a reasonable degree of accuracy. This is actually one of the first things we did when we were looking into starting Cresset – we developed a method called FieldPrint that encodes the distance matrix of field descriptors down into a fingerprint that can be used for alignment-independent similarity calculations. The concept is similar to that of GRIND 2 which was published around the same time, although the algorithmic details are somewhat different.

We put a lot of work into these techniques, but were never able to get a method that we were completely satisfied with. The problem we found is that encoding the distances in pairs/triplets of field descriptors ends up losing too much 3D information, and as a result you either end up with a slower mimic of standard 2D fingerprints, or you end up with a large false positive count. The FieldPrints have a tendency to find molecules with a similar overall pattern of positive/negative field, but can compute a very high similarity for molecules that are in reality quite dissimilar in terms of the 3D spatial arrangement of those fields. My belief now is that this is an inherent flaw of alignment-independent descriptors: they either have to be sufficiently complex that you are in effect computing an alignment, or you lose too much information and are not significantly better off than just using old-fashioned structural/pharmacophoric fingerprints.

As you move from full 3D interaction potentials to 2D correlograms to 1D fingerprints comparisons get faster but you lose information
Figure 1: As you move from full 3D interaction potentials to 2D correlograms to 1D fingerprints, comparisons get faster but you lose information

Handling conformation

A further consideration is how you handle conformation. The original GRIND papers just use a single conformer per molecule, and their validation was confined to series of rigid molecules or sets of molecules where single conformations were generated and manually adjusted to be similar. In the general case neither of these shortcuts will work. Any method that purports to be 3D but starts with a single conformation per molecule is inherently flawed: the whole point of 3D is that molecules are flexible.

There is a disturbing number of papers out there that do some sort of notionally-3D analysis on set of single CORINA-derived conformations. You can get very good enrichment factors on retrospective virtual screens doing this, but in practise the enrichments are largely bogus. CORINA is deterministic, and as a result molecules with similar structures will tend to be put into similar conformations. Combine this with the fact that many standard retrospective VS data sets have very low structural diversity, and the problem becomes apparent. The query molecule and its dozens of congeners in the “actives” data set are all placed in the same single conformation, and so application of a 3D or pseudo-3D technique can easily produce excellent-looking enrichment statistics. However, the enrichment all comes from a hidden 2D similarity.

So, single-conformation methods are a dead end and we need to consider flexibility. Once you are doing so, you need to factor in both the conformer generation time as part of the build time for the descriptor, and also factor in that your comparison speeds will now be two to four orders of magnitude slower than 2D fingerprints (assuming 100 conformers per molecule, and depending on whether you know a single bioactive conformation for one of the two molecules being compared or whether you need to compare conformer populations). 2D methods thus have an unassailable speed advantage, which is part of the reason they remain so popular.

Using FieldPrint as a filter

Our original vision for Blaze (or FieldScreen, as it was then) was that it would rely on the FieldPrints to give extremely rapid searching. You can get quite good enrichment factors from the FieldPrints in retrospective virtual screens, but when we investigated further this is largely because they act as a proxy for overall molecular size and charge. Once you control for that by more careful selection of decoys the FieldPrint performance is much less good. Analysing a molecular similarity technique through retrospective virtual screening performance is very very hard to do well, and as a result I am intrinsically wary of methods that present a set of DUD enrichments as their sole validation: FieldPrints perform quite well on DUD, but we know that they are not particularly effective in real prospective applications.

We still use the FieldPrint technology: it’s the first search stage in every Blaze run. It’s generally good enough to filter out 25-50% of decoy molecules that have no similarity to the query, but certainly not good enough to use the FieldPrint ranks directly. This is why we just use them as a pre-filter: molecules that pass that filter have much more accurate similarities computed using our alignment-based clique/simplex algorithms.

In the end, there’s no real short cut. All attempts to date to make 3D comparisons faster by simplifying descriptors and skipping the expensive alignment step just seem to leave out too much information – such techniques can be useful for cutting down the search space but if you’re going to spend CPU cycles working in 3D you might as well do it properly!

1. T. Zhou et al. / Journal of Molecular Graphics and Modelling 29 (2010) 443–449
2. Pastor, M.; Cruciani, G.; McLay, I.; Pickett, S.; Clementi, S. J. Med. Chem. 2000, 43 (17), 3233–3243.

Download an evaluation

Try our software for yourself – download a free evaluation.

Understanding SAR of circadian rhythm modulators

Anupriya KumarAnupriya Kumar, Assistant Professor, Nagoya University, summarizes the publication Development of Small-Molecule Cryptochrome Stabilizer Derivatives as Modulators of the Circadian Clock which cites how Cresset software was used to understand SAR of circadian rhythm modulators.
 
 

Introduction

Mammalian biological clocks control almost every physiological process such as metabolism, hormone secretion, body temperature and the sleep/wake cycle. Each cell in our body follows a circadian rhythm of 24 hours. Our modern lifestyle can disrupt this rhythm, for example working night shifts, or flying long distance. This in turn can lead to sleep and metabolic disorders such as obesity and type 2 diabetes. The daily oscillations of the circadian rhythm are generated by a loop which consists of clock genes such as CLOCK and BMAL1 (transcriptional activators) and two repressors: Cryptochrome (CRY1 and CRY2) and Period (PER1 and PER2).

Motivation to study the circadian rhythm modulators

The circadian rhythm can be modulated by small molecules, for example carbazole derivatives. One such molecule, KL001 was discovered by screening 60000 structurally diverse compounds in a cell-based Bmal1-dLuc reporter assay 1. KL001 lengthens the circadian rhythm by nearly 6 hours by stabilizing the clock protein cryptochrome (CRY). In the co-crystal structure of KL001 with the CRY2 protein (PDB ID. 4MLP 2), KL001 displaces the naturally binding FAD ligand (flavin adenine dinucleotide) while occupying only a small part of the binding pocket. Subsequent derivatization of KL001 resulted in a ten times more potent compound, KL044.

Building 3D-QSAR model using Forge

To understand the factors responsible for binding and stabilizing the CRY protein, we used Forge to perform 3D modeling and multivariate statistical analysis on 60 derivatives of KL001 and KL044. For generating a pharmacophore based on KL001 and KL044, we used FieldTemplater, which could reproduce the binding pose of KL001 to the CRY protein as shown in the crystal structure. We established a robust 3D-QSAR model by aligning all 60 compounds to the pharmacophore. The molecules were randomly divided into a training set of 54 and a test set of 6 compounds. Using this model, we successfully predicted activities of 5 new compounds and explained the evolution of the SAR by analyzing the contributing field points. By overlaying the field points from the 3D-QSAR model with the protein, we deciphered the important ligand-protein interactions.

Activity Cliff analysis

Activity Miner was used to study differences in the activity of the compounds by analyzing the electrostatic potential. The differences in the electrostatic potential maps (red = positive potential, cyan = negative potential) of the two matched pairs are shown below. For example, introducing a disubstituted phenyl ring and an amide group (a) enhances the activity as does the presence of a mesyl group in the same position (b).

Activity Cliff analysis

Conclusion

The current study serves as an encouragement for interdisciplinary teams to collaborate by showing that theoretical hypothesis and experimental studies complement each other. Knowledge based functional modifications can be made to KL001 and KL044 based on the activity enhancing key field points obtained from the 3D-QSAR model. Based on these findings, virtual screenings of compound databases can be performed to find novel compounds which would help in forwarding the field of drug discovery for circadian rhythm modulators which could lead to potential drugs for diabetes and other metabolic disorders.

References

1. Identification of small molecule activators of cryptochrome. Hirota T, Lee JW, St John PC, Sawa M, Iwaisako K, Noguchi T, Pongsawakul PY, Sonntag T, Welsh DK, Brenner DA, Doyle FJ 3rd, Schultz PG, Kay SA. Science. 2012 Aug 31;337(6098):1094-7.
2. Nangle S, Xing W, Zheng N. Cell Res. 2013 Dec;23(12):1417-9.

Affordable virtual screening with Blaze: Benchmarks

Introduction

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.

Conclusion

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.