Kinetics of Ligand-Protein Dissociation from All-Atom Simulations: Are

Subscriber access provided by AUSTRALIAN NATIONAL UNIV

Perspective

Kinetics of Ligand-Protein Dissociation from All-Atom Simulations: Are We There Yet? Joao Marcelo Lamim Ribeiro, Sun-Ting Tsai, Debabrata Pramanik, Yihang Wang, and Pratyush Tiwary Biochemistry, Just Accepted Manuscript • DOI: 10.1021/acs.biochem.8b00977 • Publication Date (Web): 14 Dec 2018 Downloaded from http://pubs.acs.org on December 15, 2018

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 26 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Biochemistry

Kinetics of Ligand-Protein Dissociation from All-Atom Simulations: Are We There Yet? João Marcelo Lamim Ribeiro,‡ Sun-Ting Tsai,§ Debabrata Pramanik,‡ Yihang Wang,k and Pratyush Tiwary∗,‡ ‡Department of Chemistry and Biochemistry and Institute for Physical Science and Technology, University of Maryland, College Park 20742, USA. §Department of Physics and Institute for Physical Science and Technology, University of Maryland, College Park 20742, USA. kBiophysics Program and Institute for Physical Science and Technology, University of Maryland, College Park 20742, USA. E-mail: [email protected] Phone: (301) 405-2148

1 ACS Paragon Plus Environment

Biochemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Abstract Large parallel gains in the development of both computational resources as well as sampling methods have now made it possible to simulate dissociation events in ligand-protein complexes with all–atom resolution. Such encouraging progress, together with the inherent spatiotemporal resolution associated with molecular simulations, has left their use for investigating dissociation processes brimming with potential, both in rational drug design, where it can be an invaluable tool for determining the mechanistic driving forces behind dissociation rate constants, as well as in force-field development, where it can provide a catalog of transient molecular structures on which to refine force-fields. Although much progress has been made in making force-fields more accurate, reducing their error for transient structures along a transition path could yet prove to be a critical development helping to make kinetic predictions much more accurate. In what follows we will provide a state-of-the-art compilation of the enhanced sampling methods based on molecular dynamics (MD) simulations used to investigate the kinetics and mechanisms of ligand-protein dissociation processes. Due to the timescales of such processes being slower than what is accessible using straightforward MD simulations, several ingenious schemes are being devised at a rapid rate to overcome this obstacle. Here we provide an up-to-date compendium of such methods and their achievements/shortcomings in extracting mechanistic insight into ligand-protein dissociation. We conclude with a critical and provocative appraisal attempting to answer the title of this review.

1

INTRODUCTION

Recent decades have seen enormous gains in the development of computational hardware both in generalpurpose computing resources such as graphical processing units (GPUs) and in molecular simulation-specific hardware. 1,2 Together these have pushed the timescales that are now accessible into a regime where important biochemical processes can be studied using all-atom molecular dynamics (MD) simulations. As a result, a widening catalog of biochemical processes have been studied at the atomistic level, 1–16 including the dissociation of ligand-macromolecule complexes. 17–20 These dissociation processes are relevant in various contexts, including developing a fundamental understanding for the chemical basis of life processes as well as in the rational design of drugs. 1,21 With respect to the latter, for instance, MD simulations have been found to be ever more accurate in predicting the relative binding free energies ∆∆Gb between families of ligands and/or receptors (see Table 1 for definitions of ∆∆Gb and other terms). 22–24 In particular, the determination of accurate mechanisms for ligand-protein dissociation processes, which can happen over a timespan of microseconds to several hours, 25–30 remains a classic unresolved biochemical problem. There do exist, of course, experimental techniques capable of measuring the overall rate constants


Page 2 of 26

Page 3 of 26 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Biochemistry

for ligand-protein dissociations, 25 but it nonetheless continues to be difficult for such techniques to glean an atomic residue-by-residue level understanding of the dissociative mechanism. Much of the complication – as well as richness – of the problem arise from the inherent ever-fluctuating structures, including those of the protein, ligand and the solvent medium, which are then too transient in nature to be detected in experiments. In addition, as has been highlighted in recent work, 25,26,31,32 an important aspect of ligandprotein interactions is the ligand’s residence time in the target protein, which in a number of ligand-protein complexes correlates more with mechanism and is a much stronger predictor of eventual function than binding affinities, which correlate more with structure. Often the residence time is quantified through its reciprocal, the dissociation rate constant kof f (see Table 1). MD simulations could, in principle, be ideal vehicles for investigating these dissociation processes with all-atom resolution since it tracks, per construction, the changing atomic positions as a function of time. It is rather unfortunate, then, that the relevant timescales lie far past what can be simulated even with the best computing resources if one was to use straightforward MD simulations. This is because complex biochemical processes often involve multiple metastable states in the presence of high barriers, leading to timescales that are long, while the fast vibrational motions constrain the integration time step in a MD simulation to femtosecond (fs) values. In order to tackle problems plagued with such extreme rare events, several specialized methods have been proposed and used to deal with the long timescales required to observe the desired rare event in a simulation. 33–35 The idea behind this review is to perform a critical examination of where we stand in the context of using all-atom MD simulations for investigating the dissociation process of ligand-protein complexes, which we take to include not just the calculation of dissociation rate constants but also the characterization of metastable and transition states relevant to the process, as well as various other related biochemical factors. There have been various studies in the past pointing to the fact that using dissociation kinetics as the basis of rational drug design can lead to a more efficient optimization. 25,26 Recent work from Georgi et al. 36 , meanwhile, has also suggested that low kof f values might be better at predicting clinical success than a large association rate constant, kon . That this seems to be the case has the added benefit that optimizing the dissociation kinetics of a potential drug is a problem well-suited to the medicinal chemist, who is capable of tailoring the properties of the potential drug in order to affect its rate of dissociation. 25,26 The binding of a ligand to a protein, on the other hand, is often a diffusion-limited process that is much harder to control. In this context, MD simulations could be an invaluable tool for aiding rational drug design, since it necessarily acts as in silico microscopes revealing the mechanisms giving rise to the macroscopic rate constant kof f often determined from experiments. It can provide a full catalog of intermediate and transition states that are hard to characterize in experiments due to their transient nature, revealing specific molecular attractions,


Biochemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

steric hindrances, macromolecular motions and water/solvent effects that together give rise to the measured rate constant. Hence our specific focus on, for the most part, ligand-protein dissociations. Here we provide a state-of-the-art overview for a number of methods that have been proposed to overcome the timescale limitation common to ligand-protein dissociations, making sure to highlight their main principles as well as what we believe to be their strengths and weaknesses. We, in addition, provide several specific case-studies demonstrating the use of MD methods for investigating the kinetics of ligand-protein dissociations. Our conclusion, which is half parts optimistic and half parts concerned, will then summarize the field as well as the challenges still to be surmounted. It is our belief that this review will serve as a definitive manual for a broad class of researchers, both those who wish to develop the methods further as well as those wishing to form an opinion regarding which method to use for a specific calculation/application. Table 1: Definition of important terms. Term ∆Gb ∆∆Gb kof f kon x RC

2 2.1

Definition Absolute binding free energy Relative binding free energy difference between two systems Ligand-macromolecule complex dissociation rate constant association rate constant equaling e−β∆Gb kof f Molecular configuration of ligandmacromolecule complex Reaction coordinate

METHODS Unbiased methods

As long as a sufficient number of dissociation events are observed in an unbiased MD simulation, the dissociation rate constants and mechanisms are straightforward to obtain through simple counting. Of course the timescale problem limits this approach to just the ligand-protein complexes that are fast unbinders, although with the introduction of the Anton supercomputer, 37,38 which was designed for and can run MD simulations two orders of magnitude faster than conventional supercomputers, the number of association processes that can be studied via long, brute-force unbiased simulations has expanded. 27,39–41 It remains quite difficult, however, to use long unbiased MD for ligand dissociation processes due to their much slower rates in general, 25–27 and there are but a few dissociations studied via this brute-force approach. For example, Huang and Caflisch 42 have investigated the dissociation of small ligand molecules from the FKBP protein binding sites through unbiased MD simulations in explicit water. Pan et al. 43 ran long simulations on Anton for the dissociation of several small ligands from the FKBP protein. In the remainder of this section, we summarize 4 ACS Paragon Plus Environment

Page 4 of 26

Page 5 of 26 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Biochemistry

a number of enhanced sampling methods that do not require the explicit biasing of an MD simulation. These methods, in one form or another, are designed to patch several unbiased simulations to form a picture of the overall dissociation process. In the authors’ view, it is non-trivial to assess if such approaches can lead to truly unbiased results. (see Sec. 4). 2.1.1

Markov State Models

Markov State Models (MSMs) provide a powerful framework with which to model the kinetics of biochemical processes. One of its appealing features is that it models the kinetics through much shorter MD simulations run in parallel, 44–46 which in principle helps to ameliorate the difficulties with calculating kinetic estimates of rare events with statistical confidence. Keep in mind that these short parallel MD simulations are initialized from different states, meaning that in the data sparse regime one must first generate at least an approximate set of states from simulations that can be adapted and refined in future steps. 44 MSMs, in addition, provide a natural framework with which to extract a simplified and understandable mechanistic model from a large volume of MD data, 44 albeit one should keep in mind that a trade-off exists between how simple versus how accurate the MSM will be. 46 In essence, an MSM is a parametric model, akin to a master equation, 44,47 that aims to describe the long-time behavior of a dynamical system in terms of stochastic transitions between a set of discrete states in configuration space. 44,46,47 Setting up the MSM will then involve parameterizing it, which will involve two crucial steps: (a) Defining the states relevant to describing the process of interest (b) Calculating transition probabilities for interconverting between states In practice, there is also a third critical parameter, called the lag time, which serves to ensure that the process can be modeled as a Markov chain, meaning that the transitions into future states will depend on just the current state, regardless of its past. So long as the lag time is large enough to allow thermalization, the dependence on past states is erased, although, of course, with too large a lag time the running of MD simulations will become inefficient. It is possible, with a proper choice of lag time, to estimate the transition matrix through “stitching” together several short MD simulations that are run in parallel. In general, the lag time is problem specific, but it is useful to note here that in several studies it is ∼10 ns, which means that it is often possible to build an estimate of the transition matrix via short parallel MD simulations. 44,48 Discretizing configuration space into distinct states is both a crucial and non-trivial step in building MSMs. The clustering of MD simulation data was often the standard approach towards state decomposition. 44,47 Work from Noé and Nuske 49 , however, introduced a variational approach for approximating the eigenvalues and eigenfunctions of the transition matrix, which has had a profound effect in automating MSM state decomposition. 47,49 Later, Pérez-Hernández et al. 50 showed that the time-lagged independent


Biochemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

component analysis (tICA) provides an optimal solution to the variational problem. An important point to mention is that although MSMs make it possible to fuse together several short trajectories run in parallel, it is nonetheless crucial that the model includes MD data corresponding to the rare event transition itself. It is possible, of course, that a transition is much too slow to be observed in practice during MD simulations, which can limit their practical usefulness. An extreme rare event process such as the dissociation of a ligand-protein complex whose residence time is minutes long falls outside the timescale regime for which straightforward MSM building is relevant. In such cases the use of biased simulation methods in order to observe the extreme rare event seems to be a plausible approach for constructing MSMs, although this represents a crucial open problem and remains an area of great active interest. 51,52 2.1.2

Weighted ensemble

Weighted ensemble (WE) uses the idea of resampling to speed-up the observation of a desired rare event. 53–59 In order to initiate a WE calculation, the first step is to launch an ensemble of trajectories drawn from some initial distribution in configuration space, P = P (x), which might be time-dependent or independent. 55 Once a short amount of time has passed the trajectories are re-evaluated for their progress along the binned configuration space, meaning the space of molecular configurations subdivided into a number of discrete cells, with certain configurations copied and others eliminated based on a rigorous adjustment of their weights made to ensure this new sample of configurations will conserve the current value of P (x). 55 In other words, during the resampling step that follows the short MD simulations, the set of molecular configurations is changed such that the new set remains representative of the current value of P (x) while also increasing the number of configurations visiting the regions of configuration space with low probabilities that must be traversed in order to sample the desired rare event. Using the illustrative case that the binning of configuration space is described via a binned reaction coordinate (RC), the WE method aims to increase the number of configurations further along the RC. The WE method will then re-initiate the ensemble of trajectories from the current configurations and proceed to iterate through the resampling and MD simulation steps until the desired rare event is observed. 55 One important point to consider is that the binning of configuration space in WE simulations will be specific to the problem at hand. 55 This predefinition of bins might in turn affect the ideal values of other parameters such as the time lag between resampling steps. 55 It is also worth mentioning, however, that the bins do not have to be fixed during a WE simulation. 54–56 For instance, there exist adaptive WE simulation schemes allowing the center and the size of bins to change. 54 The bins learn to maximize the distance between each other within the space that has been sampled and will spread out as more space has been sampled. 54 In general, similar to other methods discussed in this review, implementing WE in a binned configuration


Page 6 of 26

Page 7 of 26 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Biochemistry

space requires binning a low-dimensional space having sufficient overlap with the slow degree-of-freedom describing the process of interest, where having knowledge of such a low-dimensional space is far from being a trivial problem. It is also important to mention that the parallelization of trajectories in the WE method make extracting macroscopic rate constants non-trivial. 53,55 This is due to the fact that parallelization decreases the relaxation time available to each simulation in order to help the system reach steady-state. 55 Bhatt et al. 53 have devised an approach from which steady-state kinetic information such as the mean first-passage time can be calculated however. Interestingly, in recent work the WE method was used to extract kinetic and mechanistic information from a ligand-protein complex whose residence time is seconds long. 60 2.1.3

Adaptive multilevel splitting

In adaptive multilevel splitting (AMS), 61 the generic objective is to calculate the mean first passage time corresponding to a transition from an initial reactant metastable state, A, to a final product metastable state, B. For this AMS aims to estimate the committor probability, Pc , known to be the exact reaction coordinate (RC) for any arbitrary reaction, and roughly defined as the probability that a trajectory initiated in the reactant state will reach the product state before returning. 62 With the estimated Pc , the mean first passage time corresponding to a transition from A to B can be calculated. Although Pc is the perfect RC, it is important to keep in mind that it is not known a priori – it is after all the quantity that the AMS algorithm is attempting to estimate. The AMS algorithm requires that progress along the RC be measured and an approximate user-defined trial RC z different than Pc is used for this. The basic AMS algorithm begins with N initial MD trajectories launched from a point z = z0 , denoting an initial value corresponding to the reactant state A. The MD simulations are all stopped when the N replicas return to a RC value z = z0 . The next step in the AMS algorithm is to examine the progress that each replica made in sampling the RC during the simulation. For this, AMS involves systematic removal and addition of replica, until the desired rare event from A to B has been sampled and an estimate of Pc from the product of different conditional probabilities for moving between different z values can be calculated. Pc is then directly involved in calculating the rate constant through a simple expression. 61 2.1.4

Milestoning

The idea behind the milestoning method is to divide configuration space into several cells 63,64 in order to break down the dynamics of the system into the local transitions between adjacent cells. 65 In its original implementation, the division of configuration space was performed along an approximate RC, with each discretized point along this RC then defining a 3N −1 dimensional hypersurface called a milestone. 66 In more


Biochemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

recent formulations, a milestone is defined as a dividing surface between two cells. 64 The transitions between milestones represent the most critical event in this class of methods and short MD simulations in parallel are initiated at a given milestone and stopped when the simulation reaches an adjacent milestone. 64–66 These local transition events can be “stitched” together to recover the long-time kinetic properties of the rare event process. 64,65 It is worthwhile to mention that the choice of milestones have often satisfied the condition that the transitions between them are Markovian, 63,65 the consequence of this being that the first-passage-time distribution will be independent of the full history before arriving at a given milestone. A recent development is the introduction of a formally exact formulation where the mean first passage time was derived from exact statistical analysis of the individual local mean first passage times for the functional transitions between the conformational space cells. 64 We would also like to make a brief note that other methods that estimate kinetic properties based on the crossing of hypersurfaces have also been developed. 67–70 Often, these 3N − 1 dimensional hypersurfaces are defined against some order parameter. 67–70 One important class of such methods is transition interface sampling (TIS), which is based on the transition path sampling (TPS) method aiming to sample a whole ensemble of transition (dissociation) paths. 67–70 The ensemble is generated using a Monte Carlo algorithm on a reference path. 67–70 We refer the interested readers to the provided references for more detail.

2.2

Biased sampling based methods

Due to the timescale associated with ligand dissociation processes often being minutes or longer, 25,48 a large number of studies have sought to introduce biasing forces into the molecular description of the system to drive the simulation to sample the process of interest. 33–35,71 Care must be taken, however, to do so in a rigorous manner allowing the statistics of the original, unperturbed process to be recovered. 72 This is relatively straightforward to do when the aim is to recover static equilibrium properties such as free energies, but it turns out to be a much more formidable task when the goal is to recover, from a biased simulation, the unbiased kinetics, although methods to do so have been developed. 73 Unfortunately, the need to introduce bias into molecular simulations can often be severe when one is interested in capturing accurate rate constants and associated mechanisms. The reason is that several rare events must be observed in order to obtain kinetic estimates that are of statistical significance. 74 While binding free energies are state functions that require appropriate sampling of just the initial and final states 22,23 of the desired biochemical process, the kinetics are path-dependent, meaning a whole ensemble of transition paths 75 must be sampled before an averaged kinetic estimate can be produced. Here we summarize one such very popular method, namely metadynamics, that makes it possible to get unbiased rate constants and pathways from biased MD simulations. Note that there are alternative approaches to adding a biasing potential by increasing the temperature of the system, 8 ACS Paragon Plus Environment

Page 8 of 26

Page 9 of 26 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Biochemistry

either whole or along parts, for which we refer the reader to another recent thorough review of methods. 48 The metadynamics 76 method has been the most popular of these biasing approaches due to the development of the infrequent metadynamics framework, which allows unbiased kinetics to be extracted from the biased simulation using an acceleration factor. 73 The basic idea behind metadynamics is to add a historydependent bias potential to the system’s Hamiltonian using a predetermined deposition stride, and as a function of a set of collective variables (CVs) s(x) mapping the high-dimensional x into a low-dimensional representation. The bias V (s, t) is typically constructed as a sum of gaussians deposited along the trajectory in the CV space. As time goes by, this additional bias prevents the system from revisiting the same site in the CV space. The free energy surface (FES) F (s) is then obtained from the negative of the bias potential added. If the FES is in fact well-converged, including in the barrier regions, then the kinetic rate constants can also be estimated from the recovered barrier height between the metastable states using transition state theory (TST). However, a fundamental issue in using a FES derived from metadynamics, umbrella sampling or any other method for rate calculation is that of dynamical corrections to the TST rate. 77 The transition rate calculated via TST provides only a crude upper bound to the true rate, and a transmission coefficient needs to be calculated, for example, by launching trajectories from top of every barrier identified on the FES. A much simpler approach is to use the method named “infrequent metadynamics,” which was proposed as a means of recovering unbiased dynamics from the biased dynamics, and which bends around the calculation of the transmission coefficient by making the assumption that it is defined only by the properties of the barrier. 73 Thus if one could avoid adding bias on the barrier, the transmission coefficient could be kept the same between the biased and unbiased simulations. Infrequent metadynamics allows doing this assuming a decent CV is known for the process being studied. Note that this CV does not have to be the perfect RC (see Sec. 2.1.3 for descritpion of committor probability) but can be a simple combination of order parameters which can be optimized using preliminary inf requent metadynamics runs with the method “Spectral Gap Optimization of Order Parameters” (SGOOP) 78,79 Once such a CV is known, infrequent metadynamics proceeds by simply decreasing the bias addition frequency so that the time interval between two bias addition events becomes slower than the time spent in the barrier regions, thereby avoiding adding bias there. The unbiased timescale can then be obtained from the biased timescale simply through the calculation of an Rt acceleration factor, which is the running average of the exponentiated bias given by α(t) = (1/t) 0 eβV (s) dt where s is the CV being biased and t is the simulation time. The reliability of the unbiased dynamics so reconstructed can then be ascertained through a simple p-value test. 80 The idea of using an acceleration factor in metadynamics to obtain unbiased kinetic information has also been introduced in other variants,


Biochemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

which we do not discuss here. 81–84

2.3

Machine learning based methods

Recent attempts have been made to leverage machine learning to ameliorate the timescale problem in allatom MD simulation of biochemical processes. 35,85–88 Although these machine learning protocols have been used in a few different contexts, the most common approach involves using neural networks trained to learn a low-dimensional representation approximating the true RC for the process of interest. 85–87,89 Finding such low-dimensional representations is often a crucial step in both the biased and unbiased algorithms that we have described. As mentioned in Sec. 2, the binning or subdivision of configuration space into cells in several unbiased methods needs to be done, in practice, on a low-dimensional map possessing sufficient overlap with the slow process of interest, which means a good approximate RC must be known in order to produce good kinetic estimates. Two particular approaches have made direct use of the neural network based RC for either generating an MSM, 85,87 or for running a biased MD method like metadynamics 86 or umbrella sampling. 90 Meanwhile, in recent work, Mardt et al. 89 suggested the use of two parallel neural networks in conjunction with the variational principle for Markov process (VAMP) to propose VAMPnets, 89 which in principle can automate the generation of an MSM via a single end-to-end pipeline. Building an MSM often requires technical expertise on the part of the investigator and could be an important development in increasing the reach of MSM. A different approach, 35,88 which our group has proposed, uses an unsupervised machine learning method to learn both a low-dimensional approximate RC as well as a bias potential along this neural network RC, and hence is its own enhanced sampling method. One interesting feature of this approach, named Reweighted autoencoded variational Bayes for enhanced sampling (RAVE), 35,88 is that the RC and bias potential are learnt from the neural network simultaneously, which turns out to be helpful in screening through spurious machine learning solutions of the biasing parameters. The RAVE protocol involves iterating between machine learning and MD, with each iteration learning a more refined low-dimensional RC representation and bias potential. These iterations then continue until these biasing parameters are converged and sufficient to sample the desired ligand dissociation. One benefit of this approach is that a time-independent bias is produced which can be used to launch multiple independent production MD runs, which is helpful in extracting kinetics information with statistical robustness. Thus far, RAVE has been applied to the benzene-lysozyme complex in order to calculate absolute binding affinities in close agreement with other enhanced sampling methods. 88 In addition, in recent work that is soon to be published, we developed a methodological extension to RAVE and applied it to investigate the kinetics of ligand dissociation. 91


Page 10 of 26

Page 11 of 26 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Biochemistry

3 3.1

CASE STUDIES Sub-microsecond unbinders

Huang and Caflisch 42 have studied the dissociation of six small ligands from the FKBP protein using unbiased MD simulations, which was possible due to their choice of ligands whose dissociation was expected to occur in ∼20 ns, 42 based on experimental dissociation constants. 92 It was found that the kinetics for these dissociations followed single-exponential kinetics and that the unbinding times ranged between four and eighteen ns. 42 An interesting result from their work is the observation that the ligands can take various conformations at different FKBP binding site positions and that there exists a heterogeneity of dissociation paths. In addition, Pan et al. 43 performed long unbiased MD simulations that were µs in length with the Anton supercomputer on the same protein (FKBP) with various small ligands but with a different forcefield. 43 Several spontaneous association and dissociation events were observed and kon and kof f estimated. Similar to the conclusion of Huang and Caflisch 42 , ligand binding to different binding sites were observed, and multiple dissociation paths were sampled. 43 The reported residence times in 8–140 ns range 43 which interestingly is quite different from the numbers reported by some of these ligands in Ref. 42

3.2

Benzamidine-trypsin

The benzamidine-trypsin system has been one of the most popular systems for simulating the protein-ligand association/dissociation process using all-atom MD simulations, starting with the MSM based work in 2011. 93 Since then, numerous other methods have been applied in the last few years to study this system, with the only benchmark 94 that has been used in these studies being a 1970 experimental measurement of kof f =600 s−1 . MSM seems to be systematically overestimating the experimental rate constant, with predicted kof f of 95000 s−1 in 2011, 93 28000 s−1 in 2014 95 and more recently 13000 s−1 in 2015. 96 Use of WE method gives 5555 s−1 , AMS gives 240 s−1 while metadynamics gives 9 s−1 . Thus we have different MD methods giving rates that range from 160 times faster than experiments to 65 times slower than experiments – or a 4 orders of magnitude variation. This is a strikingly large range, likely arising from a combination of the use of different force-fields and different sampling methods. What is however more re-assuring that many of these different methods using different force-fields have obtained similar dominant unbinding pathway of this system, with similar metastable states as defined by the ligand position, orientation, protein conformation, and water molecule location. This more pronounced agreement between pathways than between rate constants is worthy of further discussion and is revisited in Sec. 4.


Biochemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

3.3

Protein kinases

Kinases are arguably the second most important group of drug targets, after G-protein-coupled receptors. Thus naturally there has been a strong interest in applying all-atom MD methods to study the association/dissociation of pharmaceutical ligands from protein kinases. In 2011, the special-purpose supercomputer Anton was used to study the association dynamics of the drug Dasatinib to Src kinase 39 where they were able to achieve spontaneous binding of the ligand (but not unbinding) in around 1 out of 6 trajectories. They found numerous attractor spots on the surface of the protein, and the binding process essentially involved navigating the maze of these short-lived albeit metastable states on the way to the actual bound pose. They also identified intriguing behavior of water molecules during binding, where a single layer of water molecules (lifetime 0.1 µs) was found to be one of the many impediments to binding. Later, Tiwary et al performed infrequent metadynamics simulation of the unbinding of this same system, where they achieved full unbinding in 12 out of 12 simulations, identifying metastable states and their fluxes. 97 They also found waters to play a very critical role. Namely, they found that a well-preserved Glu-Lys salt bridge had to distort before the ligand could exit the protein, aided by solvent water molecules. As we describe in Sec. 4, this finding has important repercussions for possible use of constant pH methods in the simulation of drug unbinding. The calculated association and dissociation rate constants from the Anton and metadynamics simulations respectively were in good agreement with experiments, and perhaps more significantly, with each other. In a similar vein, Casasnovas et al have applied metadynamics to study unbinding of urea-based BIRB family inhibitor from p38 kinase. 98 In this last work, they performed metadynamics using two different sets of CVs, and using the self-consistent measure of Ref., 80 they were able to select one of these as reliable. The kof f value from this agreed well with experiments, demonstrating the usefulness of the test from Ref. 80 in ascertaining the reliability of dynamics reconstructed from metadynamics.

3.4

Heat shock protein 90

Heat shock protein 90 (HSP90) is a chaperone protein generally found in all bacterial and eukaryotic cells, and inhibitors of HSP90 have strong potential to be used as anti-cancer drugs. However, the presence of various conformational states and flexibility of the HSP90 protein at the binding site makes the investigation of kinetics very challenging. 99 Mollica et al. used the smoothed-potential MD or scaled MD which is essentially an approximate and computationally efficient version of infrequent metadynamics, with the intent to focus more on relative dissociation rates rather than absolute rates. They used this to study dissociation of several ligands from HSP90 and obtained agreement with reported experimental results. In a somewhat similar vein, Random acceleration molecular dynamics (RAMD) method has also been used to study binding/unbinding


Page 12 of 26

Page 13 of 26 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Biochemistry

pathways of various proteins 100,101 and calculate relative residence times for a set of various ligands to the N-terminal domain of the HSP90, finding good agreement with reported results.

3.5

G protein–coupled receptors (GPCRs)

G protein–coupled receptors (GPCRs) are of tremendous pharmacological importance, with ∼50% of all drugs on the market designed to act through modulation of GPCRs. 102 GPCRs are trans-membrane proteins comprised of seven α-helices and have the critical cellular function of controlling the cells response to the presence of extracellular molecules. Due to both their fundamental and practical importance it is desirable to investigate the dynamics of GPCRs using all-atom MD. Using unbiased MD on the Anton supercomputer, Dror et al. studied the binding of antagonists and agonists to the β2 adrenergic receptor, 17 observing twelve binding events leading to a kon estimate of 3.1×107 M−1 s−1 , which is in good agreement with the experimental value of 1.0×107 M−1 s−1 . Notice however that the absolute kof f has not been calculated for GPCR-ligand complexes, although Mollica et al. 103 used the dissociation process in order to rank four different triazine-based ligands according to their residence times, relative to a reference. Good agreement with experimental results was reported. 103 In addition, Bortolato et al. 83 determined whether a series of ligands were fast or slow using their recent technique named adiabatic–biased metadynamics (aMetaD) in conjunction with an approximate residence time score. In more recent work, Meral et al. 104 used Metadynamics with a Maximum Caliber approach and reported that thermodynamic properties as well as kinetic properties can be achieved for GPCR activation.

3.6

HIV-1 protease

The HIV-1 protease contains two β-hairpin loops called flaps whose dynamics, according to different studies, 105,106 play an important role in the mechanisms of ligand association and dissociation from the binding pockets. For this reason the HIV-1 protease is a challenging – and interesting – system for MD simulations. Focusing on investigating the mechanism of binding for a peptide with HIV-1 protease, Pietrucci et al. 107 employed the bias-exchange metadynamics (BEMD) method and captured important features of the mechanism such as water bridges, conformational changes of the flaps, various important interactions between the ligand and protease. They calculated the kinetics from a 7-dimensional order parameter space through a weighted histogram method and the kon and kof f values were in rough agreement with the available experimental results. Sun et al. 108 studied kinetics for 6 systems including HIV protease complexed with the drug lopinavir using infrequent metadynamics, finding decent overall agreement in rate constants.


Biochemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

4

CONCLUSIONS: ARE WE THERE YET?

In this short review we have described some of the methods useful for calculating the kinetics and mechanisms of ligand-protein dissociations via all-atom MD simulations, while also highlighting some of their successes and shortcomings with case-studies. A clear message, which is embedded throughout this publication, is that the use of MD simulations in order to investigate ligand-protein dissociations is brimming with potential. This is because MD simulations are a natural complement to experimental kinetics, since it provides a wealth of mechanistic information that are difficult-to-get through experimentation. With this being said, there remain large obstacles for unleashing the full power of MD simulations to practical applications. The first obstacle is the limited timescales that can be attained in all-atom MD simulations, which tend to be quite severe for realistic ligand-protein complexes relevant for drug design. As we have shown throughout this work, a number of ingenious MD-based enhanced sampling methods have been devised to overcome this limitation. We would argue that overcoming the timescale limitation is coupled to the identification of a proper RC in both the biased as well as unbiased methods, even though it might not be obvious for the latter class of methods. There have been recent progress in solving the RC problem through the development of various automated methods, 35,78,109,110 although further developments in automated RC construction for complex biochemical problems are still needed. Meanwhile, a second obstacle is the force-field, which we have mentioned just in passing until this point, and which represents a model of the different interaction terms in the Hamiltonian. These force-fields are most often based on available structural and thermodynamic data that are much richer in describing the longer lived metastable structures. Producing kinetic estimates, however, involves sampling a full transition path such that the MD simulation will have to venture outside these data-rich regions, making it more sensitive to underlying force-field errors. The careful cataloging of the transition states describing rare event processes should find use in helping refine the available force-fields, which would lead to better estimates from MD simulations since their predictions depend on how accurate the force-field is modeling the microscopic interactions between the atoms. In a sense, these two limitations are coupled – developing force-fields that better characterize not just thermodynamic but also kinetic data essentially needs development of reliable sampling methods that can deal with long timescales in order to venture onto configurations outside the data rich regime. In fact, one crucial metric missing thus far is the degree to which the errors in sampling are decoupled from the force-field errors. Here, the benzamidine–trypsin complex serves as a perfect example, with different all-atom MD methods using different force-fields producing dissociation timescales that differ a staggering four orders of magnitude, although it is somewhat reassuring to note that a number of these studies discovered similar dissociation paths and crucial intermediate states. Both the sampling error as well as the force-field


Page 14 of 26

Page 15 of 26 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Biochemistry

error could be significant in deviations of kinetic predictions from experimental values. We would argue that a thorough investigation of how the same all-atom MD methods handle the different available force-fields for the same ligand-protein complexes would be invaluable in providing a sound estimate of the errors from force-fields relative to that of sampling, and presents an important future direction. This would help in establishing how valid is the direct comparison of kof f and kon to experiments. An important example of this is the different water models used in force-fields, which can lead to water diffusivities that are three to four times as different. 111 Similarly, yet another parameter which has not received much mention in typical MD based publications for calculating dissociation timescales and mechanisms is the friction coefficient used in the simulation in order to maintain desired temperature. This coefficient can also have profound effects on the calculations through what is known as the Kramers’ turnover phenomenon. 112,113 Such differences can significantly affect the calculated rate constants and even lead to agreement between experiments and simulations for wrong reasons. We would also like to highlight that some of the slowest dissociating ligand-protein complexes have been found, through all-atom MD simulations, to show salt-bridge distortions, which is suggestive of a protonation state dependence. 97,114,115 Such a finding hints at the usefulness of constant pH simulations for investigating ligand dissociation mechanisms, whose development is an area of active research. In conclusion, all-atom MD simulations provide more than just a rate constant – it holds promise of revealing mechanistic paths and driving forces behind experimental kinetic measurements. It can highlight the role that specific residues and individual water molecules have in driving the overall dissociation process. Although much progress has been made, there remains a long road ahead before all-atom MD simulations can be deemed truly complimentary and even predictive of ligand-protein dissociation experiments. Acknowledgments We thank Zachary Smith and Freddy Alexis Cisneros for their careful reading of this manuscript. We also thank Deepthought2, MARCC and XSEDE (projects CHE180007P and CHE180027P) for helping our group with computational resources used to develop some of the methods reported in this work. Pratyush Tiwary would like to thank the University of Maryland Graduate School for financial support through the Research and Scholarship Award (RASA). Yihang Wang would like to thank the NCI-UMD partnership for Integrative Cancer Research for financial support.

References 1. Dror, R. O., Dirks, R. M., Grossman, J., Xu, H., and Shaw, D. E. (2012) Biomolecular simulation: a computational microscope for molecular biology. Ann. Rev. Biophys. 41, 429–452. 2. Shaw, D. E., Maragakis, P., Lindorff-Larsen, K., Piana, S., Dror, R. O., Eastwood, M. P., Bank, J. A.,


Biochemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Jumper, J. M., Salmon, J. K., and Shan, Y. (2010) Atomic-level characterization of the structural dynamics of proteins. Science 330, 341–346. 3. Faraldo-Gómez, J. D., and Roux, B. (2007) On the importance of a funneled energy landscape for the assembly and regulation of multidomain Src tyrosine kinases. Proc. Natl. Acad. Sci. 104, 13643–13648. 4. Berteotti, A., Cavalli, A., Branduardi, D., Gervasio, F. L., Recanatini, M., and Parrinello, M. (2008) Protein conformational transitions: the closure mechanism of a kinase explored by atomistic simulations. J. Amer. Chem. Soc. 131, 244–250. 5. Dror, R. O., Arlow, D. H., Borhani, D. W., Jensen, M. Ø., Piana, S., and Shaw, D. E. (2009) Identification of two distinct inactive conformations of the β2-adrenergic receptor reconciles structural and biochemical observations. Proc. Natl. Acad. Sci. 106, 4689–4694. 6. Vanni, S., Neri, M., Tavernelli, I., and Rothlisberger, U. (2009) Observation of ionic lock formation in molecular dynamics simulations of wild-type β1 and β2 adrenergic receptors. Biochemistry 48, 4789– 4797. 7. Lyman, E., Higgs, C., Kim, B., Lupyan, D., Shelley, J. C., Farid, R., and Voth, G. A. (2009) A role for a specific cholesterol interaction in stabilizing the Apo configuration of the human A2A adenosine receptor. Structure 17, 1660–1668. 8. Wang, Y., Papaleo, E., and Lindorff-Larsen, K. (2016) Mapping transiently formed and sparsely populated conformations on a complex energy landscape. eLife 5, e17505. 9. Seibert, M. M., Patriksson, A., Hess, B., and Van Der Spoel, D. (2005) Reproducible polypeptide folding and structure prediction using molecular dynamics simulations. Journal of molecular biology 354, 173–183. 10. Simmerling, C., Strockbine, B., and Roitberg, A. E. (2002) All-atom structure prediction and folding simulations of a stable protein. J. Amer. Chem. Soc. 124, 11258–11259. 11. Voelz, V. A., Bowman, G. R., Beauchamp, K., and Pande, V. S. (2010) Molecular simulation of ab initio protein folding for a millisecond folder NTL9 (1- 39). J. Amer. Chem. Soc. 132, 1526–1528. 12. Arkin, I. T., Xu, H., Jensen, M. Ø., Arbely, E., Bennett, E. R., Bowers, K. J., Chow, E., Dror, R. O., Eastwood, M. P., and Flitman-Tene, R. (2007) Mechanism of Na+/H+ antiporting. Science 317, 799– 803. 16 ACS Paragon Plus Environment

Page 16 of 26

Page 17 of 26 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Biochemistry

13. Bostick, D. L., and Brooks, C. L. (2007) Selectivity in K+ channels is due to topological control of the permeant ion’s coordinated state. Proc. Natl. Acad. Sci. 104, 9260–9265. 14. Enkavi, G., and Tajkhorshid, E. (2010) Simulation of spontaneous substrate binding revealing the binding pathway and mechanism and initial conformational response of GlpT. Biochemistry 49, 1105– 1114. 15. Jensen, M. Ø., Borhani, D. W., Lindorff-Larsen, K., Maragakis, P., Jogini, V., Eastwood, M. P., Dror, R. O., and Shaw, D. E. (2010) Principles of conduction and hydrophobic gating in K+ channels. Proc. Natl. Acad. Sci. 107, 5833–5838. 16. Noskov, S. Y., Berneche, S., and Roux, B. (2004) Control of ion selectivity in potassium channels by electrostatic and dynamic properties of carbonyl ligands. Nature 431, 830. 17. Dror, R. O., Pan, A. C., Arlow, D. H., Borhani, D. W., Maragakis, P., Shan, Y., Xu, H., and Shaw, D. E. (2011) Pathway and mechanism of drug binding to G-protein-coupled receptors. Proc. Natl. Acad. Sci. 108, 13118–13123. 18. Shan, Y., Kim, E. T., Eastwood, M. P., Dror, R. O., Seeliger, M. A., and Shaw, D. E. (2011) How does a drug molecule find its target binding site? J. Amer. Chem. Soc. 133, 9181–9183. 19. Hurst, D. P., Grossfield, A., Lynch, D. L., Feller, S., Romo, T. D., Gawrisch, K., Pitman, M. C., and Reggio, P. H. (2010) A lipid pathway for ligand binding is necessary for a cannabinoid G protein-coupled receptor. J. Biol. Chem. jbc–M109. 20. Mondal, J., Ahalawat, N., Pandit, S., Kay, L. E., and Vallurupalli, P. (2018) Atomic resolution mechanism of ligand binding to a solvent inaccessible cavity in T4 lysozyme. PLoS computational biology 14, e1006180. 21. Zeevaart, J. G., Wang, L., Thakur, V. V., Leung, C. S., Tirado-Rives, J., Bailey, C. M., Domaoal, R. A., Anderson, K. S., and Jorgensen, W. L. (2008) Optimization of azoles as anti-human immunodeficiency virus agents guided by free-energy calculations. J. Amer. Chem. Soc. 130, 9492–9499. 22. Deng, Y., and Roux, B. (2006) Calculation of standard binding free energies: aromatic molecules in the T4 lysozyme L99A mutant. J. Chem. Theor. Comp. 2, 1255–1273. 23. Deng, Y., and Roux, B. (2009) Computations of standard binding free energies with molecular dynamics simulations. J. Phys. Chem. B 113, 2234–2246. 17 ACS Paragon Plus Environment

Biochemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

24. Wang, L., Wu, Y., Deng, Y., Kim, B., Pierce, L., Krilov, G., Lupyan, D., Robinson, S., Dahlgren, M. K., and Greenwood, J. (2015) Accurate and reliable prediction of relative ligand binding potency in prospective drug discovery by way of a modern free-energy calculation protocol and force field. J. Amer. Chem. Soc. 137, 2695–2703. 25. Copeland, R. A., Pompliano, D. L., and Meek, T. D. (2006) Drug–target residence time and its implications for lead optimization. Nat. Rev. Drug. Discov. 5, 730–739. 26. Copeland, R. A. (2016) The drug–target residence time model: a 10-year retrospective. Nature Reviews Drug Discovery 15, 87. 27. Pan, A. C., Borhani, D. W., Dror, R. O., and Shaw, D. E. (2013) Molecular determinants of drug– receptor binding kinetics. Drug Discov. Today 18, 667–673. 28. Tiwary, P., Mondal, J., Morrone, J. A., and Berne, B. J. (2015) Role of water and steric constraints in the kinetics of cavityligand unbinding. Proc. Natl. Acad. Sci. 112, 12015–12019. 29. Tiwary, P., and Berne, B. J. (2016) How wet should be the reaction coordinate for ligand unbinding? J. Chem. Phys. 145, 054113–054118. 30. Wang, Y., Martins, J. M., and Lindorff-Larsen, K. (2017) Biomolecular conformational changes and ligand binding: from kinetics to thermodynamics. Chemical science 8, 6466–6473. 31. Tummino, P. J., and Copeland, R. A. (2008) Residence time of receptor- ligand complexes and its effect on biological function. Biochemistry 47, 5481–5492. 32. Tiwary, P., Limongelli, V., Salvalaglio, M., and Parrinello, M. (2015) Kinetics of proteinligand unbinding: Predicting pathways, rates, and rate-limiting steps. Proc. Natl. Acad. Sci. 112, E386–E391. 33. Torrie, G. M., and Valleau, J. P. (1977) Nonphysical sampling distributions in Monte Carlo free-energy estimation: Umbrella sampling. J. Comp. Phys. 23, 187–199. 34. Laio, A., and Parrinello, M. (2002) Escaping free-energy minima. Proc Natl Acad Sci 99, 12562–12566. 35. Ribeiro, J. M. L., Bravo, P., Wang, Y., and Tiwary, P. (2018) Reweighted autoencoded variational Bayes for enhanced sampling. J. Chem. Phys. 149, 072301–072309.


Page 18 of 26

Page 19 of 26 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Biochemistry

36. Georgi, V., Schiele, F., Berger, B.-T., Steffen, A., Marin Zapata, P. A., Briem, H., Menz, S., Preuße, C., Vasta, J. D., Robers, M. B., Brands, M., Knapp, S., and Fernández-Montalván, A. E. (2018) Binding Kinetics Survey of the Drugged Kinome. Journal of the American Chemical Society 37. Shaw, D. E., Deneroff, M. M., Dror, R. O., Kuskin, J. S., Larson, R. H., Salmon, J. K., Young, C., Batson, B., Bowers, K. J., and Chao, J. C. (2008) Anton, a special-purpose machine for molecular dynamics simulation. Communications of the ACM 51, 91–97. 38. Shaw, D. E. et al. Millisecond-scale Molecular Dynamics Simulations on Anton. Proceedings of the Conference on High Performance Computing Networking, Storage and Analysis. New York, NY, USA, 2009; pp 39:1–39:11. 39. Shan, Y., Kim, E. T., Eastwood, M. P., Dror, R. O., Seeliger, M. A., and Shaw, D. E. (2011) How Does a Drug Molecule Find Its Target Binding Site? J. Amer. Chem. Soc. 133, 9181–9183. 40. Dror, R. O., Dirks, R. M., Grossman, J., Xu, H., and Shaw, D. E. (2012) Biomolecular Simulation: A Computational Microscope for Molecular Biology. Annual Review of Biophysics 41, 429–452. 41. Lindorff-Larsen, K., Trbovic, N., Maragakis, P., Piana, S., and Shaw, D. E. (2012) Structure and Dynamics of an Unfolded Protein Examined by Molecular Dynamics Simulation. J. Amer. Chem. Soc. 134, 3787–3791. 42. Huang, D., and Caflisch, A. (2011) The free energy landscape of small molecule unbinding. PLoS Comput. Biol. 7, e1002002. 43. Pan, A. C., Xu, H., Palpant, T., and Shaw, D. E. (2017) Quantitative Characterization of the Binding and Unbinding of Millimolar Drug Fragments with Molecular Dynamics Simulations. J. Chem. Theor. Comp. 13, 3372–3377. 44. Pande, V. S., Beauchamp, K., and Bowman, G. R. (2010) Everything you wanted to know about Markov State Models but were afraid to ask. Methods 52, 99–105. 45. Prinz, J.-H., Wu, H., Sarich, M., Keller, B., Senne, M., Held, M., Chodera, J. D., Schtte, C., and No, F. (2011) Markov models of molecular kinetics: Generation and validation. J. Chem. Phys. 134, 174105–174115. 46. Chodera, J. D., and Noé, F. (2014) Markov state models of biomolecular conformational dynamics. Current opinion in structural biology 25, 135–144. 19 ACS Paragon Plus Environment

Biochemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

47. Husic, B. E., and Pande, V. S. (2018) Markov state models: From an art to a science. J. Amer. Chem. Soc. 140, 2386–2396. 48. Dickson, A., Tiwary, P., and Vashisth, H. (2017) Kinetics of Ligand Binding Through Advanced Computational Approaches: A Review. Current topics in medicinal chemistry 49. Noé, F., and Nuske, F. (2013) A variational approach to modeling slow processes in stochastic dynamical systems. Multiscale Modeling & Simulation 11, 635–655. 50. Pérez-Hernández, G., Paul, F., Giorgino, T., De Fabritiis, G., and Noé, F. (2013) Identification of slow molecular order parameters for Markov model construction. J. Chem. Phys. 139, 015102. 51. Mey, A. S., Wu, H., and Noé, F. (2014) xTRAM: Estimating equilibrium expectations from timecorrelated simulation data at multiple thermodynamic states. Phys. Rev. X 4, 041018. 52. Stelzl, L. S., Kells, A., Rosta, E., and Hummer, G. (2017) Dynamic histogram analysis to determine free energies and rates from biased simulations. Journal of chemical theory and computation 13, 6328–6342. 53. Bhatt, D., Zhang, B. W., and Zuckerman, D. M. (2010) Steady-state simulations using weighted ensemble path sampling. J. Chem. Phys. 133, 014110. 54. Zhang, B. W., Jasnow, D., and Zuckerman, D. M. (2010) The weighted ensemble path sampling method is statistically exact for a broad class of stochastic processes and binning procedures. J. Chem. Phys. 132, 054107. 55. Zuckerman, D. M., and Chong, L. T. (2017) Weighted ensemble simulation: review of methodology, applications, and software. Ann. Rev. Biophys. 46, 43–57. 56. Huber, G. A., and Kim, S. (1996) Weighted-ensemble Brownian dynamics simulations for protein association reactions. Biophysical journal 70, 97. 57. Dickson, A., and Lotz, S. D. (2016) Ligand release pathways obtained with WExplore: residence times and mechanisms. The Journal of Physical Chemistry B 120, 5377–5385. 58. Dickson, A., and Lotz, S. D. (2017) Multiple ligand unbinding pathways and ligand-induced destabilization revealed by WExplore. Biophysical journal 112, 620–629. 59. Nunes-Alves, A., Zuckerman, D. M., and Arantes, G. M. (2018) Escape of a Small Molecule from Inside T4 Lysozyme by Multiple Pathways. Biophysical journal 114, 1058–1066. 20 ACS Paragon Plus Environment

Page 20 of 26

Page 21 of 26 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Biochemistry

60. Lotz, S. D., and Dickson, A. (2018) Unbiased molecular dynamics of 11 min timescale drug unbinding reveals transition state stabilizing interactions. Journal of the American Chemical Society 140, 618–628. 61. Teo, I., Mayne, C. G., Schulten, K., and Lelivre, T. (2016) Adaptive Multilevel Splitting Method for Molecular Dynamics Calculation of Benzamidine-Trypsin Dissociation Time. J. Chem. Theor. Comp. 12, 2983–2989. 62. Peters, B., and Trout, B. L. (2006) Obtaining reaction coordinates by likelihood maximization. J. Chem. Phys. 125, 054108–054114. 63. Vanden-Eijnden, E., and Venturoli, M. (2009) Markovian milestoning with Voronoi tessellations. J. Chem. Phys. 130, 194101. 64. Bello-Rivas, J. M., and Elber, R. (2015) Exact milestoning. J. Chem. Phys. 142, 03B602 1. 65. Vanden-Eijnden, E., Venturoli, M., Ciccotti, G., and Elber, R. (2008) On the assumptions underlying milestoning. J. Chem. Phys. 129, 174102. 66. Faradjian, A. K., and Elber, R. (2004) Computing time scales from reaction coordinates by milestoning. J. Chem. Phys. 120, 10880–10889. 67. Moroni, D., Bolhuis, P. G., and van Erp, T. S. (2004) Rate constants for diffusive processes by partial path sampling. The Journal of chemical physics 120, 4055–4065. 68. Allen, R. J., Frenkel, D., and ten Wolde, P. R. (2006) Forward flux sampling-type schemes for simulating rare events: Efficiency analysis. The Journal of chemical physics 124, 194111. 69. Escobedo, F. A., Borrero, E. E., and Araque, J. C. (2009) Transition path sampling and forward flux sampling. Applications to biological systems. Journal of Physics: Condensed Matter 21, 333101. 70. Elber, R. (2016) Perspective: Computer simulations of long time dynamics. The Journal of chemical physics 144, 060901. 71. Barducci, A., Bussi, G., and Parrinello, M. (2008) Well-tempered metadynamics: A smoothly converging and tunable free-energy method. Phys Rev Lett 100, 020603–020606. 72. Tiwary, P., and Parrinello, M. (2014) A time-independent free energy estimator for metadynamics. J. Phys. Chem. B 119, 736–742.


Biochemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

73. Tiwary, P., and Parrinello, M. (2013) From Metadynamics to Dynamics. Phys. Rev. Lett. 111, 230602– 230606. 74. Bowman, G. R. (2016) Accurately modeling nanosecond protein dynamics requires at least microseconds of simulation. Journal of computational chemistry 37, 558–566. 75. Bolhuis, P. G., Chandler, D., Dellago, C., and Geissler, P. L. (2002) Transition path sampling: Throwing ropes over rough mountain passes, in the dark. Ann. Rev. Phys. Chem. 53, 291–318. 76. Valsson, O., Tiwary, P., and Parrinello, M. (2016) Enhancing Important Fluctuations: Rare Events and Metadynamics from a Conceptual Viewpoint. Ann. Rev. Phys. Chem. 67, 159–184. 77. Tiwary, P., and van de Walle, A. Multiscale Materials Modeling for Nanomechanics; Springer, 2016; pp 195–221. 78. Tiwary, P., and Berne, B. J. (2016) Spectral gap optimization of order parameters for sampling complex molecular systems. Proc. Natl. Acad. Sci. 113, 2839–2844. 79. Tiwary, P. (2017) Molecular Determinants and Bottlenecks in the Dissociation Dynamics of Biotin– Streptavidin. J. Phys. Chem. B 121, 10841–10849. 80. Salvalaglio, M., Tiwary, P., and Parrinello, M. (2014) Assessing the Reliability of the Dynamics Reconstructed from Metadynamics. J. Chem. Theor. Comp. 10, 1420–1425. 81. McCarty, J., Valsson, O., Tiwary, P., and Parrinello, M. (2015) Variationally optimized free-energy flooding for rate calculation. Physical review letters 115, 070601. 82. Fu, C. D., L. Oliveira, L. F., and Pfaendtner, J. (2017) Determining energy barriers and selectivities of a multi-pathway system with infrequent metadynamics. J. Chem. Phys. 146, 014108. 83. Bortolato, A., Deflorian, F., Weiss, D. R., and Mason, J. S. (2015) Decoding the Role of Water Dynamics in LigandProtein Unbinding: CRF1R as a Test Case. J. Chem. Inf. Mod. 55, 1857–1866. 84. Wang, Y., Valsson, O., Tiwary, P., Parrinello, M., and Lindorff-Larsen, K. (2018) Frequency adaptive metadynamics for the calculation of rare-event kinetics. The Journal of Chemical Physics 149, 072309. 85. Wehmeyer, C., and Noé, F. (2017) Time-lagged autoencoders: Deep learning of slow collective variables for molecular kinetics. arXiv preprint arXiv:1710.11239


Page 22 of 26

Page 23 of 26 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Biochemistry

86. Sultan, M. M., Wayment-Steele, H. K., and Pande, V. S. (2018) Transferable neural networks for enhanced sampling of protein dynamics. arXiv preprint arXiv:1801.00636 87. Hernández, C. X., Wayment-Steele, H. K., Sultan, M. M., Husic, B. E., and Pande, V. S. (2018) Variational encoding of complex dynamics. Phys. Rev. E 97, 062412. 88. Ribeiro, J. M. L., and Tiwary, P. (2018) Achieving Reversible Ligand-Protein Unbinding with Deep Learning and Molecular Dynamics through RAVE. bioRxiv 400002. 89. Mardt, A., Pasquali, L., Wu, H., and Noé, F. (2018) VAMPnets for deep learning of molecular kinetics. Nat. Comm. 9, 5–12. 90. Chen, W., Tan, A. R., and Ferguson, A. L. (2018) Collective variable discovery and enhanced sampling using autoencoders: Innovations in network architecture and error function design. J. Chem. Phys. 149, 072312. 91. Wang, Y., Ribeiro, J. M. L., and Tiwary, P. (2018) Past–future information bottleneck based framework for simultaneous sampling of reaction coordinate, thermodynamics and kinetics in rare event systems. In preparation 92. Burkhard, P., Taylor, P., and Walkinshaw, M. D. (2000) X-ray structures of small ligand-FKBP complexes provide an estimate for hydrophobic interaction energies1. Journal of molecular biology 295, 953–962. 93. Buch, I., Giorgino, T., and De Fabritiis, G. (2011) Complete reconstruction of an enzyme-inhibitor binding process by molecular dynamics simulations. Proc. Natl. Acad. Sci. 108, 10184–10189. 94. Guillain, F., and Thusius, D. (1970) Use of proflavine as an indicator in temperature-jump studies of the binding of a competitive inhibitor to trypsin. J. Amer. Chem. Soc. 92, 5534–5536. 95. Doerr, S., and De Fabritiis, G. (2014) On-the-fly learning and sampling of ligand binding by highthroughput molecular simulations. J. Chem. Theor. Comp. 10, 2064–2069. 96. Plattner, N., and Noé, F. (2015) Protein conformational plasticity and complex ligand-binding kinetics explored by atomistic simulations and Markov models. Nat. Comm. 6 . 97. Tiwary, P., Mondal, J., and Berne, B. J. (2017) How and when does an anticancer drug leave its binding site? Science Advances 3 .


Biochemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

98. Casasnovas, R., Limongelli, V., Tiwary, P., Carloni, P., and Parrinello, M. (2017) Unbinding Kinetics of a p38 MAP Kinase Type II Inhibitor from Metadynamics Simulations. J. Amer. Chem. Soc. 139, 4780–4788. 99. Bruce, N. J., Ganotra, G. K., Kokh, D. B., Sadiq, S. K., and Wade, R. C. (2018) New approaches for computing ligandreceptor binding kinetics. Current Opinion in Structural Biology 49, 1 – 10, Theory and simulation Macromolecular assemblies. 100. Rydzewski, J., and Nowak, W. (2015) Memetic algorithms for ligand expulsion from protein cavities. J. Chem. Phys. 143, 124101. 101. Tang, Z., and Chang, C.-e. A. (2017) Systematic Dissociation Pathway Searches Guided by Principal Component Modes. J. Chem. Theor. Comp. 13, 2230–2244. 102. Garland, S. L. (2013) Are GPCRs Still a Source of New Targets? Journal of Biomolecular Screening 18, 947–966. 103. Mollica, L., Decherchi, S., Zia, S. R., Gaspari, R., Cavalli, A., and Rocchia, W. (2015) Kinetics of protein-ligand unbinding via smoothed potential molecular dynamics simulations. Scientific Reports 5, 11539. 104. Meral, D., Provasi, D., and Filizola, M. (2018) An Efficient Strategy to Estimate Thermodynamics and Kinetics of G Protein-Coupled Receptor Activation Using Metadynamics and Maximum Caliber. bioRxiv 105. Karthik, S., and Senapati, S. (2011) Dynamic flaps in HIV-1 protease adopt unique ordering at different stages in the catalytic cycle. Proteins: Structure, Function, and Bioinformatics 79, 1830–1840. 106. Huang, Y.-m. M., Raymundo, M. A. V., Chen, W., and Chang, C.-e. A. (2017) Mechanism of the Association Pathways for a Pair of Fast and Slow Binding Ligands of HIV-1 Protease. Biochemistry 56, 1311–1323. 107. Pietrucci, F., Marinelli, F., Carloni, P., and Laio, A. (2009) Substrate Binding Mechanism of HIV-1 Protease from Explicit-Solvent Atomistic Simulations. J. Amer. Chem. Soc. 131, 11811–11818. 108. Sun, H., Li, Y., Shen, M., Li, D., Kang, Y., and Hou, T. (2017) Characterizing DrugTarget Residence Time with Metadynamics: How To Achieve Dissociation Rate Efficiently without Losing Accuracy against Time-Consuming Approaches. J. Chem. Inf. Mod. 57, 1895–1906. 24 ACS Paragon Plus Environment

Page 24 of 26

Page 25 of 26 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Biochemistry

109. M. Sultan, M., and Pande, V. S. (2017) tICA-Metadynamics: Accelerating Metadynamics by Using Kinetically Selected Collective Variables. J. Chem. Theor. Comp. 13, 2440–2447, PMID: 28383914. 110. McCarty, J., and Parrinello, M. (2017) A variational conformational dynamics approach to the selection of collective variables in metadynamics. J. Chem. Phys. 147, 204109. 111. Tiwary, P., Limongelli, V., Salvalaglio, M., and Parrinello, M. (2015) Kinetics of protein–ligand unbinding: Predicting pathways, rates, and rate-limiting steps. Proc. Natl. Acad. Sci. 112, E386–E391. 112. Klimov, D., and Thirumalai, D. (1997) Viscosity dependence of the folding rates of proteins. Phys. Rev. Lett. 79, 317. 113. Tiwary, P., and Berne, B. J. (2016) Kramers turnover: From energy diffusion to spatial diffusion using metadynamics. J. Chem. Phys. 144, 134103–134106. 114. Shan, Y., Seeliger, M. A., Eastwood, M. P., Frank, F., Xu, H., Jensen, M. Ø., Dror, R. O., Kuriyan, J., and Shaw, D. E. (2008) A conserved protonation-dependent switch controls drug binding in the Abl kinase. Proc. Natl. Acad. Sci. pnas–0811223106. 115. Wallace, J. A., and Shen, J. K. (2011) Continuous constant pH molecular dynamics in explicit solvent with pH-based replica exchange. J. Chem. Theor. Comp. 7, 2617–2629.


Biochemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Graphical TOC Entry


Page 26 of 26

Kinetics of Ligand-Protein Dissociation from All-Atom Simulations: Are

Recommend Documents