Protein Flexibility and Dissociation Pathway Differentiation Can Explain Onset of Resistance Mutations in Kinases

The importance of binding kinetics and drug-target residence time in pharmacology

A dominant assumption in pharmacology throughout the 20th century has been that in vivo target occupancy-and attendant pharmacodynamics-depends on the systemic concentration of drug relative to the equilibrium dissociation constant for the drug-target complex. In turn, the duration of pharmacodynamics is temporally linked to the systemic pharmacokinetics of the drug. Yet, there are many examples of drugs for which pharmacodynamic effect endures long after the systemic concentration of a drug has waned to (equilibrium) insignificant levels. To reconcile such data, the drug-target residence time model was formulated, positing that it is the lifetime (or residence time) of the binary drug-target complex, and not its equilibrium affinity per se, that determines the extent and duration of drug pharmacodynamics. Here, we review this model, its evolution over time, and its applications to natural ligand-macromolecule biology and synthetic drug-target pharmacology.

Machine learning-augmented molecular dynamics simulations (MD) reveal insights into the disconnect between affinity and activation of ZTP riboswitch ligands

The challenge of targeting RNA with small molecules necessitates a better understanding of RNA-ligand interaction mechanisms. However, the dynamic nature of nucleic acids, their ligand-induced stabilization, and how conformational changes influence gene expression pose significant difficulties for experimental investigation. This work employs a combination of computational and experimental methods to address these challenges. By integrating structure-informed design, crystallography, and machine learning-augmented all-atom molecular dynamics simulations (MD) we synthesized, biophysically and biochemically characterized, and studied the dissociation of a library of small molecule activators of the ZTP riboswitch, a ligand-binding RNA motif that regulates bacterial gene expression. We uncovered key interaction mechanisms, revealing valuable insights into the role of ligand binding kinetics on riboswitch activation. Further, we established that ligand on-rates determine activation potency as opposed to binding affinity and elucidated RNA structural differences, which provide mechanistic insights into the interplay of RNA structure on riboswitch activation.

Calculating Protein-Ligand Residence Times through State Predictive Information Bottleneck Based Enhanced Sampling

Understanding drug residence times in target proteins is key to improving drug efficacy and understanding target recognition in biochemistry. While drug residence time is just as important as binding affinity, atomic-level understanding of drug residence times through molecular dynamics (MD) simulations has been difficult primarily due to the extremely long time scales. Recent advances in rare event sampling have allowed us to reach these time scales, yet predicting protein-ligand residence times remains a significant challenge. Here we present a semi-automated protocol to calculate the ligand residence times across 12 orders of magnitude of time scales. In our proposed framework, we integrate a deep learning-based method, the state predictive information bottleneck (SPIB), to learn an approximate reaction coordinate (RC) and use it to guide the enhanced sampling method metadynamics. We demonstrate the performance of our algorithm by applying it to six different protein-ligand complexes with available benchmark residence times, including the dissociation of the widely studied anticancer drug Imatinib (Gleevec) from both wild-type Abl kinase and drug-resistant mutants. We show how our protocol can recover quantitatively accurate residence times, potentially opening avenues for deeper insights into drug development possibilities and ligand recognition mechanisms.

Enhanced Sampling with Machine Learning

Molecular dynamics (MD) enables the study of physical systems with excellent spatiotemporal resolution but suffers from severe timescale limitations. To address this, enhanced sampling methods have been developed to improve the exploration of configurational space. However, implementing these methods is challenging and requires domain expertise. In recent years, integration of machine learning (ML) techniques into different domains has shown promise, prompting their adoption in enhanced sampling as well. Although ML is often employed in various fields primarily due to its data-driven nature, its integration with enhanced sampling is more natural with many common underlying synergies. This review explores the merging of ML and enhanced MD by presenting different shared viewpoints. It offers a comprehensive overview of this rapidly evolving field, which can be difficult to stay updated on. We highlight successful strategies such as dimensionality reduction, reinforcement learning, and flow-based methods. Finally, we discuss open problems at the exciting ML-enhanced MD interface.

Calculating Protein-Ligand Residence Times Through State Predictive Information Bottleneck based Enhanced Sampling

Understanding drug residence times in target proteins is key to improving drug efficacy and understanding target recognition in biochemistry. While drug residence time is just as important as binding affinity, atomic-level understanding of drug residence times through molecular dynamics (MD) simulations has been difficult primarily due to the extremely long timescales. Recent advances in rare event sampling have allowed us to reach these timescales, yet predicting protein-ligand residence times remains a significant challenge. Here we present a semi-automated protocol to calculate the ligand residence times across 12 orders of magnitudes of timescales. In our proposed framework, we integrate a deep learning-based method, the state predictive information bottleneck (SPIB), to learn an approximate reaction coordinate (RC) and use it to guide the enhanced sampling method metadynamics. We demonstrate the performance of our algorithm by applying it to six different protein-ligand complexes with available benchmark residence times, including the dissociation of the widely studied anti-cancer drug Imatinib (Gleevec) from both wild-type Abl kinase and drug-resistant mutants. We show how our protocol can recover quantitatively accurate residence times, potentially opening avenues for deeper insights into drug development possibilities and ligand recognition mechanisms.

Exploring Kinase Asp-Phe-Gly (DFG) Loop Conformational Stability with AlphaFold2-RAVE

Kinases compose one of the largest fractions of the human proteome, and their misfunction is implicated in many diseases, in particular, cancers. The ubiquitousness and structural similarities of kinases make specific and effective drug design difficult. In particular, conformational variability due to the evolutionarily conserved Asp-Phe-Gly (DFG) motif adopting in and out conformations and the relative stabilities thereof are key in structure-based drug design for ATP competitive drugs. These relative conformational stabilities are extremely sensitive to small changes in sequence and provide an important problem for sampling method development. Since the invention of AlphaFold2, the world of structure-based drug design has noticeably changed. In spite of it being limited to crystal-like structure prediction, several methods have also leveraged its underlying architecture to improve dynamics and enhanced sampling of conformational ensembles, including AlphaFold2-RAVE. Here, we extend AlphaFold2-RAVE and apply it to a set of kinases: the wild type DDR1 sequence and three mutants with single point mutations that are known to behave drastically differently. We show that AlphaFold2-RAVE is able to efficiently recover the changes in relative stability using transferable learned order parameters and potentials, thereby supplementing AlphaFold2 as a tool for exploration of Boltzmann-weighted protein conformations (Meller, A.; Bhakat, S.; Solieva, S.; Bowman, G. R. Accelerating Cryptic Pocket Discovery Using AlphaFold. J. Chem. Theory Comput. 2023, 19, 4355-4363).

Data-driven classification of ligand unbinding pathways

Studying the pathways of ligand-receptor binding is essential to understand the mechanism of target recognition by small molecules. The binding free energy and kinetics of protein-ligand complexes can be computed using molecular dynamics (MD) simulations, often in quantitative agreement with experiments. However, only a qualitative picture of the ligand binding/unbinding paths can be obtained through a conventional analysis of the MD trajectories. Besides, the higher degree of manual effort involved in analyzing pathways limits its applicability in large-scale drug discovery. Here, we address this limitation by introducing an automated approach for analyzing molecular transition paths with a particular focus on protein-ligand dissociation. Our method is based on the dynamic time-warping algorithm, originally designed for speech recognition. We accurately classified molecular trajectories using a very generic descriptor set of contacts or distances. Our approach outperforms manual classification by distinguishing between parallel dissociation channels, within the pathways identified by visual inspection. Most notably, we could compute exit-path-specific ligand-dissociation kinetics. The unbinding timescale along the fastest path agrees with the experimental residence time, providing a physical interpretation to our entirely data-driven protocol. In combination with appropriate enhanced sampling algorithms, this technique can be used for the initial exploration of ligand-dissociation pathways as well as for calculating path-specific thermodynamic and kinetic properties.

Allo-targeting of the kinase domain: Insights from in silico studies and comparison with experiments

The eukaryotic protein kinase domain has been a broadly explored target for drug discovery, despite limitations imposed by its high sequence conservation as a shared modular domain and the development of resistance to drugs. One way of addressing those limitations has been to target its potential allosteric sites, shortly called allo-targeting, in conjunction with, or separately from, its conserved catalytic/orthosteric site that has been widely exploited. Allosteric regulation has gained importance as an alternative to overcome the drawbacks associated with the indiscriminate effect of targeting the active site, and it turned out to be particularly useful for these highly promiscuous and broadly shared kinase domains. Yet, allo-targeting often faces challenges as the allosteric sites are not as clearly defined as its orthosteric sites, and the effect on the protein function may not be unambiguously assessed. A robust understanding of the consequence of site-specific allo-targeting on the conformational dynamics of the target protein is essential to design effective allo-targeting strategies. Recent years have seen important advances in in silico identification of druggable sites and distinguishing among them those sites expected to allosterically mediate conformational switches essential to signal transmission. The present opinion underscores the utility of such computational approaches applied to the kinase domain, with the help of comparison between computational predictions and experimental observations.

Roles of Accelerated Molecular Dynamics Simulations in Predictions of Binding Kinetic Parameters

Rational predictions on binding kinetics parameters of drugs to targets play significant roles in future drug designs. Full conformational samplings of targets are requisite for accurate predictions of binding kinetic parameters. In this review, we mainly focus on the applications of enhanced sampling technologies in calculations of binding kinetics parameters and residence time of drugs. The methods involved in molecular dynamics simulations are applied to not only probe conformational changes of targets but also reveal calculations of residence time that is significant for drug efficiency. For this review, special attention are paid to accelerated molecular dynamics (aMD) and Gaussian aMD (GaMD) simulations that have been adopted to predict the association or disassociation rate constant. We also expect that this review can provide useful information for future drug design.

Activation Loop Plasticity and Active Site Coupling in the MAP Kinase, ERK2

Previous studies of the protein kinase, ERK2, using NMR and hydrogen-exchange measurements have shown changes in dynamics accompanying its activation by phosphorylation. However, knowledge about the conformational motions involved is incomplete. Here, we examined ERK2 using long conventional molecular dynamics (MD) simulations starting from crystal structures of phosphorylated (2P) and unphosphorylated (0P) forms. Individual trajectories were run for (5 to 25) μs, totaling 727 μs. The results show unexpected flexibility of the A-loop, with multiple long-lived (>5 μs) conformational states in both 2P- and 0P-ERK2. Differential contact network and principal component analyses reveal coupling between the A-loop fold and active site dynamics, with evidence for conformational selection in the kinase core of 2P-ERK2 but not 0P-ERK2. Simulations of 2P-ERK2 show A-loop states corresponding to restrained dynamics within the N-lobe, including regions around catalytic residues. One A-loop conformer forms lasting interactions with the L16 segment, leading to reduced RMSF and greater compaction in the active site. By contrast, simulations of 0P-ERK2 reveal excursions of A-loop residues away from the C-lobe, leading to greater active site mobility. Thus, the A-loop in ERK2 switches between distinct conformations that reflect coupling with the active site, possibly via the L16 segment. Crystal packing interactions suggest that lattice contacts with the A-loop may restrain its structural variation in X-ray structures of ERK2. The novel conformational states identified by MD expand our understanding of ERK2 regulation, by linking the activated state of the kinase to reduced dynamics and greater compaction surrounding the catalytic site.

Conformational Switch of the 250s Loop Enables the Efficient Transglycosylation in GH Family 77

Amylomaltases (AMs) play important roles in glycogen and maltose metabolism. However, the molecular mechanism is elusive. Here, we investigated the conformational dynamics of the 250s loop and catalytic mechanism of Thermus aquaticus TaAM using path-metadynamics and QM/MM MD simulations. The results demonstrate that the transition of the 250s loop from an open to closed conformation promotes polysaccharide sliding, leading to the ideal positioning of the acid/base. Furthermore, the conformational dynamics can also modulate the selectivity of hydrolysis and transglycosylation. The closed conformation of the 250s loop enables the tight packing of the active site for transglycosylation, reducing the energy penalty and efficiently preventing the penetration of water into the active site. Conversely, the partially closed conformation for hydrolysis results in a loosely packed active site, destabilizing the transition state. These computational findings guide mutation experiments and enable the identification of mutants with an improved disproportionation/hydrolysis ratio. The present mechanism is in line with experimental data, highlighting the critical role of conformational dynamics in regulating the catalytic reactivity of GHs.

Kinetics from Metadynamics: Principles, Applications, and Outlook

Metadynamics is a popular enhanced sampling algorithm for computing the free energy landscape of rare events by using molecular dynamics simulation. Ten years ago, Tiwary and Parrinello introduced the infrequent metadynamics approach for calculating the kinetics of transitions across free energy barriers. Since then, metadynamics-based methods for obtaining rate constants have attracted significant attention in computational molecular science. Such methods have been applied to study a wide range of problems, including protein-ligand binding, protein folding, conformational transitions, chemical reactions, catalysis, and nucleation. Here, we review the principles of elucidating kinetics from metadynamics-like approaches, subsequent methodological developments in this area, and successful applications on chemical, biological, and material systems. We also highlight the challenges of reconstructing accurate kinetics from enhanced sampling simulations and the scope of future developments.

Exploring kinase DFG loop conformational stability with AlphaFold2-RAVE

Kinases compose one of the largest fractions of the human proteome, and their misfunction is implicated in many diseases, in particular cancers. The ubiquitousness and structural similarities of kinases makes specific and effective drug design difficult. In particular, conformational variability due to the evolutionarily conserved DFG motif adopting in and out conformations and the relative stabilities thereof are key in structure-based drug design for ATP competitive drugs. These relative conformational stabilities are extremely sensitive to small changes in sequence, and provide an important problem for sampling method development. Since the invention of AlphaFold2, the world of structure-based drug design has noticably changed. In spite of it being limited to crystal-like structure prediction, several methods have also leveraged its underlying architecture to improve dynamics and enhanced sampling of conformational ensembles, including AlphaFold2-RAVE. Here, we extend AlphaFold2-RAVE and apply it to a set of kinases: the wild type DDR1 sequence and three mutants with single point mutations that are known to behave drastically differently. We show that AlphaFold2-RAVE is able to efficiently recover the changes in relative stability using transferable learnt order parameters and potentials, thereby supplementing AlphaFold2 as a tool for exploration of Boltzmann-weighted protein conformations.

Temperature Accelerated Sliced Sampling to Probe Ligand Dissociation from Protein

Modeling ligand unbinding in proteins to estimate the free energy of binding and probing the mechanism presents several challenges. They primarily pertain to the entropic bottlenecks resulting from protein and solvent conformations. While exploring the unbinding processes using enhanced sampling techniques, very long simulations are required to sample all of the conformational states as the system gets trapped in local free energy minima along transverse coordinates. Here, we demonstrate that temperature accelerated sliced sampling (TASS) is an ideal approach to overcome some of the difficulties faced by conventional sampling methods in studying ligand unbinding. Using TASS, we study the unbinding of avibactam inhibitor molecules from the Class C β-lactamase (CBL) active site. Extracting CBL-avibactam unbinding free energetics, unbinding pathways, and identifying critical interactions from the TASS simulations are demonstrated.

Uncovering the Kinetic Characteristics and Degradation Preference of PROTAC Systems with Advanced Theoretical Analyses

Proteolysis-targeting chimeras (PROTACs), which can selectively induce the degradation of target proteins, represent an attractive technology in drug discovery. A large number of PROTACs have been reported, but due to the complicated structural and kinetic characteristics of the target-PROTAC-E3 ligase ternary interaction process, the rational design of PROTACs is still quite challenging. Here, we characterized and analyzed the kinetic mechanism of MZ1, a PROTAC that targets the bromodomain (BD) of the bromodomain and extra terminal (BET) protein (Brd2, Brd3, or Brd4) and von Hippel-Lindau E3 ligase (VHL), from the kinetic and thermodynamic perspectives of view by using enhanced sampling simulations and free energy calculations. The simulations yielded satisfactory predictions on the relative residence time and standard binding free energy (r_p > 0.9) for MZ1 in different Brd^BD-MZ1-VHL ternary complexes. Interestingly, the simulation of the PROTAC ternary complex disintegration illustrates that MZ1 tends to remain on the surface of VHL with the BD proteins dissociating alone without a specific dissociation direction, indicating that the PROTAC prefers more to bind with E3 ligase at the first step in the formation of the target-PROTAC-E3 ligase ternary complex. Further exploration of the degradation difference of MZ1 in different Brd systems shows that the PROTAC with higher degradation efficiency tends to leave more lysine exposed on the target protein, which is guaranteed by the stability (binding affinity) and durability (residence time) of the target-PROTAC-E3 ligase ternary complex. It is quite possible that the underlying binding characteristics of the Brd^BD-MZ1-VHL systems revealed by this study may be shared by different PROTAC systems as a general rule, which may accelerate rational PROTAC design with higher degradation efficiency.

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.

Selectivity and Ranking of Tight-Binding JAK-STAT Inhibitors Using Markovian Milestoning with Voronoi Tessellations

Janus kinases (JAK), a group of proteins in the nonreceptor tyrosine kinase (NRTKs) family, play a crucial role in growth, survival, and angiogenesis. They are activated by cytokines through the Janus kinase-signal transducer and activator of a transcription (JAK-STAT) signaling pathway. JAK-STAT signaling pathways have significant roles in the regulation of cell division, apoptosis, and immunity. Identification of the V617F mutation in the Janus homology 2 (JH2) domain of JAK2 leading to myeloproliferative disorders has stimulated great interest in the drug discovery community to develop JAK2-specific inhibitors. However, such inhibitors should be selective toward JAK2 over other JAKs and display an extended residence time. Recently, novel JAK2/STAT5 axis inhibitors (N-(1H-pyrazol-3-yl)pyrimidin-2-amino derivatives) have displayed extended residence times (hours or longer) on target and adequate selectivity excluding JAK3. To facilitate a deeper understanding of the kinase-inhibitor interactions and advance the development of such inhibitors, we utilize a multiscale Markovian milestoning with Voronoi tessellations (MMVT) approach within the Simulation-Enabled Estimation of Kinetic Rates v.2 (SEEKR2) program to rank order these inhibitors based on their kinetic properties and further explain the selectivity of JAK2 inhibitors over JAK3. Our approach investigates the kinetic and thermodynamic properties of JAK-inhibitor complexes in a user-friendly, fast, efficient, and accurate manner compared to other brute force and hybrid-enhanced sampling approaches.

Probing conformational landscapes and mechanisms of allosteric communication in the functional states of the ABL kinase domain using multiscale simulations and network-based mutational profiling of allosteric residue potentials

In the current study, multiscale simulation approaches and dynamic network methods are employed to examine the dynamic and energetic details of conformational landscapes and allosteric interactions in the ABL kinase domain that determine the kinase functions. Using a plethora of synergistic computational approaches, we elucidate how conformational transitions between the active and inactive ABL states can employ allosteric regulatory switches to modulate intramolecular communication networks between the ATP site, the substrate binding region, and the allosteric binding pocket. A perturbation-based network approach that implements mutational profiling of allosteric residue propensities and communications in the ABL states is proposed. Consistent with biophysical experiments, the results reveal functionally significant shifts of the allosteric interaction networks in which preferential communication paths between the ATP binding site and substrate regions in the active ABL state become suppressed in the closed inactive ABL form, which in turn features favorable allosteric coupling between the ATP site and the allosteric binding pocket. By integrating the results of atomistic simulations with dimensionality reduction methods and Markov state models, we analyze the mechanistic role of macrostates and characterize kinetic transitions between the ABL conformational states. Using network-based mutational scanning of allosteric residue propensities, this study provides a comprehensive computational analysis of long-range communications in the ABL kinase domain and identifies conserved regulatory hotspots that modulate kinase activity and allosteric crosstalk between the allosteric pocket, ATP binding site, and substrate binding regions.

Advances in computational methods for ligand binding kinetics

Binding kinetic parameters can be correlated with drug efficacy, which in recent years led to the development of various computational methods for predicting binding kinetic rates and gaining insight into protein-drug binding paths and mechanisms. In this review, we introduce and compare computational methods recently developed and applied to two systems, trypsin-benzamidine and kinase-inhibitor complexes. Methods involving enhanced sampling in molecular dynamics simulations or machine learning can be used not only to predict kinetic rates, but also to reveal factors modulating the duration of residence times, selectivity, and drug resistance to mutations. Methods which require less computational time to make predictions are highlighted, and suggestions to reduce the error of computed kinetic rates are presented.

Ligand Unbinding Pathway and Mechanism Analysis Assisted by Machine Learning and Graph Methods

We present two methods to reveal protein-ligand unbinding mechanisms in biased unbinding simulations by clustering trajectories into ensembles representing unbinding paths. The first approach is based on a contact principal component analysis for reducing the dimensionality of the input data, followed by identification of unbinding paths and training a machine learning model for trajectory clustering. The second approach clusters trajectories according to their pairwise mean Euclidean distance employing the neighbor-net algorithm, which takes into account input data bias in the distances set and is superior to dendrogram construction. Finally, we describe a more complex case where the reaction coordinate relevant for path identification is a single intraligand hydrogen bond, highlighting the challenges involved in unbinding path reaction coordinate detection.

Transition Rates and Efficiency of Collective Variables from Time-Dependent Biased Simulations

Simulations with adaptive time-dependent bias enable an efficient exploration of the conformational space of a system. However, the dynamic information is altered by the bias. Infrequent metadynamics recovers the transition rate of crossing a barrier, if the collective variables are ideal and there is no bias deposition near the transition state. Unfortunately, these conditions are not always fulfilled. To overcome these limitations, and inspired by single-molecule force spectroscopy, we use Kramers' theory for calculating the barrier-crossing rate when a time-dependent bias is added to the system. We assess the efficiency of collective variables parameter by measuring how efficiently the bias accelerates the transitions. We present approximate analytical expressions of the survival probability, reproducing the barrier-crossing time statistics and enabling the extraction of the unbiased transition rate even for challenging cases. We explore the limits of our method and provide convergence criteria to assess its validity.