Estimated reading time: 8 minutes
Introduction
Performing virtual screens can generate large numbers of hits. Herein, we discuss a workflow to triage the results from a virtual screening experiment, which culminated in using absolute free energy perturbation (ABFE) calculations to assess their binding affinity. This process, and those analogous to it, aim to identify interesting and active alternative chemical hits for a project.
This example targets the non-receptor tyrosine kinase 2 (TYK2) which is a member of the Janus kinase (JAK) family. TYK2 is associated with cytokine and growth receptors making it a target for inflammatory diseases. In this work, a protein crystal structure of TYK2 was downloaded and prepared (PDB code: 4GIH)1 using Flare™, and the ligand extracted as the query template for the Blaze™ screen, which used Cresset’s technology2 to find novel molecules with similar electrostatic and shape features to the query molecule. Initially, the top 7500 results from the virtual screen were prioritized for further triaging.
Filtering by physiochemical properties
We first used the PAINS3 filters to identify and remove compounds with substructures that potentially lead to non-specific assay readouts. The remaining compounds were then triaged using a multi-parameter optimization function (MPO) focused on desired molecular features and physicochemical properties: this was achieved by using the radial plot feature in Flare which incorporated molecular weight, SLogP, TPSA, H-bond donors, H-bond acceptors, ring count, and rotatable bonds as part of the score. Each of the properties used in the MPO was given a specific range to prioritize lead-like compounds. The MPO thus scored molecules from 0 to 1, with our criteria for passing a compound being set to an MPO score of 0.7 or higher.
Since virtual screening algorithms can often produce unreasonable torsions, we applied a filter to remove molecules with too many torsions rarely seen in the Cambridge Structural Database (CSD). Specifically, we allowed up to one low-frequency torsion (i.e. torsions rarely seen in the CSD database) and up to three medium-frequency torsions. With the scoring based on physiochemical properties complete, we created a tag to filter out molecules that failed.
Selecting by Electrostatic Complementarity
As an orthogonal scoring metric, we built another radial plot (MPO) based on the docking score, Electrostatic Complementarity™ (EC) to the protein, and the score from Blaze of the top 7500 compounds. The docking score was computed using the fixed-pose option of the ‘score only’ method, using the pose generated by Blaze. Once the radial plot score was computed, we selected molecules that scored above 0.65 for further triaging. We then repeated the initial triaging steps (PAINS passes, preferred physicochemical properties profiles, torsions filters, and interaction and ligand similarity scores) for the next 7500 compounds in the Blaze output. All passed compounds from both tranches were then merged into one dataset.
Do my molecules have potential solubility liabilities? Let’s predict before we test
Testing compounds with poor solubility can lead to unreliable results, and so the next stage of the triage was to use a QSAR model to predict which compounds might have solubility issues and to subsequently remove them from the data set. The model was built using the gradient boosting method to predict hydration free energy based on RDKit descriptors using the dataset from Mobley and Guthrie.4 We used the 513 molecules in the training set, 65 in the test set, and 64 in the prediction set to build our model in Flare. The QSAR model built provided confidence levels in the model prediction and any ligands with predicted poor solubility values with a high confidence level were removed from the dataset of virtual hits.
Finding molecules with unique protein-ligand interactions
The next filtering step used a feature in Flare to cluster molecules by their protein-ligand interaction fingerprint (PLIF). This allowed us to efficiently identify molecules having different interaction patterns with TYK2. We also clustered molecules by 2D similarity to assess the chemical diversity of our list of compounds. The representative from each cluster taken to the next stage was the molecule with the highest Radial Plot score.
Are the poses stable?
After visually inspecting the clusters of the triaged results, representatives were selected for molecular dynamics simulations to determine if the ligand poses were stable in the protein active site. We defined stability as maintaining an RMSD < 2.0 Å from the Blaze query ligand pose during a 20 ns simulation. Stable ligands were progressed to absolute free energy perturbation calculations (ABFE) to identify those with the best predicted binding energies.
ABFE results: which molecules to synthesize and modify?
Free energy perturbation (FEP) calculations allow you to accurately estimate the binding free energy between a ligand and protein of interest. ABFE calculations calculate the binding free energy independently of any other ligand in the dataset. This makes it an ideal method for our situation where we have identified examples of interesting ligands from different chemical series but want to subsequently identify those ligands that have the higher binding free energies.
In our study, the ABFE calculations indicated high predicted activity for the reference ligand, and the known active compound that we included in the virtual screen dataset. Further to these expected results, ABFE also identified a promising compound with a 4-methyl-2-pyrimidinamine group with high predicted binding energy. In a live project the interesting hits could be purchased or synthesized to provide new starting points for a project, and in this example, a compound containing a pyrimidine group was actually identified in the original work by Liang et al.1 This demonstrates how our virtual screening and triaging approach can find chemically different starting points for a project.
Conclusion
The key takeaway from this study is that we were able to start with one single ligand from a crystal structure and identify virtual hits with different classes of structure with high predicted activity. This compound could be a solid starting point for wet chemistry, whereas without doing the computational triaging work, medicinal chemists would have to potentially explore more chemistry avenues leading to increased expense through using time-consuming reactions and subsequent screening work. ABFE therefore has the potential to reduce costs and time taken to select, purchase/synthesize and test active compounds from a virtual screen.
Why Blaze, why Flare?
Using Flare in conjunction with the Blaze virtual screening platform can provide useful results to jumpstart a medicinal chemistry project. Blaze allows you to start from a single active ligand and through searching commercially available compounds, provide a diverse set of virtual hits. Flare provides a varied suite of approaches to triage these results from a virtual screen and aiding in the identification of which compounds to focus on.
References
- Liang J. et al.; Lead identification of novel and selective TYK2 inhibitors. European Journal of Medicinal Chemistry. 2013, 67, 175-187
- Cheeseright T., Mackey M., Rose S., Vinter, A.; Molecular Field Extrema as Descriptors of Biological Activity: Definition and Validation. J. Chem. Inf. Model. 2006, 46 (2), 665-676
- Baell, J.B., Nissink, J.W.M. Seven Year Itch: Pan-Assay Interference Compounds (PAINS) in 2017 – Utility and Limitations. ACS Chemical Biology, 2018 13 (1), 36-44
- Mobley, D.L., Guthrie, J.P. FreeSolv: a database of experimental and calculated hydration free energies, with input files. J Comput. Aided Mol. Des. 2014 28(7) 711-720
