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

Multiscale modeling of unfolding and bond dissociation of rubredoxin metalloprotein

Mechanical properties of proteins that have a crucial effect on their operation. This study used a molecular dynamics simulation package to investigate rubredoxin unfolding on the atomic scale. Different simulation techniques were applied, and due to the dissociation of covalent/hydrogen bonds, this protein demonstrates several intermediate states in force-extension behavior. A conceptual model based on the cohesive finite element method was developed to consider the intermediate damages that occur during unfolding. This model is based on force-displacement curves derived from molecular dynamics results. The proposed conceptual model is designed to accurately identify bond rupture points and determine the associated forces. This is achieved by conducting a thorough comparison between molecular dynamics and cohesive finite element results. The utilization of a viscoelastic cohesive zone model allows for the consideration of loading rate effects. This rate-dependent model can be further developed and integrated into the multiscale modeling of large assemblies of metalloproteins, providing a comprehensive understanding of mechanical behavior while maintaining a reduced computational cost.

Computational insights into the conformational transition of STING: Mechanistic, energetic considerations, and the influence of crucial mutations

STING (stimulator of interferon genes) is a crucial protein in the innate immune system's response to viral and bacterial infections. In this study, we investigated the mechanistic and energetic mechanism of the conformational transition process of STING activated by cGAMP binding. We found that the STING connector region undergoes an energetically unfavorable rotation during this process, which is compensated by the favorable interaction between cGAMP and the STING ligand binding domain. We further studied several disease-causing mutations and found that the V155 M mutation facilitates a smoother transition in the STING connector region. However, the V147L mutation exhibits unfavorable conformational transition energy, suggesting it may hinder STING activation pathway that relies on connector region rotation. Despite being labeled as hyperactive, the widespread prevalence of V147L/V147I mutations across species implies a neutral character, indicating complexity in its role. Overall, our analysis deepens the understanding of STING activation within the connector region, and targeting this region with compounds may provide an alternative approach to interfering with STING's function.

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.

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

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

Worth the Weight: Sub-Pocket EXplorer (SubPEx), a Weighted Ensemble Method to Enhance Binding-Pocket Conformational Sampling

Structure-based virtual screening (VS) is an effective method for identifying potential small-molecule ligands, but traditional VS approaches consider only a single binding-pocket conformation. Consequently, they struggle to identify ligands that bind to alternate conformations. Ensemble docking helps address this issue by incorporating multiple conformations into the docking process, but it depends on methods that can thoroughly explore pocket flexibility. We here introduce Sub-Pocket EXplorer (SubPEx), an approach that uses weighted ensemble (WE) path sampling to accelerate binding-pocket sampling. As proof of principle, we apply SubPEx to three proteins relevant to drug discovery: heat shock protein 90, influenza neuraminidase, and yeast hexokinase 2. SubPEx is available free of charge without registration under the terms of the open-source MIT license: http://durrantlab.com/subpex/.

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.

Dissociation Rate Calculation via Constant-Force Steered Molecular Dynamics Simulation

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

Worth the weight: Sub-Pocket EXplorer (SubPEx), a weighted-ensemble method to enhance binding-pocket conformational sampling

Structure-based virtual screening (VS) is an effective method for identifying potential small-molecule ligands, but traditional VS approaches consider only a single binding-pocket conformation. Consequently, they struggle to identify ligands that bind to alternate conformations. Ensemble docking helps address this issue by incorporating multiple conformations into the docking process, but it depends on methods that can thoroughly explore pocket flexibility. We here introduce Sub-Pocket EXplorer (SubPEx), an approach that uses weighted ensemble (WE) path sampling to accelerate binding-pocket sampling. As proof of principle, we apply SubPEx to three proteins relevant to drug discovery: heat shock protein 90, influenza neuraminidase, and yeast hexokinase 2. SubPEx is available free of charge without registration under the terms of the open-source MIT license: http://durrantlab.com/subpex/.

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

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

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

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

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.

Umbrella Sampling Simulations Measure Switch Peptide Binding and Hydrophobic Patch Opening Free Energies in Cardiac Troponin

The cardiac troponin (cTn) complex is an important regulatory protein in heart contraction. Upon binding of Ca2+, cTn undergoes a conformational shift that allows the troponin I switch peptide (cTnISP) to be released from the actin filament and bind to the troponin C hydrophobic patch (cTnCHP). Mutations and modifications to this complex can change its sensitivity to Ca2+ and alter the energetics of the transition from the Ca2+-unbound, cTnISP-unbound form to the Ca2+-bound, cTnISP-bound form. We utilized targeted molecular dynamics (TMD) to obtain a trajectory of this transition pathway, followed by umbrella sampling to estimate the free energy associated with the cTnISP-cTnCHP binding and the cTnCHP opening events for wild-type (WT) cTn. We were able to reproduce experimental values for the cTnISP-cTnCHP binding event and obtain cTnCHP opening free energies in agreement with previous computational measurements of smaller cTnC systems. This excellent agreement for WT cTn demonstrated the strength of computational methods in studying the dynamics and energetics of the cTn complex. We then introduced mutations to the cTn complex that cause cardiomyopathy or alter its Ca2+ sensitivity and observed a general decrease in the free energy of opening the cTnCHP. For these same mutations, we observed no general trend in the effect on the cTnISP-cTnCHP binding event. Our method sets the stage for future computational studies on this system that predict the consequences of yet uncharacterized mutations on cTn dynamics and energetics.

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

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

Mebendazole's Conformational Space and Its Predicted Binding to Human Heat-Shock Protein 90

Recent experimental evidence suggests that mebendazole, a popular antiparasitic drug, binds to heat shock protein 90 (Hsp90) and inhibits acute myeloid leukemia cell growth. In this study we use quantum mechanics (QM), molecular similarity, and molecular dynamics (MD) calculations to predict possible binding poses of mebendazole to the adenosine triphosphate (ATP) binding site of Hsp90. Extensive conformational searches and minimization of the five mebendazole tautomers using the MP2/aug-cc-pVTZ theory level resulted in 152 minima. Mebendazole-Hsp90 complex models were subsequently created using the QM optimized conformations and protein coordinates obtained from experimental crystal structures that were chosen through similarity calculations. Nine different poses were identified from a total of 600 ns of explicit solvent, all-atom MD simulations using two different force fields. All simulations support the hypothesis that mebendazole is able to bind to the ATP binding site of Hsp90.

Enhanced-Sampling Simulations for the Estimation of Ligand Binding Kinetics: Current Status and Perspective

The dissociation rate (k off) associated with ligand unbinding events from proteins is a parameter of fundamental importance in drug design. Here we review recent major advancements in molecular simulation methodologies for the prediction of k off. Next, we discuss the impact of the potential energy function models on the accuracy of calculated k off values. Finally, we provide a perspective from high-performance computing and machine learning which might help improve such predictions.

Determination of Multidirectional Pathways for Ligand Release from the Receptor: A New Approach Based on Differential Evolution

Steered molecular dynamics (SMD) simulation is a powerful method in computer-aided drug design as it can be used to access the relative binding affinity with high precision but with low computational cost. The success of SMD depends on the choice of the direction along which the ligand is pulled from the receptor-binding site. In most simulations, the unidirectional pathway was used, but in some cases, this choice resulted in the ligand colliding with the complex surface of the exit tunnel. To overcome this difficulty, several variants of SMD with multidirectional pulling have been proposed, but they are not completely devoid of disadvantages. Here, we have proposed to determine the direction of pulling with a simple scoring function that minimizes the receptor–ligand interaction, and an optimization algorithm called differential evolution is used for energy minimization. The effectiveness of our protocol was demonstrated by finding expulsion pathways of Huperzine A and camphor from the binding site of Torpedo California acetylcholinesterase and P450cam proteins, respectively, and comparing them with the previous results obtained using memetic sampling and random acceleration molecular dynamics. In addition, by applying this protocol to a set of ligands bound with LSD1 (lysine specific demethylase 1), we obtained a much higher correlation between the work of pulling force and experimental data on the inhibition constant IC50 compared to that obtained using the unidirectional approach based on minimal steric hindrance.

Unbinding ligands from SARS-CoV-2 Mpro via umbrella sampling simulations

The umbrella sampling (US) simulation is demonstrated to be an efficient approach for determining the unbinding pathway and binding affinity to the SARS-CoV-2 Mpro of small molecule inhibitors. The accuracy of US is in the same range as the linear interaction energy (LIE) and fast pulling of ligand (FPL) methods. In detail, the correlation coefficient between US and experiments does not differ from FPL and is slightly smaller than LIE. The root mean square error of US simulations is smaller than that of LIE. Moreover, US is better than FPL and poorer than LIE in classifying SARS-CoV-2 Mpro inhibitors owing to the reciever operating characteristic-area under the curve analysis. Furthermore, the US simulations also provide detailed insights on unbinding pathways of ligands from the binding cleft of SARS-CoV-2 Mpro. The residues Cys44, Thr45, Ser46, Leu141, Asn142, Gly143, Glu166, Leu167, Pro168, Ala191, Gln192 and Ala193 probably play an important role in the ligand dissociation. Therefore, substitutions at these points may change the mechanism of binding of inhibitors to SARS-CoV-2 Mpro.

Predicting Ion Channel Conductance via Dissipation-Corrected Targeted Molecular Dynamics and Langevin Equation Simulations

Ion channels are important proteins for physiological information transfer and functional control. To predict the microscopic origins of their voltage-conductance characteristics, here we applied dissipation-corrected targeted molecular dynamics in combination with Langevin equation simulations to potassium diffusion through the gramicidin A channel as a test system. Performing a nonequilibrium principal component analysis on backbone dihedral angles, we find coupled protein-ion dynamics to occur during ion transfer. The dissipation-corrected free energy profiles correspond well to predictions from other biased simulation methods. The incorporation of an external electric field in Langevin simulations enables the prediction of macroscopic observables in the form of I-V characteristics.

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.

Molecular design of anticancer drugs from marine fungi derivatives

Heat shock protein 90 (Hsp90) is one of the most potential targets in cancer therapy. We have demonstrated using a combination of molecular docking and fast pulling of ligand (FPL) simulations that marine fungi derivatives can be possible inhibitors, preventing the biological activity of Hsp90. The computational approaches were validated and compared with previous experiments. Based on the benchmark of available inhibitors of Hsp90, the GOLD docking package using the ChemPLP scoring function was found to be superior over both Autodock Vina and Autodock4 in the preliminary estimation of the ligand-binding affinity and binding pose with the Pearson correlation, R = -0.62. Moreover, FPL calculations were also indicated as a suitable approach to refine docking simulations with a correlation coefficient with the experimental data of R = -0.81. Therefore, the binding affinity of marine fungi derivatives to Hsp90 was evaluated. Docking and FPL calculations suggest that five compounds including 23, 40, 46, 48, and 52 are highly potent inhibitors for Hsp90. The obtained results enhance cancer therapy research.

Accurate absolute free energies for ligand-protein binding based on non-equilibrium approaches

The accurate calculation of the binding free energy for arbitrary ligand-protein pairs is a considerable challenge in computer-aided drug discovery. Recently, it has been demonstrated that current state-of-the-art molecular dynamics (MD) based methods are capable of making highly accurate predictions. Conventional MD-based approaches rely on the first principles of statistical mechanics and assume equilibrium sampling of the phase space. In the current work we demonstrate that accurate absolute binding free energies (ABFE) can also be obtained via theoretically rigorous non-equilibrium approaches. Our investigation of ligands binding to bromodomains and T4 lysozyme reveals that both equilibrium and non-equilibrium approaches converge to the same results. The non-equilibrium approach achieves the same level of accuracy and convergence as an equilibrium free energy perturbation (FEP) method enhanced by Hamiltonian replica exchange. We also compare uni- and bi-directional non-equilibrium approaches and demonstrate that considering the work distributions from both forward and reverse directions provides substantial accuracy gains. In summary, non-equilibrium ABFE calculations are shown to yield reliable and well-converged estimates of protein-ligand binding affinity.

Hierarchical dynamics in allostery following ATP hydrolysis monitored by single molecule FRET measurements and MD simulations

We report on a study that combines advanced fluorescence methods with molecular dynamics (MD) simulations to cover timescales from nanoseconds to milliseconds for a large protein. This allows us to delineate how ATP hydrolysis in a protein causes allosteric changes at a distant protein binding site, using the chaperone Hsp90 as test system. The allosteric process occurs via hierarchical dynamics involving timescales from nano- to milliseconds and length scales from Ångstroms to several nanometers. We find that hydrolysis of one ATP is coupled to a conformational change of Arg380, which in turn passes structural information via the large M-domain α-helix to the whole protein. The resulting structural asymmetry in Hsp90 leads to the collapse of a central folding substrate binding site, causing the formation of a novel collapsed state (closed state B) that we characterise structurally. We presume that similar hierarchical mechanisms are fundamental for information transfer induced by ATP hydrolysis through many other proteins.

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.

Estimating the ligand-binding affinity via λ-dependent umbrella sampling simulations

The umbrella sampling (US) approach has been demonstrated to be a very efficient method for estimating the ligand-binding affinity. However, most of the calculated values overestimate experimental ones that are probably caused by the inaccurate representation of the interaction between the ligand and the surrounding molecules. The issue can be resolved via the implementation aspects of λ-alteration simulation into the US approach, which we call the λ-dependent umbrella sampling (λUS) scheme. In particular, the electrostatic and van der Waals interactions were simultaneously changed by using the coupling parameter λ during λUS simulations. The mean value of obtained results, ∆ G US λ = 0.20 = - 11.59 ± 1.51 kcal mol^-1 , is in good fitting to the mean value of respective experiments, ∆G_EXP =  - 11.26 ± 0.89 kcal mol^-1 . Moreover, the correlation between the proposed approach and experiment is quite good with a value of R US λ = 0.20 = 0.82 ± 0.10 . The λUS scheme significantly enhances the calculated accuracy since the RMSE of the proposed scheme is smaller than traditional US simulations, RMSE US λ = 0.20 = 2.99 ± 0.82 kcal mol^-1 versus RMSE US λ = 0.00 = 5.48 ± 0.81 kcal mol^-1 . Furthermore, the precision is increased since the computed error via λUS approach, δ US λ = 0.20 = 1.51 kcal mol^-1 , was smaller than those of the US simulation, δ US λ = 0.00 = 1.78 kcal mol^-1 . Overall, the proposed approach perhaps provides an efficient way to accurately and precisely estimate the ligand-binding free energy.

Real-time observation of ligand-induced allosteric transitions in a PDZ domain

While allostery is of paramount importance for protein regulation, the underlying dynamical process of ligand (un)binding at one site, resulting time evolution of the protein structure, and change of the binding affinity at a remote site are not well understood. Here the ligand-induced conformational transition in a widely studied model system of allostery, the PDZ2 domain, is investigated by transient infrared spectroscopy accompanied by molecular dynamics simulations. To this end, an azobenzene-derived photoswitch is linked to a peptide ligand in a way that its binding affinity to the PDZ2 domain changes upon switching, thus initiating an allosteric transition in the PDZ2 domain protein. The subsequent response of the protein, covering four decades of time, ranging from ∼1 ns to ∼μs, can be rationalized by a remodeling of its rugged free-energy landscape, with very subtle shifts in the populations of a small number of structurally well-defined states. It is proposed that structurally and dynamically driven allostery, often discussed as limiting scenarios of allosteric communication, actually go hand-in-hand, allowing the protein to adapt its free-energy landscape to incoming signals.

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.

Exhaustive search of the configurational space of heat-shock protein 90 with its inhibitor by multicanonical molecular dynamics based dynamic docking

Multicanonical molecular dynamics based dynamic docking was used to exhaustively search the configurational space of an inhibitor binding to the N-terminal domain of heat-shock protein 90 (Hsp90). The obtained structures at 300 K cover a wide structural ensemble, with the top two clusters ranked by their free energy coinciding with the native binding site. The representative structure of the most stable cluster reproduced the experimental binding configuration, but an interesting conformational change in Hsp90 could be observed. The combined effects of solvation and ligand binding shift the equilibrium from a preferred loop-in conformation in the unbound state to an α-helical one in the bound state for the flexible lid region of Hsp90. Thus, our dynamic docking method is effective at predicting the native binding site while exhaustively sampling a wide configurational space, modulating the protein structure upon binding.

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.