Alchemical free energy calculations with gromacs

55 Pages·2008·13.93 MB·English

Alchemical free energy calculations with gromacs
John D. Chodera1, David L. Mobley2 , Michael R. Shirts3
1 Department of Chemistry, Stanford University
2 Department of Pharmaceutical Chemistry, University of California, San Francisco
3 Department of Chemistry, Columbia University
Building upon the work of many, many others...
gromacs workshop - Stanford, CA - Apr 7-8, 2008
1
Many questions in biology and pharmaceutical chemistry
involve free energy differences
How does compound X partition between different environments?
e.g. how well does lipitor partition between octanol and water?
How tightly does compound X bind protein Y? Protein Z?
e.g. how well does clomifene bind/discriminate ERα/ERβ
How do I affect the binding affinity if I modify compound X?
Or if I mutate or modify the protein?
e.g. what is the difference in binding between
pseudoephedrine and phenylephrine to
adrenergic receptor?
How do one or more mutations at a protein-protein
interface affect binding?
e.g. which residues at the interface contribute most to
binding, and how can mutations cause disease?
2
Free energy differences are often easier to compute
through alchemical routes
Relative free energies of binding Absolute free energies of binding
∆G
1 ∆G
P+L PL bind
1 1 P + L PL
∆G ∆G
a b
∆G
2
P+L PL
2 2 P + ø Pø
∆∆G = ∆G − ∆G
a b
3
βH(x)
Z = e dx (1)
−
!
1
∆F = β ln Z (2)
−
−
Alchemical intermediates can facilitate convergence
e βH(x)
−
β(F H(x))
P(x) = = e (3)
−
Z
P (x) exp β(F H (x))
A A A β(∆F ∆H(x))
= − = e (4)
−
Error increasesP ex(pxo)nenetxipalβly(F withH dim(xi)n)ishing phase space overlap
B B B
−
P (x)
A
= exp (β(∆F ∆E(x)) (5)
e.g. two harmonic oscillators separated by various distances
P (x) −
B
P (∆E
A
= exp (β(∆F ∆E)) (6)
) P ( ∆E) −
T
B
k −
(

s
a
i
b
β∆F β∆E(x)
P (x) = P (x)e (7)
A B −
number of samples
β(∆F ∆E(x))
1 = (e (8)
−
B
Instead, introduce inte∆rmFed=ia"teβ st1alntese toβ∆ eE#n(xs)ure a contiguous chain of go(9o)d overlap
− −
− B
" #
N 1
Z Z Z Z −
∆F = β 1 ln N = β 1 ln 2 3 N = ∆F (1Z0) = dxe−βU(x)
1 N − − n n+1 n
− Z − Z · Z · · · Z !
→ →
1 1 2 N 1
n=1
−
$
D. Wu and D. A. Kofke. "Phase-space overlap measures. I. Fail-safe bias detection...", J. Chem. Phys. 123: 054103 (2005).
4
1
Multiple ways to extract free energy differences from
simulations of intermediate states
TI (thermodynamic integration)
λ2 ∂H ∆λ ∂H ∂H
∆F = dλ! ≈ +
∂λ 2 ∂λ ∂λ
λ λ! $ λ λ %
! 1 " # " # 1 " # 2
quadrature error hard to quantify
EXP (exponential reweighting)
0.30
0.25
0.20
∆F = −β−1 ln e−β(U2−U1) = +β−1 ln e−β(U1−U2) 0.15
! " ! "
λ1 λ2
0.10
0.05
suffers from large bias and variance
0.00-5 0 5 10 15 20
BAR (Bennett acceptance ratio)
0.30
0.25
"f(U − U )#
−1 2 1 λ 0.20
∆F = −β ln 1
"f(U − U ) exp[−β(U − U )]# 0.15
1 2 2 1 λ
2
0.10
0.05
0.00 -10 -5 0 5 10 15 20
5
Multistate Bennett acceptance ratio (MBAR) estimator
extracts optimal estimate from multiple states
Statistically optimal analysis of samples from multiple equilibrium states
Michael R. Shirts (Department of Chemistry, Columbia University), John D. Chodera (Department of Chemistry, Stanford University)
http://arxiv.org/abs/0801.1426
0.025
0.020
0.015
0.010
0.005
0.000
-150 -100 -50 0 50 100 150
- optimal estimator for free energy differences from equilibrium data
- robust estimates of uncertainties
- combine data from multiple temperatures, pressures, bias potentials
- freely-available Python implementation
- provides tools to subsample correlated data to extract independent data
https://simtk.org/home/pymbar
6
Examples:
Free energies of solvation
Method Value
EXP (fwd) 3.80 ± 0.10
EXP (rev) 2.94 ± 0.14
EXP (avg) 3.37 ± 0.17
•
Example: Solvation of
3-methyl indole BAR 3.29 ± 0.06
•
5 ns explicit water
MBAR 3.23 ± 0.07
MD at 9 states
7
Optimal alchemical intermediates depend on problem
Goals in choice of intermediates:
few intermediates
good overlap between intermediates
minimize correlation times
minimize variance and size of space to be sampled
For example:
128 intermediates for extracting toluene from water
probably overkill
T. Steinbrecher, D. L. Mobley, and D. A. Case. "Non-linear scaling schemes for Lennard-Jones interactions in free energy calculations", J. Chem. Phys. 127: 214108 (2007). (DOI)
8
Hydration free energies
H
N
N S S
gas gas, neutral Ø
gas
ΔG
S S
aq aq, neutral Ø
water
M. R. Shirts, J. W. Pitera, W. C. Swope, and V. S. Pande. J. Chem. Phys. 119:5740 (2003)
D. L. Mobley, E. Dumont, J. D. Chodera and K. A. Dill. J. Phys. Chem. B. 111:2242-2254 (2007) (DOI)
9
Choice of intermediates for hydration free energy
annihilation decoupling
(no interactions within (interactions within
annihilated molecule retained) annihilated molecule retained)
potentially large ∆G smaller ∆G
electrostatics short correlation times longer correlation times due to naked charges
requires vacuum recharging calculation no vacuum recharging calculation required
3 PME evaluations per timestep, tricky to use
Lennard-Jones requires LJ annihilation vacuum simulation no vacuum simulation required
eliminates potentially unphysical conformations
Also possible to turn off both electrostatics and LJ simultaneously (e.g. pulling in 4th dimension), but much more tricky.
10