Predicting Residence Time and Drug Unbinding Pathway through Scaled Molecular Dynamics

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.

Ligand Gaussian accelerated Molecular Dynamics 3 (LiGaMD3): Improved Calculations of Binding Thermodynamics and Kinetics of Both Small Molecules and Flexible Peptides

Binding thermodynamics and kinetics play critical roles in drug design. However, it has proven challenging to efficiently predict ligand binding thermodynamics and kinetics of small molecules and flexible peptides using conventional Molecular Dynamics (cMD), due to limited simulation timescales. Based on our previously developed Ligand Gaussian accelerated Molecular Dynamics (LiGaMD) method, we present a new approach, termed "LiGaMD3", in which we introduce triple boosts into three individual energy terms that play important roles in small-molecule/peptide dissociation, rebinding and system conformational changes to improve the sampling efficiency of small-molecule/peptide interactions with target proteins. To validate the performance of LiGaMD3, MDM2 bound by a small molecule (Nutlin 3) and two highly flexible peptides (PMI and P53) were chosen as model systems. LiGaMD3 could efficiently capture repetitive small-molecule/peptide dissociation and binding events within 2 microsecond simulations. The predicted binding kinetic constant rates and free energies from LiGaMD3 agreed with available experimental values and previous simulation results. Therefore, LiGaMD3 provides a more general and efficient approach to capture dissociation and binding of both small-molecule ligand and flexible peptides, allowing for accurate prediction of their binding thermodynamics and kinetics.

It's all about making new contacts: How being metabotropic and phasicity help D1-like receptors promote LTP in the PFC

D1-like receptors have two important qualities, they are all metabotropic and they activate with phasic dopamine. After analyzing the molecular implications of each of these qualities separately and then combining them for the specific case of the prefrontal cortex, we propose a model that explains why long term potentiation in this cortical area depends on the amount of contact between D1-like receptors and dopamine. This simple model also explains why in order to promote long term potentiation, dopamine transporters should be scarce in the prefrontal cortex. Additionally, it explains why stimulants like methamphetamine could have such detrimental cognitive effects on regular substance consumers.

Simulation of the ligand-leaving process of the human heat shock protein

The human heat shock protein plays a critical role in various diseases and is an important target for pharmacological modulation. Simulation of conformational changes and free energy profiles of the human heat shock protein derived by the ligand-leaving process is a challenging issue. In this work, steered molecular dynamics simulation was adopted to simulate the ligand-leaving process. Two composite systems of heat shock protein NHSP90 and small molecules 6FJ and 6G7 are selected as research objects. The free energy during the leaving of ligand small molecules is calculated using conventional molecular dynamics simulation, steered molecular dynamics simulation (SMD), and the umbrella sampling method. We found that the a slower pulling velocity (0.001 nm ns-1) will result in 2.19 kcal mol-1, and the umbrella sampling method gives a value of 3.26 kcal mol-1 for the free energy difference for the two systems, which reasonably agrees with experimental results. A faster-pulling velocity (0.01 nm ns-1) leads to a large overestimation of free energy. At the same time, the conformational analysis indicated that the faster pulling velocity may lead to the conformational change of NHSP90, which was proved to be false by the slower pulling velocity and the umbrella sampling method.

Elucidating protein-ligand binding kinetics based on returning probability theory

The returning probability (RP) theory, a rigorous diffusion-influenced reaction theory, enables us to analyze the binding process systematically in terms of thermodynamics and kinetics using molecular dynamics (MD) simulations. Recently, the theory was extended to atomistically describe binding processes by adopting the host-guest interaction energy as the reaction coordinate. The binding rate constants can be estimated by computing the thermodynamic and kinetic properties of the reactive state existing in the binding processes. Here, we propose a methodology based on the RP theory in conjunction with the energy representation theory of solution, applicable to complex binding phenomena, such as protein-ligand binding. The derived scheme of calculating the equilibrium constant between the reactive and dissociate states, required in the RP theory, can be used for arbitrary types of reactive states. We apply the present method to the bindings of small fragment molecules [4-hydroxy-2-butanone (BUT) and methyl methylthiomethyl sulphoxide (DSS)] to FK506 binding protein (FKBP) in an aqueous solution. Estimated binding rate constants are consistent with those obtained from long-timescale MD simulations. Furthermore, by decomposing the rate constants to the thermodynamic and kinetic contributions, we clarify that the higher thermodynamic stability of the reactive state for DSS causes the faster binding kinetics compared with BUT.

15 Years of molecular simulation of drug-binding kinetics

INTRODUCTION

Drug-binding kinetics has been increasingly recognized as an important factor to be considered in drug discovery. Long residence time could prolong the action of some drugs while produce toxicity on others. Early evaluation of the binding kinetics of drug candidates could reduce attrition rate late in the drug discovery process. Computational prediction of drug-binding kinetics is useful as compounds can be evaluated even before they are made. However, simulation of drug-binding kinetics is a challenging problem because of the long-time scale involved. Nevertheless, significant progress has been made.

AREAS COVERED

This review illustrates the rapid evolution of qualitative to quantitative molecular dynamics-based methods that have been developed over the last 15 years.

EXPERT OPINION

The development of new methods based on molecular dynamics simulations now enables computation of absolute association/dissociation rate constants. Cheaper methods capable of identifying candidates with fast or slow binding kinetics, or rank-ordering rate constants are also available. Together, these methods have generated useful insights into the molecular mechanisms of drug-binding kinetics, and the design of drug candidates with therapeutically favorable kinetics. Although predicting absolute rate constants is still expensive and challenging, rapid improvement is expected in the coming years with the continuing refinement of current technologies, development of new methodologies, and the utilization of machine learning.

Structure-Kinetic Relationship for Drug Design Revealed by a PLS Model with Retrosynthesis-Based Pre-Trained Molecular Representation and Molecular Dynamics Simulation

Drug design based on kinetic properties is growing in application. Here, we applied retrosynthesis-based pre-trained molecular representation (RPM) in machine learning (ML) to train 501 inhibitors of 55 proteins and successfully predicted the dissociation rate constant (k_off) values of 38 inhibitors from an independent dataset for the N-terminal domain of heat shock protein 90α (N-HSP90). Our RPM molecular representation outperforms other pre-trained molecular representations such as GEM, MPG, and general molecular descriptors from RDKit. Furthermore, we optimized the accelerated molecular dynamics to calculate the relative retention time (RT) for the 128 inhibitors of N-HSP90 and obtained the protein-ligand interaction fingerprints (IFPs) on their dissociation pathways and their influencing weights on the k_off value. We observed a high correlation among the simulated, predicted, and experimental -log(k_off) values. Combining ML, molecular dynamics (MD) simulation, and IFPs derived from accelerated MD helps design a drug for specific kinetic properties and selectivity profiles to the target of interest. To further validate our k_off predictive ML model, we tested our model on two new N-HSP90 inhibitors, which have experimental k_off values and are not in our ML training dataset. The predicted k_off values are consistent with experimental data, and the mechanism of their kinetic properties can be explained by IFPs, which shed light on the nature of their selectivity against N-HSP90 protein. We believe that the ML model described here is transferable to predict k_off of other proteins and will enhance the kinetics-based drug design endeavor.

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.

Predicting Biomolecular Binding Kinetics: A Review

Biomolecular binding kinetics including the association (k_on) and dissociation (k_off) rates are critical parameters for therapeutic design of small-molecule drugs, peptides, and antibodies. Notably, the drug molecule residence time or dissociation rate has been shown to correlate with their efficacies better than binding affinities. A wide range of modeling approaches including quantitative structure-kinetic relationship models, Molecular Dynamics simulations, enhanced sampling, and Machine Learning has been developed to explore biomolecular binding and dissociation mechanisms and predict binding kinetic rates. Here, we review recent advances in computational modeling of biomolecular binding kinetics, with an outlook for future improvements.

Ligand Gaussian Accelerated Molecular Dynamics 2 (LiGaMD2): Improved Calculations of Ligand Binding Thermodynamics and Kinetics with Closed Protein Pocket

Ligand binding thermodynamics and kinetics are critical parameters for drug design. However, it has proven challenging to efficiently predict ligand binding thermodynamics and kinetics from molecular simulations due to limited simulation timescales. Protein dynamics, especially in the ligand binding pocket, often plays an important role in ligand binding. Based on our previously developed Ligand Gaussian accelerated molecular dynamics (LiGaMD), here we present LiGaMD2 in which a selective boost potential was applied to both the ligand and protein residues in the binding pocket to improve sampling of ligand binding and dissociation. To validate the performance of LiGaMD2, the T4 lysozyme (T4L) mutants with open and closed pockets bound by different ligands were chosen as model systems. LiGaMD2 could efficiently capture repetitive ligand dissociation and binding within microsecond simulations of all T4L systems. The obtained ligand binding kinetic rates and free energies agreed well with available experimental values and previous modeling results. Therefore, LiGaMD2 provides an improved approach to sample opening of closed protein pockets for ligand dissociation and binding, thereby allowing for efficient calculations of ligand binding thermodynamics and kinetics.

Complete reconstruction of dasatinib unbinding pathway from c-Src kinase by supervised molecular dynamics simulation method; assessing efficiency and trustworthiness of the method

Over the past years, rational drug design has gained lots of attention since employing it gave the world targeted therapy and more effective treatment solutions. Structure-based drug design (SBDD) is an excellent tool in rational drug design that takes advantage of accurate methods such as unbiased molecular dynamics (UMD) simulation for designing and optimizing molecular entities by understanding the binding and unbinding pathways of the binders. Supervised molecular dynamics (SuMD) simulation is a branch of UMD in which long-duration simulations are turned into short simulations, called replica, and a specific parameter is monitored throughout the simulation. In this work, we utilized this strategy to reconstruct the unbinding pathway of the anticancer drug dasatinib from its target protein, the c-Src kinase. Several unbinding events with valuable details were achieved. Then, to assess the efficiency and trustworthiness of the SuMD method, the unbinding pathway was also reconstructed by conventional UMD simulation, which uncovered some of the limitations of this method, such as limited sampling of the active site and finding the metastable states in the unbinding pathway. Furthermore, in times like these, when the world is desperate to find treatments for the Covid-19 disease, we think these methods are of exceptional value.Communicated by Ramaswamy H. Sarma.

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.

Unraveling the Abnormal Molecular Mechanism of Suicide Inhibition of Cytochrome P450 3A4

Suicide inhibition of the CYP3A4 enzyme by a drug inactivates the enzyme in the drug biotransformation process and often shows safety concerns about the drug. Despite extensive experimental studies, the abnormal molecular mechanism of a suicide inhibitor that forms a covalent bond with the residue far away from the catalytically active center of CYP3A4 inactivating the enzyme remains elusive. Here, the authors used molecular simulation approaches to study in detail how diquinone methide (DQR), the metabolite product of raloxifene, unbinds from CYP3A4 and inactivates the enzyme at the atomistic level. The results clearly indicate that in one of the intermediate states formed in its unbinding process, DQR covalently binds to Cys239, a residue far away from the catalytically active center of CYP3A4, and hinders the substrate from entering or leaving the enzyme. This work therefore provides an unprecedented way of clarifying the abnormal mechanism of suicide inhibition of the CYP3A4 enzyme.

Qualitative Estimation of Protein-Ligand Complex Stability through Thermal Titration Molecular Dynamics Simulations

The prediction of ligand efficacy has long been linked to thermodynamic properties such as the equilibrium dissociation constant, which considers both the association and the dissociation rates of a defined protein-ligand complex. In the last 15 years, there has been a paradigm shift, with an increased interest in the determination of kinetic properties such as the drug-target residence time since they better correlate with ligand efficacy compared to other parameters. In this article, we present thermal titration molecular dynamics (TTMD), an alternative computational method that combines a series of molecular dynamics simulations performed at progressively increasing temperatures with a scoring function based on protein-ligand interaction fingerprints for the qualitative estimation of protein-ligand-binding stability. The protocol has been applied to four different pharmaceutically relevant test cases, including protein kinase CK1δ, protein kinase CK2, pyruvate dehydrogenase kinase 2, and SARS-CoV-2 main protease, on a variety of ligands with different sizes, structures, and experimentally determined affinity values. In all four cases, TTMD was successfully able to distinguish between high-affinity compounds (low nanomolar range) and low-affinity ones (micromolar), proving to be a useful screening tool for the prioritization of compounds in a drug discovery campaign.

Insight into Thermodynamic and Kinetic Profiles in Small-Molecule Optimization

Structure-activity relationships (SARs) and structure-property relationships (SPRs) have been considered the most important factors during the drug optimization process. For medicinal chemists, improvements in the potencies and druglike properties of small molecules are regarded as their major goals. Among them, the binding affinity and selectivity of small molecules on their targets are the most important indicators. In recent years, there has been growing interest in using thermodynamic and kinetic profiles to analyze ligand-receptor interactions, which could provide not only binding affinities but also detailed binding parameters for small-molecule optimization. In this perspective, we are trying to provide an insight into thermodynamic and kinetic profiles in small-molecule optimization. Through a highlight of strategies on the small-molecule optimization with specific cases, we aim to put forward the importance of structure-thermodynamic relationships (STRs) and structure-kinetic relationships (SKRs), which could provide more guidance to find safe and effective small-molecule drugs.

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

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.

Decisive role of water and protein dynamics in residence time of p38α MAP kinase inhibitors

Target residence time plays a crucial role in the pharmacological activity of small molecule inhibitors. Little is known, however, about the underlying causes of inhibitor residence time at the molecular level, which complicates drug optimization processes. Here, we employ all-atom molecular dynamics simulations (~400 μs in total) to gain insight into the binding modes of two structurally similar p38α MAPK inhibitors (type I and type I½) with short and long residence times that otherwise show nearly identical inhibitory activities in the low nanomolar IC₅₀ range. Our results highlight the importance of protein conformational stability and solvent exposure, buried surface area of the ligand and binding site resolvation energy for residence time. These findings are further confirmed by simulations with a structurally diverse short residence time inhibitor SB203580. In summary, our data provide guidance in compound design when aiming for inhibitors with improved target residence time.

The molecular determinants of the residence time of a small molecule inhibitor at its target protein are not well understood. Here, Pantsar et al. show that the target protein’s conformational stability and solvent exposure are key factors governing the target residence time of kinase inhibitors.

Challenges and frontiers of computational modelling of biomolecular recognition

Biomolecular recognition including binding of small molecules, peptides and proteins to their target receptors plays a key role in cellular function and has been targeted for therapeutic drug design. However, the high flexibility of biomolecules and slow binding and dissociation processes have presented challenges for computational modelling. Here, we review the challenges and computational approaches developed to characterise biomolecular binding, including molecular docking, molecular dynamics simulations (especially enhanced sampling) and machine learning. Further improvements are still needed in order to accurately and efficiently characterise binding structures, mechanisms, thermodynamics and kinetics of biomolecules in the future.

Atomistic description of molecular binding processes based on returning probability theory

The efficiency of molecular binding such as host-guest binding is commonly evaluated in terms of kinetics, such as rate coefficients. In general, to compute the coefficient of the overall binding process, we need to consider both the diffusion of reactants and barrier crossing to reach the bound state. Here, we develop a methodology of quantifying the rate coefficient of binding based on molecular dynamics simulation and returning probability (RP) theory proposed by Kim and Lee [J. Chem. Phys. 131, 014503 (2009)]. RP theory provides a tractable formula of the rate coefficient in terms of the thermodynamic stability and kinetics of the intermediate state on a predefined reaction coordinate. In this study, the interaction energy between reactants is utilized as the reaction coordinate, enabling us to effectively describe the reactants' relative position and orientation on one-dimensional space. Application of this method to the host-guest binding systems, which consist of β-cyclodextrin and small guest molecules, yields the rate coefficients consistent with the experimental results.

Mechanistic Understanding from Molecular Dynamics in Pharmaceutical Research 2: Lipid Membrane in Drug Design

We review the use of molecular dynamics (MD) simulation as a drug design tool in the context of the role that the lipid membrane can play in drug action, i.e., the interaction between candidate drug molecules and lipid membranes. In the standard "lock and key" paradigm, only the interaction between the drug and a specific active site of a specific protein is considered; the environment in which the drug acts is, from a biophysical perspective, far more complex than this. The possible mechanisms though which a drug can be designed to tinker with physiological processes are significantly broader than merely fitting to a single active site of a single protein. In this paper, we focus on the role of the lipid membrane, arguably the most important element outside the proteins themselves, as a case study. We discuss work that has been carried out, using MD simulation, concerning the transfection of drugs through membranes that act as biological barriers in the path of the drugs, the behavior of drug molecules within membranes, how their collective behavior can affect the structure and properties of the membrane and, finally, the role lipid membranes, to which the vast majority of drug target proteins are associated, can play in mediating the interaction between drug and target protein. This review paper is the second in a two-part series covering MD simulation as a tool in pharmaceutical research; both are designed as pedagogical review papers aimed at both pharmaceutical scientists interested in exploring how the tool of MD simulation can be applied to their research and computational scientists interested in exploring the possibility of a pharmaceutical context for their research.

Contact Map Fingerprints of Protein-Ligand Unbinding Trajectories Reveal Mechanisms Determining Residence Times Computed from Scaled Molecular Dynamics

The binding kinetic properties of potential drugs may significantly influence their subsequent clinical efficacy. Predictions of these properties based on computer simulations provide a useful alternative to their expensive and time-consuming experimental counterparts, even at an early drug discovery stage. Herein, we perform scaled molecular dynamics (ScaledMD) simulations on a set of 27 ligands of HSP90 belonging to more than seven chemical series to estimate their relative residence times. We introduce two new techniques for the analysis and the classification of the simulated unbinding trajectories. The first technique, which helps in estimating the limits of the free energy well around the bound state, and the second one, based on a new contact map fingerprint, allow the description and the comparison of the paths that lead to unbinding. Using these analyses, we find that ScaledMD's relative residence time generally enables the identification of the slowest unbinders. We propose an explanation for the underestimation of the residence times of a subset of compounds, and we investigate how the biasing in ScaledMD can affect the mechanistic insights that can be gained from the simulations.

Comparative analysis of the unbinding pathways of antiviral drug Indinavir from HIV and HTLV1 proteases by supervised molecular dynamics simulation

Determining the unbinding pathways of potential small molecule compounds from their target proteins is of great significance for designing efficacious treatment solutions. One of these potential compounds is the approved HIV-1 protease inhibitor, Indinavir, which has a weak effect on the HTLV-1 protease. In this work, by employing the SuMD method, we reconstructed the unbinding pathways of Indinavir from HIV and HTLV-1 proteases to compare and understand the mechanism of the unbinding and to discover the reasons for the lack of inhibitory activity of Indinavir against the HTLV-1 protease. We achieved multiple unbinding events from both HIV and HTLV-1 proteases in which the RMSD values of Indinavir reached over 40 Å. Also, we found that the mobility and fluctuations of the flap region are higher in the HTLV-1 protease, making the drug less stable. We realized that critically positioned aromatic residues such as Trp98/Trp98' and Phe67/Phe67' in the HTLV-1 protease could make strong π-Stacking interactions with Indinavir in the unbinding pathway, which are unfavorable for the stability of Indinavir in the active site. The details found in this study can make a reasonable explanation for the lack of inhibitory activity of this drug against HTLV-1 protease. We believe the details discovered in this work can help design more effective and selective inhibitors for the HTLV-1 protease.

New perspectives in cancer drug development: computational advances with an eye to design

Computational chemistry has come of age in drug discovery. Indeed, most pharmaceutical development programs rely on computer-based data and results at some point. Herein, we discuss recent applications of advanced simulation techniques to difficult challenges in drug discovery. These entail the characterization of allosteric mechanisms and the identification of allosteric sites or cryptic pockets determined by protein motions, which are not immediately evident in the experimental structure of the target; the study of ligand binding mechanisms and their kinetic profiles; and the evaluation of drug-target affinities. We analyze different approaches to tackle challenging and emerging biological targets. Finally, we discuss the possible perspectives of future application of computation in drug discovery.

Kinetic and thermodynamic stability comparison for the fibrillar form of small amyloid-β(1-42) oligomers using scaled molecular dynamics

Amyloid-β (Aβ) oligomers act as intermediates for several neurodegenerative disease-relevant fibril formations. However, gaining insight into the oligomer to fibril conversion process remains a challenge due to the transient nature of small Aβ. In this study, we probe the kinetic and thermodynamic stabilities of small Aβ(1-42) oligomers in fibrillar conformations to understand from what size these aggregates start forming stable fibrils. With no definite structures available for small Aβ42 aggregates, we have started with oligomers extracted from mature fibrils having four, five, six and nine chains stacked together, and have performed order-to-disorder transition on these systems. Using scaled molecular dynamics (sMD) simulation, the timescale for breaking the native contacts of fibrils has been compared. The results indicate that the kinetic stability of oligomers increases with size, especially at the C-terminus end beyond five-chain oligomers. The free energy of breaking the contacts at the β-sheet regions in the structures has been obtained on an unscaled potential from a free energy extrapolation (FEE) approach. The values show that although stable minima are obtained for larger oligomers due to the enhanced stability of the C-terminus ends, fully stable fibril formation may require aggregates larger than the ones considered in our study. Additionally, dissimilar kinetics for the unbinding of terminal chains across all the oligomers has been observed. The interaction energy values calculated from unscaled MD simulations reveal the crucial role of water in our observations. Our work provides the application of an easy-to-deploy method that sheds light on interactions which could be significant in the early stages of Aβ42 fibril formation.

Reduced efficacy of a Src kinase inhibitor in crowded protein solution

The inside of a cell is highly crowded with proteins and other biomolecules. How proteins express their specific functions together with many off-target proteins in crowded cellular environments is largely unknown. Here, we investigate an inhibitor binding with c-Src kinase using atomistic molecular dynamics (MD) simulations in dilute as well as crowded protein solution. The populations of the inhibitor, 4-amino-5-(4-methylphenyl)-7-(t-butyl)pyrazolo[3,4-d]pyrimidine (PP1), in bulk solution and on the surface of c-Src kinase are reduced as the concentration of crowder bovine serum albumins (BSAs) increases. This observation is consistent with the reduced PP1 inhibitor efficacy in experimental c-Src kinase assays in addition with BSAs. The crowded environment changes the major binding pathway of PP1 toward c-Src kinase compared to that in dilute solution. This change is explained based on the population shift mechanism of local conformations near the inhibitor binding site in c-Src kinase.

Molecular dynamics simulations of ortho-carborane nano-diamond storage within the nonpolar channel cavities of a right-handed coiled-coil tetrabrachion nanotube

Molecular dynamics simulations have been performed on a complex in which clusters of boron in the form of molecules of the nanodiamond ortho-carborane ( C 2 B 10 H 12 ) have been inserted into the four large nonpolar cavities of a nanotube of the right-handed coiled-coil ( R H C C ) t e t r a b r a c h i o n . The techniques of multi-configurational thermodynamic integration, steered molecular dynamics and umbrella sampling have been combined to investigate the energetics of storage of ortho-carborane in the cavities and to map out the free energy landscape of the RHCC - t e t r a b r a c h i o n - o r t h o - c a r b o r a n e complex along the central channel and along directions transverse to the central channel. The purpose of the study was to explore potential pathways for the diffusion of ortho-carborane between the cavities and the solvent and to assess the stability of the complex as a possible drug delivery system for boron neutron capture therapy (BNCT). The investigation reveals a complex free energy landscape with a multitude of peaks and valleys, all of which can be related to specific architectural elements of the RHCC - n a n o t u b e , and the activation barriers for ortho-carborane capture and release support the requirements for rapid cargo uptake coupled with tight binding to the cavities.

Perspectives on High-Throughput Ligand/Protein Docking With Martini MD Simulations

Molecular docking is central to rational drug design. Current docking techniques suffer, however, from limitations in protein flexibility and solvation models and by the use of simplified scoring functions. All-atom molecular dynamics simulations, on the other hand, feature a realistic representation of protein flexibility and solvent, but require knowledge of the binding site. Recently we showed that coarse-grained molecular dynamics simulations, based on the most recent version of the Martini force field, can be used to predict protein/ligand binding sites and pathways, without requiring any a priori information, and offer a level of accuracy approaching all-atom simulations. Given the excellent computational efficiency of Martini, this opens the way to high-throughput drug screening based on dynamic docking pipelines. In this opinion article, we sketch the roadmap to achieve this goal.

Molecular Dynamics Simulations in Drug Discovery and Pharmaceutical Development

Molecular dynamics (MD) simulations have become increasingly useful in the modern drug development process. In this review, we give a broad overview of the current application possibilities of MD in drug discovery and pharmaceutical development. Starting from the target validation step of the drug development process, we give several examples of how MD studies can give important insights into the dynamics and function of identified drug targets such as sirtuins, RAS proteins, or intrinsically disordered proteins. The role of MD in antibody design is also reviewed. In the lead discovery and lead optimization phases, MD facilitates the evaluation of the binding energetics and kinetics of the ligand-receptor interactions, therefore guiding the choice of the best candidate molecules for further development. The importance of considering the biological lipid bilayer environment in the MD simulations of membrane proteins is also discussed, using G-protein coupled receptors and ion channels as well as the drug-metabolizing cytochrome P450 enzymes as relevant examples. Lastly, we discuss the emerging role of MD simulations in facilitating the pharmaceutical formulation development of drugs and candidate drugs. Specifically, we look at how MD can be used in studying the crystalline and amorphous solids, the stability of amorphous drug or drug-polymer formulations, and drug solubility. Moreover, since nanoparticle drug formulations are of great interest in the field of drug delivery research, different applications of nano-particle simulations are also briefly summarized using multiple recent studies as examples. In the future, the role of MD simulations in facilitating the drug development process is likely to grow substantially with the increasing computer power and advancements in the development of force fields and enhanced MD methodologies.

Thermodynamics and Kinetics of Drug-Target Binding by Molecular Simulation

Computational studies play an increasingly important role in chemistry and biophysics, mainly thanks to improvements in hardware and algorithms. In drug discovery and development, computational studies can reduce the costs and risks of bringing a new medicine to market. Computational simulations are mainly used to optimize promising new compounds by estimating their binding affinity to proteins. This is challenging due to the complexity of the simulated system. To assess the present and future value of simulation for drug discovery, we review key applications of advanced methods for sampling complex free-energy landscapes at near nonergodicity conditions and for estimating the rate coefficients of very slow processes of pharmacological interest. We outline the statistical mechanics and computational background behind this research, including methods such as steered molecular dynamics and metadynamics. We review recent applications to pharmacology and drug discovery and discuss possible guidelines for the practitioner. Recent trends in machine learning are also briefly discussed. Thanks to the rapid development of methods for characterizing and quantifying rare events, simulation's role in drug discovery is likely to expand, making it a valuable complement to experimental and clinical approaches.

Data-Driven Molecular Dynamics: A Multifaceted Challenge

The big data concept is currently revolutionizing several fields of science including drug discovery and development. While opening up new perspectives for better drug design and related strategies, big data analysis strongly challenges our current ability to manage and exploit an extraordinarily large and possibly diverse amount of information. The recent renewal of machine learning (ML)-based algorithms is key in providing the proper framework for addressing this issue. In this respect, the impact on the exploitation of molecular dynamics (MD) simulations, which have recently reached mainstream status in computational drug discovery, can be remarkable. Here, we review the recent progress in the use of ML methods coupled to biomolecular simulations with potentially relevant implications for drug design. Specifically, we show how different ML-based strategies can be applied to the outcome of MD simulations for gaining knowledge and enhancing sampling. Finally, we discuss how intrinsic limitations of MD in accurately modeling biomolecular systems can be alleviated by including information coming from experimental data.

Intuitive, reproducible high-throughput molecular dynamics in Galaxy: a tutorial

This paper is a tutorial developed for the data analysis platform Galaxy. The purpose of Galaxy is to make high-throughput computational data analysis, such as molecular dynamics, a structured, reproducible and transparent process. In this tutorial we focus on 3 questions: How are protein-ligand systems parameterized for molecular dynamics simulation? What kind of analysis can be carried out on molecular trajectories? How can high-throughput MD be used to study multiple ligands? After finishing you will have learned about force-fields and MD parameterization, how to conduct MD simulation and analysis for a protein-ligand system, and understand how different molecular interactions contribute to the binding affinity of ligands to the Hsp90 protein.

Revealing the Unbinding Kinetics and Mechanism of Type I and Type II Protein Kinase Inhibitors by Local-Scaled Molecular Dynamics Simulations

Protein kinase inhibitors disrupt phosphorylation of the target kinases, which are an important class of drug for treating cancer and other diseases. Conventional structure-based design methods (such as molecular docking) focus on the static binding mode of the kinase inhibitor with its target. However, dissociation kinetic properties of a drug molecule are found to correlate with its residence time in vivo and thus have drawn the attention of drug designers in recent years. In this study, we have applied the local-scaled molecular dynamics (MD) simulation enabled in GROMACS software to explore the unbinding mechanism of a total of 41 type I and type II kinase inhibitors. Our simulation considered multiple starting configurations as well as possible protonation states of kinase inhibitors. Based on our local-scaled MD results, we discovered that the integrals of the favorable binding energy during dissociation correlated well (R² = 0.64) with the experimental dissociation rate constants of those kinase inhibitors on the entire data set. Given its accuracy and technical advantage, this method may serve as a practical option for estimating this important property in reality. Our simulation also provided a reasonable explanation of the dynamic properties of kinase and its inhibitor as well as the role of relevant water molecules in dissociation.

Recent progress in molecular simulation methods for drug binding kinetics

Due to the contribution of drug-target binding kinetics to drug efficacy, there is a high level of interest in developing methods to predict drug-target binding kinetic parameters. During the review period, a wide range of enhanced sampling molecular dynamics simulation-based methods has been developed for computing drug-target binding kinetics and studying binding and unbinding mechanisms. Here, we assess the performance of these methods considering two benchmark systems in detail: mutant T4 lysozyme-ligand complexes and a large set of N-HSP90-inhibitor complexes. The results indicate that some of the simulation methods can already be usefully applied in drug discovery or lead optimization programs but that further studies on more high-quality experimental benchmark datasets are necessary to improve and validate computational methods.

Multisecond ligand dissociation dynamics from atomistic simulations

Coarse-graining of fully atomistic molecular dynamics simulations is a long-standing goal in order to allow the description of processes occurring on biologically relevant timescales. For example, the prediction of pathways, rates and rate-limiting steps in protein-ligand unbinding is crucial for modern drug discovery. To achieve the enhanced sampling, we perform dissipation-corrected targeted molecular dynamics simulations, which yield free energy and friction profiles of molecular processes under consideration. Subsequently, we use these fields to perform temperature-boosted Langevin simulations which account for the desired kinetics occurring on multisecond timescales and beyond. Adopting the dissociation of solvated sodium chloride, trypsin-benzamidine and Hsp90-inhibitor protein-ligand complexes as test problems, we reproduce rates from molecular dynamics simulation and experiments within a factor of 2–20, and dissociation constants within a factor of 1–4. Analysis of friction profiles reveals that binding and unbinding dynamics are mediated by changes of the surrounding hydration shells in all investigated systems.

Protein-ligand unbinding processes are out of reach for atomistic simulations due to time-scale involved. Here the authors demonstrate an approach relying on dissipation-corrected targeted molecular dynamics that enables to provide binding and unbinding rates with a speed-up of several orders of magnitude.

Multiscale simulation approaches to modeling drug-protein binding

Simulations can provide detailed insight into the molecular processes involved in drug action, such as protein-ligand binding, and can therefore be a valuable tool for drug design and development. Processes with a large range of length and timescales may be involved, and understanding these different scales typically requires different types of simulation methodology. Ideally, simulations should be able to connect across scales, to analyze and predict how changes at one scale can influence another. Multiscale simulation methods, which combine different levels of treatment, are an emerging frontier with great potential in this area. Here we review multiscale frameworks of various types, and selected applications to biomolecular systems with a focus on drug-ligand binding.

Transient States and Barriers from Molecular Simulations and the Milestoning Theory: Kinetics in Ligand-Protein Recognition and Compound Design

This study presents a novel computational approach to study molecular recognition and binding kinetics for drug-like compounds dissociating from a flexible protein system. The intermediates and their free energy profile during ligand association and dissociation processes control ligand-protein binding kinetics and bring a more complete picture of ligand-protein binding. The method applied the milestoning theory to extract kinetics and thermodynamics information from running short classical molecular dynamics (MD) simulations for frames from a given dissociation path. High-dimensional ligand-protein motions (3N-6 degrees of freedom) during ligand dissociation were reduced by use of principal component modes for assigning more than 100 milestones, and classical MD runs were allowed to travel multiple milestones to efficiently obtain ensemble distribution of initial structures for MD simulations and estimate the transition time and rate during ligand traveling between milestones. We used five pyrazolourea ligands and cyclin-dependent kinase 8 with cyclin C (CDK8/CycC) as our model system as well as metadynamics and a pathway search method to sample dissociation pathways. With our strategy, we constructed the free energy profile for highly mobile biomolecular systems. The computed binding free energy and residence time correctly ranked the pyrazolourea ligand series, in agreement with experimental data. Guided by a barrier of a ligand passing an αC helix and activation loop, we introduced one hydroxyl group to parent compounds to design our ligands with increased residence time and validated our prediction by experiments. This work provides a novel and robust approach to investigate dissociation kinetics of large and flexible systems for understanding unbinding mechanisms and designing new small-molecule drugs with desired binding kinetics.

Estimation of Protein-Ligand Unbinding Kinetics Using Non-Equilibrium Targeted Molecular Dynamics Simulations

We here report on nonequilibrium targeted molecular dynamics simulations as a tool for the estimation of protein-ligand unbinding kinetics. Correlating simulations with experimental data from SPR kinetics measurements and X-ray crystallography on two small molecule compound libraries bound to the N-terminal domain of the chaperone Hsp90, we show that the mean nonequilibrium work computed in an ensemble of trajectories of enforced ligand unbinding is a promising predictor for ligand unbinding rates. We furthermore investigate the molecular basis determining unbinding rates within the compound libraries. We propose ligand conformational changes and protein-ligand nonbonded interactions to impact on unbinding rates. Ligands may remain longer at the protein if they exhibit strong electrostatic and/or van der Waals interactions with the target. In the case of ligands with a rigid chemical scaffold that exhibit longer residence times, transient electrostatic interactions with the protein appear to facilitate unbinding. Our results imply that understanding the unbinding pathway and the protein-ligand interactions along this path is crucial for the prediction of small molecule ligands with defined unbinding kinetics.

Ligand Binding, Unbinding, and Allosteric Effects: Deciphering Small-Molecule Modulation of HSP90

The molecular chaperone HSP90 oversees the functional activation of a large number of client proteins. Because of its role in multiple pathways linked to cancer and neurodegeneration, drug discovery targeting HSP90 has been actively pursued. Yet, a number of inhibitors failed to meet expectations due to induced toxicity problems. In this context, allosteric perturbation has emerged as an alternative strategy for the pharmacological modulation of HSP90 functions. Specifically, novel allosteric stimulators showed the interesting capability of accelerating HSP90 closure dynamics and ATPase activities while inducing tumor cell death. Here, we gain atomistic insight into the mechanisms of allosteric ligand recognition and their consequences on the functional dynamics of HSP90, starting from the fully unbound state. We integrate advanced computational sampling methods based on FunnelMetadynamics, with the analysis of internal dynamics of the structural ensembles visited during the simulations. We observe several binding/unbinding events, and from these, we derive an accurate estimation of the absolute binding free energy. Importantly, we show that different binding poses induce different dynamics states. Our work for the first time explicitly correlates HSP90 responses to binding/unbinding of an allosteric ligand to the modulation of functionally oriented protein motions.

Accelerating Rare Dissociative Processes in Biomolecules Using Selectively Scaled MD Simulations

Molecular dynamics (MD) simulations can be a powerful tool for modeling complex dissociative processes such as ligand unbinding. However, many biologically relevant dissociative processes occur on timescales that far exceed the timescales of typical MD simulations. Here, we implement and apply an enhanced sampling method in which specific energy terms in the potential energy function are selectively "scaled" to accelerate dissociative events during simulations. Using ligand unbinding as an example of a complex dissociative process, we selectively scaled up ligand-water interactions in an attempt to increase the rate of ligand unbinding. After applying our selectively scaled MD (ssMD) approach to several cyclin-dependent kinase-inhibitor complexes, we discovered that we could accelerate ligand unbinding, thereby allowing, in some cases, unbinding events to occur within as little as 2 ns. Moreover, we found that we could make realistic estimates of the initial unbinding times (τ_unbind^sim) as well as the accompanying change in free energy (ΔG^sim) of the inhibitors from our ssMD simulation data. To accomplish this, we employed a previously described Kramers'-based rate extrapolation method and a newly described free energy extrapolation method. Because our ssMD approach is general, it should find utility as an easy-to-deploy, enhanced sampling method for modeling other dissociative processes.

Encounter complexes and hidden poses of kinase-inhibitor binding on the free-energy landscape

Modern drug discovery increasingly focuses on the drug-target binding kinetics which depend on drug (un)binding pathways. The conventional molecular dynamics simulation can observe only a few binding events even using the fastest supercomputer. Here, we develop 2D gREST/REUS simulation with enhanced flexibility of the ligand and the protein binding site. Simulation (43 μs in total) applied to an inhibitor binding to c-Src kinase covers 100 binding and unbinding events. On the statistically converged free-energy landscapes, we succeed in predicting the X-ray binding structure, including water positions. Furthermore, we characterize hidden semibound poses and transient encounter complexes on the free-energy landscapes. Regulatory residues distant from the catalytic core are responsible for the initial inhibitor uptake and regulation of subsequent bindings, which was unresolved by experiments. Stabilizing/blocking of either the semibound poses or the encounter complexes can be an effective strategy to optimize drug-target residence time.

Time-Domain Analysis of Molecular Dynamics Trajectories Using Deep Neural Networks: Application to Activity Ranking of Tankyrase Inhibitors

Molecular dynamics simulations provide valuable insights into the behavior of molecular systems. Extending the recent trend of using machine learning techniques to predict physicochemical properties from molecular dynamics data, we propose to consider the trajectories as multidimensional time series represented by 2D tensors containing the ligand-protein interaction descriptor values for each time step. Similar in structure to the time series encountered in modern approaches for signal, speech, and natural language processing, these time series can be directly analyzed using long short-term memory (LSTM) recurrent neural networks or convolutional neural networks (CNNs). The predictive regression models for the ligand-protein affinity were built for a subset of the PDBbind v.2017 database and applied to inhibitors of tankyrase, an enzyme of the poly(ADP-ribose)-polymerase (PARP) family that can be used in the treatment of colorectal cancer. As an additional test set, a subset of the Community Structure-Activity Resource (CSAR) data set was used. For comparison, the random forest and simple neural network models based on the crystal pose or the trajectory-averaged descriptors were used, as well as the commonly employed docking and molecular mechanics Poisson-Boltzmann surface area (MM-PBSA) scores. Convolutional neural networks based on the 2D tensors of ligand-protein interaction descriptors for short (2 ns) trajectories provide the best accuracy and predictive power, reaching the Spearman rank correlation coefficient of 0.73 and Pearson correlation coefficient of 0.70 for the tankyrase test set. Taking into account the recent increase in computational power of modern GPUs and relatively low computational complexity of the proposed approach, it can be used as an advanced virtual screening filter for compound prioritization.

Machine Learning Analysis of τRAMD Trajectories to Decipher Molecular Determinants of Drug-Target Residence Times

Drug-target residence times can impact drug efficacy and safety, and are therefore increasingly being considered during lead optimization. For this purpose, computational methods to predict residence times, τ, for drug-like compounds and to derive structure-kinetic relationships are desirable. A challenge for approaches based on molecular dynamics (MD) simulation is the fact that drug residence times are typically orders of magnitude longer than computationally feasible simulation times. Therefore, enhanced sampling methods are required. We recently reported one such approach: the τRAMD procedure for estimating relative residence times by performing a large number of random acceleration MD (RAMD) simulations in which ligand dissociation occurs in times of about a nanosecond due to the application of an additional randomly oriented force to the ligand. The length of the RAMD simulations is used to deduce τ. The RAMD simulations also provide information on ligand egress pathways and dissociation mechanisms. Here, we describe a machine learning approach to systematically analyze protein-ligand binding contacts in the RAMD trajectories in order to derive regression models for estimating τ and to decipher the molecular features leading to longer τ values. We demonstrate that the regression models built on the protein-ligand interaction fingerprints of the dissociation trajectories result in robust estimates of τ for a set of 94 drug-like inhibitors of heat shock protein 90 (HSP90), even for the compounds for which the length of the RAMD trajectories does not provide a good estimation of τ. Thus, we find that machine learning helps to overcome inaccuracies in the modeling of protein-ligand complexes due to incomplete sampling or force field deficiencies. Moreover, the approach facilitates the identification of features important for residence time. In particular, we observed that interactions of the ligand with the sidechain of F138, which is located on the border between the ATP binding pocket and a hydrophobic transient sub-pocket, play a key role in slowing compound dissociation. We expect that the combination of the τRAMD simulation procedure with machine learning analysis will be generally applicable as an aid to target-based lead optimization.

Energetics of Storage and Diffusion of Water and Cyclo-Octasulfur for a Nonpolar Cavity of RHCC Tetrabrachion by Molecular Dynamics Simulations

Tetrabrachion forms the key component of the S-layer of Staphylothermus marinus. Molecular dynamics simulations have been used to study the energetics of occupancy of cavity 3 of the right-handed coiled-coil stalk of tetrabrachion by both water molecules and cyclooctasulfur S₈ crowns, as well as to determine possible pathways and free energy barriers for the diffusion of both water and cyclooctasulfur through the peptide walls of RHCC tetrabrachion between cavity 3 and bulk solvent. Calculations of the transfer free energy from solvent to cavity show that clusters of six, seven and eight water molecules are marginally stable in cavity 3, but that occupancy of the cavity by a cyclooctasulfur ring is favoured significantly over water clusters of all sizes. Thermal activation simulations at T = 400K revealed that water molecules diffusing through the wall pass through a sequence of metastable configurations where they are temporarily immobilized by forming networks of hydrogen bonds with specific wall residues. Calculations of the free energy of these metastable configurations using multi-configurational thermodynamic integration yielded a free energy profile with a principal free energy maximum ∆G~50 kJ/mol and a slight activation asymmetry in favour of the direction from cavity to solvent. Potential exit pathways for cyclooctasulfur were investigated with the methods of steered molecular dynamics and umbrella sampling. The cyclooctasulfur was steered through a gap in the tetrabrachion wall along a linear path from cavity 3 into the solvent and the resulting trajectory was subdivided into 22 sampling windows. The free energy profile created for the trajectory by umbrella sampling showed a sharp principal maximum as a function of the reaction coordinate with asymmetric free energy barriers ∆G_exit~220 kJ/mol and ∆G_entrance~100 kJ/mol for cavity exit and entrance, respectively.