Erroneous Rates and False Statistical Confirmations from Infrequent

Subscriber access provided by University of South Dakota

Statistical Mechanics

Erroneous rates and false statistical confirmations from infrequent metadynamics and other equivalent violations of the hyperdynamics paradigm Bradley M Dickson J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.8b00848 • Publication Date (Web): 23 Nov 2018 Downloaded from http://pubs.acs.org on November 25, 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 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

Journal of Chemical Theory and Computation

Erroneous rates and false statistical confirmations from infrequent metadynamics and other equivalent violations of the hyperdynamics paradigm Bradley M. Dickson∗ Center for Epigenetics, Van Andel Research Institute, Grand Rapids Michigan 49503, USA. E-mail: [email protected]

Abstract Recent combinations of hyperdynamics and adaptive biasing potentials or free energy based bias potentials introduce an implicit assumption: The bias potential need only be zero on the prominent free energy barriers of interest. Herein we demonstrate that this implicit assumption is flawed. Thus, hyperdynamics in collective variables is likely to fail unless there are no hidden barriers within the initial state. Moreoever, the one-sample Kolmogorov-Smirnov test used to declare “trust” in hyperdynamics is shown unable to classify this failure. In fact, non-zero bias potential on hidden barriers emerges as a contradiction of the hypothesis that a Poisson distributed set of hyperdynamics escape times is a correctly distributed set of escape times. We demonstrate failure in the alanine dipeptide benchmark and we present an apparent failure for a protein conformation change. The standards for “trustworthy” hyperdynamics in collective variables must be raised. ∗

To whom correspondence should be addressed

1

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 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

1

Introduction

Two decades ago, hyperdynamics 1 was introduced as one way to bridge simulation and experimental timescales. The hyperdynamics method supplements the underlying physical potential energy with a bias potential. In order that hyperdynamics correctly maps the biased simulation timescale to the unbiased one, several conditions must be obeyed. Quoting directly from Reference 1, the bias potential “must be zero at all the dividing surfaces, must not block rapid ergodic sampling within a state, and must not introduce TST-violating correlations among the transitions.” These conditions on the bias potential allow a relationship between the biased and unbiased timescales that does not require knowledge of the Transition State Theory (TST) dividing surfaces that separate stable states of the potential energy landscape. The relationship between timescales is given by the hyperdynamics boost factor. Relatively recently, hyperdynamics has been combined with adaptive biasing potentials. 2–4 Adaptive biasing potentials were first introduced as sampling enhancers that could also be used for estimation of free energy. 5 Free energy is computed from a projection of the partition function onto a low-dimensional set of collective variables or order parameters. Thus, when hyperdynamics is combined with an adaptive bias potential one finds a context that was not considered when hyperdynamcis was formulated: The bias potential is built on a low dimensional free energy landscape designed to capture coarse-grained conformational features rather than features of the many dimensional potential energy landscape (e.g., the position of dividing surfaces). Thus, the results of this study also apply to any hyperdynamics bias potential that is constructed from a free energy projection. 6 The solution to the problem lies in collective variable construction. For the correct type of collective variable free energy based bias potentials will be consistent with hyperdynamics. Examples of possible collective variables are given in the Conclusion section, none of which have been applied to protein simulation. In what follows we consider specifically whether the rule “the bias potential is zero at all the dividing surfaces” applies only to the prominent barriers on the free energy landscape 2


Page 2 of 19

Page 3 of 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


that are associated with the rare event being studied or whether the rule is more strict and must be enforced on barriers orthogonal to the transition of interest. Indeed, it has been taken as an implicit assumption in all literature combining hyperdynamics and free energy that the bias must only be zero on the barriers prominent in the free energy projection. Unfortunately, the situation turns out to be rather grim and we can not advocate the general adoption of the common implicit assumption. For the trivial alanine dipeptide case we demonstrate failure of hyperdynamics that can be attributed to non-zero bias on barriers orthogonal to the prominent barrier associated with the reaction timescale being estimated. Thus, the hyperdynamics boost factor is an inaccurate map between biased and unbiased times. We also demonstrate an apparent failure for application to a protein conformation change, where 6 microseconds of dynamics were collected in search of accuracy.

2

Results

We recently introduced a variation of the mollified adaptive biasing potential 7 (mABP) that will only fill free energy minima below a preset threshold and we combined that “filllimited” bias potential with hyperdyanmics. 4 We compared our bias potential to “infrequent metadynamics” 3 for alanine dipeptide and for a ligand-receptor problem. We showed that hyperdynamics was applicable in each of these problems for specific bias types and filllimits. For the protein-ligand system, we forced the ligand to escape from the bound state directly to solvent by placing a cylindrical restraint over the binding pocket of the protein. We did this to reduce the complexity of the system and to avoid, as much as possible, “hidden” or, parallel valleys 8 in the collective variables. This requires trial and error and is unlikely a viable strategy for arbitrary systems. Our bias potential and hyperdynamics are implemented in fABMACS 9 and all input files for the following simulations can be obtained from the fABMACS page on GitHub. 10 Here we revisit alanine dipeptide to understand why hyperdynamics fails for some fill-

3



A)

B)

80

25

6

10 kJ/mol fill-level

70

20

5

A( )

4

State B

State A

50

3

2

40

30

20

kJ/mol

60

A(φ,ψ) (kJ/mol)

V b( )

15

A( )+Vb( ) 10

5

1 10

0

0

0 0

1

2

3

4

5

6

0

φ

C)

D)

30

1

2

3

4

φ

5

6

22 20

25

18 16

kJ/mol

20

kJ/mol

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 4 of 19

15 10

0

0

1

2

3

4

12 10 8

unbiased 1D bias 2D bias

5

14


6 4 5

2

6

5

5.5

6

6.5

7

φ

Figure 1: A) Alanine dipeptide free energy in the φ, ψ dihedral angles (radians). B) 1dimensional free energy with φ as the collective variable(black), bias potential Vb (φ) with 10 kJ/mol fill-level (green), A(φ) + Vb (φ) landscape (blue). C) Profile of the 2D A(φ, ψ) and A + Vb along the dashed line in panel A, where bias potentials have a 10 kJ/mol fill-level. D) Profile of the 2D A(φ, ψ) and A + Vb along the solid line in panel A, where bias potentials have a 10 kJ/mol fill-level.

4


7.5

Page 5 of 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


limits but not others. In panel A of figure 1 we show the free energy surface for rotation of the φ, ψ dihedral angles of alanine dipeptide. The initial state, state A, and final state, state B, are also define consistent with previous definitions. 4 At 300 Kelvin, the transition from state A to state B is predominantly controlled by a pathway with a barrier height of 21.5 kJ/mol. This barrier lies along the φ direction in the 2 dimensional CV space. This barrier accounts for 90% of transitions in unbiased dynamics. The prominent pathway is highlighted by a solid line in figure 1-A. This pathway is one of several that separate state B from state A, but all the pathways require a barrier crossing in the φ direction. We will call the barriers separating states A and B “φ-associated.” Within state A, there are also barriers that lie along the ψ direction. Importantly, these ψ-associated barriers separate locally metastable populations that have different φ-associated barriers (or pathways) for exiting state A. Proper sampling within state A might, therefore, be required for a trajectory to find the reaction pathway with the highest flux. Both the A and B states are “superbasins” in the sense that each state contains metastable sub-states. The free energy computed with a 1-dimensional collective variable, φ, is plotted in figure 1-B (black). Additionally, the 1D bias potential along φ is plotted (green) along with the landscape sampled by dynamics driven by this bias potential (blue), where the bias had a 10 kJ/mol fill-level. Figure 1-B clearly depicts a situation wherein the hyperdynamics potential (green) is zero on the φ-associated TST dividing surfaces separating state A from state B. In other words, the bias vanishes on the prominent free energy barriers. However, the 1D bias potential does not vanish in the dividing surfaces of ψ-associated barriers, as evidenced in figure 1-C (purple), in stark contrast to the 2D bias with the same fill-level (green). Figure 1-D demonstrates that both the 1D and 2D bias potentials are zero in φ-associated barriers. The crux of the situation is whether hyperdynamics will remain correct even when the bias potential is non-zero on the ψ-associated barriers. Figure 2 reports the results of combining hyperdynamics with the 1D and 2D bias potentials. The figure also reports the outcome of a two-sample Kolmogorov-Smirnov test (KST)

5



that compares the unbiased escape time samples to the samples from biased simulations. Asterisk indicate a failure (P-value < 0.05). These results demonstrate that a 10 kJ/mol fill-level is perfectly tolerated if the bias potential is built on the φ, ψ CV-space but is not tolerated if the bias potential is built on only the φ CV-space. Results are summarized in tables 1 and 2 of the Appendix. We note that the two-sample KST can only be used when unbiased escape samples are available and is distinct from the one-sample test advocated in Reference 11. Figure 2 also reports results from infrequent metadynamics, summarized in table 3, showing that both bias styles fail at low boosts with a 1D collective variable. The boost, reported in tables 1 and 2, is significantly larger for all fill-levels if the bias is built on the 1D CV. This follows from the fact that the 2D bias is zero in ψ-associated barriers while the 1D bias is not. The 1D and 2D bias potentials look very similar in the φ, ψ landscape for the φ-associated barrier (figure 1-D) but only the 2D bias produces correct escape time estimates. Consistent with our previous findings, when the bias potential was defined on the 2D φ, ψ-space, hyperdynamcis performs well for biasing depths that were limited to 15 kJ/mol. At a depth of 17 kJ/mol we saw a systematic elevation of the estimated time to escape state A but we could not statistically rule out the possibility that 17 kJ/mol was safe for use with hyperdynamics. 4 In this study, after increasing the number of escape time samples to 500 (see table 1), we again can not conclude that 17 kJ/mol is inconsistent with the unbiased escape time. Interestingly, the ψ-associated barriers within state A are actually around 15 to 17 kJ/mol in height, suggesting that hyperdynamics fails in a 2D collective variable when the bias is non-zero in a dividing surface — At a 17 kJ/mol fill-level, the bias partially floods over some of the tops of barriers where the free energy is less than 17 kJ/mol. (See figure 1-A and C for heights.) Once the bias potential floods higher than the barriers within the superbasin A, hyperdynamics begins to break down. The hyperdynamics escape times from the 1D bias potential could be erroneous for either of two reasons: First, the dynamical corrections could be corrupt. Second, the boost factor

6


Page 6 of 19

Page 7 of 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


could be wrong for transitions over ψ-associated barriers because the bias potential does not vanish on those barriers. During our hyperdynamics runs, all segments of the trajectories that fell outside of states A and B from figure 1-A were recorded. To investigate the structure of the trajectory near the A-to-B dividing surface, we computed the average number of crossings of the plane at φ = 0.1 radians which is close to the free energy barrier between states A and B in the most visited reaction pathway. Reactive trajectories crossed this plane an average of 4.8±0.48 times per reactive trajectory while biased in 2-dimensions with a fill-level of 10 kJ/mol, while they crossed the plane an average of 4.22±0.35 times per trajectory while being biased in 1-dimension with a fill-level of 10 kJ/mol. Unbiased dynamics crossed the plane an average of 4.49±0.46 times per reactive trajectory. The crossing structure of these trajectories is the same, regardless of which bias was used. The above analysis indicates that hyperdynamics with a 1 dimensional collective variable fails to produce correct escape time estimates because the boost factor is not an accurate map between the biased and unbiased timescales. The boost factor over estimates the acceleration of the escapes because it is non-zero in the dividing surface of barriers orthogonal to the φcollective variable. Bear in mind that the 1D collective variable does capture the barrier of interest, rotation of φ, and is zero on that dividing surface. Now we can also appreciate that the 2D bias potential introduces hyperdynamics violations when it floods over the ψassociated barriers within state A. The draw-back of the 1D bias potential is that all fill-levels produce non-zero bias on those ψ-associated barriers, while the 2D bias accommodates the hyperdynamics tenet that “the bias potential is zero at all the dividing surfaces” for several fill-levels. This problem with a 1D collective variable is independent of whether the infrequent metadynamics or the mABP scheme is used. If one wished to compute the timescale of rotating ψ within state A using hyperdynamics, then obviously the bias potential should vanish on the ψ-associated barrier(s). One could therefore say that the timescale of motions inside state A require the bias potential to vanish

7



on the ψ-associated barrier(s). One might then ask, could the timescale of exiting state A be estimated correctly without accurate knowledge of the timescale of motions within state A? The answer, based on the above findings, is no. Going forward with hyperdynamics in collective variable spaces, this indicates that our attention should focus on resolving metastability within “state A.” Generally, any barrier large enough to introduce a metastability that gives rise to a well-defined TST rate must be respected by the bias potential. The problem of hidden barriers should be exacerbated as the dimensionality of configuration space and collective variable space diverge. For example, one should expect a significant number of hidden barriers when protein (un)folding is described by a single collective variable. Hyperdynamics, in the form of infrequent metadynamics, has been seen to fail for exactly this case. 12 Notice that the one-sample KST for trustworthiness 11 fails to identify spoiled hyperdynamics data in tables 1-3 and in Reference 12. The escape times from the 1D bias potential are exponentially distributed, but the boost factor is incorrect due to hidden barriers. The hidden barrier problem (and the non-zero bias in those dividing surfaces) demonstrate that the one-sample KST can not indicate whether hyperdynamics has failed in general. In order to demonstrate how hidden barriers might play out in an application to a protein system, we sought to compute the escape time for a conformation change in the protein UHRF1 (Ubiquitin Like with PHD and Ring Finger domain-containing protein 1). This protein is thought to undergo a conformational change in response to certain ligands. 13–15 The histone binding cassette of UHRF1 is comprised of two domains, the tandem tudor (TTD) and plant homeodomain (PHD), that are connected by a disordered linker. 16,17 Arginine 296 of the linker is found buried in the TTD, resulting in coupled binding action of the TTD and PHD. The histone binding cassette is shown in figure 3-A. Figure 3-B illustrates the initial state for our simulations, where hyperdynamics is used to compute the average residence time of the linker-R296 bound state. The product state for this calculation is represented by any state wherein the distance between the central C atom

8


Page 8 of 19

Page 9 of 19

20

*

18

Escape time (ns)

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


mABP-1D mABP-2D IM-1D

16

14

* 12

*

10

**

p=0.055 8

unbiased

*

*

p=0.175

6

4 1

10

100

Boost Figure 2: Escape time estimates from unbiased dynamics (black), mABP with 1D bias (purple), mABP with 2D bias (green), and infrequent metadynamics (blue). Asterisks indicate inconsistency between sample data and unbiased data in a two-sample KST at 95% confidence interval. All P-values are reported in tables 1-3. of the R296 guanidine group is at least 9 angstrom from the Cβ atom of D142 for at least 50000 force calls. One such state is shown in figure 3-C. The equilibrated GROMACS files, topology, and fABMACS inputs are available within the RErun directory of the fABMACS project. 10 For the hyperdynamics bias we used two collective variables. The first is an RMSD between the guanidine group of R296 referenced to the crystallographic pose. The second collective variable is an RMSD between the backbone C-N-Cα atoms of the R296 residue and its crystallographic pose. As indicated in figure 3-B and C, we removed the PHD domain to reduce computational burden and complexity of these simulations. This transition is similar to an event along an unfolding pathway. Figure 4 summarizes our hyperdynamics escape time estimates for the escape of R296 from the TTD. Table 4 reports the fill-levels, average boosts, number of samples, and onesample KST results. Importantly, the one-sample KST that is widely used to indicate trustworthy dynamics does not indicate that hyperdynamics is not convergent for the set of tested fill-levels. In fact, having done any one of these calculations would suffice as 9



A) TTD

B)

Figure 3: A) The histone binding cassette of UHRF1. The TTD and PHD in white surfaces, the interdomain linker in red, the histone N-terminus in green. B) The initial state of hyperdynamics simulations where R296 and D142 are shown as sticks. C) An example of the product state for hyperdynamics where R296 is ejected from the binding pocket.

10


Page 10 of 19

Page 11 of 19

“trustworthy” given the current standard of querying only the one-sample KST. (See table 4 for P-values and boosts.) It is only after comparing calculations at different fill-levels, with boosts different by an order of magnitude, that we see the lack of convergence in hyperdynamics times. In direct simulation time, figure 4 represents 6.2 microseconds of combined trajectories, with 1.9 microseconds being collected at fill-level 52 kJ/mol alone. A great deal of computational resources were used to generate figure 4 as the intention was to obtain accurate estimates of the timescale for this conformation change. This apparent, and costly, failure of hyperdynamics is what lead to the analysis of the hyperdynamics assumptions presented above. 100

Escape time (s)

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


10

1

0.1 50

52

level (kJ/mol) Figure 4: Escape times estimated by hyperdynamics at different fill-levels for R296 escape.

One is inclined to think that at some small fill-limit hyperdynamics should converge, allowing correct estimation of the timescale for linker unbinding. However, there are several rearrangements of hydrogen bonding networks among the residues immediately adjacent to R296 that represent hidden barriers analogous to the ψ-associated barriers of alanine dipeptide. In order to get the state-to-state dynamics correct, which is constituted by changes 11



in hydrogen bonding partners during linker dissociation, the bias needs to vanish in the dividing surface of each change. There are, in fact, several hydrogen bonding interactions between R295 of the linker and D190, E153, and D189 of the TTD within the R296-bound state. In the many dimensional configuration space the interaction partners of R295 (and others) resolve the R296 unbinding pathway into different saddle points. Two pathways of linker unbinding (at 52 kJ/mol filllevel) are shown in figure 5. Simulations started from the crystallographic configuration quickly settle into a state where the R296 backbone amine hydrogen bonds with E153 and the guanidine of R295 hydrogen bonds with D190. From this state, R295 can move to engage either E153 (figure 5-(2)) or D189 (figure 5-(2’)). The next transition bears out a major divergence of mechanisms. When figure 5-(3) and figure 5-(3’) are compared, the most notable difference is the position of R295. In figure 5-(3), R295 is interacting with D190, where R295 is solvated in figure 5-(3’). For the lower pathway in figure 5 there is only one remaining transition, dissociation of the R295-D142 salt-bridge, before the linker is free from the TTD. For the upper pathway, the sidechain of R296 rotates so that the guanidine and D142 are not directly interacting. Finally, for the upper pathway, R295 looses contact with D190 and the linker escapes the TTD. Figure 5 shows two pathways of escape that are separated in the configuration space by barriers that are hidden from the collective variable. The lack of convergence in figure 4 can be understood as resulting from violations of the core assumptions of hyperdynamics. Namely, the bias potential “must be zero at all the dividing surfaces.” We stress that at a 95% confidence interval, all of the TTD-linker data in figure 4 passes the self-reference one-sample KST 11 that is used to determine trust in hyperdyanmics times. Without understanding that non-zero bias on hidden barrier tops is a violation of the basic assumptions of hyperdynamics, it would be possible to claim a linker escape time at roughly one escape per 7 seconds at a fill-level of 60 kJ/mol. This is clearly incorrect in light of the above discussion.

12


Page 12 of 19

Page 13 of 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


2) R

3

D

R

4

D

R

D

E D

D

1)

R

R

R

D

R

D E

2

3 D

D R

D R

E E

D

R

D

R R

Figure 5: Two pathways from hyperdynamics simulations with a fill-level of 52 kJ/mol. The two paths demonstrate very different hydrogen bonding networks. Each change in the network represents crossing a hidden barrier, a barrier that has non-zero bias on the dividing surface.

3

Conclusion

In summary, the hyperdynamics bias potential must be zero in all dividing surfaces to produce accurate dynamics. This requirement is particularly difficult to guarantee when the bias potential is built on a low-dimensional collective variable space. Thus, a new “pitfall” of collective variable selection emerges in the context of hyperdyanmics. Prior to this work it has been acceptable to implicitly assume hyperdynamics is valid so long as the bias potential vanished on the prominent barriers of free energy. We have shown this not to be the case. We have also shown that there is no statistical test in use that can correctly deem a single set of hyperdynamics escape times as trustworthy. The problem of hidden barriers is, in fact, a logical argument against the hypothesis that a one-sample KST would be a means of testing hyperdynamics. It is, therefore, still true that generally applicable hyperdynamics bias potentials are nontrivial to construct. 18 One of the main challenges to applying hyperdynamics in collective variable spaces will be to resolve all metastability in the state to be escaped. A collective variable that returns the smallest distance to a dividing surface would isolate the current basin of attraction and

13



eliminate accedental construction of superbasins. Machine learning has been used to estimate the location of the dividing surface, allowing a hyperdynamics bias potential that is defined on the smallest distance to the ridge. 19 This type of collective variable has not been used in the study of proteins and while it would produce a bias potential consistent with the assumptions of hyperdynamics, there is no guarantee for efficiently accelerated sampling. A path-based collective variable has also been introduced in hyperdynamics 20 and might also be useful for isolating superbasins. The local-hyperdynamics method 21 will also likely be a key ingredient for protein systems because it allows different zones of the physical system to be decoupled in the bias. Several alternative techniques exist for computing rates of transition. 22–26 Until a hyperdynamics-consistent collective variable is adapted for application to protein systems, these alternative techniques should be employed when a rate must be computed.

4

Acknowledgements

All computations were performed on the Van Andel Research Institute compute cluster.

References (1) Voter, A. F. Hyperdynamics: Accelerated molecular dynamics of infrequent events. Phys. Rev. Lett. 1997, 78, 3908. (2) Ishii, A.; Ogata, S.; Kimizuka, H.; Li, J. Adaptive-boost molecular dynamics simulation of carbon diffusion in iron. Phys. Rev. B 2012, 85, 064303. (3) Tiwary, P.; Parrinello, M. From metadynamics to dynamics. Phys. Rev. Lett. 2013, 111, 230602. (4) Dickson, B. M. Overfill protection and hyperdynamics in adaptively biased simulations. J. Chem. Theory Comput. 2017, 13, 5925–5932. (5) Laio, A.; Parrinello, M. Escaping free-energy minima. PNAS 2002, 99, 12562–12566. 14


Page 14 of 19

Page 15 of 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


(6) McCarty, J.; Valsson, O.; Tiwary, P.; Parrinello, M. Variationally optimized free-energy flooding for rate calculation. Phys. Rev. Lett. 2015, 115, 070601. (7) Dickson, B. M. Approaching a parameter-free metadynamics. Phys. Rev. E 2011, 84, 037701–037705. (8) Minoukadeh, K.; Chipot, C.; Lelievre, T. Potential of mean force calculations: a multiple-walker adaptive biasing force approach. J. Chem. Theory Comput. 2010, 6, 1008–1017. (9) Dickson, B. M.; de Waal, P. W.; Ramjan, Z. H.; Xu, H. E.; Rothbart, S. B. A fast, open source implementation of adaptive biasing potentials uncovers a ligand design strategy for the chromatin regulator BRD4. J. Chem. Phys. 2016, 145, 154113. (10) Dickson, B. Author’s github page. https://github.com/BradleyDickson/fABMACS 2016, (11) Salvalaglio, M.; Tiwary, P.; Parrinello, M. Assessing the reliability of the dynamics reconstructed from metadynamics. J. Chem. Theory Comput. 2014, 10, 1420–1425. (12) Tung, H.-J.; Pfaendtner, J. Kinetics and mechanism of ionic-liquid induced protein unfolding: application to the model protein HP35. Mol. Syst. Des. Eng. 2016, 1, 382– 390. (13) Gelato, K. A.; Tauber, M.; Ong, M. S.; Winter, S.; Hiragami-Hamada, K.; Sindlinger, J.; Lemak, A.; Bultsma, Y.; Houliston, S.; Schwarzer, D. with others. Accessibility of different histone H3-binding domains of UHRF1 is allosterically regulated by phosphatidylinositol 5-phosphate. Mol. Cell 2014, 54, 905–919. (14) Houliston, R. S.; Lemak, A.; Iqbal, A.; Ivanochko, D.; Duan, S.; Kaustov, L.; Ong, M. S.; Fan, L.; Senisterra, G.; Brown, P. J. Conformational dynamics of the

15



TTD-PHD histone reader module of UHRF1 reveals multiple histone binding states, allosteric regulation and druggability. J. Biol. Chem. 2017, jbc–M117. (15) Ferry, L.; Fournier, A.; Tsusaka, T.; Adelmant, G.; Shimazu, T.; Matano, S.; Kirsh, O.; Amouroux, R.; Dohmae, N.; Suzuki, T. Methylation of DNA ligase 1 by G9a/GLP recruits UHRF1 to replicating DNA and regulates DNA methylation. Mol. Cell 2017, 67, 550–565. (16) Arita, K.; Isogai, S.; Oda, T.; Unoki, M.; Sugita, K.; Sekiyama, N.; Kuwata, K.; Hamamoto, R.; Tochio, H.; Sato, M. Recognition of modification status on a histone H3 tail by linked histone reader modules of the epigenetic regulator UHRF1. Proc. Natl Acad. Sci. U.S.A. 2012, 109, 12950–12955. (17) Rothbart, S. B.; Dickson, B. M.; Ong, M. S.; Krajewski, K.; Houliston, S.; Kireev, D. B.; Arrowsmith, C. H.; Strahl, B. D. Multivalent histone engagement by the linked tandem Tudor and PHD domains of UHRF1 is required for the epigenetic inheritance of DNA methylation. Genes Dev. 2013, 27, 1288–1298. (18) Henkelman, G.; J´onsson, H. Long time scale kinetic Monte Carlo simulations without lattice approximation and predefined event table. J. Chem. Phys. 2001, 115, 9657– 9666. (19) Xiao, P.; Duncan, J.; Zhang, L.; Henkelman, G. Ridge-based bias potentials to accelerate molecular dynamics. J. Chem. Phys. 2015, 143, 244104. (20) Huang, C.; Perez, D.; Voter, A. F. Hyperdynamics boost factor achievable with an ideal bias potential. J. Chem. Phys. 2015, 143, 074113. (21) Kim, S. Y.; Perez, D.; Voter, A. F. Local hyperdynamics. J. Chem. Phys. 2013, 139, 144110.

16


Page 16 of 19

Page 17 of 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


(22) Weinan, E.; Ren, W.; Vanden-Eijnden, E. Finite temperature string method for the study of rare events. J. Phys. Chem. B 2005, 109, 6688–6693. (23) Vanden-Eijnden, E.; Tal, F. A. Transition state theory: Variational formulation, dynamical corrections, and error estimates. J. Chem. Phys. 2005, 123, 184103. (24) Peters, B. Reaction coordinates and mechanistic hypothesis tests. Annu. Rev. Phys. Chem. 2016, 67, 669–690. (25) Elber, R. A new paradigm for atomically detailed simulations of kinetics in biophysical systems. Q. Rev. Biophys. 2017, 50 . (26) Peters, B. Reaction rate theory and rare events; Elsevier, 2017.

5

Appendix: Tables of hyperdynamics data

Table 1: Table of 2D CV results: P-values using either a self-referenced KST or using a twosample KST against the unbiased sample set. flim = 0 corresponds to unbiased dynamics. flim

self-KST

2-sample KST

Boost

Samples

0

0.805

1

1

375

10

0.471

0.271

5.01

300

15

0.700

0.964

18.29

300

17

0.699

0.175

31.48

500

18

0.993

0.024

41.18

300

20

0.873

0.00

73.69

200

17



Table 2: Table of 1D CV results: P-values using either a self-referenced KST or using a two-sample KST against the unbiased sample set. flim

self-KST

2-sample KST

Boost

Samples

5

0.495

0.723

2.85

352

10

0.687

0.038

12.43

390

11

0.820

0.00

16.63

465

12

0.561

0.00

22.24

500

15

0.937

0.00

54.98

100

Table 3: Table of 1D CV results for “infrequent metadynamics”: P-values using either a self-referenced KST or using a two-sample KST against the unbiased sample set. “Steps” indicates the number of force-calls between bias updates. Steps

self-KST

2-sample KST

Boost

Samples

1000

0.440

0.001

3.71

300

5000

0.112

0.055

10.8

428

Table 4: Table of TTDlinker unbinding results: Fill-level (kJ/mol), P-values using selfreferenced one-sample KST, boost, and sample number. flim

self-KST

Boost

Samples

52

0.748

9 × 106

20

54

0.853

22 × 106

20

56

0.654

43 × 106

17

58

0.683

96 × 106

20

60

0.297

202 × 106

20

18


Page 18 of 19

Page 19 of 19

Dashed path

Solid path

30

80

22

6

20 25

18

60

3

2

40

30

20

16

20

kJ/mol

State A

50

kJ/mol

4

State B

uninteresting motion

70 5

A( ψ) (kJ/mol)

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


15 10

10

0

0

0 0

1

2

3

4

5

6

0

1

interesting motion

2

3

4

uninteresting motion

Hyperdynamics boost = biased =

KTST biased

12 10 8


5

1

14


6 4 5

6

2

5

5.5

6

Hyperdynamics boost = biased ≠

KTST

Figure 6: TOC image

19


6.5

interesting motion

KTST biased KTST

7

7.5

Erroneous Rates and False Statistical Confirmations from Infrequent

Recommend Documents