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

Learning Protein-Ligand Unbinding Pathways via Single-Parameter Community Detection

Understanding the dynamics of biomolecular complexes, e.g., of protein-ligand (un)binding, requires the comprehension of paths such systems take between metastable states. In MD simulations, paths are usually not observable per se, but they need to be inferred from simulation trajectories. Here, we present a novel approach to cluster trajectories based on a community detection algorithm that necessitates only the definition of a single parameter. The unbinding of the streptavidin-biotin complex is used as a benchmark system and the A2a adenosine receptor in complex with the inhibitor ZM241385 as an elaborate application. We demonstrate how such clusters of trajectories correspond to pathways and how the approach helps in the identification of reaction coordinates for a considered (un)binding process.

Increased throughput in methods for simulating protein ligand binding and unbinding

By incorporating full flexibility and enabling the quantification of crucial parameters such as binding free energies and residence times, methods for investigating protein-ligand binding and unbinding via molecular dynamics provide details on the involved mechanisms at the molecular level. While these advancements hold promise for impacting drug discovery, a notable drawback persists: their relatively time-consuming nature limits throughput. Herein, we survey recent implementations which, employing a blend of enhanced sampling techniques, a clever choice of collective variables, and often machine learning, strive to enhance the efficiency of new and previously reported methods without compromising accuracy. Particularly noteworthy is the validation of these methods that was often performed on systems mirroring real-world drug discovery scenarios.

How exascale computing can shape drug design: A perspective from multiscale QM/MM molecular dynamics simulations and machine learning-aided enhanced sampling algorithms

Molecular simulations are an essential asset in the first steps of drug design campaigns. However, the requirement of high-throughput limits applications mainly to qualitative approaches with low computational cost, but also low accuracy. Unlocking the potential of more rigorous quantum mechanical/molecular mechanics (QM/MM) models combined with molecular dynamics-based free energy techniques could have a tremendous impact. Indeed, these two relatively old techniques are emerging as promising methods in the field. This has been favored by the exponential growth of computer power and the proliferation of powerful data-driven methods. Here, we briefly review recent advances and applications, and give our perspective on the impact that QM/MM and free-energy methods combined with machine learning-aided algorithms can have on drug design.

Surface-Accelerated String Method for Locating Minimum Free Energy Paths

We present a surface-accelerated string method (SASM) to efficiently optimize low-dimensional reaction pathways from the sampling performed with expensive quantum mechanical/molecular mechanical (QM/MM) Hamiltonians. The SASM accelerates the convergence of the path using the aggregate sampling obtained from the current and previous string iterations, whereas approaches like the string method in collective variables (SMCV) or the modified string method in collective variables (MSMCV) update the path only from the sampling obtained from the current iteration. Furthermore, the SASM decouples the number of images used to perform sampling from the number of synthetic images used to represent the path. The path is optimized on the current best estimate of the free energy surface obtained from all available sampling, and the proposed set of new simulations is not restricted to being located along the optimized path. Instead, the umbrella potential placement is chosen to extend the range of the free energy surface and improve the quality of the free energy estimates near the path. In this manner, the SASM is shown to improve the exploration for a minimum free energy pathway in regions where the free energy surface is relatively flat. Furthermore, it improves the quality of the free energy profile when the string is discretized with too few images. We compare the SASM, SMCV, and MSMCV using 3 QM/MM applications: a ribozyme methyltransferase reaction using 2 reaction coordinates, the 2'-O-transphosphorylation reaction of Hammerhead ribozyme using 3 reaction coordinates, and a tautomeric reaction in B-DNA using 5 reaction coordinates. We show that SASM converges the paths using roughly 3 times less sampling than the SMCV and MSMCV methods. All three algorithms have been implemented in the FE-ToolKit package made freely available.

Thermodynamically Optimized Machine-Learned Reaction Coordinates for Hydrophobic Ligand Dissociation

Ligand unbinding is mediated by its free energy change, which has intertwined contributions from both energy and entropy. It is important, but not easy, to quantify their individual contributions to the free energy profile. We model hydrophobic ligand unbinding for two systems, a methane particle and a C60 fullerene, both unbinding from hydrophobic pockets in all-atom water. Using a modified deep learning framework, we learn a thermodynamically optimized reaction coordinate to describe the hydrophobic ligand dissociation for both systems. Interpretation of these reaction coordinates reveals the roles of entropic and enthalpic forces as the ligand and pocket sizes change. In both cases, we observe that the free-energy barrier to unbinding is dominated by entropy considerations. Furthermore, the process of methane unbinding is driven by methane solvation, while fullerene unbinding is driven first by pocket wetting and then fullerene wetting. For both solutes, the direct importance of the distance from the binding pocket to the learned reaction coordinate is present, but low. Our framework and subsequent feature important analysis thus give useful thermodynamic insight into hydrophobic ligand dissociation problems that are otherwise difficult to glean.

Finite temperature string by K-means clustering sampling with order parameters as collective variables for molecular crystals: application to polymorphic transformation between β-CL-20 and ε-CL-20

Polymorphic transformation of molecular crystals is a fundamental phase transition process, and it is important practically in the chemical, material, biopharmaceutical, and energy storage industries. However, understanding of the transformation mechanism at the molecular level is poor due to the extreme simulating challenges in enhanced sampling and formulating order parameters (OPs) as the collective variables that can distinguish polymorphs with quite similar and complicated structures so as to describe the reaction coordinate. In this work, two kinds of OPs for CL-20 were constructed by the bond distances, bond orientations and relative orientations. A K-means clustering algorithm based on the Euclidean distance and sample weight was used to smooth the initial finite temperature string (FTS), and the minimum free energy path connecting β-CL-20 and ε-CL-20 was sketched by the string method in collective variables, and the free energy profile along the path and the nucleation kinetics were obtained by Markovian milestoning with Voronoi tessellations. In comparison with the average-based sampling, the K-means clustering algorithm provided an improved convergence rate of FTS. The simulation of transformation was independent of OP types but was affected greatly by finite-size effects. A surface-mediated local nucleation mechanism was confirmed and the configuration located at the shoulder of potential of mean force, rather than overall maximum, was confirmed to be the critical nucleus formed by the cooperative effect of the intermolecular interactions. This work provides an effective way to explore the polymorphic transformation of caged molecular crystals at the molecular level.

Standard Binding Free-Energy Calculations: How Far Are We from Automation?

Recent success stories suggest that in silico protein-ligand binding free-energy calculations are approaching chemical accuracy. However, their widespread application remains limited by the extensive human intervention required, posing challenges for the neophyte. As such, it is critical to develop automated workflows for estimating protein-ligand binding affinities with minimum personal involvement. Key human efforts include setting up and tuning enhanced-sampling or alchemical-transformation algorithms as a preamble to computational binding free-energy estimations. Additionally, preparing input files, bookkeeping, and postprocessing represent nontrivial tasks. In this Perspective, we discuss recent progress in automating standard binding free-energy calculations, featuring the development of adaptive or parameter-free algorithms, standardization of binding free-energy calculation workflows, and the implementation of user-friendly software. We also assess the current state of automated standard binding free-energy calculations and evaluate the limitations of existing methods. Last, we outline the requirements for future algorithms and workflows to facilitate automated free-energy calculations for diverse protein-ligand complexes.

How Robust Is the Ligand Binding Transition State?

For many drug targets, it has been shown that the kinetics of drug binding (e.g., on rate and off rate) is more predictive of drug efficacy than thermodynamic quantities alone. This motivates the development of predictive computational models that can be used to optimize compounds on the basis of their kinetics. The structural details underpinning these computational models are found not only in the bound state but also in the short-lived ligand binding transition states. Although transition states cannot be directly observed experimentally due to their extremely short lifetimes, recent successes have demonstrated that modeling the ligand binding transition state is possible with the help of enhanced sampling molecular dynamics methods. Previously, we generated unbinding paths for an inhibitor of soluble epoxide hydrolase (sEH) with a residence time of 11 min. Here, we computationally modeled unbinding events with the weighted ensemble method REVO (resampling of ensembles by variation optimization) for five additional inhibitors of sEH with residence times ranging from 14.25 to 31.75 min, with average prediction accuracy within an order of magnitude. The unbinding ensembles are analyzed in detail, focusing on features of the ligand binding transition state ensembles (TSEs). We find that ligands with similar bound poses can show significant differences in their ligand binding TSEs, in terms of their spatial distribution and protein-ligand interactions. However, we also find similarities across the TSEs when examining more general features such as ligand degrees of freedom. Together these findings show significant challenges for rational, kinetics-based drug design.

Perspectives on Ligand/Protein Binding Kinetics Simulations: Force Fields, Machine Learning, Sampling, and User-Friendliness

Computational techniques applied to drug discovery have gained considerable popularity for their ability to filter potentially active drugs from inactive ones, reducing the time scale and costs of preclinical investigations. The main focus of these studies has historically been the search for compounds endowed with high affinity for a specific molecular target to ensure the formation of stable and long-lasting complexes. Recent evidence has also correlated the in vivo drug efficacy with its binding kinetics, thus opening new fascinating scenarios for ligand/protein binding kinetic simulations in drug discovery. The present article examines the state of the art in the field, providing a brief summary of the most popular and advanced ligand/protein binding kinetics techniques and evaluating their current limitations and the potential solutions to reach more accurate kinetic models. Particular emphasis is put on the need for a paradigm change in the present methodologies toward ligand and protein parametrization, the force field problem, characterization of the transition states, the sampling issue, and algorithms' performance, user-friendliness, and data openness.

Investigating the Unbinding of Muscarinic Antagonists from the Muscarinic 3 Receptor

Patient symptom relief is often heavily influenced by the residence time of the inhibitor-target complex. For the human muscarinic receptor 3 (hMR3), tiotropium is a long-acting bronchodilator used in conditions such as asthma or chronic obstructive pulmonary disease (COPD). The mechanistic insights into this inhibitor remain unclear; specifically, the elucidation of the main factors determining the unbinding rates could help develop the next generation of antimuscarinic agents. Using our novel unbinding algorithm, we were able to investigate ligand dissociation from hMR3. The unbinding paths of tiotropium and two of its analogues, N-methylscopolamin and homatropine methylbromide, show a consistent qualitative mechanism and allow us to identify the structural bottleneck of the process. Furthermore, our machine learning-based analysis identified key roles of the ECL2/TM5 junction involved in the transition state. Additionally, our results point to relevant changes at the intracellular end of the TM6 helix leading to the ICL3 kinase domain, highlighting the closest residue L482. This residue is located right between two main protein binding sites involved in signal transduction for hMR3's activation and regulation. We also highlight key pharmacophores of tiotropium that play determining roles in the unbinding kinetics and could aid toward drug design and lead optimization.

Dissociation Rate Calculation via Constant-Force Steered Molecular Dynamics Simulation

Steered molecular dynamics (SMD) simulations are used to study molecular dissociation events by applying a harmonic force to the molecules and pulling them at a constant velocity. Instead of constant-velocity pulling, we use a constant force: the constant-force SMD (CF-SMD) simulation. The CF-SMD simulation employs a constant force to reduce the activation barrier of molecular dissociation, thereby enhancing the dissociation event. Here, we present the capability of the CF-SMD simulation to estimate the dissociation time at equilibrium. We performed all-atom CF-SMD simulations for NaCl and protein-ligand systems, producing dissociation time at various forces. We extrapolated these values to the dissociation rate without a constant force using Bell's model or the Dudko-Hummer-Szabo model. We demonstrate that the CF-SMD simulations with the models predicted the dissociation time in equilibrium. A CF-SMD simulation is a powerful tool for estimating the dissociation rate in a direct and computationally efficient manner.

Predicting Protein-Ligand Binding and Unbinding Kinetics with Biased MD Simulations and Coarse-Graining of Dynamics: Current State and Challenges

The prediction of drug-target binding and unbinding kinetics that occur on time scales between milliseconds and several hours is a prime challenge for biased molecular dynamics simulation approaches. This Perspective gives a concise summary of the theory and the current state-of-the-art of such predictions via biased simulations, of insights into the molecular mechanisms defining binding and unbinding kinetics as well as of the extraordinary challenges predictions of ligand kinetics pose in comparison to binding free energy predictions.

Path separation of dissipation-corrected targeted molecular dynamics simulations of protein-ligand unbinding

Protein-ligand (un)binding simulations are a recent focus of biased molecular dynamics simulations. Such binding and unbinding can occur via different pathways in and out of a binding site. Here, we present a theoretical framework on how to compute kinetics along separate paths and on how to combine the path-specific rates into global binding and unbinding rates for comparison with experimental results. Using dissipation-corrected targeted molecular dynamics in combination with temperature-boosted Langevin equation simulations [S. Wolf et al., Nat. Commun. 11, 2918 (2020)] applied to a two-dimensional model and the trypsin-benzamidine complex as test systems, we assess the robustness of the procedure and discuss the aspects of its practical applicability to predict multisecond kinetics of complex biomolecular systems.

Critical comparison of general-purpose collective variables for crystal nucleation

The nucleation of crystals is a prominent phenomenon in science and technology that still lacks a full atomic-scale understanding. Much work has been devoted to identifying order parameters able to track the process, from the inception of early nuclei to their maturing to critical size until growth of an extended crystal. We critically assess and compare two powerful distance-based collective variables, an effective entropy derived from liquid state theory and the path variable based on permutation invariant vectors using the Kob-Andersen binary mixture and a combination of enhanced-sampling techniques. Our findings reveal a comparable ability to drive nucleation when a bias potential is applied, and comparable free-energy barriers and structural features. Yet, we also found an imperfect correlation with the committor probability on the barrier top which was bypassed by changing the order parameter definition.

Simulation of ligand dissociation kinetics from the protein kinase PYK2

Early-stage drug discovery projects often focus on equilibrium binding affinity to the target alongside selectivity and other pharmaceutical properties. The kinetics of drug binding are ignored but can have significant influence on drug efficacy. Therefore, increasing attention has been paid on evaluating drug-binding kinetics early in a drug discovery process. Simulating drug-binding kinetics at the atomic level is challenging for the long time scale involved. Here, we used the transition-based reweighting analysis method (TRAM) with the Markov state model to study the dissociation of a ligand from the protein kinase PYK2. TRAM combines biased and unbiased simulations to reduce computational costs. This work used the umbrella sampling technique for the biased simulations. Although using the potential of mean force from umbrella sampling simulations with the transition-state theory over-estimated the dissociation rate by three orders of magnitude, TRAM gave a dissociation rate within an order of magnitude of the experimental value.