Outline Simulation of biomolecular systems ! Basic TPS - PowerPoint PPT Presentation

A dvanced path sampling 
 of rare events in 
 biomolecular systems Peter Bolhuis van ‘t Hoff institute for Molecular Sciences University of Amsterdam, The Netherlands

Outline • Simulation of biomolecular systems ! • Basic TPS – shooting algorithms – stable state definitions – example Photoactive Yellow Protein – reaction coordinate analysis ! • Advanced path sampling – rates by Transition Interface Sampling (TIS) – replica exchange TIS and multiple state TIS – single replica multiple state TIS – example Trp-cage folding network ! • Conclusion

Protein self-assembly ! understanding the cellular processes – folding – structure formation (cytoskeleton) – complex formation (regulation) – neurodegenerative/genetic diseases – ….. ! novel self assembling biomaterials – artificial tissue – smart packaging – self healing coatings – sensors Challenge understanding and predicting protein assembly with advanced molecular simulation

All-atom force fields for biomolecules • Potential energy for protein 1 1 2,3 θ 3 φ r 2 4 4 vdW interactions only between non-bonded |i-j|>4

Currently available empirical force fields • CHARMm (MacKerrel et 96) • AMBER (Cornell et al. 95) • GROMOS (Berendsen et al 87) • OPLS-AA (Jorgensen et al 95) • ENCAD (Levitt et al 83) • ….. ! ! • Subtle differences in improper torsions, scale factors 1-4 bonds, united atom rep. • Partial charges based on empirical fits to small molecular systems • Amber & Charmm also include ab-initio calculations • Not clear which FF is best : top 4 mostly used ! • Water models also included in description – TIP3P , TIP4P – SPC/E • Current limit: 10 6 atoms, microseconds ( with Anton ms)

Molecular Dynamics of proteins MD yields equilibrium statistics: free energy landscapes, • stable structures, transition states, … kinetics: rates, mechanisms, transport • Free energy properties, … FE landscape : high dimensional trajectory can be projected onto collective variable λ e − β F ( λ ) = h δ ( λ ( r ) � λ ) i Conformational order parameter λ

Timescales in proteins Reactions Local flexibility Collective motions Bond Methyl 
 Loop 
 Larger domain 
 vibration rotation motion motions Side-chain rotamer μ s to ms Free energy ms s fs ps ns μ s ns ps • Straightforward MD inefficient Β ! • enhanced sampling: thermodynamic integration, Α umbrella sampling, hyper dynamics, adaptive biasing force, metadynamics, …. 
 Conformational space makes exponential barrier problem linear 
 requires good reaction coordinate

Transition path sampling Importance sampling of the rare event path ensemble: all dynamical trajectories that lead over (high) barrier and connect stable states. Why TPS? ! -selects important rare paths -yields unbiased molecular dynamics PGB, D. Chandler, 
 C. Dellago, P.L. Geissler, 
 -no reaction coordinate needed Annu. Rev. Phys. Chem 2002 -reaction coordinate from committor C. Dellago, PGB, Adv Polym Sci, 2009 -rate constants

Shooting moves accept reject

Flexible one way shooting (TIS) Biased shooting L ∑ S = b ( x i ) i = 0 ⎛ ( o ) ⎞ ( n ) )min 1, L ⎛ ⎞ acc [ x ( o ) → x ( n ) ] = h A ( x 0 ( n ) )min 1, S ( o ) ( n ) ) h B ( x L acc [ x ( o ) → x ( n ) ] = h A ( x 0 P ( n ) ) h B ( x L ⎜ ⎟ P ⎜ ⎟ ( n ) L S ( n ) ⎝ ⎠ ⎝ ⎠ Aimless shooting Spring shooting Δτ P sp acc [ τ → τ 0 ) = min[1 , e sk ∆ τ ] acc [ x ( o ) → x ( n ) ] = h A ( x 0 ( n ) ) h B ( x L ( n ) ) P

Flexible one way shooting ! • Algorithm – Choose new shooting point randomly from old path p sel = 1/L – Do not alter momenta – Integrate in one direction until one stable states is reached – keep old partial path, accept new partial with probability PGB 2003, Juraszek & PGB 2006) ! ! ⎛ ⎞ ( o ) ( n ) )min 1, L acc [ x ( o ) → x ( n ) ] = h A ( x 0 ! ( n ) ) h B ( x L P ⎜ ⎟ ( n ) L ⎝ ⎠ ! densities.30 10b 11b ! 12f 13b 14b 15b 16b 17b ! 18b 19f 20f 21f 22f • higher acceptance, better convergence 23f 24b 25f 26b 27b for diffusive transitions and long pathways 28b 29f 30f 31f 32b 33f • requires some stochastic dynamics, e.g. 34f 35b 36f 37b thermostat 38f 39f 40f 41f 42f 43f • needs check for decorrelation of paths 44b 45f 46f 47b 48b 49f • useful for diffusive (bio)systems 50b 51f 52f 53b

Path sampling indicators Least changed path Path length distribution

Spring shooting for asymmetric barriers • uniform one way shoot has bad decorrelation ! ! 20 ! 10 ! 0 V � x � ! � 10 ! ! bad decorrelation � 20 • spring shooting algorithm: � 5 0 5 x Δτ P sp acc [ τ → τ 0 ) = min[1 , e sk ∆ τ ] pro: much better decorrelation con: need optimisation of k and Δτ max good decorrelation Z.F. Brotzakis, PGB, JCP, in press (2016)

Protein association/dissociation System • beta-lac dimer (2AKQ) • ~65000 atoms • AMBER99SB Τ =300K, P=1 atm • dissociation Δ G=30 kJ/mol • TPS of extremely asymmetric barrier • spring shooting algorithm • spring const 5 kT, dmax = 200 • 50 ps between frames faststore.00 • max path length 50 ns 1b 2b 3b 4b 5f • acceptance 30% 6b 7b 8f 9f 10f 11b Z.F. Brotzakis, PGB, JCP, (2016)

Definition of stable states a) b) basin A basin A basin B B basin B A A B c) d) basin A basin A basin B basin B A B A B

catalysis reactions crystallisation enzyme reactions complex fluids solvent effects folding & binding isomerization

Photoactive yellow protein Ground state pG signal ms-sec fs-ns signal transduction response proteins DNA membrane μ s-ms Signalling state pB Question: What is the mechanism for amplifying signal? Tyr42 ! Glu46 We studied 2 steps: 1) proton transfer chromophore 2) partial unfolding Cys69 Thr50 Photoactive yellow protein

TPS of proton transfer • 28244 atoms • CPMD/QMMM • BLYP functional • Electronic mass 750 au • QM region: pCA, Glu46,Tyr42, 
 Thr50, Arg52 • Gromos96 force field ! • TPS: two way shooting, 
 perturbation temp 35 K • 160 paths/ 50% acceptance • average path length 0.5-1.5 ps • reaction time microseconds stable states pR (reaction) pB’ (product) pCA-Glu46(H) > 1.60 A < 0.98 A OX2-Tyr42 > 3.70 A < 1.80 A OX1-Tyr42 > 5.30 A < 1.80 A

Transition path sampling of partial unfolding Table 1. Statistics of the TPS ensembles. The average path length is a weighted average over the whole ensemble. Decorrelated pathways have lost the memory of the previous decorrelated pathway. The aggregate time is the ensemble aggregate length pB 0 − I α U α − S E U α − S X S E − pB acceptance 41% 25% 38% 44% avg. path length 105 ps 1.8 ns 1.5 ns 1.7 ns accepted paths 3847 305 584 311 decorr. paths 180 18 7 29 aggregate time ( μ s) 1.0 2.3 2.3 1.2 Vreede, Juraszek, PGB, PNAS 2010

Transition states by committor probability that a trajectory initiated at r relaxes into B r B A r is a transition state (TS) if p B (r) = p A (r) =0.5 B A TSE: Intersections of transition pathways with the L. Onsager, Phys. Rev. 54 , 554 (1938). M. M. Klosek, B. J. Matkowsky, Z. Schuss, Ber. Bunsenges. Phys. p B =1/2 surface Chem. 95 , 331 (1991) V. Pande, A. Y. Grosberg, T. Tanaka, E. I. Shaknovich, J. Chem. Phys. 108 , 334 (1998) W.E, E. Vanden-Eijnden, J. Stat.Phys, 123 503 (2006)

Reaction coordinate analysis • Each TPS shot is a committor attempt. 
 Use this information to optimise model of reaction coordinate r ! • The probability that structure x with rc r is on a transition path (for diffusive dynamics) Peters & Trout, JCP 125 054108(2006) ! p ( TP | r ) = 2 p B ( r )(1 − p B ( r )) see also Best & Hummer PNAS (2005) ! • Assume committor function to be 1 p B (r) ! p B ( x ) = 1 2 + 1 2 tanh [ r ( q ( x )] ! 0.8 • parametrize r as linear combination of q 0.6 ! � 0.4 r ( x ) = α i q ( x ) + α 0 p(TP|r) ! i 0.2 • best r is maximizing likelihood 0 N B N A -3 -2 -1 0 1 2 3 p B ( r ( q ( x ( B ) (1 − p B ( r ( q ( x ( B ) � � L ( α ) = )) ))) i i i =1 i =1

Included order parameters


 Reaction coordinate of helix α 3 unfolding Reaction coordinate by likelihood maximization (Peters & Trout, JCP 2006) ! Order Parameters involved (out of 78): RMSD α nwY42 : water molecules around Tyr42 dPA : distance Ala44(N) - Pro54(C γ ) dhb2 : distance Ala44(O) - Asp48(H) δ Lmin = 4.17 n ln L RC 1 -2117 3.89–29.10 × rmsd α 2 -2098 3.88–26.35 × rmsd α − 0.19 × nwY42 3 -2085 5.11–16.81 × rmsd α − 4.68 × dhb2 − 2.55 × dPA

Reaction coordinate pB’ → I α r =5.11 -16.28 rmsd α 3 -4.68 dhb2- 2.55 d PA Committor check: 30% of predicted TSE is a true transition Vreede, Juraszek, PGB, PNAS 2010 state.

Solvent exposure transitions rc = − 2.03 + 2.70 dXE rate limiting step 16 k B T: k ≈ 1 ms -1 TS U α -S E TS U α -S X rc = − 5.05 + 5.02 dXYcom − 2.51 dXEcom + 4.30 dXE

pB’ (a) I α S X U α (c) (b) S E pB Juraszek , Vreede, PGB, Chem Phys 2011