Computational Experiments in Coarse-Graining Atomistic Simulations

1. IMA, Minneapolis, January 2008 Computational Experiments in coarse graining atomistic simulations (Equation-free-& Variable Free- Framework For Complex/Multiscale Systems) I.G. Kevrekidis -C. W. Gear, G. Hummer, R. Coifman and several other good people- Department of Chemical Engineering, PACM & Mathematics Princeton University, Princeton, NJ 08544

2. SIAM– July, 2004 Clustering and stirring in a plankton model Young, Roberts and Stuhne, Nature 2001

3. Dynamics of System with convection

4. Simulation Method • Random (equal) birth and death, probability: l = m. • Brownian motion. • Advective stirring. (j, qare random phases) • IC: 20000 particles randomly placed in 1*1 box • Analytical Equation for G(r):

5. Stirring by a random field (color = y)

6. Dynamics of System with convection

8. Projective Integration: From t=2,3,4,5 to 10

9. RESTRICTION - a many-one mapping from a high-dimensional description (such as a collection of particles in Monte Carlo simulations) to a low-dimensional description - such as a finite element approximation to a distribution of the particles. LIFTING - a one-many mapping from low- to high-dimensional descriptions. We do the step-by-step simulation in the high-dimensional description. We do the macroscopic tasks in the low-dimensional description.

10. So, the main points: • You have a “microscopic code” • Somebody tells you what are good coarse variable(s) • Somebody tells you what KIND of equation this variable satisfies (deterministic, stochastic…) but NOT what the equation looks like. • Then you can use the IDEA that such an equation exists and closes to accelerate the simulation/ extraction of information. • KEY POINT – make micro states consistent with macro states

11. Application to Micelle Formation • Hydrophobic tail (T) • Hydrophilic head (H)

12. Lattice Model (Larson et al., 1985) • Surfactant = chain of H and T beads (H4T4) • No explicit solvent • Hydrophobic-hydrophilic interactions • modeled as direct interaction between • H and T beads • Pairwise interactions with nearest sites:

13. Snapshot of Micellar System T = 7.0, µ = - 47.40

14. Kinetic Approach to Rare Events(Hummer and Kevrekidis, JCP 118, 10762 (2003)) • Evolution of coarse variables is slow • micelle birth/death are rare events • Reconstruction of free energy surface: • long equilibrium simulation • series of short nonequilibrium simulations

15. Reconstruction of Free Energy Surface Obtain G() and () from short-scale nonequilibrium simulations

16. Reconstructed free energy curve

17. Results: Micelle Destruction Rate Kramers’ formula • Result of nonequilibrium simulation: k = 5.58 x 10-9 • Equilibrium result: k = 7.70 x 10-9 • CPU time required: less than 7% of equilibrium simulation • Extension to multi-dimensional systems (Hummer and Kevrekidis, 2003) • Chapman-Kolmogorov equation

18. Reverse Projective Integration – a sequence of outer integration steps backward; based on forward steps + estimation 1 2 3 Reverse Integration: a little forward, and then a lot backward ! We are studying the accuracy and stability of these methods

19. Reverse coarse integration from both sides

20. Reconstructed free energy curve

21. Details of Multidimensional Dynamics Small clusters: 2d dynamics Larger clusters: 1d dynamics

22. Multidimensional Dynamics (2nd variable = E)

23. w/ Gerhard Hummer, NIDDK / J.Chem.Phys. 03 Alanine Dipeptide In 700 tip3p waters The waters The dipeptide and the Ramachandran plot

24. Start with constrained MD • Let 50 configurations free • Estimate d/dt of average • Perform projective step

25. Alanine Dipeptide Energy Landscape  G.Hummer, I.G.Kevrekidis. Coarse molecular dynamics of a peptide fragment: Free energy, kinetics, and long-time dynamics computations J.Chem. Phys. 118 (2003).

27. Protocols for Coarse MD (CMD) using Reverse Ring Integration MD Simulations run FORWARD in Time Initialize ring nodes Step Ring BACKWARD in Energy

28. Fokker-Planck Equation (2D)for distribution P(x1,x2) Probability Currents 2D FPE: Diffusion Coefficients Drift Coefficient

29. Drift and Diffusion Coefficients   MD Simulations run FORWARD in Time Ring node ICs multiple replicas per node compute drift (v) and diffusion (D) coefficients P Dihedral angles use (v, D) estimates to check for existence of potential () and compute potential values at each node

30. Ring Stepping in Generalized Potential Ring at + Ring at  Goal: Compute potential  associated with ring nodes (ring “height” above x1-x2 plane)

31. Right-handed -helical minimum

32. ( D = scalar) Case I: Diffusion matrix is proportional to unit matrix: Probability Current: Probability Current vanishes Potential Conditions I Drift Coefficient Potential Conditions Conditions for existence of generalized potential : Path Path

33. Short bursts of MD simulation initialized at each node in the ring Data analysis of forward in time MD provides gradient information small step FORWARD in time large step BACKWARD in energy CMD Exploration of Alanine Dipeptideusing drift coefficients Ring nodes  Reverse ring integration stagnates at saddle points on coarse free energy landscape ONLY drifts estimated here ring step size is scaled by unknown diffusivity D  nodal “heights” unknown stepsize ~

34. UNFOLDED TRANSITION STATE FOLDED Fraction non-native contacts formed Fraction native contacts formed Src homology 3 (SH3) domain Protein modeled using off-lattice C representation associated with simplified minimally frustrated Hamiltonian Small Modular Domain 55-75 amino acids long blue (N-terminus) to red (C-terminus) n-Src Diverging Turn C-terminus Distal Loop Protein folding intrinsically low-dimensional “Collective” (coarse, slow) coordinates fully describe long time system dynamics N-terminus Characteristic fold consisting of five or six β-strands arranged as two tightly packed anti-parallel β sheets P. Das, S. Matysiak, and C. Clementi, Proc. Natl. Acad. Sci. USA 102, 10141 (2005).

35. Reverse Ring Integration and MDCoarse MD (CMD) Transition state identification using CMD

36. So, again, the main points • Somebody needs to tell you what the coarse variables are • And what TYPE of equation they satisfy • And then you can use this information to bias the simulations “intelligently” accelerating the extraction of information (also need hypothesis testing).

37. and now for something completely different: Little stars ! (well…. think fishes)

39. Fish Schooling Models Initial State Position, Direction, Speed INFORMED UNINFORMED Compute Desired Direction Zone of Attraction Rij< Zone of Deflection Rij< Normalize  Update Direction for Informed Individuals ONLY Couzin, Krause, Franks & Levin (2005) Nature (433) 513 Update Positions

40. INFORMED DIRN STICK STATES STUCK ~ typically around rxn coordinate value of about 0.5 INFORMED individual close to front of group away from centroid

41. INFORMED DIRN SLIP STATES SLIP ~ wider range of rxn coordinate values for slip 00.35 INFORMED individual close to group centroid

42. Effective Fokker-Planck Equation FPE: Diffusion Coefficient Drift Coefficient Potential:

43. Coarse Free Energy Calculation Estimate Drift and Diffusion coefficients numerically from simulation “bursts” X0 X0 X0 Improved estimates using Maximum Likelihood Estimation (MLE) Y. Aït-Sahalia. Maximum Likelihood Estimation of Discretely Sampled Diffusions: A Close-Form Approximation Approach. Econometrica 70 (2002). Kopelevich, Panagiotopoulos & Kevrekidis J Chem Phys 122 (2005) Hummer & Kevrekidis J Chem Phys 118 (2003)

44. Energy Landscape – Fish Swarming Problem CENTROID Informed DIRN SLIP SLIP STUCK

45. Using the computer to select variable Rationale: Straight Line Distance between locations NOT representative of actual transition difficulty\distance Lake Carnegie, Princeton, NJ

46. Lake Carnegie, Princeton, NJ Using the computer to select variable Rationale: Straight Line Distance Curved Transition Distance Actual transition difficulty represented by curved path

47. Using the computer to select good variables Rationale: Straight Line DistanceIS representative of actual transition difficulty\distance in small LOCAL patches Patch size related to problem “geography” Lake Carnegie, Princeton, NJ

48. Z X Y Using the computer to select variable Rationale: 3D Dataset with 2D manifold Selected Datapoint Euclidean Distance Euclidean distance in input space may be weak indicator of INTRINSIC similarity of datapoints Geodesic distance is good for this dataset

49. Dataset as Weighted Graph edges 3 W23 weights 2 W12 1 vertices parameter t

50. Unequal separation (Euclidean distance) between IC ( ) and limits of random walk ( , ) Multiple random walks through simulation data initialized at