The ability to predict the mechanisms and the associated rate constants of protein-ligand unbinding is of great practical importance in drug design. In this work we demonstrate how a recently introduced metadynamics-based approach allows exploration of the unbinding pathways, estimation of the rates, and determination of the rate-limiting steps in the paradigmatic case of the trypsin-benzamidine system. Protein, ligand, and solvent are described with full atomic resolution. Using metadynamics, multiple unbinding trajectories that start with the ligand in the crystallographic binding pose and end with the ligand in the fully solvated state are generated. The unbinding rate k off is computed from the mean residence time of the ligand. Using our previously computed binding affinity we also obtain the binding rate k on . Both rates are in agreement with reported experimental values. We uncover the complex pathways of unbinding trajectories and describe the critical ratelimiting steps with unprecedented detail. Our findings illuminate the role played by the coupling between subtle protein backbone fluctuations and the solvation by water molecules that enter the binding pocket and assist in the breaking of the shielded hydrogen bonds. We expect our approach to be useful in calculating rates for general protein-ligand systems and a valid support for drug design. protein-ligand unbinding | kinetics | enhanced sampling | drug design U nderstanding the thermodynamics and kinetics of proteinligand interactions is of paramount relevance in the early stages of drug discovery (1-3). So far the major emphasis has been placed on predicting the most likely binding pose as determined by the highest binding affinity (4, 5). In contrast, it has not been possible to predict the pathways for unbinding and the associated rates. However, it is by now well-recognized that one of the most pertinent factors for sustained drug efficacy and safety is not just its affinity, but possibly even more so, the mean lifetime of the protein-ligand complex (1-3). The latter property is strictly related to the time during which the ligand remains in the binding site (1, 2), and is typically expressed by its inverse, the dissociation rate k off (2). In principle k off should be amenable to calculations through all-atom molecular dynamics (MD) simulations. These simulations could give detailed and useful insights into the atomic interactions at work during unbinding, especially in the ephemeral but kinetically most relevant transition state ensemble (TSE) (6, 7). Such information is of great value in designing modifications of the ligand that might improve its pharmaceutical properties. However, despite the potential of MD simulations no such calculation has yet been reported. This is a consequence of the limited timescales of MD simulations. Even with the most modern purpose-built supercomputers or massive distributed computing, one can barely reach the timescale of milliseconds (3). Unfortunately most of the reported ligand-protein dissociation times far exceed this timescale (2). These timescales can be reached either by transition path sampling methods (8, 9), quasiclassical approximations (10), by the construction of Markov state models (11, 12), or through carefully designed enhanced sampling methods (8, 13-30) that make accessible the timescale of seconds and beyond in a controlled and accurate way. The enhanced sampling method we use in this work is based on metadynamics (13-15), which has been widely and successfully applied to a variety of systems including complex protein-ligand systems (25-30), and has been rigorously proven to converge to the correct free-energy surface (31, 32) . Recently, we have extended the scope of metadynamics by showing that it can also be used to recover kinetic information (15). Furthermore, we showed that by using an a posteriori statistical analysis (33) one can also establish the reliability of the kinetics thus generated. The use of metadynamics for obtaining kinetic information is still in its infancy, however its usefulness has been tested by us and other groups in a range of systems (15, (33) (34) (35) (36) . In this work, we demonstrate that the scope of the method reported in ref. 15 can be extended to study protein-ligand dissociation pathways and to determine in an accurate way the ligand unbinding rates. We reach well into the hundreds of milliseconds regime and longer, maintaining at the same time full atomic resolution for protein, ligand, and solvent. Specifically, we study the unbinding of the inhibitor benzamidine from trypsin, a serine protease protein (27, 37, 38) using classical force fields (39, 40). Using our acceleration method (15, 33) we are able to harness 21 independent successful unbinding trajectories in which the ligand goes from the bound to the fully unbound state. We find that one of the most distinctive features of the unbinding process is the role played by the water molecules (41, 42). In particular, the solvent promotes unbinding by assisting in the breakage of shielded hydrogen bonds through the formation of water bridge interactions (41) . From the analysis of the unbinding trajectories we find that along the unbinding pathways the ligand rests for times ranging from nanoseconds to milliseconds in a number of intermediate structures. We calculate the rates for all possible transitions Significance A crucial factor for drug efficacy is not just the binding affinity, but also the mean residence time in the binding pocket, usually quantified by its inverse, k off . This is an important parameter that regulates the time during which the drug is active. Whereas the calculation of the binding affinity is by now routine, the calculation of k off has proven more challenging because the timescales involved far exceed the limits of standard molecular dynamics simulation. We propose a metadynamics-based strategy that allows reaching timescales of seconds, and estimate k off along with unbinding pathways and associated dynamical bottlenecks. The protocol is exemplified for trypsin-benzamidine unbinding. This work is a step towards a more effective computer-based drug design. Author contributions: P.T
The metadynamics can be used to pass a ligand through a tunnel in a protein REF and evaluate the thermodynamics and kinetics of the process.
Kinetics of protein–ligand unbinding: Predicting pathways, rates, and rate-limiting steps
Proceedings of the National Academy of Sciences of the United States of America
Proceedings of the National Academy of Sciences of the United States of America
Chemistry Medicine