High-throughput protein-ligand binding free energy calculations with automated molecular dynamics thermodynamic integration protocol

Accurate estimation of protein-ligand binding free energy (ABFE) is of primary importance for computer-aided drug design. Molecular dynamics (MD) based alchemical free energy simulation methods such as thermodynamic integration (TI) allow for rigorous ABFE calculation. Recently, with the development of GPU-computing, ABFE simulations are increasingly used in drug discovery pipelines. Here, we describe an automated TI protocol for ABFE calculations utilizing GPU-accelerated MD simulations and its applications to two promising drug targets.
The first system is WDR domain of leucine-rich repeat kinase 2 (LRRK2), a protein considered a promising target for Parkinson disease treatment. Critical assessment of computational hit-finding experiments (CACHE) challenge aimed to find ligands targeting the central cavity of the WDR domain. For that, we performed ABFE calculations for approximately 1000 diverse small-molecule compounds from Enamine REAL and Mcule databases and identified 43 molecules with nanomolar predicted binding affinity.
The second system is AMPA-subtype ionotropic glutamate receptor (AMPAR) which plays a key role in epileptogenesis and represents a promising target for antiepileptic drugs. One of the small-molecule non-competitive AMPAR antagonists is GYKI 53655 (GYKI), which binds to an allosteric site of the AMPAR located at the interface of trans-membrane domain (TMD) and linkers connecting it to the ligand binding domains (LBD). Our previous docking and MD simulation studies demonstrated stability of the GYKI crystal and flipped orientation in the binding pocket as well as flexibility of TMD-LBD linkers. To account for protein flexibility, we calculated AMPAR-GYKI ensemble-averaged ABFE for both ligand orientations by averaging ABFEs computed for structures selected equidistantly from 500 ns AMPAR-GYKI MD trajectories. In total, we calculated ABFE for 800 individual AMPAR-GYKI structures which amounted to approximately 38 µs of MD sampling.