Protein-ligand (un)binding kinetics as a new paradigm for drug discovery at the crossroad between experiments and modelling

Pseudoirreversible inhibition elicits persistent efficacy of a sphingosine 1-phosphate receptor 1 antagonist

Sphingosine 1-phosphate receptor 1 (S1PR1), a G protein-coupled receptor, is required for lymphocyte trafficking, and is a promising therapeutic target in inflammatory diseases. Here, we synthesize a competitive S1PR1 antagonist, KSI-6666, that effectively suppresses pathogenic inflammation. Metadynamics simulations suggest that the interaction of KSI-6666 with a methionine residue Met124 in the ligand-binding pocket of S1PR1 may inhibit the dissociation of KSI-6666 from S1PR1. Consistently, in vitro functional and mutational analyses reveal that KSI-6666 causes pseudoirreversible inhibition of S1PR1, dependent on the Met124 of the protein and substituents on the distal benzene ring of KSI-6666. Moreover, in vivo study suggests that this pseudoirreversible inhibition is responsible for the persistent activity of KSI-6666.

MolDy: molecular dynamics simulation made easy

MOTIVATION

Molecular dynamics (MD) is a computational experiment that is crucial for understanding the structure of biological macro and micro molecules, their folding, and the inter-molecular interactions. Accurate knowledge of these structural features is the cornerstone in drug development and elucidating macromolecules functions. The open-source GROMACS biomolecular MD simulation program is recognized as a reliable and frequently used simulation program for its precision. However, the user requires expertise, and scripting skills to carrying out MD simulations.

RESULTS

We have developed an end-to-end interactive MD simulation application, MolDy for Gromacs. This front-end application provides a customizable user interface integrated with the Python and Perl-based logical backend connecting the Linux shell and Gromacs software. The tool performs analysis and provides the user with simulation trajectories and graphical representations of relevant biophysical parameters. The advantages of MolDy are (i) user-friendly, does not requiring the researcher to have prior knowledge of Linux; (ii) easy installation by a single command; (iii) freely available for academic research; (iv) can run with minimum configuration of operating systems; (v) has valid default prefilled parameters for beginners, and at the same time provides scope for modifications for expert users.

AVAILABILITY AND IMPLEMENTATION

MolDy is available freely as compressed source code files with user manual for installation and operation on GitHub: https://github.com/AIBResearchMolDy/Moldyv01.git and on https://aibresearch.com/innovations.

Structural and mechanistic insights into Quinolone Synthase to address its functional promiscuity

Quinolone synthase from Aegle marmelos (AmQNS) is a type III polyketide synthase that yields therapeutically effective quinolone and acridone compounds. Addressing the structural and molecular underpinnings of AmQNS and its substrate interaction in terms of its high selectivity and specificity can aid in the development of numerous novel compounds. This paper presents a high-resolution AmQNS crystal structure and explains its mechanistic role in synthetic selectivity. Additionally, we provide a model framework to comprehend structural constraints on ketide insertion and postulate that AmQNS's steric and electrostatic selectivity plays a role in its ability to bind to various core substrates, resulting in its synthetic diversity. AmQNS prefers quinolone synthesis and can accommodate large substrates because of its wide active site entrance. However, our research suggests that acridone is exclusively synthesized in the presence of high malonyl-CoA concentrations. Potential implications of functionally relevant residue mutations were also investigated, which will assist in harnessing the benefits of mutations for targeted polyketide production. The pharmaceutical industry stands to gain from these findings as they expand the pool of potential drug candidates, and these methodologies can also be applied to additional promising enzymes.

In silico investigation of falcipain-2 inhibition by hybrid benzimidazole-thiosemicarbazone antiplasmodial agents: A molecular docking, molecular dynamics simulation, and kinetics study

The emergence of artemisinin-resistant variants of Plasmodium falciparum necessitates the urgent search for novel antimalarial drugs. In this regard, an in silico study to screen antimalarial drug candidates from a series of benzimidazole-thiosemicarbazone hybrid molecules with interesting antiplasmodial properties and explore their falcipain-2 (FP2) inhibitory potentials has been undertaken herein. FP2 is a key cysteine protease that degrades hemoglobin in Plasmodium falciparum and is an important biomolecular target in the development of antimalarial drugs. Pharmacokinetic properties, ADMET profiles, MM/GBSA-based binding free energies, reaction mechanisms, and associated barrier heights have been investigated. DFT, molecular dynamics simulation, molecular docking, and ONIOM methods were used. From the results obtained, four 4N-substituted derivatives of the hybrid molecule (E)-2-(1-(5-chloro-1H-benzo[d]imidazol-2-yl)ethylidene)hydrazine-1-carbothioamide (1A) denoted 1B, 1C, 1D, and 1E are drug-like and promising inhibitors of FP2, exhibiting remarkably small inhibitory constants (5.94 × 10-14 - 2.59 × 10-04 n M) and favorable binding free energies (-30.32 to -17.17 kcal/mol). Moreover, the ONIOM results have revealed that 1B and possibly 1C and 1D may act as covalent inhibitors of FP2. The rate-determining step of the thermodynamically favorable covalent binding mechanism occurs across a surmountable barrier height of 24.18 kcal/mol in water and 28.42 kcal/mol in diethyl ether. Our findings are useful for further experimental investigations on the antimalarial activities of the hybrid molecules studied.

Fluorescence-Based Protein Stability Monitoring-A Review

Proteins are large biomolecules with a specific structure that is composed of one or more long amino acid chains. Correct protein structures are directly linked to their correct function, and many environmental factors can have either positive or negative effects on this structure. Thus, there is a clear need for methods enabling the study of proteins, their correct folding, and components affecting protein stability. There is a significant number of label-free methods to study protein stability. In this review, we provide a general overview of these methods, but the main focus is on fluorescence-based low-instrument and -expertise-demand techniques. Different aspects related to thermal shift assays (TSAs), also called differential scanning fluorimetry (DSF) or ThermoFluor, are introduced and compared to isothermal chemical denaturation (ICD). Finally, we discuss the challenges and comparative aspects related to these methods, as well as future opportunities and assay development directions.

Computational Workflow for Refining AlphaFold Models in Drug Design Using Kinetic and Thermodynamic Binding Calculations: A Case Study for the Unresolved Inactive Human Adenosine A3 Receptor

A structure-based drug design pipeline that considers both thermodynamic and kinetic binding data of ligands against a receptor will enable the computational design of improved drug molecules. For unresolved GPCR-ligand complexes, a workflow that can apply both thermodynamic and kinetic binding data in combination with alpha-fold (AF)-derived or other homology models and experimentally resolved binding modes of relevant ligands in GPCR-homologs needs to be tested. Here, as test case, we studied a congeneric set of ligands that bind to a structurally unresolved G protein-coupled receptor (GPCR), the inactive human adenosine A3 receptor (hA3R). We tested three available homology models from which two have been generated from experimental structures of hA1R or hA2AR and one model was a multistate alphafold 2 (AF2)-derived model. We applied alchemical calculations with thermodynamic integration coupled with molecular dynamics (TI/MD) simulations to calculate the experimental relative binding free energies and residence time (τ)-random accelerated MD (τ-RAMD) simulations to calculate the relative residence times (RTs) for antagonists. While the TI/MD calculations produced, for the three homology models, good Pearson correlation coefficients, correspondingly, r = 0.74, 0.62, and 0.67 and mean unsigned error (mue) values of 0.94, 1.31, and 0.81 kcal mol-1, the τ-RAMD method showed r = 0.92 and 0.52 for the first two models but failed to produce accurate results for the multistate AF2-derived model. With subsequent optimization of the AF2-derived model by reorientation of the side chain of R1735.34 located in the extracellular loop 2 (EL2) that blocked ligand's unbinding, the computational model showed r = 0.84 for kinetic data and improved performance for thermodynamic data (r = 0.81, mue = 0.56 kcal mol-1). Overall, after refining the multistate AF2 model with physics-based tools, we were able to show a strong correlation between predicted and experimental ligand relative residence times and affinities, achieving a level of accuracy comparable to an experimental structure. The computational workflow used can be applied to other receptors, helping to rank candidate drugs in a congeneric series and enabling the prioritization of leads with stronger binding affinities and longer residence times.

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.

Uncovering the Structural and Binding Insights of Dual Inhibitors Simultaneously Targeting Two Distinct Sites on EGFR Kinase

Epidermal growth factor receptor (EGFR) is the first growth factor receptor identified in normal cells that is related to the receptor tyrosine kinase, which causes regular cell division. A point mutation in EGFR intracellular kinase domain forces the abnormal cell divisions throughout time, leading to non-small cell lung cancer (NSCLC) transformation. Thus, competitive inhibitors that bind to the ATP binding pocket have been developed as a targeted therapy for NSCLC. The third-generation kinase inhibitor Osimertinib is currently playing a very vital role in the treatment of NSCLC. However, it is not effective toward the C797S kinase domain mutation. For this reason, fourth-generation kinase noncompetitive inhibitors are introduced which work through binding to an allosteric pocket near the ATP binding region and act as a better binding agent for this mutated kinase domain. However, the problem is that these single fourth-generation kinase inhibitors may not be as effective as a single agent. The aim of this work was to apply combinations of these two inhibitors together in different binding regions of EGFR without overlapping the resistance mechanism to obtain the key direct and indirect interactions occurring between them. Moreover, the free energy of dissociation of an inhibitor from its binding sites in the presence of a second inhibitor immobilized in another binding site was also the focus of the study. To realize this aim, we performed conventional molecular dynamics simulations and principal component analysis and dynamic cross-correlation matrices along with umbrella sampling. Our results demonstrated that binding of dual inhibitors triggered conformational changes of the protein more toward the inactive state. Furthermore, allosteric inhibitors bound more strongly to protein kinase EGFR than the orthosteric inhibitors in the presence of dual inhibitors. Finally, the binding mechanism and important hydrogen-bonding residues during unbinding of the inhibitors were fully elucidated. This study provides insight into the binding of the receptor-orthosteric inhibitor-allosteric inhibitor, which can be helpful for further design of novel inhibitors that have a better inhibitory action.

Electrochemical Nanoreactor Provides a Comprehensive View of Isocitrate Dehydrogenase Cancer-drug Kinetics

The ability to control enzyme cascades entrapped in a nanoporous electrode material (the "Electrochemical Leaf", e-Leaf) has been exploited to gain detailed kinetic insight into the mechanism of an anti-cancer drug. Ivosidenib, used to treat acute myeloid leukemia, acts on a common cancer-linked variant of isocitrate dehydrogenase 1 (IDH1 R132H) inhibiting its "gain-of-function" activity-the undesired reduction of 2-oxoglutarate (2OG) to the oncometabolite 2-hydroxyglutarate (2HG). The e-Leaf quantifies the kinetics of IDH1 R132H inhibition across a wide and continuous range of conditions, efficiently revealing factors underlying the inhibitor residence time. Selective inhibition of IDH1 R132H by Ivosidenib and another inhibitor, Novartis 224, is readily resolved as a two-stage process whereby initial rapid non-inhibitory binding is followed by a slower step to give the inhibitory complex. These kinetic features are likely present in other allosteric inhibitors of IDH1/2. Such details, essential for understanding inhibition mechanisms, are not readily resolved in conventional steady-state kinetics or by techniques that rely only on measuring binding. Extending the new method and analytical framework presented here to other enzyme systems will be straightforward and should rapidly reveal insight that is difficult or often impossible to obtain using other methods.

Electrochemical Nanoreactor Provides a Comprehensive View of Isocitrate Dehydrogenase Cancer-drug Kinetics

The ability to control enzyme cascades entrapped in a nanoporous electrode material (the "Electrochemical Leaf", e-Leaf) has been exploited to gain detailed kinetic insight into the mechanism of an anti-cancer drug. Ivosidenib, used to treat acute myeloid leukemia, acts on a common cancer-linked variant of isocitrate dehydrogenase 1 (IDH1 R132H) inhibiting its "gain-of-function" activity-the undesired reduction of 2-oxoglutarate (2OG) to the oncometabolite 2-hydroxyglutarate (2HG). The e-Leaf quantifies the kinetics of IDH1 R132H inhibition across a wide and continuous range of conditions, efficiently revealing factors underlying the inhibitor residence time. Selective inhibition of IDH1 R132H by Ivosidenib and another inhibitor, Novartis 224, is readily resolved as a two-stage process whereby initial rapid non-inhibitory binding is followed by a slower step to give the inhibitory complex. These kinetic features are likely present in other allosteric inhibitors of IDH1/2. Such details, essential for understanding inhibition mechanisms, are not readily resolved in conventional steady-state kinetics or by techniques that rely only on measuring binding. Extending the new method and analytical framework presented here to other enzyme systems will be straightforward and should rapidly reveal insight that is difficult or often impossible to obtain using other methods.

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.

Physiologically Based Pharmacokinetic Model of Brain Delivery of Plasma Protein Bound Drugs

INTRODUCTION

A physiologically based pharmacokinetic (PBPK) model is developed that focuses on the kinetic parameters of drug association and dissociation with albumin, alpha-1 acid glycoprotein (AGP), and brain tissue proteins, as well as drug permeability at the blood-brain barrier, drug metabolism, and brain blood flow.

GOAL

The model evaluates the extent to which plasma protein-mediated uptake (PMU) of drugs by brain influences the concentration of free drug both within the brain capillary compartment in vivo and the brain compartment. The model also studies the effect of drug binding to brain tissue proteins on the concentration of free drug in brain.

METHODS

The steady state and non-steady state PBPK models are comprised of 11-12 variables, and 18-23 parameters, respectively. Two model drugs are analyzed: propranolol, which undergoes modest PMU from the AGP-bound pool, and imipramine, which undergoes a high degree of PMU from both the albumin-bound and AGP-bound pools in plasma.

RESULTS

The free propranolol concentration in brain is under-estimated 2- to fourfold by in vitro measurements of free plasma propranolol, and the free imipramine concentration in brain is under-estimated by 18- to 31-fold by in vitro measurements of free imipramine in plasma. The free drug concentration in brain in vivo is independent of drug binding to brain tissue proteins.

CONCLUSIONS

In vitro measurement of free drug concentration in plasma under-estimates the free drug in brain in vivo if PMU in vivo from either the albumin and/or the AGP pools in plasma takes place at the BBB surface.

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.

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.

Estimation of Drug-Target Residence Time by Targeted Molecular Dynamics Simulations

Drug-target residence time has emerged as a key selection factor in drug discovery since the binding duration of a drug molecule to its protein target can significantly impact its in vivo efficacy. The challenge in studying the residence time, in early drug discovery stages, lies in how to cost-effectively determine the residence time for the systematic assessment of compounds. Currently, there is still a lack of computational protocols to quickly estimate such a measure, particularly for large and flexible protein targets and drugs. Here, we report an efficient computational protocol, based on targeted molecular dynamics, to rank drug candidates by their residence time and to obtain insights into ligand-target dissociation mechanisms. The method was assessed on a dataset of 10 arylpyrazole inhibitors of CDK8, a large, flexible, and clinically important target, for which the experimental residence time of the inhibitors ranges from minutes to hours. The compounds were correctly ranked according to their estimated residence time scores compared to their experimental values. The analysis of protein-ligand interactions along the dissociation trajectories highlighted the favorable contribution of hydrophobic contacts to residence time and revealed key residues that strongly affect compound residence time.

Computationally Elucidating the Binding Kinetics for Different AChE Inhibitors to Access the Rationale for Improving the Drug Efficacy

Traditional drug discovery is based on a binding affinity (thermodynamics)-driven paradigm. Numerous examples, however, demonstrated that drug efficacy does not always depend only on binding affinity but positively correlates with binding kinetics, that is, the dissociation rate constant (k_off). Binding free energy landscape (BFEL) constructor is a computational binding kinetics prediction method, previously developed by us, that estimates the binding kinetics for ligand-protein based on their constructed binding free energy landscape, but it also reveals the detailed molecular mechanism of the binding event, hence, providing the position of transition states at the molecular level to modify/improve the binding kinetics. Acetylcholinesterase (AChE) is a well-known Alzheimer's disease (AD) target for which there is still not an ideal drug on the market. Therefore, to improve the drug design strategy for AD, the binding kinetics and binding molecular mechanisms of the four inhibitors of AChE, that is, E2020 (Aricept), HupA, Rivastigmine, and Galantamine, were studied. Also, the differentiation of the binding kinetics between mAChE and TcAChE was studied to evaluate the sensitiveness of BFEL constructor. The flexibility of molecules has a noticeable effect on the nature of BFEL. To the same target, flexible molecules (i.e., E2020 and Rivastigmine) which contain more rotatable bonds tend to have more complicated BFELs reflecting more complicated molecular action mechanisms than the rigid ones (i.e., HupA and Galantamine), which therefore could be more challenging to be optimized. The binding kinetics is highly dependent on the structure of the molecules, such as the length and the functional groups. Therefore, E2020 presents better binding kinetic and thermodynamic properties with either TcAChE or mAChE. Therefore, it is the most promising lead drug for binding kinetics-based drug design. In addition, the binding kinetics of a drug may present different values in the proteins of different organisms because the residue compositions of the binding gorges of the targets are variant, that is, E2020 shows lower binding affinity and association energy barrier in binding with mAChE than TcAChE. However, HupA presents a better binding property with TcAChE than mAChE.

LINES: Log-Probability Estimation via Invertible Neural Networks for Enhanced Sampling

It is very challenging to sample a molecular process with large activation energies using molecular dynamics simulations. Current enhanced sampling methodologies, such as umbrella sampling and metadynamics, rely on the identification of appropriate reaction coordinates for a system. In this paper, we developed a method for log-probability estimation via invertible neural networks for enhanced sampling (LINES). This iterative scheme utilizes a normalizing flow machine learning model to learn the underlying free energy surface (FES) of a system as a function of molecular coordinates and then applies a gradient-based optimization method to the learned normalizing flow to identify reaction coordinates. A biasing potential is then evaluated over a tabulated grid of the reaction coordinate values, which can be applied to the next round of simulations for enhanced sampling, resulting in more efficient sampling. We tested the accuracy and efficiency of the LINES method in sampling the FES using the alanine dipeptide system. We also demonstrated the effectiveness of identification of reaction coordinates through simulation of cyclobutanol unbinding from β-cyclodextrin and the folding/unfolding of CLN025─a variant of the peptide Chignolin. The LINES method can be extended to the study of large-scale protein systems with complex nonlinear reaction pathways.

Application of molecular dynamics simulation in biomedicine

Molecular dynamics (MD) simulation has been widely used in the field of biomedicine to study the conformational transition of proteins caused by mutation or ligand binding/unbinding. It provides some perspectives those are difficult to find in traditional biochemical or pathological experiments, for example, detailed effects of mutations on protein structure and protein-protein/ligand interaction at the atomic level. In this review, a broad overview on conformation changes and drug discovery by MD simulation is given. We first discuss the preparation of protein structure for MD simulation, which is a key step that determines the accuracy of the simulation. Then, we summarize the applications of commonly used force fields and MD simulations in scientific research. Finally, enhanced sampling methods and common applications of these methods are introduced. In brief, MD simulation is a powerful tool and it can be used to guide experimental study. The combination of MD simulation and experimental techniques is an a priori means to solve the biomedical problems and give a deep understanding on the relationship between protein structure and function.

RNA-ligand molecular docking: advances and challenges

With rapid advances in computer algorithms and hardware, fast and accurate virtual screening has led to a drastic acceleration in selecting potent small molecules as drug candidates. Computational modeling of RNA-small molecule interactions has become an indispensable tool for RNA-targeted drug discovery. The current models for RNA-ligand binding have mainly focused on the docking-and-scoring method. Accurate docking and scoring should tackle four crucial problems: (1) conformational flexibility of ligand, (2) conformational flexibility of RNA, (3) efficient sampling of binding sites and binding poses, and (4) accurate scoring of different binding modes. Moreover, compared with the problem of protein-ligand docking, predicting ligand binding to RNA, a negatively charged polymer, is further complicated by additional effects such as metal ion effects. Thermodynamic models based on physics-based and knowledge-based scoring functions have shown highly encouraging success in predicting ligand binding poses and binding affinities. Recently, kinetic models for ligand binding have further suggested that including dissociation kinetics (residence time) in ligand docking would result in improved performance in estimating in vivo drug efficacy. More recently, the rise of deep-learning approaches has led to new tools for predicting RNA-small molecule binding. In this review, we present an overview of the recently developed computational methods for RNA-ligand docking and their advantages and disadvantages.

Practical Protocols for Efficient Sampling of Kinase-Inhibitor Binding Pathways Using Two-Dimensional Replica-Exchange Molecular Dynamics

Molecular dynamics (MD) simulations are increasingly used to study various biological processes such as protein folding, conformational changes, and ligand binding. These processes generally involve slow dynamics that occur on the millisecond or longer timescale, which are difficult to simulate by conventional atomistic MD. Recently, we applied a two-dimensional (2D) replica-exchange MD (REMD) method, which combines the generalized replica exchange with solute tempering (gREST) with the replica-exchange umbrella sampling (REUS) in kinase-inhibitor binding simulations, and successfully observed multiple ligand binding/unbinding events. To efficiently apply the gREST/REUS method to other kinase-inhibitor systems, we establish modified, practical protocols with non-trivial simulation parameter tuning. The current gREST/REUS simulation protocols are tested for three kinase-inhibitor systems: c-Src kinase with PP1, c-Src kinase with Dasatinib, and c-Abl kinase with Imatinib. We optimized the definition of kinase-ligand distance as a collective variable (CV), the solute temperatures in gREST, and replica distributions and umbrella forces in the REUS simulations. Also, the initial structures of each replica in the 2D replica space were prepared carefully by pulling each ligand from and toward the protein binding sites for keeping stable kinase conformations. These optimizations were carried out individually in multiple short MD simulations. The current gREST/REUS simulation protocol ensures good random walks in 2D replica spaces, which are required for enhanced sampling of inhibitor dynamics around a target kinase.

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.

The prediction of protein-ligand unbinding for modern drug discovery

INTRODUCTION

Drug-target thermodynamic and kinetic information have perennially important roles in drug design. The prediction of protein-ligand unbinding, which can provide important kinetic information, in experiments continues to face great challenges. Uncovering protein-ligand unbinding through molecular dynamics simulations has become efficient and inexpensive with the progress and enhancement of computing power and sampling methods.

AREAS COVERED

In this review, various sampling methods for protein-ligand unbinding and their basic principles are firstly briefly introduced. Then, their applications in predicting aspects of protein-ligand unbinding, including unbinding pathways, dissociation rate constants, residence time and binding affinity, are discussed.

EXPERT OPINION

Although various sampling methods have been successfully applied in numerous systems, they still have shortcomings and deficiencies. Most enhanced sampling methods require researchers to possess a wealth of prior knowledge of collective variables or reaction coordinates. In addition, most systems studied at present are relatively simple, and the study of complex systems in real drug research remains greatly challenging. Through the combination of machine learning and enhanced sampling methods, prediction accuracy can be further improved, and some problems encountered in complex systems also may be solved.

Xeno-Nucleic Acid (XNA) 2'-Fluoro-Arabino Nucleic Acid (FANA) Aptamers to the Receptor-Binding Domain of SARS-CoV-2 S Protein Block ACE2 Binding

The causative agent of COVID-19, SARS-CoV-2, gains access to cells through interactions of the receptor-binding domain (RBD) on the viral S protein with angiotensin-converting enzyme 2 (ACE2) on the surface of human host cells. Systematic evolution of ligands by exponential enrichment (SELEX) was used to generate aptamers (nucleic acids selected for high binding affinity to a target) to the RBD made from 2ʹ-fluoro-arabinonucleic acid (FANA). The best selected ~79 nucleotide aptamers bound the RBD (Arg319-Phe541) and the larger S1 domain (Val16-Arg685) of the 1272 amino acid S protein with equilibrium dissociation constants (K_D,app) of ~10–20 nM, and binding half-life for the RBD, S1 domain, and full trimeric S protein of 53 ± 18, 76 ± 5, and 127 ± 7 min, respectively. Aptamers inhibited the binding of the RBD to ACE2 in an ELISA assay. Inhibition, on a per weight basis, was similar to neutralizing antibodies that were specific for RBD. Aptamers demonstrated high specificity, binding with about 10-fold lower affinity to the related S1 domain from the original SARS virus, which also binds to ACE2. Overall, FANA aptamers show affinities comparable to previous DNA aptamers to RBD and S1 protein and directly block receptor interactions while using an alternative Xeno-nucleic acid (XNA) platform.

A novel assay based on pre-equilibrium titration curves for the determination of enzyme inhibitor binding kinetics

Selection of pharmacological agents based on potency measurements performed at equilibrium fail to incorporate the kinetic aspects of the drug–target interaction. Here we describe a method for screening or characterization of enzyme inhibitors that allows the concomitant determination of the equilibrium inhibition constant in unison with rates of complex formation and dissociation. The assay is distinct from conventional enzymatic assays and is based on the analysis of inhibition curves recorded prior to full equilibration of the system. The methodology is illustrated using bicyclic peptide inhibitors of the serine protease plasma kallikrein.

Identifying requirements for RSK2 specific inhibitors

Identifying isoform-specific inhibitors for closely related kinase family members remains a substantial challenge. The necessity for achieving this specificity is exemplified by the RSK family, downstream effectors of ERK1/2, which have divergent physiological effects. The natural product, SL0101, a flavonoid glycoside, binds specifically to RSK1/2 through a binding pocket generated by an extensive conformational rearrangement within the RSK N-terminal kinase domain (NTKD). In modelling experiments a single amino acid that is divergent in RSK3/4 most likely prevents the required conformational rearrangement necessary for SL0101 binding. Kinetic analysis of RSK2 association with SL0101 and its derivatives identified that regions outside of the NTKD contribute to stable inhibitor binding. An analogue with an n-propyl-carbamate at the 4” position on the rhamnose moiety was identified that forms a highly stable inhibitor complex with RSK2 but not with RSK1. These results identify a SL0101 modification that will aid the identification of RSK2 specific inhibitors.

Binding kinetics of liposome conjugated E-selectin and P-selectin glycoprotein ligand-1 measured with atomic force microscopy

Various ligand-functionalized liposomes have been developed for targeted therapies. Typically, the binding properties of the ligands and targeted proteins are measured with surface plasmon resonance (SPR), where the proteins are immobilized on a rigid surface. However, the difference of protein-ligand binding kinetics between liposome-conjugated protein and rigid surface-conjugated protein is not fully understood. In this work, the binding kinetics of P-selectin glycoprotein ligand-1 (PSGL-1) and E-selectin conjugated on liposome and on rigid surfaces are investigated with Atomic Force Microscopy (AFM). The results suggest that protein orientation and diffusion on liposomal membrane can alter the binding kinetics of the protein-ligand interaction. Specifically, the association and dissociation rate constant of AFM probe-conjugated E-selectin and glass-conjugated PSGL-1 are measured as 9.32 × 10⁴ M^-1s^-1 and 1.54 s^-1, respectively. While for the liposome-conjugated E-selectin and glass-conjugated PSGL-1, the kinetic constants are measured as 5.00 × 10⁷ M^-1s^-1 and 2.76 s^-1, respectively. Thus, there is an order's magnitude increase of binding affinity (from k_d = 16.51 μM to k_d = 0.06 μM) when protein is attached to liposome compared to attached to a rigid surface. The results might provide better understanding and pave the way for the future design of the ligand-targeted liposomes.

Combining Well-Tempered Metadynamics Simulation and SPR Assays to Characterize the Binding Mechanism of the Universal T-Lymphocyte Tetanus Toxin Epitope TT830-843

Peptide TT830-843 from the tetanus toxin is a universal T-cell epitope. It helps in vaccination and induces T-cell activation. However, the fine molecular interaction between this antigen and the major histocompatibility complex (MHC) remains unknown. Molecular analysis of its interaction with murine MHC (H-2) was proposed to explore its immune response efficiency. Molecular dynamics simulations are important mechanisms for understanding the basis of protein-ligand interactions, and metadynamics is a useful technique for enhancing sampling in molecular dynamics. SPR (surface plasmon resonance) assays were used to validate whether the metadynamics results are in accordance with the experimental results. The peptide TT830-843 unbinding process was simulated, and the free energy surface reconstruction revealed a detailed conformational landscape. The simulation described the exiting path as a stepwise mechanism between progressive detachment states. We pointed out how the terminus regions act as anchors for binding and how the detachment mechanism includes the opening of α-helices to permit the peptide's central region dissociation. The results indicated the peptide/H-2 receptor encounter occurs within a distance lesser than 27.5 Å, and the encounter can evolve to form a stable complex. SPR assays confirmed the complex peptide/H-2 as a thermodynamically stable system, exhibiting enough free energy to interact with TCR on the antigen-presenting cell surface. Therefore, combining in silico and in vitro assays provided significant evidence to support the peptide/H-2 complex formation.

Prediction of the Drug-Target Binding Kinetics for Flexible Proteins by Comparative Binding Energy Analysis

There is growing consensus that the optimization of the kinetic parameters for drug-protein binding leads to improved drug efficacy. Therefore, computational methods have been developed to predict kinetic rates and to derive quantitative structure-kinetic relationships (QSKRs). Many of these methods are based on crystal structures of ligand-protein complexes. However, a drawback is that each ligand-protein complex is usually treated as having a single structure. Here, we present a modification of COMparative BINding Energy (COMBINE) analysis, which uses the structures of ligand-protein complexes to predict binding parameters. We introduce the option of using multiple structures to describe each ligand-protein complex in COMBINE analysis and apply this to study the effects of protein flexibility on the derivation of dissociation rate constants (k_off) for inhibitors of p38 mitogen-activated protein (MAP) kinase, which has a flexible binding site. Multiple structures were obtained for each ligand-protein complex by performing docking to an ensemble of protein configurations obtained from molecular dynamics simulations. Coefficients to scale ligand-protein interaction energies determined from energy-minimized structures of ligand-protein complexes were obtained by partial least squares regression, and they allowed for the computation of k_off values. The QSKR model obtained using single, energy-minimized crystal structures for each ligand-protein complex had higher predictive power than the QSKR model obtained with multiple structures from ensemble docking. However, incorporation of ligand-protein flexibility helped to highlight additional ligand-protein interactions that lead to longer residence times, such as interactions with residues Arg67 and Asp168, which are close to the ligand in many crystal structures. These results show that COMBINE analysis is a promising method to guide the design of compounds that bind to flexible proteins with improved binding kinetics.

Modifying the catalytic preference of alpha-amylase toward n-alkanes for bioremediation purposes using in silico strategies

Since the beginning of oil exploration, whole ecosystems have been affected by accidents and bad practices involving petroleum compounds. In this sense, bioremediation stands out as the cheapest and most eco-friendly alternatives to reverse the damage done in oil-impacted areas. However, more efforts must be made to engineer enzymes that could be used in the bioremediation process. Interestingly, a recent work described that α-amylase, one of the most evolutionary conserved enzymes, was able to promiscuously degrade n-alkanes, a class of molecules abundant in the petroleum admixture. Considering that α-amylase is expressed in almost all known organisms, and employed in numerous biotechnological processes, using it can be a great leap toward more efficient applications of enzyme or microorganism-consortia bioremediation approaches. In this work, we employed a strict computational approach to design new α-amylase mutants with potentially enhanced catalytic efficiency toward n-alkanes. Using in silico techniques, such as molecular docking, molecular dynamics, metadynamics, and residue-residue interaction networks, we generated mutants potentially more efficient for degrading n-alkanes, L183Y, and N314A. Our results indicate that the new mutants have an increased binding rate for tetradecane, the longest n-alkane previously tested, which can reside in the catalytic center for more extended periods. Additionally, molecular dynamics and network analysis showed that the new mutations have no negative impact on protein structure than the WT. Our results aid in solidifying this enzyme as one more tool in the petroleum bioremediation toolbox.

Impact of the Coverage of Aptamers on a Nanoparticle on the Binding Equilibrium and Kinetics between Aptamer and Protein

Knowledge of the interaction between aptamer and protein is integral to the design and development of aptamer-based biosensors. Nanoparticles functionalized with aptamers are commonly used in these kinds of sensors. As such, studies into how the number of aptamers on the nanoparticle surface influence both kinetics and thermodynamics of the binding interaction are required. In this study, aptamers specific for interferon gamma (IFN-γ) were immobilized on the surface of gold nanoparticles (AuNPs), and the effect of surface coverage of aptamer on the binding interaction with its target was investigated using fluorescence spectroscopy. The number of aptamers were adjusted from an average of 9.6 to 258 per particle. The binding isotherm between AuNPs-aptamer conjugate and protein was modeled with the Hill-Langmuir equation, and the determined equilibrium dissociation constant (K'_D) decreased 10-fold when increasing the coverage of aptamer. The kinetics of the reaction as a function of coverage of aptamer were also investigated, including the association rate constant (k_on) and the dissociation rate constant (k_off). The AuNPs-aptamer conjugate with 258 aptamers per particle had the highest k_on, while the k_off was similar for AuNPs-aptamer conjugates with different surface coverages. Therefore, the surface coverage of aptamers on AuNPs affects both the thermodynamics and the kinetics of the binding. The AuNPs-aptamer conjugate with the highest surface coverage is the most favorable in biosensors considering the limit of detection, sensitivity, and response time of the assay. These findings deepen our understanding of the interaction between aptamer and target protein on the particle surface, which is important to both improve the scientific design and increase the application of aptamer-nanoparticle based biosensor.

Kinetics-Driven Drug Design Strategy for Next-Generation Acetylcholinesterase Inhibitors to Clinical Candidate

The acetylcholinesterase (AChE) inhibitors remain key therapeutic drugs for the treatment of Alzheimer's disease (AD). However, the low-safety window limits their maximum therapeutic benefits. Here, a novel kinetics-driven drug design strategy was employed to discover new-generation AChE inhibitors that possess a longer drug-target residence time and exhibit a larger safety window. After detailed investigations, compound 12 was identified as a highly potent, highly selective, orally bioavailable, and brain preferentially distributed AChE inhibitor. Moreover, it significantly ameliorated cognitive impairments in different mouse models with a lower effective dose than donepezil. The X-ray structure of the cocrystal complex provided a precise binding mode between 12 and AChE. Besides, the data from the phase I trials demonstrated that 12 had good safety, tolerance, and pharmacokinetic profiles at all preset doses in healthy volunteers, providing a solid basis for its further investigation in phase II trials for the treatment of AD.

Revealing the cholinergic inhibition mechanism of Alzheimer's by galantamine: a metadynamics simulation study

Galantamine is one of the approved drugs based on the cholinergic hypothesis for the symptomatic treatment of mild to moderate Alzheimer's disease (AD). The etiology of AD is not fully known; however, the reported cholinergic hypothesis suggests the inadequate synthesis of the neurotransmitter acetylcholine (ACh) is responsible for this disease. The crystal structure of galantamine bound human acetylcholinesterase (hAChE) has been reported; however, the inhibition mechanism of hAChE by galantamine is not well understood. A Well-tempered metadynamics (WTMtD) simulation study has been performed with the crystal structure of galantamine bound hAChE. The reported mechanism for the degradation of ACh is suggested through a proton transfer process from a carboxylic group of Glu334 to the hydroxyl group of Ser203, which attacks ACh for the degradation to acetic acid and choline. Such proton transfer process is lowered in the presence of galantamine due to the separation of catalytic triad inside the gorge of AChE as observed with WTMtD. A docking study has been performed to examine the ACh's binding with the catalytic triad of galantamine bound hAChE. The docking results reveal that the approach of ACh to the catalytic triad is interrupted due to the galantamine's presence in the gorge of the enzyme.

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.

Influence of ligand's directional configuration, chrysenes as model compounds, on the binding activity with aryl hydrocarbon receptor

Understanding what and how physico-chemical factors of a ligand configure conditions for ligand-receptor binding is a key to accurate assessment of toxic potencies of environmental pollutants. We investigated influences of the dipole-driven orientation and resulting directional configuration of ligands on receptor binding activities. Using physico-chemical properties calculated by ab initio density functional theory, directional reactivity factors (DRF) were devised as main indicators of toxic potencies, linking molecular ligand-receptor binding to in vitro responses. The directional reactive model was applied to predict variation of aryl hydrocarbon receptor-mediated toxic potencies among homologues of chrysene with structural modifications such as the numbers of constituent benzene rings, methylation and hydroxylation. Results of predictive models were consistent with empirical potencies determined by use of the H4IIE-luc transactivation bioassay. The experiment-free approach based on first principles provides an analytical framework for estimating molecular bioactivity in silico and complements conventional empirical approaches to studying molecular initiating events in adverse outcome pathways.

On Calculating Free Energy Differences Using Ensembles of Transition Paths

The free energy of a process is the fundamental quantity that determines its spontaneity or propensity at a given temperature. In particular, the binding free energy of a drug candidate to its biomolecular target is used as an objective quantity in drug design. Recently, binding kinetics—rates of association (k_on) and dissociation (k_off)—have also demonstrated utility for their ability to predict efficacy and in some cases have been shown to be more predictive than the binding free energy alone. Some methods exist to calculate binding kinetics from molecular simulations, although these are typically more difficult to calculate than the binding affinity as they depend on details of the transition path ensemble. Assessing these rate constants can be difficult, due to uncertainty in the definition of the bound and unbound states, large error bars and the lack of experimental data. As an additional consistency check, rate constants from simulation can be used to calculate free energies (using the log of their ratio) which can then be compared to free energies obtained experimentally or using alchemical free energy perturbation. However, in this calculation it is not straightforward to account for common, practical details such as the finite simulation volume or the particular definition of the “bound” and “unbound” states. Here we derive a set of correction terms that can be applied to calculations of binding free energies using full reactive trajectories. We apply these correction terms to revisit the calculation of binding free energies from rate constants for a host-guest system that was part of a blind prediction challenge, where significant deviations were observed between free energies calculated with rate ratios and those calculated from alchemical perturbation. The correction terms combine to significantly decrease the error with respect to computational benchmarks, from 3.4 to 0.76 kcal/mol. Although these terms were derived with weighted ensemble simulations in mind, some of the correction terms are generally applicable to free energies calculated using physical pathways via methods such as Markov state modeling, metadynamics, milestoning, or umbrella sampling.

Rapid Quantification of Protein-Ligand Binding via 19F NMR Lineshape Analysis

Fluorine incorporation is ideally suited to many NMR techniques, and incorporation of fluorine into proteins and fragment libraries for drug discovery has become increasingly common. Here, we use one-dimensional 19F NMR lineshape analysis to quantify the kinetics and equilibrium thermodynamics for the binding of a fluorine-labeled Src homology 3 (SH3) protein domain to four proline-rich peptides. SH3 domains are one of the largest and most well-characterized families of protein recognition domains and have a multitude of functions in eukaryotic cell signaling. First, we showe that fluorine incorporation into SH3 causes only minor structural changes to both the free and bound states using amide proton temperature coefficients. We then compare the results from lineshape analysis of one-dimensional 19F spectra to those from two-dimensional 1H-15N heteronuclear single quantum coherence spectra. Their agreement demonstrates that one-dimensional 19F lineshape analysis is a robust, low-cost, and fast alternative to traditional heteronuclear single quantum coherence-based experiments. The data show that binding is diffusion limited and indicate that the transition state is highly similar to the free state. We also measured binding as a function of temperature. At equilibrium, binding is enthalpically driven and arises from a highly positive activation enthalpy for association with small entropic contributions. Our results agree with those from studies using different techniques, providing additional evidence for the utility of 19F NMR lineshape analysis, and we anticipate that this analysis will be an effective tool for rapidly characterizing the energetics of protein interactions.

Teaching indicators to unravel the kinetic features of host–guest inclusion complexes

Three new, practically convenient methods are introduced for measuring kinetic parameters of supramolecular host–guest and protein–ligand complexes. Combined with thermodynamic data, this allows for an in-depth of the binding mechanism.

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.

Variational implicit-solvent predictions of the dry-wet transition pathways for ligand-receptor binding and unbinding kinetics

Ligand-receptor binding and unbinding are fundamental biomolecular processes and particularly essential to drug efficacy. Environmental water fluctuations, however, impact the corresponding thermodynamics and kinetics and thereby challenge theoretical descriptions. Here, we devise a holistic, implicit-solvent, multimethod approach to predict the (un)binding kinetics for a generic ligand-pocket model. We use the variational implicit-solvent model (VISM) to calculate the solute-solvent interfacial structures and the corresponding free energies, and combine the VISM with the string method to obtain the minimum energy paths and transition states between the various metastable ("dry" and "wet") hydration states. The resulting dry-wet transition rates are then used in a spatially dependent multistate continuous-time Markov chain Brownian dynamics simulation and the related Fokker-Planck equation calculations of the ligand stochastic motion, providing the mean first-passage times for binding and unbinding. We find the hydration transitions to significantly slow down the binding process, in semiquantitative agreement with existing explicit-water simulations, but significantly accelerate the unbinding process. Moreover, our methods allow the characterization of nonequilibrium hydration states of pocket and ligand during the ligand movement, for which we find substantial memory and hysteresis effects for binding vs. unbinding. Our study thus provides a significant step forward toward efficient, physics-based interpretation and predictions of the complex kinetics in realistic ligand-receptor systems.

Kinetics of Drug Binding and Residence Time

The kinetics of drug binding and unbinding is assuming an increasingly crucial role in the long, costly process of bringing a new medicine to patients. For example, the time a drug spends in contact with its biological target is known as residence time (the inverse of the kinetic constant of the drug-target unbinding, 1/ koff). Recent reports suggest that residence time could predict drug efficacy in vivo, perhaps even more effectively than conventional thermodynamic parameters (free energy, enthalpy, entropy). There are many experimental and computational methods for predicting drug-target residence time at an early stage of drug discovery programs. Here, we review and discuss the methodological approaches to estimating drug binding kinetics and residence time. We first introduce the theoretical background of drug binding kinetics from a physicochemical standpoint. We then analyze the recent literature in the field, starting from the experimental methodologies and applications thereof and moving to theoretical and computational approaches to the kinetics of drug binding and unbinding. We acknowledge the central role of molecular dynamics and related methods, which comprise a great number of the computational methods and applications reviewed here. However, we also consider kinetic Monte Carlo. We conclude with the outlook that drug (un)binding kinetics may soon become a go/no go step in the discovery and development of new medicines.

Multiscale simulations on human Frizzled and Taste2 GPCRs

Recently, molecular dynamics simulations, from all atom and coarse grained to hybrid methods bridging the two scales, have provided exciting functional insights into class F (Frizzled and Taste2) GPCRs (about 40 members in humans). Findings include: (i) The activation of one member of the Frizzled receptors (FZD4) involves a bending of transmembrane helix TM7 far larger than that in class A GPCRs. (ii) The affinity of an anticancer drug targeting another member (Smoothened receptor) decreases in a specific drug-resistant variant, because the mutation ultimately disrupts the binding cavity and affects TM6. (iii) A novel two-state recognition mechanism explains the very large agonist diversity for at least one member of the Taste2 GPCRs, hTAS2R46.

Protein-ligand interaction fingerprints for accurate prediction of dissociation rates of p38 MAPK Type II inhibitors

Binding/unbinding kinetics are key determinants of drug potencies. However, there are still a lot of challenges in predicting kinetic properties during early-stage drug development. In this work, position-restrained molecular dynamics simulations combined with energy decomposition were applied to extract protein-ligand interaction (PLI) fingerprints along the unbinding pathway of 20 p38 mitogen-activated protein kinase (p38 MAPK) Type II inhibitors. The results showed that the electrostatic and/or van der Waals interaction fingerprints at three key positions can be used for accurate prediction of the dissociation rate constants (koff) of p38 MAPK Type II inhibitors. The strategy proposed in this paper can provide not only an efficient method of predicting the dissociation rates of the p38 MAPK Type II inhibitors, but also the atom-level mechanism of enthalpy-driven unbinding process.

Van der Waals Potential in Protein Complexes

Van der Waals forces are determinants of the formation of protein-ligand complexes. Physical models based on the Lennard-Jones potential can estimate van der Waals interactions with considerable accuracy and with a computational complexity that allows its application to molecular docking simulations and virtual screening of large databases of small organic molecules. Several empirical scoring functions used to evaluate protein-ligand interactions approximate van der Waals interactions with the Lennard-Jones potential. In this chapter, we present the main concepts necessary to understand van der Waals interactions relevant to molecular recognition of a ligand by the binding pocket of a protein target. We describe the Lennard-Jones potential and its application to calculate potential energy for an ensemble of structures to highlight the main features related to the importance of this interaction for binding affinity.