In silico screening of LRRK2 WDR domain inhibitors using deep docking and free energy simulations

Abstract

The Critical Assessment of Computational Hit-Finding Experiments (CACHE) Challenge series is focused on identifying small molecule inhibitors of protein targets using computational methods. Each challenge contains two phases, hit-finding and follow-up optimization, each of which is followed by experimental validation of the computational predictions. For the CACHE Challenge #1, the Leucine-Rich Repeat Kinase 2 (LRRK2) WD40 Repeat (WDR) domain was selected as the target for in silico hit-finding and optimization. Mutations in LRRK2 are the most common genetic cause of the familial form of Parkinson's disease. The LRRK2 WDR domain is an understudied drug target with no known molecular inhibitors. Herein we detail the first phase of our winning submission to the CACHE Challenge #1. We developed a framework for the high-throughput structure-based virtual screening of a chemically diverse small molecule space. Hit identification was performed using the large-scale Deep Docking (DD) protocol followed by absolute binding free energy (ABFE) simulations. ABFEs were computed using an automated molecular dynamics (MD)-based thermodynamic integration (TI) approach. 4.1 billion ligands from Enamine REAL were screened with DD followed by ABFEs computed by MD TI for 793 ligands. 76 ligands were prioritized for experimental validation, with 59 compounds successfully synthesized and 5 compounds identified as hits, yielding a 8.5% hit rate. Our results demonstrate the efficacy of the combined DD and ABFE approaches for hit identification for a target with no previously known hits. This approach is widely applicable for the efficient screening of ultra-large chemical libraries as well as rigorous protein-ligand binding affinity estimation leveraging modern computational resources.

Fig. 1. Thermodynamic cycle for alchemical ABFE simulations. The ligand (L) in its interacting and non-interacting form is shown by a gray and white hexagon correspondingly, and the protein (P) is shown by green. ΔG⁰_bind corresponds to the ABFE. ΔG^wat_int and ΔG^prot_int correspond to the free energy of turning off interactions of the ligand with the environment in solvent and in complex with the protein, respectively. ΔG^prot_+VB and ΔG^prot_−VB correspond to the addition and removal of the virtual bond restraints.

Fig. 2. Overview of hit identification pipeline. (A) Broad overview of the computational pipeline. Each number within the funnel relates to the number of molecules after the corresponding step has been taken. (B) In-depth breakdown of the computational pipeline: (1) Deep Docking virtual screening of the initial 4.1B library leading to the focused library for conventional docking, (2) conventional docking via multiple docking programs, (3) filtering by consensus of 2 and 3 docking programs, (4) in-parallel to step 3—Expert selection of prioritized molecules, (5) poverty & diversity filtering and docking-score based prioritization, (6) assessment of absolute binding free energy via molecular dynamics simulations. See Methods section for detailed descriptions of each step.

Fig. 3. (A) Protein–ligand interaction maps for hit compounds. Protein residues are shown as circles. Protein–ligand interactions are shown as lines. The residues and interactions are colored depending on the interaction type according to the legend. (B) Simulated representative structure (ligand shown in opaque yellow) of five hit compounds superimposed on the initial docked structure using protein C_α atoms. The initial position of the protein is not shown. The initial ligand position is shown in transparent yellow. The simulated protein is shown by a transparent cartoon representation. The protein residues that interact with the ligand in the simulated structure are shown in green stick representation. Protein–ligand hydrogen bonds are shown as dotted lines.

Fig. 4. (A) The binding site of WDR domain of LRRK2 predicted by MOE Site finder. (B) Putatively important hydrophobic residues for ligand binding, in complex with one of the virtual hits.

Fig. 5. (A) Distribution of ligand RMSD of the MD representative poses with respect to the docked pose. (B–D) Structure of ligands with various ligand RMSDs in the binding pocket of the WDR domain. The docked pose of the ligand is shown by transparent blue. The MD representative pose of the ligand is shown by yellow. The protein is presented by the gray cartoon. The ligand RMSDs are displayed in the lower left corner of each panel.

Fig. 6. (A) Distribution of ABFE calculated with MD TI for all molecules for which free energy simulations were performed. The Expert set was selected by visual inspection of docked ligands. The Consensus set was selected by consensus docking score. There were no overlaps between the Expert and Consensus sets. (B) Scatter plot of computed ABFE and consensus docking scores for molecules from the Consensus set. (C) Scatter plot of computed ABFE and ligand RMSD of the MD representative poses with respect to the docked pose.