Exploring landscapes . . .

1. 700 K replica energy 200 K replica “important coordinates” Exploring landscapes . . .

2. Exploring landscapes for protein folding and binding using replica exchange simulations Effective potential Important coordinates • The AGBNP all atom effective solvation potential & REMD • Peptide free energy surfaces & folding pathways from all atom simulations and network models • Temp. dependence of folding: physical kinetics and replica exchange kinetics using network models • Replica exchange on a 2-d continuous potential with an entropic barrier to folding

3. AGBNP effective solvation potential(Analytical Generalized Born + Non Polar) • OPLS-AA AGBNP effective potential, an all atom model • Novel pairwise descreening Generalized Born model. • Separate terms for cavity free energy and solute-solvent van der Waals interaction energy. • Fully analytical. • Applicable to small molecules and macromolecules. Generalized Born Surface area model Born radius-based estimator E. Gallicchio, and R.M. Levy, JCC, 25, 479 (2004)

4. AGBNP: Pairwise Descreening Scheme Born radii: rescaled pairwise descreening approximation: j i Rescale according to self-volume of j: Self-volume of j (Poincarè formula, ca. 1880): Hawkins, Cramer, and Truhlar, JPC 1996 Schaefer and Karplus, JPC 1996 Qiu, Shenkin, Hollinger, and Still, JPC 1997 E. Gallicchio, R. Levy, J. Comp. Chem. (2004)

5. Non-Polar Hydration Free Energy Non-polar hydration free energy estimator: : Surface area of atom i : Estimator based on Born radius : Surface tension and van der Waals adjustable parameters R.M. Levy, L. Y. Zhang, E. Gallicchio, and A.K. Felts, JACS, 125, 9523 (2003) (proteins in water) E. Gallicchio, M. Kubo, and R.M. Levy, JPC, 104, 6271 (2000) (hydrocarbons in water)

6. Enthalpy-Entropy and Cavity Decomposition of Alkane Hydration Free Energies: Numerical Results and Implications for Theories of Hydrophobic Solvation E. Gallicchio, M. Kubo, R. M. Levy, J. Phys. Chem., 104, 6271 (2000)

7. The replica exchange method for structural biology problems • Has been successfully applied to protein and peptide folding, ligand binding, and NMR structure determination • Questions have been raised about the efficiency of the algorithm relative to MD e.g. Nymeyer, Gnanakaran & García (2004) Meth. Enz. 383: 119 Ravindranathan, Levy, et al. (2006) JACS 128: 5786 Chen, Brooks, et al. (2005) J. Biomol. NMR 31: 59 Beck, White & Daggett (2007) J. Struct. Biol. 157: 514 Zuckerman & Lyman (2006) JCTC 2: 1200 (with erratum)

8. 200 K Replica exchange molecular dynamics rough energy landscapes and distributed computing MD MD MD MD MD 700 K 450 K 320 K energy Y. Sugita, Y. Okamoto Chem. Phys. Let., 314, 261 (1999) “important coordinates”

9. 700 K replica walker 1 walker 2 energy walker 3 walker 4 200 K replica “important coordinates” Replica exchange molecular dynamics rough energy landscapes and distributed computing MD MD MD MD MD MD MD 700 K 450 K replica 320 K 200 K Y. Sugita, Y. Okamoto (1999) Chem. Phys. Let., 314:261

10. F2 U2 F1 U1 Protein folding: REM and kinetic network models • free energy surfaces of the GB1 peptide from • REM and comparison with experiment • kinetic network model of folding pathways • for GB1 Andrec M, Felts AK, Gallicchio E, Levy RM.. PNAS (2005) 102:6801. • kinetic network model of REMD • (simulations of simulations) Zheng W, Andrec M, Gallicchio E, Levy RM. PNAS (2007) 104:15340.

11. The -Hairpin of B1 Domain of Protein G Folding nucleus of the B1 domain Blanco, Serrano. Eur. J. Biochem. 1995, 230, 634. Kobayashi, Honda, Yoshii, Munekata. Biochemistry 2000, 39, 6564. Features of a small protein: stabilized by 1) formation of secondary structure 2) association of hydrophobic residues Munoz, Thompson, Hofrichter, Eaton. Nature 1997, 390, 196. Computational studies using Explicit and Implicit solvent models Pande, PNAS 1999 Dinner,Lazaridis,Karplus,PNAS,1999 Ma & Nussinov, JMB, 2000 Pande, et al., JMB, 2001 Garcia & Sanbonmatsu, Proteins, 2001 Zhou & Berne, PNAS, 2002

12. The b-Hairpin of B1 Domain of Protein G Simple (surf area) nonpolar model OPLS/AGBNP -hairpin > 90% -helix < 10% G ~ 2 kcal/mol The potential of mean force of the capped peptide. A Felts, Y. Harano, E. Gallicchio, and R. Levy, Proteins, 56, 310 (2004)

13. Kinetic network models for folding Andrec, Felts, Gallicchio & Levy (2005) PNAS, 102, 6801 Network nodes are snapshots from multiple temperatures of a replica exchange simulation. Tcold Thot 800,000 nodes 7.4 billion edges Dynamical/kinetic considerations: Transition rates (edges) are motivated by Kramers theory: transitions are allowed if there is sufficient structural similarity, and forbidden otherwise. Simulations are performed using the Gillespie algorithm for simulating Markov processes on discrete states: • Waiting time in a state is an exponential random variable with mean = 1/(Sj kij) • Next state is chosen with probability proportional to kij Equilibrium considerations: Sufficiently long trajectories must reproduce WHAM results.

14. Connection between kinetic model and equilibrium populations Equilibrium populations for temperature T0 are preserved if for each pair of nodes (i, j) the ratio of transition rates follows WHAM weighting: node j from temperature TB having energy Ej node i from temperature TAhaving energy Ei where fA(0) and fB(0) are free energy weights for the TA and TB simulations at reference temperature T0 These weights are order-parameter independent and will give correct PMFs for any projection. T-WHAM PMF at low temperature contains information from high temperature simulations

15. The majority of beta-hairpin folding trajectories pass through alpha helical intermediate states b a t = 9 units ≈ 180 ns t = 2500 units ≈ 50 µs b a Fraction of hairpin conformation averaged over 4000 stochastic trajectories run at 300 K and begun from an initial state ensemble equilibrated at 700 K. 91% of 4000 temperature-quenched stochastic trajectories begun from high-energy coil states pass through states with a-helical content Andrec, Felts, Gallicchio & Levy (2005) PNAS, 102, 6801

16. Evidence for a-helical intermediates in b-sheet folding and misfolding • Non-native helices have been observed in b-lactoglobulin folding • Rapid formation of a structure • Can exist as a stable thermodynamic species and as intermediates • May be important in protecting exposed ends of b-sheet from intermolecular interactions Forge, Hoshino, Kuwata, Arai, Kuwajima, Batt & Goto (2000) JMB 296:1039 Kuwata, Shastry, Cheng, Hoshino, Batt, Goto & Roder (2001) Nat. Struct. Biol. 8:151 • Amyloid b-sheets can form from a-helical precursors • Myoglobin and coiled-coil proteins can form amyloid fibrils Fändrich, Forge, Buder, Kittler, Dobson & Diekmann (2003) PNAS 100:15463 Kammerer, Dobson, Steinmetz et al. (2004) PNAS 101: 4435 • Fibril formation in amyloid b-protein may occur via a helical intermediate Kirkitadze, Condron & Teplow (2001) JMB 312:1103 Fezoui & Teplow (2002) JBC 277: 36948 • Computational and theoretical evidence • Helical structures have been observed in G-peptide simulations García & Sanbonmatsu (2001) Proteins 42:345 Zagrovic, Sorin & Pande (2001) JMB 313:151 Wei, Mousseau & Derreumaux (2004) Proteins 56:464 • Entropy-stabilized helical intermediates may be generic in b-sheet protein folding landscapes Chikenji & Kikuchi (2000) PNAS 97:14273

17. Exploring landscapes for protein folding and binding using replica exchange simulations Effective potential Important coordinates • The AGBNP all atom effective solvation potential & REMD • Peptide free energy surfaces & folding pathways from all atom simulations and network models • Temp. dependence of folding: physical kinetics and replica exchange kinetics with a network model • Replica exchange on a 2-d continuous potential with an entropic barrier to folding

18. F2 U2 F1 U1 F2U1 U2U1 U2F1 F2F1 F1U2 U1U2 F1F2 U1F2 One walker Two walkers N walkers Network models of Replica Exchange ku2 ku kuN F U kf2 FN UN kRE kRE kf kfN ku1 kf1 ku and kf: physical kinetics kRE: replica exchange “kinetics” ku2 F2 U2 kf2 kRE ku1 F1 U1 kf1 5 walkers: 3840 states N walkers: 2N N! states 2 walkers: 8 states Convergence at low temperature depends on the number of F1 to U1 to F1 “transition events” Gillespie “simulation of protein folding simulations”

19. † † Speed limit for replica exchange efficiency The number of transition events at low temperature is approximately equal to the average of the harmonic means of the rate constants at all temperatures: Results for 2 walkers: Non-Arrhenius case (∆Cp† < 0)

20. Replica exchange convergence is dependent on the physical kinetics of the system Zheng W, Andrec M, Gallicchio E, Levy RM. PNAS (2007) 104:15340. • The number of transition events depends on the average of the harmonic mean rates, and sets a “speed limit” for efficiency • Maximizing the rate of temperature diffusion is appropriate if the underlying kinetics is Arrhenius • For non-Arrhenius kinetics, an optimal temperature exists which maximizes the number of transition events and convergence • “Training” simulations (like those used for the multicanonical method) may be useful to locate optimal maximal temperatures

21. Replica exchange on a 2-d continuous potential with an entropic barrier to folding 2-d continuous potential Potential energy along x F U Simple Continuous and Discrete Models for Simulating Replica Exchange Simulations of Protein Folding W. Zheng, M. Andrec, E. Gallicchio, R. M. Levy,J. Phys. Chem., in press s

22. Rate constants extracted from (Uncoupled) simulations on the 2-d continuous potential ku kf

23. Uncoupled Reverse-engineering rates fex= 10-4 fex= 10-2 kf(T1) 6.1 6.4 6.3 ku(T1) 0.0036 0.0038 0.0037 kf(T2) 0.30 0.29 0.31 ku(T2) 0.42 0.42 0.43 Reverse-Engineering rates from the trajectory on the continuous potential using lifetime & branching ratios T1=296K T2=474K

24. RE on the continuous potential vs RE on the kinetic network The faster the replica exchange rate, the bigger the discrepancy. Kinetic network Continuous potl Total # of transitions fex = 5·10-2 fex = 5·10-3 fex = 1·10-3 Infinitely fast exchange limit* *Calculated using harmonic mean of rate constants

25. Non-Markovian effects -- History dependence Probability Calculated from the continuous traj. at different exchange rates fex =5·10-3 fex =10-4 P(U2F1U1F2) 0.849 0.103 P(F2F1U2F1U1F2) 0.521 0.094 P(U2F1F2F1) 0.150 0.886 P(F2F1U2F1F2F1) 0.477 0.895

26. Summary • Non-Markovian effects are observed in Replica Exchange simulations on the continuous potential • When the frequency of replica exchange exceeds the time scale for relaxation in the F and U macrostates, the convergence rate slows • The efficiency of RE in more complex systems is fundamentally limited by the time scale of conformational diffusion within the free energy basins.

27. Exploring landscapes for protein folding and binding using replica exchange simulations Effective potential Important coordinates • The AGBNP all atom effective solvation potential & REMD • Emilio Gallicchio • Peptide free energy surfaces & folding pathways from all atom simulations and network models • Tony Felts, Zenmei Ohkubo, and Michael Andrec • Temp. dependence of folding: physical kinetics and replica exchange kinetics • Weihua Zheng, Michael Andrec, Emilio Gallicchio • Replica exchange on a 2-d continuous potential with an entropic barrier to folding Weihua Zheng, Michael Andrec, Emilio Gallicchio

28. Protein Folding with All Atom Potentials Insights using Replica Exchange and Network Models Effective potential Important coordinates • The AGBNP all atom effective solvation potential • Emilio Gallicchio, Tony Felts • Peptide free energy surfaces and folding pathways • Tony Felts, Zenmei Ohkubo, and Michael Andrec • Network models and kinetics in the replica exchange ensemble • Michael Andrec, Emilio Gallicchio

29. Potential of mean force (PMF) along x F U

30. Replica exchange convergence is dependent on the physical kinetics of the system • The number of transition events depends on the average of the harmonic mean rates, and sets a “speed limit” for efficiency • Maximizing the rate of temperature diffusion is appropriate if the underlying kinetics is Arrhenius • For non-Arrhenius kinetics, an optimal temperature exists which maximizes the number of transition events and convergence • “Training” simulations (like those used for the multicanonical method) may be useful to locate optimal maximal temperatures

31. Replica Exchange and Ligand Binding Binding free energy landscape contains multiple minima Effect of binding & temperature is to shift distribution of conformations Replica Exchange addresses the sampling problem while providing estimates of populations Reduced Coordinate Folding Landscape Binding Landscape

32. The P450 puzzle • Cytochrome P450s metabolize many aliphatic molecules and 90% of pharmaceutical ligands NPG P450BM-3/NPG Phe87 Heme • Several X-ray crystal structures of P450s show substrate distant from active site • Hypothesized a conformational equilibrium between productive and unproductive conformational states

33. Experimental and Modeling Clues • UV-VIS and SSNMR experiments indicate temperature-dependent equilibrium between Fe-bound and un-bound species. Induced Fit docking finds a Fe-bound conformation of higher energy than the X-ray conformation ω1-Fe X-ray (Distal) Induced Fit Model* (Proximal) • Questions: • Do the Xray and Induced Fit structures correspond to the low and high temperature conformations? • Are there other states? • What’s the mechanism of interconversion between states? *Jovanovic, T.; Farid, R.; Friesner, R. A.; McDermott, A. E. J. Am. Chem. Soc. 2005, 127, 13548.

34. REMD of P450 NPG Complex • Provides populations of conformational states (canonical sampling) as a function of temperature • Can be used to construct free energy landscape • Model ligand and 120 active site residues • 24 replicas between 260 and 463 K • 72 ns total aggregate simulation time Ravindranathan, K.P., E. Gallicchio, R.A. Friesner, A.E. McDermott, and R.M. Levy. J. Am. Chem. Soc.,128, 5786-5791 (2006).

35. Population of proximal conformations Free Energy Landscape Phe87 1 Proximal ligand-free • New proximal ligand-free state, most populated at physiological temperature.Entropically stabilized • Conversion from distal state goes through proximal ligand-locked conformation • Barrier from proximal to distal is about 4 Kcal/mol. T-WHAM used to resolve barrier region Phe87 2 Proximal ligand-locked 1-Fe Distance [Å] Distal

36. Conclusions • REMD shows the conformational transition and supports thermal activation hypothesis • Proximal state stabilized by conformational entropy • Conformational states exist at all temperatures: relative populations change with temperature