Combined Free-Energy Calculation and Machine Learning Methods for Understanding Ligand Unbinding Kinetics

Abstract

The determination of drug residence times, which define the time an inhibitor is in complex with its target, is a fundamental part of the drug discovery process. Synthesis and experimental measurements of kinetic rate constants are, however, expensive and time consuming. In this work, we aimed to obtain drug residence times computationally. Furthermore, we propose a novel algorithm to identify molecular design objectives based on ligand unbinding kinetics. We designed an enhanced sampling technique to accurately predict the free-energy profiles of the ligand unbinding process, focusing on the free-energy barrier for unbinding. Our method first identifies unbinding paths determining a corresponding set of internal coordinates (ICs) that form contacts between the protein and the ligand; it then iteratively updates these interactions during a series of biased molecular dynamics (MD) simulations to reveal the ICs that are important for the whole of the unbinding process. Subsequently, we performed finite-temperature string simulations to obtain the free-energy barrier for unbinding using the set of ICs as a complex reaction coordinate. Importantly, we also aimed to enable the further design of drugs focusing on improved residence times. To this end, we developed a supervised machine learning (ML) approach with inputs from unbiased “downhill” trajectories initiated near the transition state (TS) ensemble of the string unbinding path. We demonstrate that our ML method can identify key ligand–protein interactions driving the system through the TS. Some of the most important drugs for cancer treatment are kinase inhibitors. One of these kinase targets is cyclin-dependent kinase 2 (CDK2), an appealing target for anticancer drug development. Here, we tested our method using two different CDK2 inhibitors for the potential further development of these compounds. We compared the free-energy barriers obtained from our calculations with those observed in available experimental data. We highlighted important interactions at the distal ends of the ligands that can be targeted for improved residence times. Our method provides a new tool to determine unbinding rates and to identify key structural features of the inhibitors that can be used as starting points for novel design strategies in drug discovery.

Simulation of ligand dissociation kinetics from the protein kinase PYK2

“…As Badaoui et al 70 used umbrella sampling and string methods to study the dissociation of several ligands from another protein kinase, cyclin dependent kinase 2 (CDK2), it is useful to compare their dissociation mechanisms to ours. In their method, an initial dissociation pathway is constructed by performing biased simulations using a biasing potential based on a periodically updated list of interatomic distances.…”

Investigating the Unbinding of Muscarinic Antagonists from the Muscarinic 3 Receptor

“…Based on the sum of these interacting distances, a collective variable (CV) is defined and restrained harmonically. 29 During an iterative process, subsequent simulations of 10 ns use this biasing CV with a force constant of 10 kcal mol -1 A -2 . The constraint position (i.e., the length) of the CV is monotonically increased.…”

“…To aid the identification of the main CVs driving the system across the TS and to pinpoint novel descriptors that determine the fate of a binding/unbinding events, we used our MLTSA 29 . In this approach, we train an ML model to predict the outcome of downhill simulations with data close to the TS.…”

“…28 Given the homology between hMR2 and hMR3 similar limitations are expected to arise, which have been considered for this study. In this work, we applied our recently developed unbinding algorithm 29 to hMR3, to investigate the dissociation of tiotropium (1) and two structurally similar ligands, Nmethylscopolamin (2) and homatropine methylbromide (3) (Figure 1). The obtained unbinding pathways were refined using an adaptation of the finite temperature string method.…”

A critical comparison of general-purpose collective variables for crystal nucleation

“…This result points to the relevance of developing algorithms combin- ing CV-optimization and sampling acceleration in order to obtain accurate barriers. [25,26,62] To shed light on the issue related to the non-peaked distribution of CPA, we inspected the size of the largest crystalline cluster, N crys , by computing the value of the Steinhardt's bond-orientational order parameter averaged over the first neighbor shell, and defined ordered atoms as having q 6 larger than 0.25 [14,15]. In order to identify the shortcomings in the employed CVs, we focused on structures that were selected in Fig.…”