Xin Zhou xzhou@apctp.org

1. Structuring and Sampling in Complex Conformational Space Weighted Ensemble Dynamics Simulation XinZhou xzhou@apctp.org Asia Pacific Center for Theoretical Physics, Dep. of Phys., POSTECH, Pohang, Korea Oct 12, 2009 Beijing

2. Independent Junior Research Group Multiscale Modeling & Simulations in soft materials Members: X.Z. (Leader, 2008.6-) PakpoomReunchan(Postdoctor, 2009.10-) Shun Xu(Ph.D candidate, 2008.11-) Linchen Gong (Ph.D candidate, 2008.12-) Shijing Lu (Ph.D candidate, 2009.1-) One available position for postdoctor/exchange Ph.D student xzhou@apctp.org http://www.apctp.org/jrg/members/xzhou

3. Multiscale simulations: Extend spatial and temporal scales but keep necessary details Improve efficiencies of simulations coarse-graining, enhanced sampling, accelerated slow-dynamics Understanding results of simulations traditional: project to low-dimensional (reaction coordinates) space new : kinetic transition network 1.More sufficient simulation provides more complete understanding 2.The understanding of systems is helpful to design more efficient simulation algorithm

4. multiple scales Vibration of bonds: 10-14 second Protein folding > 10-6second There are coupling among different scales !

5. Barriers Entropic barrier Energetic barrier Due to high free energy (energy and/or entropy) barriers, standard MC/MD simulations need very long time to reach equilibrium Current advanced simulation techniques are not very helpful in overcoming entropic barriers

6. Ensemble Dynamics A. F. Voter (1998) V. S. Pande (2000) n trajectories: transition rate Independently generate multiple short trajectories Statistically analyze slow transition dynamics

7. Weighted Ensemble Dynamics Linchen Gong & X.Z. (2009) arbitrarily select initial conformations Independently generate multiple short trajectories weight the trajectories Statistically analyze slow transition dynamics analyze state structure and equilibrium properties

8. Weighted Ensemble Dynamics A single t-length MD trajectory is not sufficient to reach global equilibrium Multiple t-length trajectories can be used to reproduce global equilibrium properties by reweighting the trajectories Each trajectory has an unique weight (wi) in contributing to equilibrium properties The weight of trajectory is only dependent on its initial conformation

9. Weighted Ensemble Dynamics Usually impractical The initial distribution might be unknown The fluctuation of weights might be too huge to be practice in reproducing equilibrium properties {Wi} satisfies a self-consistent equation for any selected initial conformations

10. Expansion of Probability Density

11. Theory of WED Self-consistent equation: the (short) initial segments of trajectories replace the initial configurations

12. Theory of WED A symmetric linear homogeneous equation:

13. The ground state of H (eigenvector with zero eigenvalue) gives weights of trajectories If the ground state of H is non-degenerate, a unique w is obtained, the equilibrium distribution is reproduced

14. Equilibrium Criterion In principle not practice In practice for any A(x) for complete independent basis functions WED: Judge if simulations reach equilibrium Reweighting trajectories to reach equilibrium distribution parallel simulation from any initial conformations

15. Meta-stable states If the ground state is degenerate, the trajectories are limited in different conformational regions, which are separated each other within the scale of total simulation time: meta-stable states The number of degenerated ground states equals to the number of meta-stable states in the total simulation time scale Simulation trajectories visit in a few completely separated conformation regions, the relative weights of the regions are unknown

16. States and eigenvalues of H Eigenvalue = 0 : separated states in the time scale Eigenvalue = 1 : trajectories in a same state 0<eigenvalue<1 : partially separated states in the time scale a (small) fraction of trajectories happen transitions between states

17. Weights and eigenvector Trajectories are grouped into states Transition trajectories slightly split the degenerate ground states The weights of trajectories in the same state are almost constant

18. Projection in ground states S1 (1.75) : 77 S2 (0.75) : 92 S3 (-0.25) : 37 S4 (-1.25) : 66 S4-S3-S2 : 9 S4-S3 : 119 Non-transition trajectories transition trajectories 1. Non-transition trajectories inside a state are mapped to the same point 2. Transition trajectories between two states are mapped to the line connected by the states 3. Transition trajectories among three states are mapped to the plane of the states

19. Occupation fraction vs projection Multi-time transition trajectories single-time transition trajectories The occupation fraction of a trajectory in states is linearly related to its projection

20. Transition time vs projection single-time transition trajectories Transition state ensemble Without requiring knowledge of states and transitions

21. Free energy reconstruction Two different initial distributions are re-weighted to the accurate free energy profile Weights of trajectories started from the same state are almost same

22. Transition network in 2D multi-well potential Topology of transition network is kept

23. Mexico-hat: entropy effects Topology of transition network is kept Eigenvalues gradually increase from 0 to 1

24. Alaninedipeptide in waters An alaninedipeptide solvated in 522 TIP3P water molecules: 1588 atoms 500 initial conformations, generated from a 10 ns simulation at T=600K 500 WED trajectories (600 ps each)

25. Eigenvalue of H

26. 450K and 300K Started from Ceq7 Started from Alpha7 Started from Cax7 Initial psi

27. Modified potential at 300K projection Free energy reconstruction

28. Occupation fraction vs projection (300K mod) Single-time transition trajectories

29. Transition time vs projection (300K mod) Single-time transition trajectories Real transition time by checking along the single-time transition trajectories

30. Dipeptide at 150K Eigenvalues of H continuously increase from zero to unity at 150K: Inter-trajectory difference due to entropy effects makes multiple eigenvalues be smaller than (but close to) unity

31. 150K More dispersive projection Count of traj. 300K modified

32. Diffusive dynamics at 150K Histogram at transition regions is significant Diffusively cross the transition regions Trajectories do not sufficiently cover whole the state at the low temperature Statistical difference between trajectories is large

33. Solvent effects Include solvent-related functions in expansion

34. Completeness of Basis functions It does not require the expansion is accurate at everywhere, but distinguish conformational regions S quickly reaches saturation while the number of basis functions is far smaller than the size of sample

35. More complex cases While there are n small eigenvalues, trajectories should be projected to n-1 dimensional space Cluster analysis is required Trajectories might need to be split into multiple shorter segments to distinguish transition and non-transition trajectories Meta-stable states can be clustered in different time scales

36. Generalization i=1,…,n trajectories generated from the same potential but different initial configurations, the distributions are denoted as Pi(x) The overlap matrix of trajectories Each trajectory (set of conformations) is mapped to a point

37. mapping 0 1 If two samples come from the same distribution, their mapped points locate at the same position The error satisfies a Gaussian distribution (the center limit theorem)

38. Trajectory mapping t-length MD trajectory: Inside a state and reach local equilibrium Transition among a few states and reach local equilibrium in each of the states, but not reach the inter-state equilibrium Inside a state (conformational region) but not reach local equilibrium 0 1 Equilibrium distribution of meta-stable states

39. Trajectory clustering t-length MD trajectory: Concentrated points (clusters): non-transition 2. Points on lines connected with clusters: transition 3. Diffusion dynamics in entropy-dominated regions

40. Dimension of manifold nd is the dimension of basis functions trajectories ns is number of meta-stable states It is an equality while the set of applied basis functions is sufficient

41. Hierarchic kinetic network Split trajectory into short segments: detect kinetic network in shorter time scales Hierarchic meta-stable state structure

42. Understanding of system The trajectory ensemble is mapped as a trajectory Fraction in states: correlation: Form hierarchic Kinetic network which involving complete equilibrium and transition kinetic/dynamical properties 2. Calculate weights of trajectories and and correlation in diffusion regions

43. Example 74 atoms, charged terminals; Implicit solvent simulation: Generalized Born; 1000 trajectories; 20 ns and 40000 conformationsper trajectory.

44. Application 12-alanine peptide 1000 20-ns trajectories 172 basis functions from torsion angles

45. Dimension of manifold Principle Component Analysis Total dimension: 172

46. clustering minimal spanning tree clustering algorithm 2

47. sub-trajectory clustering (1 ns)

48. G1, 1ns clustering

49. 2D might be insufficient G1, 2D free energy profile

50. G11, sub-trajectory clustering