Pratyush Tiwary

From metadynamics to dynamics.

Metadynamics is a commonly used and successful enhanced sampling method. By the introduction of a history dependent bias which depends on a restricted number of collective variables it can explore complex free energy surfaces characterized by several metastable states separated by large free energy barriers. Here we extend its scope by introducing a simple yet powerful method for calculating the rates of transition between different metastable states. The method does not rely on a previous knowledge of the transition states or reaction coordinates, as long as collective variables are known that can distinguish between the various stable minima in free energy space. We demonstrate that our method recovers the correct escape rates out of these stable states and also preserves the correct sequence of state-to-state transitions, with minimal extra computational effort needed over ordinary metadynamics. We apply the formalism to three different problems and in each case find excellent agreement with the results of long unbiased molecular dynamics runs.

Kinetics of protein–ligand unbinding: Predicting pathways, rates, and rate-limiting steps

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, koff. 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 koff 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 koff 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. 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 koff is computed from the mean residence time of the ligand. Using our previously computed binding affinity we also obtain the binding rate kon. Both rates are in agreement with reported experimental values. We uncover the complex pathways of unbinding trajectories and describe the critical rate-limiting 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.

A Time-Independent Free Energy Estimator for Metadynamics

Metadynamics is a powerful and well-established enhanced sampling method for exploring and quantifying free energy surfaces of complex systems as a function of appropriately chosen variables. In the limit of long simulation time, metadynamics converges to the exact free energy surface plus a time-dependent constant. In this article, we analyze in detail this time-dependent constant. We show an easy way to calculate it, and by explicitly calculating the time dependence of this constant, we are able to derive a time-independent and locally convergent free energy estimator for metadynamics. We also derive an alternate procedure for obtaining the full unbiased distributions of generic operators from biased metadynamics simulations and explicitly test its usefulness.

Enhancing Important Fluctuations: Rare Events and Metadynamics from a Conceptual Viewpoint

Atomistic simulations play a central role in many fields of science. However, their usefulness is often limited by the fact that many systems are characterized by several metastable states separated by high barriers, leading to kinetic bottlenecks. Transitions between metastable states are thus rare events that occur on significantly longer timescales than one can simulate in practice. Numerous enhanced sampling methods have been introduced to alleviate this timescale problem, including methods based on identifying a few crucial order parameters or collective variables and enhancing the sampling of these variables. Metadynamics is one such method that has proven successful in a great variety of fields. Here we review the conceptual and theoretical foundations of metadynamics. As demonstrated, metadynamics is not just a practical tool but can also be considered an important development in the theory of statistical mechanics.

Assessing the reliability of the dynamics reconstructed from metadynamics

Sampling a molecular process characterized by an activation free energy significantly larger than kBT is a well-known challenge in molecular dynamics simulations. In a recent work [Tiwary and Parrinello, Phys. Rev. Lett. 2013, 111, 230602], we have demonstrated that the transition times of activated molecular transformations can be computed from well-tempered metadynamics provided that no bias is deposited in the transition state region and that the set of collective variables chosen to enhance sampling does not display hysteresis. Ensuring though that these two criteria are met may not always be simple. Here we build on the fact that the times of escape from a long-lived metastable state obey Poisson statistics. This allows us to identify quantitative measures of trustworthiness of our calculation. We test our method on a few paradigmatic examples.

Spectral gap optimization of order parameters for sampling complex molecular systems

Significance Molecular-dynamics (MD) simulations have become a versatile tool for exploration of complex molecular systems. However, they are limited in the timescales that can be reached. Thus, over the years, a suite of enhanced-sampling algorithms have been proposed that assist MD to transcend the timescale limitation, with diverse applications across physical and life sciences. A continuing grand challenge in the success of many such sampling methods pertains to a judicious choice of order parameters. In this work, we propose a new method for designing order parameters that minimizes the role played by human intuition and makes the progress significantly more automated than before. We expect this algorithm to be of great use in furthering the success of enhanced sampling. In modern-day simulations of many-body systems, much of the computational complexity is shifted to the identification of slowly changing molecular order parameters called collective variables (CVs) or reaction coordinates. A vast array of enhanced-sampling methods are based on the identification and biasing of these low-dimensional order parameters, whose fluctuations are important in driving rare events of interest. Here, we describe a new algorithm for finding optimal low-dimensional CVs for use in enhanced-sampling biasing methods like umbrella sampling, metadynamics, and related methods, when limited prior static and dynamic information is known about the system, and a much larger set of candidate CVs is specified. The algorithm involves estimating the best combination of these candidate CVs, as quantified by a maximum path entropy estimate of the spectral gap for dynamics viewed as a function of that CV. The algorithm is called spectral gap optimization of order parameters (SGOOP). Through multiple practical examples, we show how this postprocessing procedure can lead to optimization of CV and several orders of magnitude improvement in the convergence of the free energy calculated through metadynamics, essentially giving the ability to extract useful information even from unsuccessful metadynamics runs.

Prediction of Protein–Ligand Binding Poses via a Combination of Induced Fit Docking and Metadynamics Simulations

Ligand docking is a widely used tool for lead discovery and binding mode prediction based drug discovery. The greatest challenges in docking occur when the receptor significantly reorganizes upon small molecule binding, thereby requiring an induced fit docking (IFD) approach in which the receptor is allowed to move in order to bind to the ligand optimally. IFD methods have had some success but suffer from a lack of reliability. Complementing IFD with all-atom molecular dynamics (MD) is a straightforward solution in principle but not in practice due to the severe time scale limitations of MD. Here we introduce a metadynamics plus IFD strategy for accurate and reliable prediction of the structures of protein-ligand complexes at a practically useful computational cost. Our strategy allows treating this problem in full atomistic detail and in a computationally efficient manner and enhances the predictive power of IFD methods. We significantly increase the accuracy of the underlying IFD protocol across a large data set comprising 42 different ligand-receptor systems. We expect this approach to be of significant value in computationally driven drug design.

Role of water and steric constraints in the kinetics of cavity–ligand unbinding

Significance The unbinding of ligand–substrate systems in molecular water is a problem of great theoretical and practical interest. To understand the dynamical nature of the unbinding, it is desirable to use atomistic techniques like molecular dynamics (MD). However, the associated timescales are typically too long for MD to be applicable. In this work, we apply a recent metadynamics scheme that allows the use of MD to get information on thermodynamics and kinetics of systems that are virtually impossible to treat with unbiased MD. We calculate unbinding pathways and timescales, and show that solvent molecules play a crucial role in cooperation with steric effects. The approaches used here demonstrate a broadly applicable methodology for studying ligand unbinding. A key factor influencing a drug’s efficacy is its residence time in the binding pocket of the host protein. Using atomistic computer simulation to predict this residence time and the associated dissociation process is a desirable but extremely difficult task due to the long timescales involved. This gets further complicated by the presence of biophysical factors such as steric and solvation effects. In this work, we perform molecular dynamics (MD) simulations of the unbinding of a popular prototypical hydrophobic cavity–ligand system using a metadynamics-based approach that allows direct assessment of kinetic pathways and parameters. When constrained to move in an axial manner, the unbinding time is found to be on the order of 4,000 s. In accordance with previous studies, we find that the cavity must pass through a region of sharp wetting transition manifested by sudden and high fluctuations in solvent density. When we remove the steric constraints on ligand, the unbinding happens predominantly by an alternate pathway, where the unbinding becomes 20 times faster, and the sharp wetting transition instead becomes continuous. We validate the unbinding timescales from metadynamics through a Poisson analysis, and by comparison through detailed balance to binding timescale estimates from unbiased MD. This work demonstrates that enhanced sampling can be used to perform explicit solvent MD studies at timescales previously unattainable, to our knowledge, obtaining direct and reliable pictures of the underlying physiochemical factors including free energies and rate constants.

Unbinding Kinetics of a p38 MAP Kinase Type II Inhibitor from Metadynamics Simulations

Understanding the structural and energetic requisites of ligand binding toward its molecular target is of paramount relevance in drug design. In recent years, atomistic free energy calculations have proven to be a valid tool to complement experiments in characterizing the thermodynamic and kinetic properties of protein/ligand interaction. Here, we investigate, through a recently developed metadynamics-based protocol, the unbinding mechanism of an inhibitor of the pharmacologically relevant target p38 MAP kinase. We provide a thorough description of the ligand unbinding pathway identifying the most stable binding mode and other thermodynamically relevant poses. From our simulations, we estimated the unbinding rate as koff = 0.020 ± 0.011 s-1. This is in good agreement with the experimental value (koff = 0.14 s-1). Next, we developed a Markov state model that allowed identifying the rate-limiting step of the ligand unbinding process. Our calculations further show that the solvation of the ligand and that of the active site play crucial roles in the unbinding process. This study paves the way to investigations on the unbinding dynamics of more complex p38 inhibitors and other pharmacologically relevant inhibitors in general, demonstrating that metadynamics can be a powerful tool in designing new drugs with engineered binding/unbinding kinetics.

Hybrid deterministic and stochastic approach for efficient atomistic simulations at long time scales

We propose a hybrid deterministic and stochastic approach to achieve extended time scales in atomistic simulations that combines the strengths of molecular dynamics (MD) and Monte Carlo (MC) simulations in an easy-to-implement way. The method exploits the rare event nature of the dynamics similar to most current accelerated MD approaches but goes beyond them by providing, without any further computational overhead, (a) rapid thermalization between infrequent events, thereby minimizing spurious correlations, and (b) control over accuracy of time-scale correction, while still providing similar or higher boosts in computational efficiency. We present two applications of the method: (a) Vacancy-mediated diffusion in Fe yields correct diffusivities over a wide range of temperatures and (b) source-controlled plasticity and deformation behavior in Au nanopillars at realistic strain rates (10^4/s and lower), with excellent agreement with previous theoretical predictions and in situ high-resolution transmission electron microscopy observations. The method gives several orders-of-magnitude improvements in computational efficiency relative to standard MD and good scalability with the size of the system.