Free energy simulations: theory and applications O. Michielin (1,2,4) - - PowerPoint PPT Presentation

Free energy simulations: theory and applications

• O. Michielin(1,2,4)

(1) Ludwig Institute for Cancer Research Epalinges, Switzerland (2) Swiss Institute for Bioinformatics Epalinges, Switzerland (4) Centre Pluridiscipinaire d'oncologie CHUV, Lausanne, Switzerland

Dynamical aspects of molecular recognition

Free energy: classical definition

+

Enthalpic Entropic

Hydrogen bonds Polar interactions Van der Waals interactions ... Loss of degrees of freedom Gain of vibrational modes Loss of solvent/protein structure ...

Theoretical Predictions: Approximate: empirical formula for all contributions

Exact: using statistical physics definition of G

G = -KBT ln(Z)

G = H TS

The free energy is the energy left for once you paid the tax to entropy:

Examples of factors determining the binding free energy

Electrostatic interactions

• Strength depends on microscopic environment ()
• Case of hydrogen bonds
• 2.4 to -4.8 kcal/mol

Charge assisted :

• 1.2 ± 0.6 kcal/mol

Neutral :

O H O H H O H H O H H O H H O H H O H O H S

solvent complex EH-bond (solv.) - EH-bond (comp.) determines if H-bonds contributes to affinity or not Unpaired polar groups upon binding are detrimental Strong directional nature Specificity of molecular recognition

Free energy: statistical mechanics definition

G = kBT ln(Z) Z = eEi

i

• where

is the partition function Free energy differences between 2 states (bound/unbound, …) are, therefore, ratios of partition functions G = GA GB = kBT ln ZA ZB

• Free energy simulations aim at computing this ratio using various

techniques.

Relation with chemical equilibrium

A + B

• A’B’

A + B

• A’B’

Kb : binding constant, Kd : dissociation constant, Ki : inhibition constant KD (mol/l) Gbinding (kcal/mol)

• 2
• 4
• 6
• 8
• 10
• 12
• 14
• 16

10 -12 10 -9 10 -6 10 -3

Weak asso. Strong asso.

KD = Ki = A

[ ] B [ ]

A'B'

[ ]

KA= Kb = A'B'

[ ]

A

[ ] B [ ]

Gbinding = RTlnKA = RTlnKD = H TS

Connection micro/macroscopic: thermodynamics and kinetics

Free Energy Association Constant

e - G/RT = KA

Microscopic Structure Biological function

Relative binding free energies: G KA’ / KA Absolute binding free energies: G KA Binding free energy profiles: G() KA, Kon, Koff

The free energy is the main function behind all process

A) Chemical equilibrium B) Conformational changes C) Ligand binding D) …

+

Gbinding = RTlnKA Gconf = kBTln P

Conf 1

P

Conf 2

Gbinding = kBTln P

Bound

P

Unbound

KA = AB

[ ]

A

[ ] B [ ]

A B AB KD = 1/ KA

R = kBNA

Free energy: computational approaches

G = GA GB = kBT ln ZA ZB

• Free energy simulations techniques aim at computing ratios of

partition functions using various techniques. Z = eEi

i

• Sampling of important

microstates of the system (MD, MC, GA, …) Computation of energy

• f each microstate

(force fields, QM, CP, …)

Connection micro/macroscopic: intuitive view

E1, P1 ~ e-E1 E2, P2 ~ e-E2 E3, P3 ~ e-E3 E4, P4 ~ e-E4 E5, P5 ~ e-E5

Where is the partition function Expectation value

O = 1 Z OieEi

i

• Z =

eEi

i

Central Role of the Partition Function

G = -kBT ln(Z)

. . .

Expectation Value Internal Energy Pressure Gibbs free energy

Z = eEi

i

• O = 1

Z OieEi

i

• E =
• ln(Z) = U

p = kBT ln(Z) V

• N,T

Molecular Modeling Principles

1) Modeling of molecular interactions 2) Simulation of time evolution (Newton) 3) Computation of average values

O = < O >Ensemble = < O >Temps (Ergodicity) Macroscopic value Average simulation value

Connection microscopic/ macroscopic Free energy landscape

Electrostatics Van der Waals Covalent bonds Solvent

The CHARMM Force Field

V = Kb

Bonds

• (b b0)2 +

K

Angles

• ( 0)2

+ K

Impropers

• ( 0)2

+ K

Dihedrals

• 1 cos(n )

[ ]

+ qiq j 4

i> j

• 1

r

i, j

+ 4ij ( ij /r

ij)12 ( ij /r ij)6

[ ]

i> j

• E

3N Spatial coordinates

Ergodic Hypothesis

NVT simulation

?

MD Trajectory

NVE simulation

“Alanine” Protein

O Ensemble = 1 Z O(, )eE(, )

• dd = 1
• O(t)dt
• = O Time

Free energy calculation: Main approaches

Free Energy Perturbation (FEP) Thermodynamical Integration (TI)

CPU Time

Linear Interaction Energy (LIE) Molecular Mechanics/Poisson- Boltzmann/Surface area (MM-PBSA) Quantitative Structure Activity Relationship (QSAR) Non Equilibrium Statistical Mechanics (Jarzynski)

Sampling, Exact Sampling, Approx. Approx.

G k 0 k iX i

(X is a descriptor)

G = F(X)

Free energy calculation: Main approaches

Free Energy Perturbation (FEP) Thermodynamical Integration (TI)

CPU Time

Linear Interaction Energy (LIE) Molecular Mechanics/Poisson- Boltzmann/Surface area (MM-PBSA) Quantitative Structure Activity Relationship (QSAR) Non Equilibrium Statistical Mechanics (Jarzynski)

Sampling, Exact Sampling, Approx. Approx.

G k 0 k iX i

(X is a descriptor)

G = F(X)

Theoretical approaches for the estimation of binding affinities

Without 3D structure of the complex

• 2D - QSAR

With 3D structure of the complex

• Binding free energy decomposition (MM-PBSA, MM-GBSA)
• Knowledge based approaches
• regression based methods
• potential of mean force
• Free energy simulation
• Linear interaction energy
• 3D - QSAR

Quantitative Structure Activity Relationships (QSAR)

Advantage

No need for structural information about the target

Assumption

Affinity is a function of the ligand physico-chemical properties

Requirements (drawbacks)

Known affinities for a series of ligands Structurally related ligands or similar binding modes Chemical similarity of ligands Similarity of biological response

2D - QSAR

N N O R1 R2 R3 O R4

• +

=

• i

i bind

X k k G

n structurally related molecules Quantitative description Measured activities

( )

i bind

X network neural G =

• Descriptors (Xi):
• Substituents: surface, volume,

electrostatics (), hydrophobicity (), partial charges, etc...

• Molecule : volume, MR, HOMO,

dipole moment, etc...

• C. Hansch and T. Fujita, JACS, 1964, 86, 1616

S.S. So and M. Karplus, J. Med. Chem., 1996, 39, 1521 S.S. So and M. Karplus, J. Med. Chem., 1996, 39, 5246

2D - QSAR: limitations

Needs the experimental activity of a series ligands Not for ab initio studies Use limited to the descriptor’s range in the training set :

N N O R1 R2 R3 O R4

If only hydrophobic groups at R1 in the training set Influence of a hydrophilic group at R1 ? If only methyl, ethyl, propyl, butyl at R1 in the training set Contribution of pentyl, hexyl, etc... ? Limited to structurally related molecules Overfitting

• Method for selecting the descriptor (genetic algorithm)
• Estimation of the predictive ability

(test set, randomization test, Leave-One-Out method, ...)

3D - QSAR

Example : Comparative molecular field analysis (CoMFA)

R.D. Cramer et al., JACS, 1988, 110, 5959

• +

+ =

• i

i bind

E S k G

i i

• Steric field

Electrostatic field

(x1,y1,z1) (x2,y2,z2) (x1,y1,z1) (x2,y2,z2)

Ligands mutually aligned Common 3D lattice PLS

(xi,yi,zi) (xj,yj,zj)

3D - QSAR: limitations

Needs the experimental activity of a series of ligands Not limited to structurally related molecules Alignment of the molecules in their bioactive conformation. Use of:

• known structure of a complex (QSAR?)
• conformationally rigid example in the dataset
• functional groups in agreement with a

pharmacophore hypothesis Others : CoMSIA, HASL, Compass, APEX-3D, YAK, ...

Knowledge based approaches. Regression based methods

Example : LUDI score

H.J. Bohm, J. Comput.-Aided Mol. Des., 1994, 8, 623 H.J. Bohm, J. Comput.-Aided Mol. Des., 1998, 12, 309

Trained using a 82 protein-ligand complexes dataset with known experimental Gbind

Polar interactions

Ghb = -0.81, Gion = -1.41 and Gesrep = +0.10 kcal/mol

Apolar interactions

Glipo = -0.81 and Garo = -0.62 kcal/mol

Desolvation effect

Active site filled with water molecules MD Unbound water molecules

Glipo water = -0.33 kcal/mol

Ligand flexibility

Grot = +0.26 kcal/mol Nrot : number of rotable bonds (acyclic sp3-sp3, sp3-sp2)

Gbind = G0 + Gpolar + Gapolar + Gsolv + Gflexi

Gpolar = Ghb f R,

( ) f Nneighb

( ) fpcs

hb

• +Gion

f R,

( ) f Nneighb

( ) fpcs

ion

• +GesrepNesrep

Gapolar = GlipoAlipo + Garo f (R)

aro

• Gsolv = Glipo wat

unbound water

• Gflex = GrotNrot

Knowledge based approaches. Regression based methods

Example : LUDI score SD ~2 kcal/mol Advantages : Drawbacks :

• Rapid estimation of the affinities
• Structurally different ligands
• Allows identification of high affinity ligands
• “Universal” (different proteins)
• Method biased:
• Somewhat large errors
• certain type of proteins
• only good complementarity

protein/ligand

• Some interactions ignored (cation – )

82 complexes of the training set Others : ChemScore, VALIDATE

Knowledge based approaches. Potential of mean force

Example : PMF score

• I. Muegge et al., J. Med. Chem., 1999, 42, 791
• I. Muegge et al., Persp. In Drug Disc. And Des., 2000, 20, 99

Trained using 697 complexes from the PDB. No need for experimental Gbind

atom type i for protein and j for ligand ij(r) : number density of atom pairs of type ij at distance r ij

bulk(r) : number density of atom pairs of type ij in a sphere with radius R

f j

volv_corr : ligand volume correction

16 protein atom types, 34 ligand atom types

PMF (kcal/mol)

Atom pair distance (Å) Atom pair distance (Å) Atom pair distance (Å)

NC positively charged nitrogen ND nitrogen as hydrogen bond donor OC negatively charged oxygen OD oxygen as hydrogen bond donnor OA oxygen as hydrogen bond acceptor

Aij r

( ) = kbTln fvol_corr

j

(r) ij(r) bulk

ij

• PMF score =

Aij(r)

kl of type ij r<r

• ij

Knowledge based approaches. Potential of mean force

Example : PMF score

PMF score

log Ki (experimental)

77 complexes, 5 different proteins SD ~2 kcal/mol Advantages : Drawbacks :

• Rapid estimation of the affinities
• Structurally different ligands
• Allows identification of high affinity ligands
• “Universal” (different ligand and protein types)
• Somewhat large errors
• No measure of directionality of H-bonds
• No fitting parameters to measured Gbind

Others : SMoG-Score, BLEEP, DrugScore

Free energy calculation: Main approaches

Free Energy Perturbation (FEP) Thermodynamical Integration (TI)

CPU Time

Linear Interaction Energy (LIE) Molecular Mechanics/Poisson- Boltzmann/Surface area (MM-PBSA) Quantitative Structure Activity Relationship (QSAR) Non Equilibrium Statistical Mechanics (Jarzynski)

Sampling, Exact Sampling, Approx. Approx.

G k 0 k iX i

(X is a descriptor)

G = F(X)

Linear interaction energy (LIE)

• J. Åqvist, J. Phys. Chem., 1994, 98, 8253

Free state “solvent” = water Bound state “solvent” = water + protein

Two MD runs : free state and bound state

• J. Åqvist, J. Phys. Chem., 1994, 98, 8253

=0.165 and =0.5

• T. Hansson et al., J. Comp.-Aided Molec. Design, 1998, 12, 27

=0.181 and =0.5, 0.43, 0.37, 0.33

• W. Wang, Proteins, 1999, 34, 395

function of binding site hydrophobicity

Gbind = Els

vdw bound Els vdw free

( ) +

Els

elec bound Els elec free

( )

Linear interaction energy (LIE)

Advantages : Drawbacks :

• Faster than free energy simulation
• More structurally different ligands

than for free energy simulation. But generally restricted to rather similar ligands.

• Slower than scores based on a

single conformation (LUDI, PMF, ...)

• Not really universal

( and system dependent) Modifications :

• Additional term proportional to buried surface upon complexation
• Use of continuum solvent model instead of explicit solvent
• R. Zhou and W.L. Jorgensen et al., J. Phys. Chem., 2001, 105, 10388

D.K. Jones-Hertzog and W.L. Jorgensen, J. Med. Chem., 1997, 40, 1539

• Need experimental binding

affinities of known complexes

elec lig prot

E

• desolv

G

• buried

SASA

Electrostatic interactions between protein and ligand

Electrostatic contribution to desolvation energy (continuum model)

Solvent accessible surface of protein and ligand buried upon complexation

New LIE-like method

• V. Zoete, O. Michielin and M. Karplus, J. Comput.-Aided Molec. Design, 2004, in press

MD simulation of the complex without explicit solvent , and fitted using a training set

Gbind = Eprotlig

elec

+ Gdesolv + SASAburied

New LIE-like method

Training set :

Trained with a set of 16 known inhibitors of the HIV-1 protease belonging to different chemical families Used to rank new HIV-1 protease ligands

New LIE-like method

SD ~ 0.85 kcal/mol Training set New ligands Advantages : Drawbacks :

• Faster than free energy simulation

LIE

• Ligands with very different chemical

structures

• Slower than scores based on a

single conformation

• Not universal
• Need experimental binding

affinities of known complexes

• Restricted to ligands uncharged

when bound

• 17
• 15
• 13
• 11
• 9
• 17
• 15
• 13
• 11
• 9
• Exp. Gbind (kcal/mol)
• Calc. Gbind (kcal/mol)
• 14
• 13
• 12
• 11
• 10
• 9

4 5 6 7 8 9

• Exp. pIC50
• Calc. Gbind (kcal/mol)

Free energy calculation: Main approaches

Free Energy Perturbation (FEP) Thermodynamical Integration (TI)

CPU Time

Linear Interaction Energy (LIE) Molecular Mechanics/Poisson- Boltzmann/Surface area (MM-PBSA) Quantitative Structure Activity Relationship (QSAR) Non Equilibrium Statistical Mechanics (Jarzynski)

Sampling, Exact Sampling, Approx. Approx.

G k 0 k iX i

(X is a descriptor)

G = F(X)

Binding free energy decomposition: MM-PBSA, MM-GBSA

Lig + Prot Lig:Prot Lig:Prot

Gbind

Lig + Prot

Gaz Sol

Averaged over an MD simulation trajectory

• f the complex (and isolated parts)
• B. Tidor and M. Karplus, J. Mol. Biol., 1994, 238, 405

Molecular mechanics – Poisson-Boltzmann Surface Area (MM- PBSA) Molecular mechanics – Generalized Born Surface Area (MM- GBSA)

• J. Srinivasan, P.A. Kollmann et al., J. Am. Chem. Soc., 1998, 120, 9401
• H. Gohlke, C. Kiel and D.A. Case, J. Mol. Biol., 2003, 330, 891

Depending on the way Gsolv,elec is calculated:

Gbind = Egaz + Gdesolv T S Egaz = Eelec + Evdw + Eint ra Gdesolv = Gsolv

comp Gsolv lig + Gsolv prot

( )

TS = T(Scomp (S prot + Slig))

Gsolv

lig

Gsolv

prot

Gsolv

comp

S = Strans + Srot + Svib

Gsolv = Gsolv,elec + Gsolv,np

Gdesolv = Gsolv,elec

comp

Gsolv,elec

lig

+ Gsolv,elec

prot

( ) + SASAcomp SASAlig + SASA prot ( )

( )

Egaz

Expérimental Théorique

MM-GBSA Method: application to TCR-p-MHC

Gbind = Egaz + Gdesolv T S

MM-GBSA Method: application to TCR-p-MHC

Gbind = Egaz + Gdesolv T S

Examples of TCR optimization: 2C TCR

Improvement of favorable TCR residues Replacement of unfavorable TCR residues

Examples of TCR optimization: 2C TCR

Binding free energy decomposition

Advantages : Drawbacks :

• Used for ligand:protein and protein:protein

complexes

• Could be applied to structurally different

ligands (but in fact applied to similar

• nes)
• “Universal” (no parameter to be fitted)
• Rather time consuming
• In some cases, found unable to

rank ligands

• TS is necessary to find the order of magnitude
• f the absolute binding free energies but, in some

cases, it is not necessary to estimate relative binding free energies

• MM-GBSA allows a per-atom decomposition
• f Gbind (e.g. contribution of side chains)

MM- PBSA, MM-GBSA

• H. Gohlke, C. Kiel and D.A. Case, J. Mol. Biol., 2003, 330, 891
• W. Wang and P.A. Kollman, J. Mol. Biol., 2000, 303, 567

Free energy calculation: Main approaches

Free Energy Perturbation (FEP) Thermodynamical Integration (TI)

CPU Time

Linear Interaction Energy (LIE) Molecular Mechanics/Poisson- Boltzmann/Surface area (MM-PBSA) Quantitative Structure Activity Relationship (QSAR) Non Equilibrium Statistical Mechanics (Jarzynski)

Sampling, Exact Sampling, Approx. Approx.

G k 0 k iX i

(X is a descriptor)

G = F(X)

Free energy calculation: Main approaches

Free Energy Perturbation (FEP) Thermodynamical Integration (TI)

CPU Time

Linear Interaction Energy (LIE) Molecular Mechanics/Poisson- Boltzmann/Surface area (MM-PBSA) Quantitative Structure Activity Relationship (QSAR) Non Equilibrium Statistical Mechanics (Jarzynski)

Sampling, Exact Sampling, Approx. Approx.

G k 0 k iX i

(X is a descriptor)

G = F(X)

“Alchemical” Free Energy Calculations

Alchemical free energy formalism

Hybrid Side Chain for P A Mutation

Free energy formalism

From the statistical definition of the free energy,

Results of the Free Energy Simulations

Total (Path Independent)

Experimental: 2.9 (0.2) kcal/mol Theoretical: 2.9 (1.1) kcal/mol

Components (Path Dependent)

Experimental:

• Theoretical:

TCR 25% Solvent 20% HLA A2 40% Peptide 15%

} }

45% 55%

Free energy derivative

100% Pro 100% Ala

Free energy derivative components

100% Pro 100% Ala

Solvent Contribution

Solvent contribution to A:

• vdW:

+1.4 kcal/mol

• Elec:
• 0.7 kcal/mol
• Total:

+0.7 kcal/mol The solvent favors association with Tax(P6)-A2 compared to the A6 mutant that is more so- luble.

p-MHC conformational change contribution

p-MHC conformational change observed upon TCR binding is indirectly computed along the horizontal paths:

• the free energy cost of the PA mutation is higher in the bound

simulation because of the numerous favorable interaction of P6 with HLA A2. This brings a net addition of +1.2 kcal/mol to A.

Combined approach for peptide design

p-MHC affinity TCR-(p-MHC) affinity

Check on repertoire selection Select modifications that increase MHC affinity

Concluding Remarks

1) Free energy simulations can reproduce accurately experimental changes in association constant between to closely related protein systems if detailed structural knowledge is available (X-ray, NMR or model) 2) The formalism is exact from a statistical physics stand point and accurate treatment of entropic terms, solvent effect or conformational changes can be obtained 3) Convergence of the free energy derivative is still problematic. The situation should improve with new methodological enhancements as well as longer simulation time 4) Absolute free energies can also be computed but the convergence is even more difficult 5) Much details about the specificity of the association can be gained using component analysis, opening the door to rational peptide or protein design

Free energy calculation: Main approaches

Free Energy Perturbation (FEP) Thermodynamical Integration (TI)

CPU Time

Linear Interaction Energy (LIE) Molecular Mechanics/Poisson- Boltzmann/Surface area (MM-PBSA) Quantitative Structure Activity Relationship (QSAR) Non Equilibrium Statistical Mechanics (Jarzynski)

Sampling, Exact Sampling, Approx. Approx.

G k 0 k iX i

(X is a descriptor)

G = F(X)

time reaction

Independent starting pts (canonical ensemble) reference trajectory

Wadia = G

KA

(Infinitely slow) W = Wadia + Wdiss (Finite rate) Let G be the free energy and W the work,

W

C

V

Pulld

Pull

Computation of absolute TCR binding free energy

(Jarzynsky)

eW = eG

“Proof” of the Jarzynski equation: Park & al. 2004

Suppose a system S in equilibrium with a thermal bath B at temperature T, and

composite system SB is thermally isolated interaction energy of SB is negligible:

where and denotes the phases (p,q) of S and B, resp. then

H

SB(,) = H S() + H B()

Z

SB =

ddexp H

SB(,)

[ ]

• =

dexp H

S()

[ ]

• dexp H B()

[ ]

• = Z

SZ B

Let’s now consider a process for which the (time dependent) parameter goes from 0 to and drives the system from (0,0) to (,). The work done during the process is

W = H

SB( , ) H0 SB(0,0)

Since SB is an isolated system, it follows a micorcanonical distribution. However, for large system, one can use the thermodynamical limit and use a canonical distribution:

1 Z0

SB exp H0 SB(0,0)

[ ]

“Proof” of the Jarzynski equation: Park & al. 2004

The expectation value for the exponential work is

W = H

SB( , ) H0 SB(0,0)

eW = d0

• d0

1 Z0

SB exp H0 SB(0,0)

[ ]

x exp H

SB( , ) H0 SB(0,0)

[ ] { }

= d0

• d0

1 Z0

SB exp H SB( , )

[ ]

The last step consists in transforming the integration variable from (0,0) to (,). According to Liouville’s theorem (cave!) the Jacobian of the transformation is unity and d0d0 = dd

eW = d

• d

1 Z0

SB exp H SB( , )

[ ]

= Z

SB

Z0

SB

Using the negligible interaction between S and B (see above)

eW = Z

SB

Z0

SB

= Z

SZ B

Z0

SZ B

= Z

S

Z0

S

= exp(GS)

eW = exp(GS)

i.e

Simulation setup

• Gromos96 Force Field
• Gromacs Engine
• Particle Mesh Ewald (PME)
• Periodic boundary conditions
• Box: 80x80x150 A
• 26000 Water molecules
• 85000 Atoms
• Hydrogen shaken
• 2 fs timestep
• 0.5 ns / 24h on 4 alpha CPU

Absolute free energy results

Pull

Binding free energy profile:

eW = eG

8.3 13 10.6 11 19 24 Exp. Jarzynski

Agreement with experimental values

Free energy simulation: conclusion

Advantages : Drawbacks :

• Rigorous
• Estimates influence of small

modifications

• Partitioning of the free energy (TI)
• Restricted to small mutations of

ligand or protein

• Time consuming
• S. Boresch et al., Proteins, 1994, 20, 25
• S. Boresch and M. Karplus, J. Mol. Biol., 1995, 254, 801

... ) ( H ) ( H G

1 angles 1 bonds bind

+

• +
• =
• =

= = =

• d

d ... G G G

angles bond bind

+

• +
• =
• No parameter to be fitted
• Most often: relative Gbind