Absolute binding free energies

1. ...of ligands from umbrella sampling and steered MD simulations; ...applications to ions, small molecules and toxin peptides. Absolute binding free energies

2. Contents of this presentation I Introduction • General principles of free energy calculations • Illustrations by simple examples • Toy models: ions in small water boxes • Theoretical Underpinnings • Ways to affirm and void the fidelity of simulation models to reality • Advanced examples • Ligands in gramicidin-A channel, K+ channels

3. Contents of this presentation II Introduction • Full analysis of an absolute free energy of binding calculation • Using toxin peptides (my work) • Summary and notes • Implications of various assumptions • Warning signs to watch out for ...let’s go!

4. Goal of FE calculations Introduction • Accuracy • Chief motivation of conducting calculations is to predict/replicate experimental values • We must not lose sight of this • Therefore, it is worthwhile to test and verify our foundations • Especially for calculations that last over a month

5. Mechanics of FE calculations Introduction • Umbrella sampling and steered MD are path-dependentcalculations e.g. the ligand must physically move out of the receptor • A well-constructed simulation can: • Take information from known reaction mechanisms • Provide information on how those mechanisms occure

6. Foundation of FE calculations Introduction F = ma • All simulations are models • Classical approximation of QM phenomena • System proceeds via forces imposed on atoms • Energy calculations occur via manipulation of the these forces • Small loss in accuracy must be granted ...motivated?

7. Steered molecular dynamics A simple manipulation F= ½ k ( r – r0 + tv ) • SMD is a “moving spring” • Force is imposed on an atom or set of atoms • Equilibrium position moves to create a path • Work done by this spring is used to find the local free energy surface.

8. A simple model Steered Molecular Dynamics • A small box, containing: • ~400 water atoms, • one neon atom. • Suppose that we create an artificial gaussian barrier (of ~5 kT) for the Ne atom. • How do I measure the free energy surface of that barrier, using the Ne as a probe?

9. Simulation design Steered Molecular Dynamics • Constrain position of Ne in (x,y,z) to a starting location • Define path along water box • Measure the total amount of work done W = ∫ F(z) dz • W → manipulations → free energy surface

10. Work done(?) Steered Molecular Dynamics • Resulting potential of mean force • Includes environmental influences • i.e. action of water around the atom • NB: notice that PMF does not return to zero readily

11. Averaging Steered Molecular Dynamics • A single trajectory does not “sample” the reaction sufficiently • i.e. A free energy surface requires comprehensive knowledge of the sub-states, ala possible trajectories • Trajectories contain influence from the SMD pulling itself • Some of the work is directed at displacing the solvent around Ne • This needs to be accounted for

12. Jarzynski’s Equality Steered Molecular Dynamics • .˙. Take multiple trajectories • The average of multiple trajectories should even out random influences on the PMF < e-W/kT > ≈ Σi=1ne-Wi / kT / n • .˙. Apply Jarzynski’s Equality (given certain assumptions) e-ΔG/kT = < e-W/kT >

13. Average work Steered Molecular Dynamics • The Boltzmann average of 10 calculations • v = 10 Å ns-1 • A very fast calculation, converges because Neon does not interact with environment ...next: theory

14. Criteria of SMD-based FE calculations SMD Assumptions • Reaction coordinate must be well chosen • This coordinate measures all contributions to the real ΔG, and only those contributions. • The simulation must be near-equilibrium • Jarzynski's Equality holds when no dissipative work is done by the pulling

15. Reaction coordinates SMD Assumptions In SMD, the dimension(s) controlled by the constraining potential is areaction coordinate • The reaction coordinate can be: • Distance • Center of mass (collective variable) • RMSD to target structure • The path traced through the SMD simulation is a reaction path

16. Theoretical interpretation SMD Assumptions • Our systems are canonical ensembles: Z(n,P,T) • We wish to measure particular sub-states of the system • e.g. ligand-bound and ligand-unbound. • Thus, steered molecular dynamics (SMD): • Constrains the system to certain sub-states according to some coordinate • Drives the system along this coordinate to new sub-states by application of forces.

17. Implications on path design SMD Assumptions • Both the reaction-path and the simulations must sufficiently sample the required sub-space • May not be complete, but must be representative • In a ligand-unbinding process • Translation must be primary • Rotation may be secondary Phase space Bulk Bound

18. Implications on path design SMD Assumptions • Sampling and stability of a system suffers where large barriers must be crossed • Therefore, choose path of least resistance Simply using the distance between the ligand and the whole receptor may not be a good choice…

19. Criteria of SMD-based FE calculations SMD Assumptions • Reaction coordinate must be well chosen • This coordinate measures all contributions to the real ΔG, and only those contributions. • The simulation must be conducted in a near-equilibrium state • Jarzynski's Equality holds when no dissipative work is done by the pulling

20. Theoretical implications SMD Assumptions SMD –(Jarzynski's Equality)→ Free energy • JE comes with certain qualifications: • Dissipative work can be done on the system • Pulling velocity i.e. System must equilibrate around SMD perturbation, else perturbation will also be measured. (We will show this later.) ...next: 

21. Umbrella Sampling Collecting local information F= ½ k ( r – ri), i along path • US is a static potential • Force is also imposed on an atom or set of atoms • Multiple overlapping states are constructed to cover reaction path. • Each state provides information about local surface • Link to derive complete surface ...

22. Second toy model Umbrella sampling • A box containing two ions • sodium and chloride • Solution known • FE surface related to radial distribution function • This was done ab-initio yesterday

23. Second toy model Umbrella sampling • Reaction coordinate • Na – Cl separation • Umbrella potential • 1 Å apart, 2.5-9.5 Å • k = 10 kcal mol-1 Å-2 • Derive original by WHAM analysis

24. PMF convergence Umbrella sampling Phase space Bulk Bound

25. PMF Umbrella sampling Phase space Bulk Bound ...pretty...

26. Criteria of US-based FE calculations US assumptions • Reaction coordinate must be well chosen • This coordinate measures all contributions to the real ΔG, and only those contributions. • Sufficient sampling over the entire path: • Convergence of PMF curve means that environmental variables are well sampled. • Overlap between adjacent windows.

27. Reaction coordinate US assumptions • Same arguments as for SMD (underlying physics identical) • Umbrella sampling is capable of treating two/three dimensions • As long as all dimensions are properly sampled

28. Criteria of US-based FE calculations US assumptions • Reaction coordinate must be well chosen • This coordinate measures all contributions to the real ΔG, and only those contributions. • Sufficient sampling over the entire path: • Convergence of PMF curve means that environmental variables are well sampled • Overlap between adjacent windows

29. Environmental variables Bulk US assumptions • All coordinates not included in your reaction coordinate(s) must be well sampled • This means that simulation has visited dimensions perpendicular to reaction path that may contribute to FEbind Bound

30. Window Overlap Bulk US assumptions • Require enough sampling to accurately interpolate between windows Bound

31. Window-overlap US assumptions • Define measure of overlap between two distributions • When underlying surface is flat, harmonic potential produces gaussian distributions • Theoretical overlap: • Ω = [ 1- erf(d/8σ)]

32. Ω = [ 1- erf(d/8σ)] US assumptions • In practice, minimum overlap is ~2% • Overlap should agree with theoretical value when in bulkk = 20 kcal/mol/Å-2d = 0.5 ÅΩ = 15%k = 40 kcal/mol/Å-2d = 0.5 ÅΩ = 4%

33. Summary Steered MD –vs– Umbrella Sampling • Constructing FE surfaces via SMD • Straightforward construction • Relies on JE conditions: Difficult to achieve in practice • Constructing FE surfaces via US • Additional checks and balances • Both dependent on sufficient sampling ...take a break?

34. JE/US comparisons J Chem. Phys. 128:155104 (2008) • Using several testcases of ions andmolecules inchannel systems • Ion transit through membrane • Ion binding to gramicidin exterior • Organic-cation binding to gramicidin-A

35. JE/US comparisons J Chem. Phys. 128:155104 (2008) • Comparing PMFs • Tests balanced byequalising the totalsimulation time SMD setup @ v=5 Å ns-1 ~ US setup • SMD: Also use different pulling velocities to test reversibility of JE

36. Ion transit (nanotube) J Chem. Phys. 128:155104(2008) • Results equivalentbetween two methods • Energy surface not equal at both openings • resulting from system setup • JE valid

37. Ion transit (gA) J Chem. Phys. 128:155104(2008) • Umbrella sampling, not SMD, gives symmetric surface • Pulling at different velocities do not seem to help • Can potentially usev < 1 Å ns-1 • But more time consuming than equivalent US setup

38. JE: Practical problems? J Chem. Phys. 128:155104(2008) • Equilibration time sharply increases for peptide environments • Nanotube highly ordered .˙. Fast dissipation • Reliance on sampling “negative work” trajectories • High v: low probability for the environment to push SMD particle

39. K+-binding to gA J Chem. Phys. 128:155104(2008) ( next test case: ) • gA has a weak ion binding site at the entrances • Smaller potentials • Perhaps using a smaller k will reduce the perturbations v =2.5 Å ns-1

40. K+-binding to gA J Chem. Phys. 128:155104(2008) • What about using different force constants? • No significant help in repairing JE assumptions • Using small k may reduce perturbations on system • However, binding site shape lost • k must be greater than binding well ‘potential’ v =2.5 Å ns-1

41. K+-binding to gA J Chem. Phys. 128:155104(2008) • There are hard limits to varying parameters k = 2 kcal mol-1Å-2 k = 20 kcal mol-1Å-2

42. EA and TEA binding J Chem. Phys. 128:155104(2008) • Ethylammonium (EA) and tetra-ethylammonium (TEA) bind weakly to gA • v = 2.5 A /ns • k = 20 kcal mol-1 Å-2 • Reducing barrier height produces no difference here

43. CnErg1 Toy comparison J Chem. Phys. 128:155104(2008) • If it doesn’t work for small cations, it won’t work for a peptide • Test for a purported binding of CnERG1 toxin to hERG channel since rate of dissipation to environment is less than rate of work done… • SMD-PMFs essentially measures work done to move solvent

44. Mechanisms J Chem. Phys. 128:155104 (2008) • Input: work done on system • Carried out by imposed SMD forces (irreversible) • Contribution from underlying FE surface (reversible) • Output: dissipation to heat-bath (NPT systems) • Equilibration occurs by two means: • Forces bleed out to atoms far from SMD location • Temperature coupling to thermostat • JE maintained in O < I conditions (only in nanotube)

45. Mechanisms J Chem. Phys. 128:155104 (2008) • The time for equilibration issuch that v << 2.5 Ang/ns is required for JE condition to hold • This velocity requirement become more stringent as: • ligand size increases • interactions increase • Not as efficient as umbrella sampling. ...paper finished.

46. Organic cations-gA J. Phys. Chem. B 111:11303 (2007) • Block of GA by various small molecular ligands • Energetics dependent on ligand size and partial charges • Comparison between ligands, and with extant experimental data

47. Organic cations J. Phys. Chem. B 111:11303 (2007) • Use Autodock3 to find potential binding sites of molecules • MD, umbrella sampling to find PMF and free energy of binding • COM-coordinates of ligands 10 kcal mol-1Å-2 0.5 Å ns-1

48. Molecule list J. Phys. Chem. B 111:11303 (2007) • Use of six different molecules • Varying sizes and polarity • Determines strength of binding, and whether molecules can permeate through gA

49. MA and EA J. Phys. Chem. B 111:11303 (2007)

50. FMI and GNI J. Phys. Chem. B 111:11303 (2007)