Reaction Coordinate-Free Approach to Recovering Kinetics from

Subscriber access provided by UCL Library Services

Article

A Reaction Coordinate-Free Approach to Recovering Kinetics from Potential-Scaled Simulations: Application of Kramers’ Rate Theory Aaron Terrence Frank, and Ioan Andricioaei J. Phys. Chem. B, Just Accepted Manuscript • DOI: 10.1021/acs.jpcb.6b02654 • Publication Date (Web): 24 May 2016 Downloaded from http://pubs.acs.org on May 25, 2016

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 free 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 accessible to all readers and 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.

The Journal of Physical Chemistry B 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 23

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

The Journal of Physical Chemistry

A Reaction Coordinate-Free Approach to Recovering Kinetics from Potential-Scaled Simulations: Application of Kramers’ Rate Theory Aaron T. Frank∗,†,‡ and Ioan Andricioaei∗,† †Department of Chemistry, The University of California, Irvine, 4212 Natural Sciences 1, Irvine, CA 92697, USA ‡Current address: Departments of Chemistry and Biophysics, University of Michigan, 930 North University Avenue, Ann Arbor, Michigan 48109, USA E-mail: [email protected]; [email protected] Phone: (949) 824-3569. Fax: (949) 824-9920

1

ACS Paragon Plus Environment


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 Enhanced sampling techniques are used to increase the frequency of “rare events” during computer simulations of complex molecules. Although methods exist that allow accurate thermodynamics to be recovered from enhanced simulations, recovering kinetics proves to be more challenging. Here we present an extrapolation approach that allows reliable kinetics to be recovered from potential-scaled MD simulations. The approach, based on Kramers’ rate theory, is simple and computationally efficient, and allows kinetics to be recovered without defining reaction coordinates. To test our approach, we use it to determine the kinetics of barrier crossing between two metastable states on the 2D-M¨ uller potential and the C7eq to αR transition in alanine dipeptide. The mean first passage time estimates obtained are in excellent agreement with reference values obtained from direct simulations on the unscaled potentials performed over times that are orders of magnitude longer.

Introduction Computing the kinetics and relaxation of complex physico-chemical systems presents significant difficulties 1 . A particularly outstanding challenge is the task of estimating rates of long-time conformational transitions in biomolecules. Examples include protein folding, transitions between high populated and transient substates in proteins and nucleic acids, or the association and dissociation of a ligand to and from a receptor. In biomolecules, the high-dimensional energy landscape is rugged, with functionally important conformational transitions frequently occurring between states that are separated by barriers that are ≫ kB T , and that are either energetic or entropic in origin. Molecular dynamics (MD) simulation 2 is usually the tool of choice for studying the exact time-dependendence of complex conformational transitions that occur in biomolecules – MD generates atomic-level trajectories of a system, from which the kinetics can be extracted based on the number of conformational transitions that occur along these trajectories. This 2


Page 2 of 23

Page 3 of 23

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


direct approach, however, is rarely useful – typical simulations are too short to observe a statistically relevant number of transitions between vastly different conformational substates. This “mismatch” between simulation time, typically ranging from nanoseconds to microseconds, and the time required to observe a statistically relevant number of functionally important transitions, usually on the order of milliseconds and beyond, precludes the extraction of reliable kinetics from brute-force simulations. Recent advances in computing architecture and networks have gone a long way toward addressing the timescale-mismatch issue of molecular dynamics simulations: Shaw and coworkers, for example, recently developed special-purpose MD simulation machines 3 that enable relatively long (i.e., microsecond to millisecond) simulations to be carried out 4,5 . Another example involves distributed computing networks 6 which enable one to carry out a large number of nanosecond to microsecond simulations. Subsequently, using a Markov modeling approach, kinetics can be extracted from these simulations 7–9 . In similar work, a distributed computing network based on GPUs (GPUgrid) has been used to explore long time scale processes like ligand association and dissociation 10,11 . Exciting and worth pursuing as these developments are, complementary approaches are still needed to enable kinetics to be extracted from MD simulations executed on readily available computational resources and for processes that span beyond the millisecond time window, a time domain during which much molecular biology occurs. Alternative strategies to those highlighted above use enhanced sampling methods that circumvent the timescale-mismatch issue by maximizing sampling efficiency. These approaches maximize efficiency by increasing the frequency of so called rare events (i.e. transitions between metastable states) 12–18 ; unfortunately they also “distort” time such that kinetics cannot be extracted directly from the resulting trajectories. A number of methods have been developed that allow “true” kinetics to be recovered; these include (but are not limited to): the transition-state theory based time-renormalization method 13 ; the Onsanger action based reweighting method 19,20 ; milestoning 21,22 ; transition path sampling 23,24 ; the

3



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 ACS Paragon Plus Environment

Page 4 of 23

Page 5 of 23

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


accelerated, but here we assume, for purposes of presenting the proof of concept, without loss of generality, uniform scaling. In what follows we describe the theoretical underpinnings of the extrapolation method, and then show the results of its application to conformational transitions in two model systems, namely transitions between two local minima on the M¨ uller potential, and the C7eq to αR transition in alanine dipeptide. Our results indicate that the method can accurately and efficiently recover the kinetics of such rare events. We also compare and contrast our approach to the approach of Hamelberg and coworkers, which is also based on Kramers’ rate theory, but was developed to recover kinetics from the accelerated MD simulation framework developed by McCammon 36–48 and his coworkers. 49–51

Method The starting point for our method is the Kramers’ rate equation 52 for a transition from a “reactant” region (F) to a “product” region (U ): 1 ∼ ωmin ωmax ∆U † kK = = . exp − τ 2πζ kB T

(1)

Here kK is the Kramers’ rate constant, τ is the mean first passage time of the transitions between F and U , ωmin and ωmax are the curvature of the effective potential (in effect, a free energy profile) near the minima and maxima respectively, ζ is the frictional constant, kB is the Boltzmann constant, T is temperature, and ∆U † is the free energy barrier for the transition. This equation is valid in the “spatial-diffusive regime” (i.e., in the high friction limit and when barrier heights ≫ kB T). Unfortunately, direct application of Eq. (1) requires some knowledge of the topological features of the free energy landscape. The potential of mean force, i.e., the free energy along a reaction coordinate z (a 1-d manifold chosen such that dynamics along it is much slower than dynamics in the perpendicular

5



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

Page 6 of 23

directions) can be expressed as

exp (−βU (z)) =

Z

exp (−βV (x)) δ(x − z)dx

(2)

where V (x) is the potential energy as a function of the coordinates, x. From this equation, it is evident that the scaling of the potential energy V (x) → αV (x), with 0 < α < 1, amounts to scaling the free energy U by α. More precisely, for any z, Z

(1) (α) = exp −β ′ Uβ ′ exp (−βαV (x)) dx = exp −βUβ

(3)

where, by definition,

(α)

Uβ

(1)

= αUβ ′

(4)

and β ′ = βα.

(5)

Using Eq. 1 and Eq. 4, the Kramers’ predicted rate at β and on the unscaled potential, τβ (1), can be expressed as

τβ (1) =

2πζ β β ωmin (1)ωmax (1)

†(1) exp β∆Uβ

(6)

the Kramers’ predicted rate at β ′ and on the unscaled potential, τβ′ (1), as τβ′ (1) =

2πζ β β ωmin (α)ωmax (α) ′

′

†(α) exp β ′ ∆Uβ ′

(7)

and, the Kramers’ predicted rate at β and on the potential scaled by α, τβ , (α), as

τβ (α) =

2πζ β β ωmin (α)ωmax (α)

6

†(α) exp β∆Uβ


(8)

Page 7 of 23

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


Next, for relatively low temperatures, given that ωmin ∼ β and ωmax ∼ β and that

=

s

d2 Uβ ′ = dz 2

s

1 d2 Uβ α dz 2

ωmax (1) =

s

(1) d2 Uβ ′ dz 2

s

1 d2 Uβ α dz 2

β′ ωmin (1)

β′

(1)

=

(α)

=

r ! 1 β ωmin (α), α

(9)

=

r ! 1 β ωmax (α), α

(10)

(α)

and finally, that, with a slight abuse of notation, β∆Uβ† = β∆H † − ∆S †

(11)

1 τβ (α) = τβ ′ (1) α

(12)

it is easy to show that

and that τβ (α) 1 = 2 exp β∆H † (α − 1) . τβ (1) α

(13)

By determining τ as function of α and then fitting to Eq. 13, τ (1) and ∆H † can been determined. Eq. 13 therefore provides a simple means of recovering reaction rates from simulations on scaled potentials. For more complicated temperatures dependences of the curvature of the minima and maxima, it can be shown that the polynomial prefactor 1/α2 is generalized to 1/αm , with 1 < m < 2 (for example, m = 3/2 corresponds to a temperature insensitive barrier top, while m = 1 corresponds to a constant entropy (temperature independent) free energy profile. We however stress that other than the assumptions needed for the applicability of Kramers’ rate equation, no additional knowledge about the topology of the energy surface is assumed. All that is required is an ability to partition conformational space into reactant (F) and product (U ) regions. This Kramers’-based extrapolation approach can therefore be used to recover kinetics from a series of potential-scaled MD simulations, without defining a reaction coordinates (although an effective one is assumed for

7



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

Page 8 of 23

validity of the Kramers’ theory) that describes the transition between F and U ; in this sense, ours is a reaction-coordinate free approach.

Computational Details To test our approach to recovering kinetics from potential-scaled MD, we examined the kinetics of barrier crossing on the two-dimensional M¨ uller potential and the kinetics of a transition between two metastable states in alanine dipeptide. In each case, the kinetic parameters calculated using our method were compared to those obtained from reference simulations carried out on the untransformed potentials.

2D M¨ uller potential As an initial test of the applicability of our approach, the mean first passage time τ for the transition between two local minima on the M¨ uller potential 53 was determined. The 2d M¨ uller potential is given by:

U (x, y) =

4 X

Ai exp[ai (x − x0,i )2 + bi (x − x0,i )(y − y0,i ) + ci (y − y0,i )2 ]

(14)

i=1

where A ∈ {0,-10,-10,0.5}, a ∈ {-1,-1,-6.5,0.7}, b ∈ {0,0,11,0.6}, c ∈ {-10,-10,-6.5,0.7}, x0 ∈ {1,0,-0.5,-1} and y0 ∈ {0,0.5,1.5,1}. We partitioned the space into F and U regions using the dividing line y = 1.15x+0.53. Region F was defined as the basin to the left of the dividing line and region U defined the as the basin to right of the dividing line (Fig. 2A). To apply Eq. 13, τ (α) was calculated from a series of potential-scaled MD simulations on the M¨ uller surface with α = 0.5, 0.6, 0.7, 0.8, and 0.9. To make comparison with τ (1.0) obtained from fitting to Eq. 13, τ was also determined on the untransformed M¨ uller potential (i.e. with α=1.0). For a given α, the potential was scaled and 1500 canonically distributed states were prepared in F. These initial states were then propagated using

8


Page 9 of 23

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


Langevin dynamics with γ= 100 (in arbitrary units) and kB T =0.5. The trajectories were terminated the first instance they arrived in U by crossing the dividing line from F. The transition times was recorded, and then τ (α) was calculated as ensemble averages over the initial states. All simulations were carried out using in-house code.

Alanine dipeptide To test the applicability of our approach on a more biological relevant system, we studied kinetics of transition between two metastable states in alanine dipeptide. Alanine dipeptide is a widely used test case because it contains the complexities of typical proteins, yet it is small enough that accurate reference calculations can be carried out. Alanine dipeptide has a number of metastable states, but we focus on C7eq and αR . C7eq was defined as residing in the region where ψ ≥ 39◦ and ψ ≤ -117◦ , and αR was defined as residing in the region where -117◦ < ψ < 39◦ . In both cases, -180◦ ≤ φ ≤ 0◦ (Fig. 3A). To apply Eq. 13, τ (α) was calculated from a series of potential-scaled MD simulations with α = 0.5, 0.6, 0.7, 0.8, and 0.9. To make comparison with τ (1.0) obtained from fitting to Eq. 13, τ was also determined on the untransformed potential. A set of 100 canonically distributed states were prepared in the C7eq basin on the untransformed potential. For a given α, the potential was scaled and the initial states propagated using Langevin dynamics. For each trajectory at a given α, the transition times from C7eq to αR were determined. τ (α) was then calculated as an average over the total number of transitions observed during all 100 simulations. τ (α) was fitted to Eq. 13 and τ (1) and ∆H † determined. For comparison, adaptive umbrella sampling 54 was used to calculate a 2D ψ-φ potential of mean force (PMF) on the untransformed potential. For alanine dipeptide, all simulations of were carried out with CHARMM macromolecular package 55 using the united atom (PARAM19) force field. Initial coordinates were generated using the IC build function in CHARMM and solvent effects were accounted for using the ACE implicit solvation model 56 . For Langevin dynamics the integration timestep was 0.002 9




Page 10 of 23

Page 11 of 23




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

scaled MD simulations (α = 0.5, 0.6, 0.7, 0.8, 0.9). Shown in Fig. 3B is the plot of τ (α) and the fits to the Kramers’ rate equation (Eq. 13). The fit obtained with α = 0.5, 0.6, 0.7 and 0.8 predicts τ (1.0) for the C7eq to αR transition to be 274.1 ps, in excellent agreement with the value of 272.1 ps calculated from reference simulations on the untransformed potential (Fig. 3B). Additionally, the enthalpic component of free energy barrier (∆H † ) associated with the transition was predicted to be 4.57 kcal/mol. By comparison the total free energy barrier (∆U † ) is estimated to be about 3.7 kcal/mol (Fig. 3A). Similar estimates for τ (1.0) and ∆H † , 315.0 ps and 4.8 kcal/mol respectively, were obtained when Eq. 13 was fitted using only τ (0.5) and τ (0.6)(Fig. 3B).

Discussion and Summary We have introduced a simple and computationally efficient method, which enables kinetics to be recovered from potential-scaled enhanced MD simulations. The approach, which is based on a modified Kramers’ rate equation, is particularly attractive because it obviates the need to define reaction coordinates. We applied our approach to two model systems. In the case of the 2D M¨ uller potential, the predicted mean first passage time, τ , which characterized the rate of barrier crossing between two metastable states (Fig. 2A), is in excellent agreement with the reference value calculated from direct simulations on the untransformed potential (Fig. 2B). Likewise, for the alanine dipeptide test case, the predicted τ , which characterized the rate of transition between C7eq and αR , agrees well with the value calculated directly from simulations on the untransformed potential (Fig. 3B). Encouragingly, even when applying our approach using data from potential-scaled MD at a couple low values of α (0.5 and 0.6), we found that reliable estimates of the τ can still be obtained. This is significant because τ at these low α values were at least an order of magnitude shorter than τ on the untransformed potential, indicating that our approach is not only accurate, but also computationally efficient as it

12


Page 12 of 23

Page 13 of 23

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


allows short simulations to be used to study longer timescale processes. Taken together, these results demonstrate the effectiveness of our Kramers-based extrapolation approach in recovering kinetics from potential-scaled MD. The validity of our extrapolation approach hinges on the assumption that motion corresponding to the conformational transition under study is spatially diffusive (i.e., in the strong damping limit). For biomolecules, this is a reasonable assumption since motion tends to be strongly coupled to the “environment”. This observation motivates the application of the high friction, high barrier Kramers’ rate theory with the additional assumption that the rate of conformational transition is determined by a single dominant barrier, i.e. we assume first-order kinetics. Fortunately, many kinetic problems can be formulated using “pseudo” first-order rate equations and so the method can be applied to these problems as well. We note of the fact that our approach to recovering microscopic rates from potentialscaled MD is similar to the approach used by Hamelberg and coworkers to recover rates from accelerated MD – like that used here, their approach is based on Kramers’ rate theory 50,57 . In accelerated MD, a boost potential is added to the local minima regions of the potential 13 , or in the case of a recently describe approach, the boost is added only to regions which fall below a predefined energy 36 . This boosting has the effect of reducing energy barriers, which in turn increases the frequency of transitions between metastable states. To recover kinetics from accelerated MD simulations, Kramers’ rate equation is expressed as a function of the boost energy q

Fo − α+b ωmin ωmax D kK = exp − 2πkB T kB T

!

(15)

where Fo is the free energy barrier height, a, b and α are parameters that determine the magnitude of the boost potential. The magnitude of the boost energy dictates the degree to which sampling is accelerated. In their approach, ωb and ωo are calculated by reweighting the free energy profile along a suitable reaction coordinate. As such, knowledge of topological feature of the free energy landscape is needed. In our approach, the only requirement is a means of partitioning conformational space into reactant and product regions. 13



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

A potential application of our approach is its use in determining drug binding kinetics. Up to this point, the focus of drug discovery efforts tends to be on optimizing the affinity of drugs to their target. Interestingly however, a growing body of evidence now points toward kinetic parameters, in particular “residence time”, i.e. the time a drug stays bound to its target, as being a key parameter that significantly impacts drug efficacy and safety 58–64 . The approach presented here can provide a simple and computationally efficient framework for determining drug residence time from simulations: a series of short potential-scaled simulations can be carried out in which the interaction potential between the drug and target is scaled. Using the Kramers’-based extrapolation method (Eq. 13), these potential-scaled simulations can be used to recover τ , which is equivalent to the residence time for the drug. Future efforts will explore this application further.

Acknowledgement ATF was supported by the National Science Foundation Graduate Fellowship program. IA acknowledges support from the National Science Foundation Career award (CHE-0918817).

References (1) Klippenstein, S. J.; Pande, V. S.; Truhlar, D. G. Chemical Kinetics and Mechanisms of Complex Systems: A Perspective on Recent Theoretical Advances. J. Am. Chem. Soc. 2014, 136, 528–546. (2) McCammon, J. A.; Gelin, B. R.; Karplus, M. Dynamics of Folded Proteins. Nature 1977, 267, 585–590. (3) Deneroff, M. M.; Dror, R. O.; Kuskin, J. S.; Larson, R. H.; Salmon, J. K.; Young, C.; Batson, B.; Bowers, K. J.; Chao, J. C.; Shaw, D. E. Anton, a Special-Purpose Machine for Molecular Dynamics Simulation. Commun. ACM 2008, 51, 91–97. 14


Page 14 of 23

Page 15 of 23

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) Maragakis, P.; Lindorff-Larsen, K.; Piana, S.; Dror, R. O.; Eastwood, M. P.; Bank, J. A.; Jumper, J. M.; Salmon, J. K.; Shan, Y.; Shaw, D. E. Atomic-Level Characterization of the Structural Dynamics of Proteins. Science 2010, 330, 341–346. (5) Piana, S.; Lindorff-Larsen, K.; Shaw, D. E. Protein Folding Kinetics and Thermodynamics from Atomistic Simulation. Proc. Natl. Acad. Sci. USA 2012, 109, 17845–17850. (6) Beberg, A. L.; Ensign, D. L.; Jayachandran, G.; Khaliq, S.; Pande, V. S. Folding@Home: Lessons from Eight Years of Volunteer Distributed Computing. Parallel & Distributed Processing, 2009. IPDPS 2009. IEEE International Symposium on. 2009; pp 1–8. (7) Beauchamp, K. A.; Bowman, G. R.; Lane, T. J.; Maibaum, L.; Haque, I. S.; Pande, V. S. MSMBuilder2: Modeling Conformational Dynamics on the Picosecond to Millisecond Scale. J. Chem. Theory Comput. 2011, 7, 3412–3419. (8) Cronkite-Ratcliff, B.; Pande, V. MSMExplorer: Visualizing Markov State Models for Biomolecule Folding Simulations. Bioinformatics 2013, 29, 950–952. (9) Lane, T. J.; Shukla, D.; Beauchamp, K. A.; Pande, V. S. To Milliseconds and Beyond: Challenges in the Simulation of Protein Folding. Curr. Opin. Struct. Biol. 2013, 23, 58–65. (10) Buch, I.; Harvey, M. J.; Giorgino, T.; Anderson, D. P.; De Fabritiis, G. HighThroughput All-Atom Molecular Dynamics Simulations Using Distributed Computing. J. Chem. Inf. Model. 2010, 50, 397–403. (11) Buch, I.; Giorgino, T.; De Fabritiis, G. Complete Reconstruction of an EnzymeInhibitor Binding Process by Molecular Dynamics Simulations. Proc. Natl. Acad. Sci. USA 2011, 108, 10184–10189. (12) Kumar, S.; Rosenberg, J.; Bouzida, D.; Swendsen, R.; Kollman, P. The Weighted His-

15



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

togram Analysis Method for Free-Energy Calculations on Biomolecules. I. The Method. J. Comput. Chem. 1992, 13, 1011–1021. (13) Voter, A. A Method for Accelerating the Molecular Dynamics Simulation of Infrequent Events. J. Chem. Phys. 1997, 106, 4665–4677. (14) Huber, T.; van Gunsteren, W. SWARM-MD: Searching Conformational Space by Cooperative Molecular Dynamics. J. Phys. Chem. A 1998, 102, 5937–5943. (15) Sugita, Y.; Okamoto, Y. Replica-Exchange Molecular Dynamics Method for Protein Folding. Chem. Phys. Lett. 1999, 314, 141–151. (16) Wu, X.; Brooks, B. Self-Guided Langevin Dynamics Simulation Method. Chem. Phys. Lett. 2003, 381, 512–518. (17) Zheng, W.; Rohrdanz, M. A.; Clementi, C. Rapid Exploration of Configuration Space with Diffusion-Map-Directed Molecular Dynamics. J. Phys. Chem. B 2013, 117, 12769– 12776. (18) Dickson, A.; Brooks III, C. L. WExplore: Hierarchical Exploration of High-Dimensional Spaces Using the Weighted Ensemble Algorithm. J. Phys. Chem. B 2014, 118, 3532– 3542. (19) Xing, C.; Andricioaei, I. On the Calculation of Time Correlation Functions by Potential Scaling. J. Chem. Phys. 2006, 124, 034110. (20) Nummela, J.; Andricioaei, I. Exact Low-Force Kinetics from High-Force Single-Molecule Unfolding Events. Biophys. J. 2007, 93, 3373–3381. (21) Faradjian, A.; Elber, R. Computing Time Scales From Reaction Coordinates by Milestoning. J. Chem. Phys. 2004, 120, 10880.

16


Page 16 of 23

Page 17 of 23

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


(22) Cardenas, A. E.; Elber, R. Modeling Kinetics and Equilibrium of Membranes with Fields: Milestoning Analysis and Implication to Permeation. J. Chem. Phys. 2014, 141, 054101. (23) Dellago, C.; Bolhuis, P.; Csajka, F.; Chandler, D. Transition Path Sampling and the Calculation of Rate Constants. J. Chem. Phys. 1998, 108, 1964. (24) Bolhuis, P. G.; Chandler, D.; Dellago, C.; Geissler, P. L. Transition Path Sampling: Throwing Ropes over Rough Mountain Passes, in the Dark. Annu. Rev. Phys. Chem. 2002, 53, 291–318. (25) MacFadyen, J.; Andricioaei, I. A Skewed-Momenta Method to Efficiently Generate Conformational-Transition Trajectories. J. Chem. Phys. 2005, 123, 074107. (26) Sangha, A. K.; Keyes, T. Proteins Fold by Subdiffusion of the Order Parameter. J. Phys. Chem. B 2009, 113, 15886–15894. (27) Tiwary, P.; Parrinello, M. From Metadynamics to Dynamics. Phys. Rev. Lett. 2013, 111, 230602. (28) Daru, J.; Stirling, A. Divided Saddle Theory: A New Idea for Rate Constant Calculation. J. Chem. Theory Comput. 2014, 10, 1121–1127. (29) Dixit, P. D.; Dill, K. A. Inferring Microscopic Kinetic Rates from Stationary State Distributions. J. Chem. Theory Comput. 2014, 10, 3002–3005. (30) Mark, A. E.; van Gunsteren, W. F.; Berendsen, H. J. Calculation of Relative Free Energy via Indirect Pathways. J. Chem. Phys. 1991, 94, 3808–3816. (31) Tsujishita, H.; Moriguchi, I.; Hirono, S. Potential-Scaled Molecular Dynamics and Potential Annealing: Effective Conformational Search Techniques for Biomolecules. J. Phys. Chem. 1993, 97, 4416–4420.

17



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

(32) Torrie, G. M.; Valleau, J. P. Nonphysical Sampling Distributions in Monte Carlo FreeEnergy Estimation: Umbrella Sampling. J. Comput. Phys. 1977, 23, 187–199. (33) Sinko, W.; Miao, Y.; de Oliveira, C. A. F.; McCammon, J. A. Population Based Reweighting of Scaled Molecular Dynamics. J. Phys. Chem. B 2013, 117, 12759–12768. (34) Kappel, K.; Miao, Y.; McCammon, J. A. Accelerated Molecular Dynamics Simulations of Ligand Binding to a Muscarinic G-Protein-Coupled Receptor. Q. Rev. Biophys. 2015, 48, 479–487. (35) Miao, Y.; Feher, V. A.; McCammon, J. A. Gaussian Accelerated Molecular Dynamics: Unconstrained Enhanced Sampling and Free Energy Calculation. J. Chem. Theory Comput. 2015, 11, 3584–3595. (36) Hamelberg, D.; Mongan, J.; McCammon, J. Accelerated Molecular Dynamics: A Promising and Efficient Simulation Method for Biomolecules. J. Chem. Phys. 2004, 120, 11919. (37) Hamelberg, D.; Shen, T.; McCammon, J. A. Phosphorylation Effects on Cis/Trans Isomerization and the Backbone Conformation of Serine-Proline Motifs: Accelerated Molecular Dynamics Analysis. J. Am. Chem. Soc. 2005, 127, 1969–1974. (38) Hamelberg, D.; Shen, T.; McCammon, J. A. Relating Kinetic Rates and Local Energetic Roughness by Accelerated Molecular-Dynamics Simulations. J. Chem. Phys. 2005, 122, 241103. (39) De Oliveira, C. A. F.; Hamelberg, D.; McCammon, J. A. Estimating Kinetic Rates from Accelerated Molecular Dynamics Simulations: Alanine Dipeptide in Explicit Solvent as a Case Study. J. Chem. Phys. 2007, 127, 175105. (40) Minh, D. D.; Hamelberg, D.; McCammon, J. A. Accelerated Entropy Estimates with Accelerated Dynamics. J. Chem. Phys. 2007, 127, 154105. 18


Page 18 of 23

Page 19 of 23

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


(41) Fajer, M.; Hamelberg, D.; McCammon, J. A. Replica-Exchange Accelerated Molecular Dynamics (REXAMD) Applied to Thermodynamic Integration. J. Chem. Theory Comput. 2008, 4, 1565–1569. (42) Grant, B. J.; Gorfe, A. A.; McCammon, J. A. Ras Conformational Switching: Simulating Nucleotide-Dependent Conformational Transitions with Accelerated Molecular Dynamics. PLoS Comput. Biol. 2009, 5, e1000325. (43) Williams, S. L.; De Oliveira, C. A. F.; McCammon, J. A. Coupling Constant pH Molecular Dynamics with Accelerated Molecular Dynamics. J. Chem. Theory Comput. 2010, 6, 560–568. (44) Wereszczynski, J.; McCammon, J. A. Using Selectively Applied Accelerated Molecular Dynamics to Enhance Free Energy Calculations. J. Chem. Theory Comput. 2010, 6, 3285–3292. (45) Markwick, P. R.; Pierce, L. C.; Goodin, D. B.; McCammon, J. A. Adaptive Accelerated Molecular Dynamics (Ad-AMD) Revealing the Molecular Plasticity of P450cam. J. Phys. Chem. Letters 2011, 2, 158–164. (46) Sinko, W.; de Oliveira, C. A. F.; Pierce, L. C.; McCammon, J. A. Protecting High Energy Barriers: A New Equation to Regulate Boost Energy in Accelerated Molecular Dynamics Simulations. J. Chem. Theory Comput. 2011, 8, 17–23. (47) Ortiz-S´anchez, J. M.; Bucher, D.; Pierce, L. C.; Markwick, P. R.; McCammon, J. A. Exploring the Photophysical Properties of Molecular Systems Using Excited State Accelerated Ab Initio Molecular Dynamics. J. Chem. Theory Comput. 2012, 8, 2752–2761. (48) Lindert, S.; Bucher, D.; Eastman, P.; Pande, V.; McCammon, J. A. Accelerated Molecular Dynamics Simulations with the AMOEBA Polarizable Force Field on Graphics Processing Units. J. Chem. Theory Comput. 2013, 9, 4684–4691.

19



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

(49) Doshi, U.; Hamelberg, D. Reoptimization of the AMBER Force Field Parameters for Peptide Bond (Omega) Torsions Using Accelerated Molecular Dynamics. J. Phys. Chem. B 2009, 113, 16590–16595. (50) Doshi, U.; Hamelberg, D. Extracting Realistic Kinetics of Rare Activated Processes from Accelerated Molecular Dynamics Using Kramers’ Theory. J. Chem. Theory Comput. 2011, 7, 575–581. (51) Doshi, U.; Hamelberg, D. Improved Statistical Sampling and Accuracy with Accelerated Molecular Dynamics on Rotatable Torsions. J. Chem. Theory Comput. 2012, 8, 4004– 4012. (52) Kramers, H. Brownian Motion in a Field of Force and the Diffusion Model of Chemical Reactions. Physica 1940, 7, 284–304. (53) Muller, K.; Brown, L. Location of Saddle Points and Minimum Energy Paths by a Constrained Simplex Optimization Procedure. Theor. Chem. Acc. 1979, 53, 75–93. (54) Bartels, C.; Karplus, M. Multidimensional Adaptive Umbrella Sampling: Applications to Main Chain and Side Chain Peptide Conformations. J. Comput. Chem. 1997, 18, 1450–1462. (55) Brooks, B.; Bruccoleri, R.; Olafson, B.; Karplus, M. CHARMM: A Program for Macromolecular Energy, Minimization, and Dynamics Calculations. J. Comput. Chem. 1983, 4, 187–217. (56) Schaefer, M.; Karplus, M. A Comprehensive Analytical Treatment of Continuum Electrostatics. J. Phys. Chem. 1996, 100, 1578–1599. (57) Xin, Y.; Doshi, U.; Hamelberg, D. Examining the Limits of Time Reweighting and Kramers’ Rate Theory to Obtain Correct Kinetics from Accelerated Molecular Dynamics. J. Chem. Phys. 2010, 132, 224101. 20


Page 20 of 23

Page 21 of 23

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


(58) Copeland, R. A.; Pompliano, D. L.; Meek, T. D. Drug–Target Residence Time and Its Implications for Lead Optimization. Nat. Rev. Drug Discovery 2006, 5, 730–739. (59) Zhang, R.; Monsma, F. Binding Kinetics and Mechanism of Action: Toward the Discovery and Development of Better and Best in Class Drugs. Expert Opin. Drug Discov. 2010, 5, 1023–1029. (60) Holdgate, G. A.; Gill, A. L. Kinetic Efficiency: The Missing Metric for Enhancing Compound Quality? Drug Discov. Today 2011, 16, 910–913. (61) Vilums, M.; Zweemer, A. J. M.; Yu, Z.; de Vries, H.; Hillger, J. M.; Wapenaar, H.; Bollen, I. A. E.; Barmare, F.; Gross, R.; Clemens, J. et al. Structure-Kinetic Relationships–An Overlooked Parameter in Hit-to-Lead Optimization: A Case of Cyclopentylamines as Chemokine Receptor 2 Antagonists. J. Med. Chem. 2013, 56, 7706– 7714. (62) Pan, A. C.; Borhani, D. W.; Dror, R. O.; Shaw, D. E. Molecular Determinants of Drug–Receptor Binding Kinetics. Drug Discov. Today 2013, 18, 667–673. (63) Guo, D.; Hillger, J. M.; IJzerman, A. P.; Heitman, L. H. Drug-Target Residence Time–A Case for G Protein-Coupled Receptors. Med. Res. Rev. 2014, 34, 856–892. (64) Lee, K. S. S.; Liu, J.-Y.; Wagner, K. M.; Pakhomova, S.; Dong, H.; Morisseau, C.; Fu, S. H.; Yang, J.; Wang, P.; Ulu, A. et al. Optimized Inhibitors of Soluble Epoxide Hydrolase Improve in Vitro Target Residence Time and in Vivo Efficacy. J. Med. Chem. 2014, 57, 7016–7030.

21




Page 22 of 23

check$V2[seq(1, 2e+05, 100)]

check$V2[seq(1, 2e+05, 100)]

α-scaled potential surface untransformed potential surface (”short” simulations) (“long” simulations) how that tb (a) F1 F = tb 0 (1) a U U 0

20

40

60

80 100

2 1 0 -1

1 0 -1

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23


2

Page 23 of 23

0 ⇣ t (a) 1 Indexb = 2 exp b DH † (a tb (1) a ACS Paragon Plus Environment

100

⌘200 1) Index .

KRAMERS’ RATE EQUATION

300

400

†

Reaction Coordinate-Free Approach to Recovering Kinetics from

Recommend Documents