530 likes | 695 Views
Structuring and Sampling in Complex Conformational Space Weighted Ensemble Dynamics Simulation . Xin Zhou xzhou@apctp.org. Asia Pacific Center for Theoretical Physics, Dep. o f Phys., POSTECH , Pohang , Korea . Oct 12, 2009 Beijing. Independent Junior Research Group .
E N D
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
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
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
multiple scales Vibration of bonds: 10-14 second Protein folding > 10-6second There are coupling among different scales !
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
Ensemble Dynamics A. F. Voter (1998) V. S. Pande (2000) n trajectories: transition rate Independently generate multiple short trajectories Statistically analyze slow transition dynamics
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
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
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
Theory of WED Self-consistent equation: the (short) initial segments of trajectories replace the initial configurations
Theory of WED A symmetric linear homogeneous equation:
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
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
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
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
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
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
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
Transition time vs projection single-time transition trajectories Transition state ensemble Without requiring knowledge of states and transitions
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
Transition network in 2D multi-well potential Topology of transition network is kept
Mexico-hat: entropy effects Topology of transition network is kept Eigenvalues gradually increase from 0 to 1
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)
450K and 300K Started from Ceq7 Started from Alpha7 Started from Cax7 Initial psi
Modified potential at 300K projection Free energy reconstruction
Occupation fraction vs projection (300K mod) Single-time transition trajectories
Transition time vs projection (300K mod) Single-time transition trajectories Real transition time by checking along the single-time transition trajectories
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
150K More dispersive projection Count of traj. 300K modified
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
Solvent effects Include solvent-related functions in expansion
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
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
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
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)
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
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
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
Hierarchic kinetic network Split trajectory into short segments: detect kinetic network in shorter time scales Hierarchic meta-stable state structure
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
Example 74 atoms, charged terminals; Implicit solvent simulation: Generalized Born; 1000 trajectories; 20 ns and 40000 conformationsper trajectory.
Application 12-alanine peptide 1000 20-ns trajectories 172 basis functions from torsion angles
Dimension of manifold Principle Component Analysis Total dimension: 172
clustering minimal spanning tree clustering algorithm 2
2D might be insufficient G1, 2D free energy profile