Molecular Motion Pathways: Computation of Ensemble Properties with Probabilistic Roadmaps

1. Molecular Motion Pathways: Computation of Ensemble Properties with Probabilistic Roadmaps • A.P. Singh, J.C. Latombe, and D.L. Brutlag. A Motion Planning Approach to Flexible Ligand Binding. Proc. 7th Int. Conf. on Intelligent Systems for Molecular Biology (ISMB), AAAI Press, Menlo Park, CA, pp. 252-261, 1999. • N.M. Amato, K.A. Dill, and G. Song. Using Motion Planning to Map Protein Folding Landscapes and Analyze Folding Kinetics of Known Native Structures. J. Comp. Biology, 10(2):239-255, 2003. • M.S. Apaydin, D.L. Brutlag, C. Guestrin, D. Hsu, J.C. Latombe, and C. Varma. Stochastic Roadmap Simulation: An Efficient Representation and Algorithm for Analyzing Molecular Motion. J. Comp. Biology, 10(3-4):257-281, 2003. • N. Singhal, C.D. Snow, and V.S. Pande. Using Path Sampling to Build Better Markovian State Models: Predicting the Folding Rate and Mechanism of a Tryptophan Zipper Beta Hairpin, J. Chemical Physics, 121(1):415-425, 2004. • J. Cortés, T. Siméon, M. Renaud-Siméon, and V. Tran. Geometric Algorithms for the Conformational Analysis of Long Protein Loops. J. Comp. Chemistry, 25:956-967, 2004.

2. Molecular motion is an essential process of life Mad cow disease is caused by misfolding Drug molecules act bybinding to proteins

3. Stanford BioX cluster NMR spectrometer So, studying molecular motion is of critical importance in molecular biology However, few tools are available • Computer simulation: • Monte Carlo simulation • Molecular Dynamics

4. Two Major Drawbacks of MD and MC Simulation • Each simulation run yields a single pathway, while molecules tend to move along many different pathways  Interest in ensemble properties

5. 1- pfold pfold Example of Ensemble Property: Probability of Folding pfold Measure kinetic distance to folded state Du, Pande, Grosberg, Tanaka, and Shakhnovich. On the Transition Coordinate for Protein Folding. Journal of Chemical Physics (1998). Unfolded state Folded state

6. Other Examples of Ensemble Properties • Folding: • Order of formation of SSE’s • Folding rate / Mean first passage time • Key intermediates • Binding: • Average time to escape from active site • Average energy barrier

7. Two Major Drawbacks ofMD and MC Simulation • Each simulation run yields a single pathway, while molecules tend to move along many different pathways • Each simulation run tends to waste much time in local minima

8.  Roadmap-Based Representation • Compact representation of many motion pathways • Coarse resolution relative to MC and MD simulation • Efficient algorithms for analyzing multiple pathways

9. free space [Kavraki, Svetska, Latombe,Overmars, 96] Roadmaps for Robot Motion Planning

10. Initial WorkA.P. Singh, J.C. Latombe, and D.L. Brutlag. A Motion Planning Approach to Flexible Ligand Binding. Proc. 7th ISMB, pp. 252-261, 1999 • Study of ligand-protein binding • The ligand is a small flexible molecule, but the protein is assumed rigid • A fixed coordinate system P is attached to the protein and a moving coordinate system L is defined using three bonded atoms in the ligand • A conformation of the ligand is defined by the position and orientation of L relative to P and the torsional angles of the ligand

11. Roadmap Construction (Node Generation) • The nodes of the roadmap are generated by sampling conformations of the ligand uniformly at random in the parameter space (around the protein) • The energy E at each sampled conformation is computed: E = Einteraction + Einternal Einteraction = electrostatic + van der Waals potential Einternal = Snon-bonded pairs of atoms electrostatic + van der Waals

12. P = if Emin E  Emax Roadmap Construction (Node Generation) • The nodes of the roadmap are generated by sampling conformations of the ligand uniformly at random in the parameter space (around the protein) • The energy E at each sampled conformation is computed: E = Einteraction + Einternal Einteraction = electrostatic + van der Waals potential Einternal = Snon-bonded pairs of atoms electrostatic + van der Waals • A sampled conformation is retained as a node of the roadmap with probability:0 if E > Emax Emax-E Emax-Emin 1 if E < Emin  Denser distribution of nodes in low-energy regions of conformational space

13. qi qi+1 q q’ Emax Roadmap Construction (Edge Generation) • Each node is connected to its closest neighbors by straight edges • Each edge is discretized so that between qi and qi+1 no atom moves by more than some ε (= 1Å) • If any E(qi) > Emax , then the edge is rejected E

14. qi qi+1 q q’ Heuristic measureof energetic difficultyor moving from q to q’ Roadmap Construction (Edge Generation) • Any two nodes closer apart than some threshold distance are connected by a straight edge • Each edge is discretized so that between qi and qi+1 no atom moves by more than some ε (= 1Å) • If all E(qi)  Emax , then the edge is retained and is assigned two weights w(qq’) and w(q’q) where: (probability that the ligand moves from qi to qi+1 when it is constrained to move along the edge)

15. Querying the Roadmap • For a given goal node qg (e.g., binding conformation), the Dijkstra’s single-source algorithm computes the lowest-weight paths from qg to each node (in either direction) in O(NlogN) time, where N = number of nodes • Various quantities can then be easily computed in O(N) time, e.g., average weights of all paths entering qg and of all paths leaving qg(~ binding and dissociation rates Kon and Koff) Protein: Lactate dehydrogenase Ligand: Oxamate (7 degrees of freedom)

16. Experiments on 3 Complexes • PDB ID: 1ldm Receptor: Lactate Dehydrogenase (2386 atoms, 309 residues) Ligand: Oxamate (6 atoms, 7 dofs) • PDB ID: 4ts1 Receptor: Mutant of tyrosyl-transfer-RNA synthetase (2423 atoms, 319 residues) Ligand: L- leucyl-hydroxylamine (13 atoms, 9 dofs) • PDB ID: 1stp Receptor: Streptavidin (901 atoms, 121 residues) Ligand: Biotin (16 atoms, 11 dofs)

17. active site Computation of Potential Binding Conformations • Sample many (several 1000’s) ligand’s conformations at random around protein • Repeat several times: • Select lowest-energy conformations that are close to protein surface • Resample around them • Retain k (~10) lowest-energy conformations whose centers of mass are at least 5Å apart lactate dehydrogenase

18. Results for 1ldm • Some potential binding sites have slightly lower energy than the active site  Energy is not a discriminating factor • Average path weights (energetic difficulty) to enter and leave binding site are significantly greater for the active site  Indicates that the active site is surrounded by an energy barrier that “traps” the ligand

19. Energy Conformation Potential bindingsite Potential bindingsite Active site

20. Application of Roadmaps to Protein FoldingN.M. Amato, K.A. Dill, and G. Song. Using Motion Planning to Map Protein Folding Landscapes and Analyze Folding Kinetics of Known Native Structures. J. Comp. Biology, 10(2):239-255, 2003 • Known native state • Degrees of freedom: φ-ψ angles • Energy: van der Waals, hydrogen bonds, hydrophobic effect • New idea: Sampling strategy • Application: Finding order of SSE formation

21. Sampling Strategy(Node Generation) • High dimensionality  non-uniform sampling • Conformations are sampled using Gaussian distribution around native state • Conformations are sorted into bins by number of native contacts (pairs of C atoms that are closeapart in native structure) • Sampling ends when all bins have minimum number of conformations  “good” coverage of conformational space

22. Application: Order of Formation of Secondary Structures • The lowest-weight path is extracted from each denatured conformation to the folded one • The order of formation of SSE’s is computed along each path • The formation order that appears the most often over all paths is considered the SSE formation order of the protein

23. Method • The contact matrix showing the time step when each native contact appears is built

24. Protein CI2 (1a + 4 b)

25. 60 5 The native contact between residues 5 and 60 appears at step 216 Protein CI2 (1a + 4 b)

26. Method • The contact matrix showing the time step when each native contact appears is built • The time step at which a structure appears is approximated as the average of the appearance time steps of its contacts

27. a forms at time step 122 (II) b3 and b4 come together at 187 (V) b2 and b3 come together at 210 (IV) b1 and b4 come together at 214 (I) a and b4 come together at 214 (III) Protein CI2 (1a + 4 b)

28. Method • The contact matrix showing the time step when each native contact appears is built • The time step at which a structure appears is approximated as the average of the appearance time steps of its contacts

29. SSE’s roadmap size 1a+4b 5126, 70k 5471, 104k 3a 7975, 104k 1a+4b 8357, 119k 1a+5b CI2 Comparison with Experimental Data

30. vi Pij vj Stochastic RoadmapsM.S. Apaydin, D.L. Brutlag, C. Guestrin, D. Hsu, J.C. Latombe and C. Varma. Stochastic Roadmap Simulation: An Efficient Representation and Algorithm for Analyzing Molecular Motion. J. Comp. Biol., 10(3-4):257-281, 2003 New Idea:Capture the stochastic nature of molecular motion by assigning probabilities to edges

31. Pii vi Pij Edge probabilities Follow Metropolis criteria: Self-transition probability: vj [Roadmap nodes are sampled uniformly at random and energy profilealong edges is not considered]

32. Stochastic roadmap simulation and Monte Carlo simulation converge to the Boltzmann distribution, i.e., the number of times SRS is at a node in V converges towardwhen the number of nodes grows (and they are uniformly distributed) Stochastic Roadmap Simulation V Pij

33. Roadmap as Markov Chain j Pij i • Transition probability Pij depends only on i and j

34. 1- pfold pfold Example #1: Probability of Folding pfold Unfolded state Folded state

35. U: Unfolded state F: Folded state =1 =1 First-Step Analysis • One linear equation per node • Solution gives pfold for all nodes • No explicit simulation run • All pathways are taken into account • Sparse linear system l k j Pik Pil Pij m Pim i Pii Let fi = pfold(i) After one step: fi = Pii fi + Pij fj + Pik fk + Pil fl + Pim fm

36. Number of Self-Avoiding Walks on a 2D Grid 1, 2, 12, 184, 8512, 1262816, 575780564, 789360053252, 3266598486981642, (10x10) 41044208702632496804, (11x11) 1568758030464750013214100, (12x12) 182413291514248049241470885236 > 1028 http://mathworld.wolfram.com/Self-AvoidingWalk.html

37. In contrast … Computing pfold with MC simulation requires: Foreveryconformation q of interest • Perform many MC simulation runs from q • Count number of times F is attained first

38. Computational Tests • 1HDD (Engrailed homeodomain) • 3 a helices • 12 DOF • 1ROP (repressor of primer) • 2 a helices • 6 DOF H-P energy model with steric clash exclusion [Sun et al., 95]

39. Correlation with MC Approach 1ROP

40. pfold for ß hairpin Immunoglobin binding protein (Protein G) Last 16 amino acids Cα based representation Go model energy function 42 DOFs [Zhou and Karplus, `99]

41. Computation Times (ß hairpin) Monte Carlo (30 simulations): Over 107energy computations ~10 hours of computer time 1 conformation Roadmap: ~50,000energy computations 23 seconds of computer time 2000 conformations ~6 orders of magnitude speedup!

42. Example #2: Ligand-Protein Interaction Computation of escape timefrom funnels of attraction around potential binding sites Funnel of attraction = ball of 10Å rmsd around bound state [Camacho and Vajda, 01]

43. Computation Through Simulation [Sept, Elcock and McCammon `99] 10K to 30K independent simulations

44. l k Pil Pik m Pij j Pim i Pii Funnel of Attraction Computing Escape Time with Roadmap ti = 1 + Piiti + Pijtj+ Piktk + Piltl + Pimtm (escape time is measured as number of stepsof stochastic simulation) = 0

45. Distinguishing Active Site Given several potential binding sites,which one is the active one? Energy: electrostatic + van der Waals + solvation free energy terms

46. Complexes Studied

47. Distinction Using Escape Time (# steps)

48. Using Path Sampling to Construct RoadmapsN. Singhal, C.D. Snow, and V.S. Pande. Using Path Sampling to Build Better Markovian State Models: Predicting the Folding Rate and Mechanism of a Tryptophan Zipper Beta Hairpin, J. Chemical Physics, 121(1):415-425, 2004 New idea: Paths computed with Molecular Dynamics simulation techniques are used to create the nodes of the roadmap More pertinent/better distributed nodes  Edges are labeled with the time needed to traverse them

49. ~dt Sampling Nodes from Computed Paths (Path Shooting) F U

50. j tij pij i Example: Langevin dynamics equation of motion is where R is a Gaussian random force Sampling Nodes from Computed Paths (Path Shooting) F U