Our proposed hit-identification workflow extends pipelines developed by our team during  CACHE Challenge 2. The core methodology consists of high-throughput docking followed by binding affinity estimation using Molecular Mechanics Poisson-Boltzman Surface Area (MMPBSA) on multiple poses drawn from a molecular dynamics (MD) run of the protein-ligand complex. Given the computational cost of running MD and MMPBSA, we developed two novel strategies to ensure that we are able to explore most of the molecular space of Enamine REAL while only doing the expensive evaluations on the most promising molecules. The first approach uses pharmacophore matching on fragments followed by fragment combination. The second approach trains a Bayesian graph neural network to predict the binding affinity scores derived from MMPBSA.

 

Pharmacophore-based fragment matching and combining: we leverage a fragment combination approach to search a vast combinatorial space of reagent combinations while reducing the need for computationally expensive methods. We construct several pharmacophore models based on experimentally observed binders at the NSP3 target site. Each pharmacophore is divided into non-overlapping sub-components. Fragment-sized compounds are matched against each pharmacophore subcomponent. Hits from complementary pharmacophore sets are fed into a reaction predictor which determines whether the fragments can be combined, based on both chemical and geometrical constraints. If the combined compound maintains compatibility with the full pharmacophore this compound is passed for scoring. 

 

We use Enamine (building blocks) fragment libraries, as these are combined by Enamine to derive the full Enamine REAL database. This way we only yield compounds that are available from Enamine. Furthermore, this approach enables us to effectively explore the full 5.5bn compound space while only having to run pharmacophore matching on ~300,000 fragment compounds. This approach explores a combinatorial reaction product space that grows quadratically (in the case of two-component reactions) or cubically (in the case of three-component reactions) while computational costs are linear with the number of fragments. Using this approach in CACHE Challenge 2 we were able to generate leads that shared chemical properties with experimentally observed fragment hits and we found that more than 90% of the compounds generated using this approach were available from Enamine. 

 

Fragment combinations that pass the pharmacophore model are further scored by consensus docking [Autodock, rDock, DOCK]. For the most promising pool of candidates, we compute the absolute free energy of binding using MMPBSA [GROMACS and gmx-MMPBSA], as well as QM/MMPBSA using semi-empirical quantum mechanical calculations. Free energy perturbation calculations will be employed to choose between top-ranking ligands with similar structures and to evaluate minor modifications to the final structures [GROMACS PyAutoFEP]. 

 

Bayesian Graph Neural Network (BGNN) binding affinity predictor: we create a training dataset by running MMPBSA on a few hundred molecules to generate binding affinity estimates. We then train a BGNN using the Chemprop library to predict binding affinity. This model is used to generate predictions on the Enamine REAL database and those molecules which have the highest upper-confidence bound are selected to pass through the docking + MMPBSA pipeline. After obtaining more MMPBSA scores the model can be re-trained and the cycle repeated. In CACHE 2 we were able to generate predictions for 1bn molecules per day and we found good concordance between the MMPBSA predicted score and the BGNN score.

 

A key aspect of our strategy is ensuring efficient sampling of the 5.5bn compound space of the Enamine REAL database, while also ensuring that the most computationally intensive methods (MMPBSA & FEP) are reserved for the most promising candidate molecules. Using cluster compute resources we are able to dock several million molecules: this includes carrying out ensemble docking against multiple conformations of the receptor.

 

The ~3 million curated set of molecules that we dock and score will be composed of the following:

 

• Pharmacophore matches generated from combining Enamine building blocks
• Molecules from Enamine REAL with high score from the BGNN
• ~2.1m molecules from the Enamine high-throughput screening library
• ~460k molecules from the Enamine hit locator library
• ~200k molecules from the Enamine building blocks library

Using this combination we will probe most of the relevant molecular space of the 5.5bn compound Enamine REAL database while reducing docking computations by multiple orders of magnitude.

 

Ligands are initially screened for molecular weight, solubility and exclusion of carboxylic acids. Final selections are screened against PAINS [BadApples webserver] to exclude ligands likely to be toxic or promiscuous binders. To ensure a diversity of structures in the final submission, we will cluster compounds based on structure and submit one sample from each cluster. FEP will be used to decide between structurally similar molecules. 

 

Our Reaction informed fragment fuser (RiFF) uses pharmacophore fragment matching followed by fragment combination. This allows us to reduce pharmacophore embedding and searching from billions of molecules to a few hundred thousand and yields compounds that share structural similarities with the experimentally observed binders. Second, we optimise the use of computationally expensive MD-based scoring by iteratively training a BGNN classifier to predict these scores and guide which molecules are evaluated next using a bayesian optimization strategy. This allows us to quickly screen Enamine REAL. Each of these methods complements our virtual screening funnel which uses docking followed by binding affinity estimation using MMPBSA.

Reaction informed fragment fuser (RiFF)

None

AutoDock Vina

AutoDock4

rDock

GROMACS

DOCK v6 & v3 (Kuntz Group UCSF)

PyAutoFEP

 

Badaoui, Magd, Pedro J. Buigues, Dénes Berta, Gaurav M. Mandana, Hankang Gu, Tamás Földes, Callum J. Dickson, et al. ‘Combined Free-Energy Calculation and Machine Learning Methods for Understanding Ligand Unbinding Kinetics’. Journal of Chemical Theory and Computation 18, no. 4 (12 April 2022): 2543–55. 

https://doi.org/10.1021/acs.jctc.1c00924

 

Berta, Dénes, Magd Badaoui, Sam Alexander Martino, Pedro J. Buigues, Andrei V. Pisliakov, Nadia Elghobashi-Meinhardt, Geoff Wells, Sarah A. Harris, Elisa Frezza, and Edina Rosta. ‘Modelling the Active SARS-CoV-2 Helicase Complex as a Basis for Structure-Based Inhibitor Design’. Chemical Science 12, no. 40 (2021): 13492–505. https://doi.org/10.1039/D1SC02775A

 

Bradshaw, John, Brooks Paige, Matt J. Kusner, Marwin H. S. Segler, and José Miguel Hernández-Lobato. ‘Barking up the Right Tree: An Approach to Search over Molecule Synthesis DAGs’. arXiv, 21 December 2020. 

http://arxiv.org/abs/2012.11522

 

Cook, Nicola J, Wen Li, Dénes Berta, Magd Badaoui, Allison Ballandras-Colas, Andrea Nans, Abhay Kotecha, Edina Rosta, Alan N Engelman, Peter Cherepanov. ‘Structural basis of second-generation HIV integrase inhibitor action and viral resistance’. Science 367 (6479), 806-810.

https://www.science.org/doi/abs/10.1126/science.aay4919  

 

Karlova, Andrea, Wim Dehaen, and Andrei Penciu. ‘How to Reward Your Drug Agent?’, NeurIPS Workshop 2021. 

 

Our approach to predicting actives from the pooled candidate molecules will use the physics-based scoring techniques described in the Hit Optimization stage. Our approach generates 3d coordinates of the ligands [OpenBabel], uses docking to get an initial binding location [AutoDock Vina], followed by MD-based methods: MMPBSA [GROMACS, gmx-MMPBSA] and FEP [GROMACS, PyAutoFEP] to estimate absolute and relative binding affinity.

Not given

See section 1

See section 1

Our optimisation strategy consists of two components. Firstly, searching for molecular optimisations to the experimental hits via a modify and test cycle. Proposed modifications to functional groups and atoms will be scored using free energy perturbation and accepted with a probability in proportion to the change in binding free energy (subject to a simulated annealing regime). Secondly, Experimental hits (and to some extent negatives) will be used to update the pharmacophore and re-train the QSAR model. Updating these methods at the top of the Hit Identification pipeline will yield new candidates that will be scored in the same way as described in sections 1 & 2. Additional candidates to put through the scoring pipeline will be generated from a structure similarity search of the experimental hits against the Enamine database. Our BGNN binding affinity predictor allows us to rapidly score millions of structural analogs of the experimental hits. This allows us to efficiently select which compounds should be evaluated with more computationally expensive methods such as MM-PBSA and FEP. Rapid prediction with the BGNN means that we can enumerate a larger space of molecules that encompass multi-step modifications to the lead molecules which would not be possible if we relied solely on MD-based binding affinity estimations.

Our team includes researchers experienced in drug discovery from multiple disciplines including physics, pharmacy, structural molecular biology and computer science with additional industry experience in medicinal chemistry. Our approach reflects this diverse range of experience, combining expertise in state-of-the-art machine learning architectures, pharmacophore design, and advanced free energy calculations using classical and QM-based MD simulations.

Not given

See section 1

See section 1

See section 1