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
+ The free energy is the energy left for once you paid the tax to entropy:
!G = !H " T!S Enthalpic ! ! ! !
Hydrogen bonds Polar interactions Van der Waals interactions ...
Entropic ! ! ! !
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)
Examples of factors determining the binding free energy Electrostatic interactions - Strength depends on microscopic environment (!) - Case of hydrogen bonds Neutral :
-1.2 ± 0.6 kcal/mol
Charge assisted :
-2.4 to -4.8 kcal/mol
EH-bond (solv.) - EH-bond (comp.) determines if H-bonds contributes to affinity or not
H H
H
O
O H O
H H O
H O
H O
H H
H O H
O
S
H
solvent
complex
Unpaired polar groups upon binding are detrimental Strong directional nature Specificity of molecular recognition
Free energy: statistical mechanics definition
G = !k BT ln(Z )
where
Z = # i e! " Ei
is the partition function Free energy differences between 2 states (bound/unbound, …) are, therefore, ratios of partition functions
# ZA & !G = GA " GB = "kBT ln % ( $Z ' B
Free energy simulations aim at computing this ratio using various techniques.
Relation with chemical equilibrium A’B’
A'B'] [ KA = Kb = [ A ] [ B]
A’B’ # " A+B
A ] [ B] [ KD = Ki = [ A'B']
A+B # "
Kb : binding constant, Kd : dissociation constant, Ki : inhibition constant
!G binding = "RTln K A = RTln K D = !H " T!S "Gbinding (kcal/mol) Weak asso.
KD (mol/l)
-2
-4
-6 10 -3
-8
-10 10 -6
-12
-14
10 -9
-16
Strong asso.
10 -12
Connection micro/macroscopic: thermodynamics and kinetics
e - RT!G Free Energy
Microscopic Structure
=
Absolute binding free energies: !G " KA
KA Association Constant
Biological function
Relative binding free energies: !!G " KA’ / KA Binding free energy profiles: !G(#) " KA, Kon, Koff
The free energy is the main function behind all process KA =
A) Chemical equilibrium
!G binding = RTln K A
+ A
B) Conformational changes
!G conf = kB Tln
PConf 1 PConf 2
C) Ligand binding
!G binding = kB Tln D) …
PUnbound PBound
R = kB N A
[ AB] [ A ] [ B]
B KD = 1 / KA
AB
Free energy: computational approaches # ZA & !G = GA " GB = "kBT ln % ( $Z ' B
Free energy simulations techniques aim at computing ratios of partition functions using various techniques.
Z = # i e! " Ei Sampling of important microstates of the system (MD, MC, GA, …)
Computation of energy of each microstate (force fields, QM, CP)
Connection micro/macroscopic: intuitive view E1, P1 ~ e-$E1
Expectation value
O =
1 Oie! "Ei # Z i
#
E2, P2 ~ e-$E2
E3, P3 ~ e-$E3
! "Ei
Where Z = e i is the partition function
E4, P4 ~ e-$E4
E5, P5 ~ e-$E5
Central Role of the Partition Function
Z = # i e! " Ei
O =
1 ! " Ei O e # i Z i
...
Expectation Value
E =
! ln(Z) = U !"
Internal Energy
" ! ln(Z ) % p = k BT $ # !V '& N ,T Pressure
G = -kBT ln(Z) Gibbs free energy
Molecular Modeling Principles 1) Modeling of molecular interactions Electrostatics
Van der Waals
Covalent bonds
Solvent
Free energy landscape
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
The CHARMM Force Field V=
!K
2 (b " b ) + b 0
Bonds
"K
+
!K
2 ( # " # ) # 0
Angles
2 ( ! # ! ) ! 0
Impropers
+
" K [1# cos(n ! # $ )] !
!
!
Dihedrals
qiq j 1 i> j 4 !" ri, j
+#
+$ 4!ij [(" ij /rij )12 # (" ij /rij ) 6 ] i> j
Ergodic Hypothesis MD Trajectory E
$ Dialanine
O
Ensemble
#
NVT simulation
NVE simulation 3N Spatial coordinates
Protein
? 1& 1 = % O(! ," )e# $ E (! ," )d! d" = % O(t)dt = O Z & 0
Time
Approx.
Free Energy Perturbation (FEP)
Thermodynamical Integration (TI)
Non Equilibrium Statistical Mechanics (Jarzynski)
Linear Interaction Energy (LIE)
Molecular Mechanics/PoissonBoltzmann/Surface area (MM-PBSA)
Quantitative Structure Activity Relationship (QSAR)
G k0
kiXi
!G = F(X)
(X is a descriptor)
CPU Time
Sampling, Approx.
Sampling, Exact
Free energy calculation: Main approaches
Approx.
Free Energy Perturbation (FEP)
Thermodynamical Integration (TI)
Non Equilibrium Statistical Mechanics (Jarzynski)
Linear Interaction Energy (LIE)
Molecular Mechanics/PoissonBoltzmann/Surface area (MM-PBSA)
Quantitative Structure Activity Relationship (QSAR)
G k0
kiXi
!G = F(X)
(X is a descriptor)
CPU Time
Sampling, Approx.
Sampling, Exact
Free energy calculation: Main approaches
Theoretical approaches for the estimation of binding affinities Without 3D structure of the complex - 2D - QSAR - 3D - QSAR
With 3D structure of the complex - Knowledge based approaches • regression based methods • potential of mean force - Free energy simulation - Linear interaction energy - Binding free energy decomposition (MM-PBSA, MM-GBSA)
Quantitative Structure Activity Relationships (QSAR) (Drug design) Assumption Chemical similarity of ligands
Similarity of biological response
Affinity is a function of the ligand physico-chemical properties
Advantage No need for structural information about the target
Requirements (drawbacks) Known affinities for a series of ligands Structurally related ligands or similar binding modes
2D - QSAR n structurally related molecules O R4
R1
N
O
N
R2
R3
Quantitative description
Measured activities
Descriptors (Xi): - Substituents: surface, volume, electrostatics (%), hydrophobicity (&), partial charges, etc... - Molecule : volume, MR, HOMO, dipole moment, etc...
X 1
"G bind = k 0 + ! k i X i !G bind = neural network (X i )
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 Limited to structurally related molecules Needs the experimental activity of a series ligands Not for ab initio studies Overfitting - Method for selecting the descriptor (genetic algorithm) - Estimation of the predictive ability (test set, randomization test, Leave-One-Out method, ...) Use limited to the descriptor’s range in the training set : 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... ?
O R4 O
R1
N N R3
R2
3D - QSAR Example : Comparative molecular field analysis (CoMFA) R.D. Cramer et al., JACS, 1988, 110, 5959
(x1,y1,z1)
(x2,y2,z2)
S
S
1
2
1
(x1,y1,z1)
. . .
Ligands mutually aligned Common 3D lattice
(x2,y2,z2)
E
E
1
2
. . .
2
" G b i n d
3 . . .n
Steric field
PLS
Electrostatic field
"G bind = k 0 + ! $ iSi + ! # i E i
(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
!G bind = !G 0 + !G polar + !G apolar + !G solv + !G flexi Trained using a 82 protein-ligand complexes dataset with known experimental "Gbind
Polar interactions
Desolvation effect
(
)
!G polar = !G hb $ f ( !R, !" ) # f N neighb # fpcs hb
(
)
+!G ion $ f ( !R, !" ) # f N neighb # fpcs ion
+!G esrep N esrep
Active site filled with water molecules
!G solv = !G lipo wat
MD
Unbound water molecules
" unbound water
"Glipo water = -0.33 kcal/mol
"Ghb = -0.81, "Gion = -1.41 and "Gesrep = +0.10 kcal/mol
Ligand flexibility
Apolar interactions !G apolar = !G lipo A lipo + !G aro " f (R) aro
"Glipo = -0.81 and "Garo = -0.62 kcal/mol
!G flex = !G rot N rot "Grot = +0.26 kcal/mol Nrot : number of rotable bonds (acyclic sp3-sp3, sp3-sp2)
Knowledge based approaches. Regression based methods Example : LUDI score 82 complexes of the training set SD ~2 kcal/mol
Advantages :
Drawbacks :
- Allows identification of high affinity ligands
- Somewhat large errors
- Rapid estimation of the affinities
- Method biased: • certain type of proteins • only good complementarity protein/ligand - Some interactions ignored (cation – &)
- Structurally different ligands - “Universal” (different proteins)
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
# j " ij (r) & A ij ( r ) = !k b Tln % fvol_corr (r) ij ( "bulk ' $
%ij(r) : number density of atom pairs of type ij at distance r %ijbulk(r) : number density of atom pairs of type ij in a sphere with radius R f jvolv_corr : ligand volume correction
16 protein atom types, 34 ligand atom types PMF (kcal/mol)
NC positively charged nitrogen OC negatively charged oxygen ND nitrogen as hydrogen bond donor OA oxygen as hydrogen bond acceptor OD oxygen as hydrogen bond donnor Atom pair distance (Å)
Atom pair distance (Å)
PMF score = !
Atom pair distance (Å)
!
ij kl of type ij r
