An Introduction to Free Energy Calculations

1
An Introduction to Free Energy Calculations
Alan E. Mark
School of Molecular and Microbial Sciences
(SMMS)Chemistry Building (68)University of
QueenslandBrisbane, QLD 4072,AustraliaEmail
a.e.mark_at_uq.edu.auPhone 61-7-3365 4180 FAX
61-7-3365 3872
2
Free Energy The Holy Grail of Computational
Chemistry
Almost all measurable properties depend (directly
or indirectly) on the free energy of the system
Hamiltonian (energy function)
Free Energy
Probability of finding the system in a given
state
How the system preferentially evolves in time.
- rates of reaction - transition
paths (non-equilibrium processes)
- conformational equilibria - association
constants (equilibrium properties)
3

MD Simulations
Table 1 Classification of Properties
Class of Properties Examples Experimental Quantity
Instantaneous properties Only current state -potential energy -kinetic energy -configuration ... (easy) -heat of vaporization -temperature -pressure -structure
Time dependent properties Current state previous states -position, f(t) -velocity, f(t) -orientation, f(t) -dipole moment, f(t) ... (depends on time scale) -diffusion -order parameters -dielectric properties
Global properties All possible states -free energy -entropy (hard, depends on the extent of the available phase space) -binding constants -relative stability -barrier heights -heat capacity -compressibility -most things experimentally measurable
4
Thermodynamic Relations
5
Classical Statistical Mechanics
For a system of N particles
masses Cartesian
coordinates Conjugate momenta
Interaction function
Hamiltonian
kinetic
potential
6
Basic statistical mechanics
For a system in equilibrium which can be in M
energy states, the probability Pn of the system
having an energy En is given by k
Boltzmann constant T temperature Z
canonical partition function
Z the canonical partition function The sum
or integral over all possible states of the
system (normalization).
Link to Thermodynamics Helmholtz free energy
7
(No Transcript)
8
(No Transcript)
9
ltAgt ensemble average of A
10
Classical Statistical Mechanics
h Plancks constant k Boltzmanns constant T
temperature
Partition Function
F kT lnZ F ideal gas Fconfigurational
Free energy
11
Methods to Compute Free Energies (cont.)
Direct Determination of the Absolute Free Energy

• Problems
• Never sample all configurational space.
• Each addition part of phase space sampled makes a
  negative (favorable) contribution to the free
  energy

Integrand always gt o
i) Absolute free energy ? only if complete
partition function can be enumerated. ii)
Simulation techniques ? absolute free energy
systematically over-estimated.
12
Methods to Compute Free Energies (cont.)
Absolute Free Energy as Ensemble Average (1)
Q. Can one estimate the absolute free energy as
an ensemble average?
excess free energy
13
Methods to Compute Free Energies (cont.)
Absolute Free Energy as Ensemble Average (2)
Absolute free energy as an ensemble average
NOTE Exponent in ve. Low energy
states dominate ensemble BUT high energy states
dominate the free energy (entropy is
important). Probability of sampling states
contributing little (or nothing) to F ? high
Probability of sampling states contributing a
lot to F ? low. Very poor
convergence!!
Q. Can one estimate the absolute free energy as
an ensemble average? NO
14
Methods to Compute Free Energies (cont.)
Absolute Free Energy Simulation ?
no Experiment ? no
Experimentally only measure Relative Free
Energies i.e. difference in free energy ?F
between two (or more) states of a system.
?FAB FB FA
15
Methods to Compute Free Energies (cont.)
Relative free energies Examples
?Fbinding Fbound Funbound ?Fstability
F folded Funfolded ?Fpartition F solvent
B Fsolvent A
16
Methods to Compute Free Energies ?G from
experiment

• Two Possibilities
• Relative probability of finding the system in one
  of two states.
• 2. Work required to go from an initial to a
  final state via a reversible path.

Example 1. Relative Probability Free energy of
binding, ?Gbinding
KA
ligand acceptor ? ligandacceptor
complex
Association constant, KA relative probability
of bound/free
17
Methods to Compute Free Energies ?G from
experiment

• Two Possibilities
• Relative probability of finding the system in one
  of two states.
• 2. Work required to go from an initial to a
  final state via a reversible path.

Example 2. Reversible work Change of state.
Using
V(x) potential
constant T
18
Methods to Compute Free Energies ?G from
simulation

• Two Possibilities
• Relative probability of finding the system in one
  of two states.
• 2. Work required to go from an initial to a final
  state via a reversible path.

Same as experiment
A. In terms of a probability function, P(r)
B. As reversible work along coordinate
work integral of force (derivative of the
potential w.r.t a given coordinate)
V(r) potential
Choice of coordinate system is arbitrary
19
Methods to Compute Free Energies (cont.)
Temperature Integration
Method Perform constant V,N simulations at
different temperatures to obtain average of
E. Integrate numerically
Reference State Can use low temperature harmonic
state as a reference but phase transitions are a
problem.
20
Methods to Compute Free Energies (cont.)
Pressure Integration
As for temperature integration.
Method Perform constant N,T simulations at
different volumes to compute pressure. Integrate
numerically
Reference State A state with a large volume
(treat as ideal gas). Problems with phase
transitions.
21
Free Energy as a Function of a Spatial
Coordinate, ?.
1. Work along a reversible path
Integration along a path
Derivative of free energy as function of ?
22
Free Energy as a Function of a Spatial
Coordinate, ?.
1. Work along a reversible path
Derivative of free energy as function of ?
? is a spatial coordinate therefore
(cartesian or internal )
Free energy profile (PMF) (average force
acting along, ?
)
PMF Potential of mean force
Method Simulate at given values of ?, average
force acting along ?, integrate.

23
Free energy barrier estimation of unfolding the
a-helical surfactant-associated polypeptide C
Compare the stability of lung epithelium
surfactant protein C in water and methanol
SP-C in water.
SP-C (LRIPCCPVNLKRLLVVVVVVVLVVVVIVGALLMGL)
Free energy unwinding two helix turns (Val25 to
Leu32)
water
methanol
Zangi, R., et al. Proteins 43 (2001) 395-402
24
Free Energy as a Function of a Spatial
Coordinate, ?.
2. Probability Distribution
Estimate F(?) from probability distribution along
?, P(?)
?, 1-dimension of multi-dimensional
space. sampling possible limited range (e.g. -?
lt ? lt ? for dihedral) sampling possible
25
Free Energy as a Function of a Spatial
Coordinate, ?.
2. Probability Distribution
2. Probability Distribution (cont.)
Estimate F(?) from probability distribution
along ?, P(?)
Example PMF from radial distribution function
(rdf).
PMF average force from other H2O molecules Use
in implicit solvent model to reduce vacuum
effects.
26
Free Energy as a Function of a Spatial
Coordinate, ?.
2. Probability Distribution (cont.)
Estimate F(?) from probability distribution
along ?, P(?)
Problems Sampling must be continuous and
sufficient for all ?.
27
Methods to Compute Free Energies Umbrella
Sampling.
Aim Bias or focus sampling in important regions
of phase space.
Method 1. Modify energy function V(r) to
restrict sampling only to region(s) of interest.
weighting function
28
Methods to Compute Free Energies Umbrella
Sampling.
Aim Bias or focus sampling in important regions
of phase space.
Method 1. Modify energy function V(r) to
restrict sampling only to region(s) of interest.
weighting function
e.g. harmonic restraining potential
2. Correct for Umbrella (weighting function).
Non-Boltzmann weights.
Correction Unbiased ensemble average from a
biased ensemble.
Umbrella Sampling Widely used general method to
improve sampling.
29
Methods to Compute Free Energies Umbrella
Sampling.
3. Umbrella Sampling. A. Modify potential energy
function to focus (bias) sampling along ?.
weighting function
e.g. harmonic restraining potential
B. Umbrella potential ? biased sampling
probability distribution P(?)U.
C. Determine F(?) (unbiased) from biased
distribution.
umbrella potential
biased distribution
unknown offset work to apply umbrella potential
30
Methods to Compute Free Energies Umbrella
Sampling.
3. Umbrella Sampling.
C. Determine F(?) (unbiased) from biased
distribution.
umbrella potential
biased distribution
unknown offset work to apply umbrella potential
31
Interaction of the cell-penetrating peptides with
lipid bilayers
Pulling of a single penetratin molecule through
a DPPC bilayer.
Semen Yesylevskyy
32
Thermodynamic Cycles Basic principles.
1. Free Energy is a State Function Difference in
free energy between two states is independent of
the path used to go between them.
33
Calculation of relative free energy via indirect
pathways
Mark, et al. (1991) J. Chem. Phys. 94, 3808-3816.
Optimal choice of path maximizes the relaxation
of environment
34
Thermodynamic Cycles Basic principles.
2. ?G for Any Closed Path 0 Thermodynamic
Cycle ? a closed path for which ?G 0
35
Thermodynamic Cycles Relative free energies
1. Hydration Free Energy Difference between A
and B
36
Thermodynamic Cycles (cont.).
2. Binding Free Energy Modification of a ligand
37
Methods to Compute Free Energies Coupling
Parameter Approach
Aim General expression for the difference in
free energy between two arbitrary states A and B.
Example mass of atom
Note A and B effectively become two states of a
single system
Hamiltonian ? function of ?, H(?) Partition
function ? function of ?, Z(?) Free Energy ?
function of ?, F(?)
38
Methods to Compute Free Energies Coupling
Parameter Approach
1. Perturbation Formula Express ?FAB as the
ratio of the probability of state A and state B.
Estimate ?FAB from the probability of finding a
configuration appropriate to state B H(?B) in
an ensemble of states generated at state A
H(?A) or visa versa
39
Methods to Compute Free Energies Coupling
Parameter Approach
1. Perturbation Formula (cont.)
In principle the perturbation formula is correct
for a mutation between any two states A and
B. ltgt refers to sampling over all possible
states. In practice the perturbation approach
only converges to the correct answer if there is
a strong overlap of configurational space (low
energy states) in state A and B. (i.e if the
difference between A and B is small).
Expressed generally as an integral over a series
of small changes.
perturbation formula
40
Methods to Compute Free Energies Coupling
Parameter Approach
2. Integration Formula Reversible work along a
path (defined by ? dependence of Hamiltonian).
Analytical derivative for F(?)
41
Methods to Compute Free Energies Coupling
Parameter Approach
2. Integration Formula Reversible work along a
path (defined by ? dependence of Hamiltonian).
Analytical derivative for F(?)
integration formula
Integral of the average derivative of the
potential w.r.t. an arbitrary coordinate
42
Methods to Compute Free Energies Coupling
Parameter Approach
1. Perturbation Formula
perturbation formula
Requirements 1. Equilibrium 2. Convergence of
the ensemble average
Note Convergence depends on overlap of low
energy regions of configurational space at ? and
? ??. ?? must have small effect on
configurations sampled
43
Methods to Compute Free Energies Coupling
Parameter Approach
integration formula

• Perform integration by
• 1. Making ? function of time and slowly change
  during simulation.
• (single-configurational TI or slow growth)
• 2. Simulating at fixed ? values and integrate
  numerically
• (multi-configurational TI or numerical
  quadrature)

44
Methods to Compute Free Energies Coupling
Parameter Approach
Example Difference in free energy of hydration
between water and methanol.
-0.548
-0.824
O
O
0.398
0.412
0.150
0.412
CH3
H
H
H
water
methanol
Gradually change water into methanol by changing
the charge, LJ and angle parameters during the
course of a simulation. Calculated the work done.
methanol
free energy
water
coupling parameter
45
Perturbation or Integration
V
Ensembles overlap Use perturbation
Derivative of free energy changes slowly with X
Use integration
C
C
x0
x0
X
46
Perturbation Phase space overlap
efficient
efficient
in efficient
in efficient
in efficient
47
Methods to Compute Free Energies Entropy and
Enthalpy
Q. Why not calculate entropy and enthalpy?
Free Energy
and
combine
Entropy
Derivative of S w.r.t. ?
Entropy is expressed in terms of the correlations
between
48
Methods to Compute Free Energies Entropy and
Enthalpy (cont.)
Thermodynamic Integration Formula for ?SAB
Enthalpy
Enthalpy difference, ?EAB is given by the
difference between the two end states.
NOTE ?F expressed only in terms of lt?H/??gt? i.e.
interactions which change as a function of ?
(local interactions). ?S and ?E depend on
ltHgt (all the interactions). Accuracy of ?S and
?E ltlt ?F (1-2 orders of magnitude)
49
Table 2 Summary of Methods to Calculate Free
Energies
50
Methods to Compute Free Energies Applications
Aim Estimate excess chemical potential,
µ. excess chemical potential excess free energy
per particle
Methods 1. Switch off all intermolecular
interactions (N,V,T). (i.e turn system into
an ideal gas and estimate work). 2. Add or
remove 1 particle (V,T or P,T).
excess free energy free energy in excess of
ideal gas (configurational part)
51
(Widom) Particle Insertion (Test Particle
Insertion)
Particle Insertion A perturbation approach

• For an equilibrium configuration (MD or MC) of N
  particles.
• Randomly attempt to add in 1 additional particle
  (single step).
• Estimate free energy using the perturbation
  formula.

?V change in V for N1 particles
State A N particles State B N1 particles
Randomly insert a test (ghost) particle M times)
52
(Widom) Particle Insertion (Test Particle
Insertion)
Dense systems
Low energy (good) configurations are never
sampled. No spontaneously formed cavities large
enough to accept test (ghost particle)
Particle insertion fails.
Solution Slowly grow particles, allow the
system to adapt.
53
(Widom) Particle Insertion (or Test Particle
Insertion)
Particle insertion works well for low density
systems but not for high density systems.
Question For what type of system would you
expect particle deletion to work well? Hint
For particle insertion you must sample locations
for which it is favourable to add particle (i.e.
a cavity). For particle deletion you need to
sample locations where it is unfavourable to have
a particle (i.e. where particles overlap).
54
(Widom) Particle Insertion (or Test Particle
Insertion)
Example Creation of a cavity (purely repulsive)
in SPC water
SPC Simple Point Charge (water model)
repulsive part of LJ
Particle insertion Works well if many
appropriately sized cavities are sampled. If no
appropriate cavities method fails
Thermodynamic integration Larger cavities can
be created because the system can relax.
55
Methods to Compute Free Energies
Example Creation of a cavity (purely repulsive)
in SPC water
Simulate at a number of different cavity sizes
then integrate.
Thermodynamic integration
Relative accuracy of the excess Gibbs free
energy (derivative) for different cavity sizes.
Depending on how the water molecules pack around
the cavity the free energy of certain cavity
sizes is harder to predict.
56
Binding of p-substituted phenols to ?-Cyclodextrin
Relatively small host-guest system
?-cyclodextrin 6 sugar (glucose) units (cyclic)
p-methoxyphenol showing the orientation of the
guest when inserted.
?-cyclodextrin,
System ?-cyclodextrin 500 H2O Constant
T,P Analyze a) Orientation and motion of ligand
in pocket b) Interaction energies c) Binding
free energies
57
Cyclodextrin
VDW surface of ?-cyclodextrin front cut away
?-cyclodextrin with p-chloro-phenol. Note, the
tight fit of the ring and the distortion of the
cyclodextrin
58
Motion of Guest in ?-Cyclodextrin
Side view of ?-CD with p-Cl-phenol showing motion
of the guest inside the cavity (water not shown).
?CD with p-Cl-phenol showing van der Waals
contacts.
59
Studies on Structure and Binding of Cyclodextrin
Motion of p-Cl-phenol in ?-cyclodextrin ring
(rotation of guest)
Relaxation times (guest) a) Rotational
averaging 20-40ps b) In/out motion 60-80ps c)
Tilt averaging gt80ps
60
Studies on Structure and Binding of Cyclodextrin
Part C. Estimation of Binding Energies of
p-Substituted Phenols to ?-CD
Comparison of host-guest interaction energies and
experimental binding free energies. (no entropic
contributions)
Average error between interaction energy and
binding free energy gt 12kJ/mol (4.6 kT or 2
orders of magnitude in association constants)
61
Studies on Structure and Binding of Cyclodextrin
Part D. Binding Free Energies of p-Substituted
Phenols to ?-CD Free Energy Calculations
Example Difference in Binding Free Energy
?Gbinding p-methylphenol ? p-chorophenol to
?-cyclodextrin

Coupling Parameter Approach
Express the Hamiltonian as a function
of a coupling parameter ? such that
Calculate the work to go from the initial to the
final state over a reversible path
integration formula
62
Studies on Structure and Binding of Cyclodextrin
Part D. Binding Free Energies of p-Substituted
Phenols to ?-CD Free Energy Calculations
63
Studies on Structure and Binding of Cyclodextrin
Part D. Binding Free Energies of p-Substituted
Phenols to ?-CD Free Energy Calculations
Difference in Binding Free Energy ??Gbinding
p-methyl phenol ? p-chloro phenol bound to
?-cyclodextrin (p-CH3)
(p-Cl)
Thermodynamic Cycle Express difference in
experimental binding free energies ?G1 and ?G2 in
terms of non-physical mutations ?G3 and ?G4
Free energy is a state function ? ?G is
independent of path. ?G (any cycle) 0
64
Methods To Calculate Free Energy Make
Hamiltonian dependent on ?.
? dependency determines pathway from A to B ?
0 state A ? 1 state B
65
Methods To Calculate Free Energy Slow growth
Integration Formula

Slow Growth 1. Make ? function of time and
slowly change during simulation
Approximate integral by sum (errors accumulate)
Compare ?Gforward and ?Greverse Estimate error
from hysteresis.
66
Methods To Calculate Free Energy Slow growth
Requirements (for each ?) A. System in
equilibrium ?equil. gt ?system B. Sufficient
sampling ?sample gtgt ?system ? relaxation
time
67
Studies on Structure and Binding of Cyclodextrin
Part D. Binding Free Energies of p-Substituted
Phenols to ?-CD Free Energy Calculations
Slow Growth Slowly change ? as a function of
time during simulation
68
Studies on Structure and Binding of Cyclodextrin
Part D. Binding Free Energies of p-Substituted
Phenols to ?-CD Free Energy Calculations
Slow Growth Slowly change ? as a function of
time during simulation
Hysteresis first increases then decreases with
longer sampling.
69
Studies on Structure and Binding of Cyclodextrin
Motion of p-Cl-phenol in ?-cyclodextrin ring
(rotation of guest)
Relaxation times (guest) a) Rotational
averaging 20-40ps b) In/out motion 60-80ps c)
Tilt averaging gt80ps
70
Methods To Calculate Free Energy Integration
Formula
integration formula
Numerical Quadrature
1. Simulate at fixed ? values and integrate
numerically
71
Studies on Structure and Binding of Cyclodextrin
Part D. Binding Free Energies of p-Substituted
Phenols to ?-CD Free Energy Calculations
Ensemble average of must converge.
72
Studies on Structure and Binding of Cyclodextrin
Part D. Binding Free Energies of p-Substituted
Phenols to ?-CD Free Energy Calculations
Ensemble average of must converge.
Change in free energy (kJ mol-1) as a function of
sampling time in water. The integral from ?
0 to 1 9 equally spaced points, Simpsons rule
approximation.
73
Studies on Structure and Binding of Cyclodextrin
Part D. Binding Free Energies of p-Substituted
Phenols to ?-CD Free Energy Calculations
Ensemble average of must converge.
Change in free energy (kJ mol-1) as a function of
sampling time in cyclodextrin. The integral
from ? 0 to 1 9 equally spaced points,
Simpsons rule approximation.
74
Methods To Calculate Free Energy
Example Binding Free Energies of p-Substituted
Phenols to ?-CD
Effect of Increasing the number of points used to
perform the integration.
F(?) must be a smooth function
75
Studies on Structure and Binding of Cyclodextrin
Part D. Binding Free Energies of p-Substituted
Phenols to ?-CD Free Energy Calculations
Comparison of Experimental and Calculated Binding
Free Energies.
Average difference between the experimental and
calculated binding free energies including
entropic contributions lt 3 kJ/mol (approx. 1kT or
factor of 2 in association constant)
76
Example Glycogen phosphorylase (GP)
Anti-diabetic target glucose-analog active site
inhibitors. Large number of high resolution
GP-inhibitor complex crystal structures (validate
calculations)
R Ki (mM)
1 H 3.1
2 NH2 146
3 OH 39.3
4 CH3 1200
5 NHCOCH3 550
R R2 Ki (mM)
6 NHCO2CH3 H 86
7 NHCOCH3 CONH2 310
8 NHCO2CH3 CONH2 16.5
9 H CONH2 370
10 NHCOCH3 H 32
77
TI cycle in water
TI cycle in protein
Cycles well converged in water and
protein. Accurately predicts structures of
derivatives. Force Field leads to deviations
5-10 kJ/mole from experimental binding
constants.
78
Multi-configurational TI (MCTI) Vs.
Single-configurational TI (slow-growth)
DG7,8 7 to 8 in water 8 to 7 in water SG
(20 ns) -81.10 80.54 hysteresis
0.56 MCTI -79.39 78.68 (15 l points)
The value of dG(l)/dl as a function of l for the
forward (black) and reverse (red) MCTI mutation
of 7 to 8. Error bars are also shown.
79
Practical Challenges
Avoiding singularities and numerical
instabilities when growing or deleting atoms.
Lennard-Jones potential Singularity as r-gt 0
Derivative has a singularity at rij 0
Beutler, T. C. et al. (1994) Chem. Phys. Lett.
222, 529-539.
80
Scale 0.01
81
(No Transcript)
82
Free Energies of Transfer of Trp Analogs from
Chloroform to Water
Example Simultaneously grow multiple atoms
Ac-Trp to Trp-NMe. Solid line water Broken line
methanol
J. Am. Chem. Soc. 1996, 118, 6285-6294
83
Non-Equilibrium Methods Jarzynski Equality
Given a sufficiently large sample of W
Described as remarkable, amazing,
unexpected.
Gives equilibrium work from non-equilibrium
simulations (or experiments) (closely related to
fluctuation dissipation theorem and reaction path
sampling)
Can be considered as coordinate transformation
where time is treated as a constraint (importance
sampling in trajectory space)
84
Very simple test case for non-equilibrium
pulling in a dissipative environment.
Move particle from location A to location B in
box of water.
?GA-gtB 0
B
A
85
Method 1. Reversible work

1. Very slowly pull particle from A to B (system
   always in equilibrium).
2. Determine average force.

x
trivial
B
A
Move particle from A to B
?GA-gtB 0
86
Non-equilibrium pulling
many paths w gt 0 some paths w 0 must find
paths for which w lt 0
x
longer the path faster I pull more likely w gt 0
must find path with w ltlt0
B
A
Move particle from A to B
?GA-gtB 0
87
Non-equilibrium pulling
Do any paths give w 0
x
Yes but rare
B
A
Move particle from A to B
?GA-gtB 0
88
Non-equilibrium pulling
Which if any paths give w lt 0
x
Particle must be pushed from A to B by the
environment
B
A
Move particle from A to B
?GA-gtB 0
89
Equilibrium versus non-equilibrium work Moving
particles in a condensed phase
Methane in water
90
Non-Equilibrium Methods Jarzynski Equality
Does is work YES

• Advantages
• Can be only option (AFM experiments).
• All of the tricks developed for equilibrium free
  energy calculations can be applied.
• Unified free energy calculations and keeps many
  theoreticians off the streets.

• Disadvantages
• Suffers same convergence problems as perturbation
  methods.
• Number of trials grows exponentially with the
  average dissipative work.
• Changes can only correspond to observed
  fluctuations in the system at equilibrium.
• Equilibrium methods generally more efficient

91
Hamiltonian (energy function)
Free Energy
Probability of finding the system in a given
state
How the system preferentially evolves in time.
If we know one we can derive the others Choice
coordinate system is arbitrary (cartesian,
parameter, trajectory) Must sample a
representative ensemble. Secret is to transform
the system or bias the sampling to improve
sampling.