Combined QMMM studies of enzymes

1
Combined QM/MM studies of enzymes
Walter Thiel Max-Planck-Institut für
Kohlenforschung, Mülheim

• Introduction
• p-Hydroxybenzoate hydroxylase
• Cytochrome P450

2
QM/MM approach General overview
QM ab initio, DFT, semiempirical MM standard
force field
Inner subsystem
QM MM interactions electronic embedding

• Border region
• hydrogen link atoms L
• charge shift for q(M1)

Codes ChemShell as control module Interfaces to
standard QM and MM codes
H. M. Senn and W. Thiel, Top. Curr. Chem. 268,
173-290 (2007).
3
ChemShell A modular QM/MM package
Chemshell

CHARMM27academic
GAUSSIAN98
Tcl scripts
TURBOMOLE
CHARMm26MSI
Integratedroutines
GAMESS-UK
datamanagement
GROMOS96
geometryoptimisation
DL_POLY
moleculardynamics
GULP
genericforce fields
QM/MMcoupling
QM codes
MM codes
P. Sherwood et al, J. Mol. Struct. Theochem 632,
1-28 (2003).
4
Exploring potential surfaces of complex systems
Goal Compute barriers for enzymatic
reactions. Molecular dynamics - thermodynamic
integration Determine free energy barriers by
performing molecular dynamics simulations along
an assumed reaction path and integrating over the
resulting constraint forces. Geometry
optimization Determine energy barriers by
locating representative transition states and the
associate reactants and products.
S. R. Billeter, A. J. Turner, and W. Thiel, Phys.
Chem. Chem. Phys. 2, 2177 (2000).H. M. Senn, S.
Thiel, and W. Thiel, J. Chem. Theory Comput. 1,
494 (2005).
5
PHBH p-hydroxybenzoate hydroxylase
6
Aromatic hydroxylation of p-hydroxybenzoate

• rate-determining step oxygen transfer
• from cofactor FADHOOH to p-OHB
• (FAD flavin adenine dinucleotide)
• electrophilic aromatic substitution with
• heterolytic cleavage of the peroxide bond
• activation energy 12 kcal/mol

7
PHBH General setup
8
PHBH Motion in transition state (40 ps snapshot)
9
PHBH Active site
10
PHBH Comparing different snapshots
80 ps
1.63
2.06
1.57
1.49
1.83
1.86
165.2
1.28
2.72
1.50
1.50
2.92
2.20
92.1
23.3
11
Thermodynamic integration
12
(No Transcript)
13
PHBH Role of Pro293
Distance / Å
dRC / a0
14
PHBH DE versus DA for different snapshots
30.0
25.0
20.0
DE,DA / kcal/mol
15.0
10.0
5.0
0.0
40 ps
80 ps
120 ps
160 ps
200 ps
15
Umbrella sampling biased MD

• MD with a restraint (bias) on the reaction
  coordinate
• Windows with different ?i are sampled
• The biased distribution of ?, Pib(?), is sampled

bias potential
potential energy surface
G. M. Torrie and J. P. Valleau, Chem. Phys. Lett.
28, 578 (1974).
16
Umbrella sampling unbiasing

• Unbiasing provides the free energy for each
  window
• Unknown constant Fi for each window
• Combination of the windows (weighted average
  and estimation of Fi)
• Weighted histogram analysis method (WHAM), or
• Umbrella integration

WHAM M. A. Ferrenberg and R. H. Swendsen, Phys.
Rev. Lett. 61, 2635 (1988).
17
Umbrella integration method

• Conceptual combination of thermodynamic
  integration and umbrella sampling
• Analysis of umbrella sampling data
• Calculate the mean force for each window
• Combine the windows by a weighted sum
• Integrate to obtain A(?)

J. Kästner and W. Thiel, J. Chem. Phys. 123,
144104 (2005).
18
Umbrella integration weighted average

• Reaction coordinate is divided into bins of
  uniform width
• The unbiased mean forces of the windows are
  averaged on the grid provided by the bins
• Weight with Ni being the number of MD
  steps for window i
• Numerical integration yields DA

19
Umbrella integration normal distribution

• Full distribution Pib(?)
• Normal distribution of Pib(?) through
• Truncation of Ai(?) after the quadratic term in
  ?
• Truncation of a cumulant expansion of Pib(?)
• Results depend only on the mean and
  the variance for each window

20
Power series truncation
analytic example

• Noise reduction
• Linear and quadratic contributions contain
  relevant information
• Higher terms (residuum) predominantly contain
  noise
• The central region contributes mainly (gt 50 ) to
  A(?)

21
Umbrella integration analytic potential

• 2-dimensional function
• Barrier between two minima
• Monte Carlo sampling, T 300 K
• 20,000 steps in each of 40 windows
• Results
• Umbrella integration converges with the number of
  bins. Errors in barrier heights 0.097 and -0.136
  kJ/mol.With better sampling (80 windows, 80,000
  steps) -0.013 and -0.035 kJ/mol
• WHAM does not converge with the number of bins.
  Errors with 4500 bins 1.040 and 0.831 kJ/mol

22
Application PHBH - results

• Snapshot after 40 ps
• Molecular dynamics, T300K
• 8000 steps in each of 38 windows
• Tests for equilibration of and (sib)2
• Activation barrier
• Umbrella integration 101.5 kJ/mol
• WHAM 100.1 102.3 kJ/mol, depending on
  the number of bins
• Thermodynamic integration 1012 kJ/mol

?
J. Kästner and W. Thiel, J. Chem. Phys. 123,
144104 (2005).
23
Error analysis summary

• Data collection in each window
• Combining the windows
• Integration
• Confidence interval (95)
• This estimate only covers the statistical error,
  not the systematic error.

J. Kästner and W. Thiel, J. Chem. Phys. 124,
234106 (2006).
24
Umbrella integration summary

• Combines thermodynamic integration and umbrella
  sampling
• Advantages over WHAM analysis
• Enables control of the equilibration of the
  system
• Analysis independent of the bin width
• Error bar estimate is available
• Advantages over thermodynamic integration
• Metric tensor correction is avoided
• Easier implementation of new types of reaction
  coordinates

J. Kästner and W. Thiel, J. Chem. Phys. 123,
144104 (2005).
25
QM/MM free-energy perturbation (FEP)

• Reaction profile
• Full QM/MM calculations
• QM and MM atoms optimized

• Sampling
• Frozen QM part
• Density replaced by ESP charges

Y. Zhang, H. Liu, W. Yang J. Chem. Phys. 112,
3483 (2000)
26
FEP applied to PHBH
?A 101 2 108.2 1.0 112.3
?rA 212 2 198.6 1.3 184.2

• Good agreement between FEPand termodynamic
  integration

J. Kästner, H.-M. Senn, S. Thiel, N. Otte, and W.
Thiel, J. Chem. Theory Comput. 2, 452 (2006).
27
PHBH B3LYP/GROMOS results (TZVP basis)
a) Snapshots labeled TI taken from MD of
thermodynamic integration.
28
SP LMP2/GROMOS barriers (TZ basis)
Activation energy DE / kcal/mol
LMP2/TZ B3LYP/TZ
Snapshot
29
SP LMP2/GROMOS and LCCSD(T0)/GROMOS barriers (TZ
basis)
Activation energy DE / kcal/mol
LCCSD(T0)/TZDZ LMP2/TZ B3LYP/TZ
Snapshot
30
PHBH Comparison of barriers
QM/GROMOS results (kcal/mol) at B3LYP geometries
QM method Range Average rms B3LYP (a)
5.2 - 9.6 7.9 1.3 DF-LMP2 (b)
9.0 - 13.6 12.0 1.3 DF-LCCSD(T0) (b) 11.1
- 16.6 14.6 1.6 a) TZVP basis. b) cc-pVTZ
basis in general, aug-cc-pVTZ for O.
Experimentally derived enthalpy of activation
12 kcal/mol 1, from temperature-dependent
measurements of the overall rate. Experimentally
derived free enthalpy of activation 14 - 15
kcal/mol 2, from measured individual and
overall rate constants. Estimate for the
zero-point vibrational and thermal enthalpic
corrections to barrier from AM1 gas-phase
calculations of 102-atom QM region -1.3
kcal/mol (at 300 K) Resulting LCCSD(T0)-based
prediction of activation enthalpy 13.3 kcal/mol
(at 300 K) Average entropic contribution to
barrier from QM/MM TI runs 0.4 kcal/mol (at
300 K) Best prediction of free energy
barrier 13.7 kcal/mol (at 300 K)
1 W. J. H. van Berkel, F. Müller, Eur. J.
Biochem. 179, 307 (1989). 2 B. Entsch, B. A.
Palfey, D. P. Ballou, V. Massey, J. Biol. Chem.
266, 17341 (1991).
31
PHBH and CM Comparison of barriers

• PHBH p-hydroxybenzoate hydroxylase,
  electrophilic aromatic substitution
• CM chorismate mutase, pericyclic Claisen
  rearrangement
• Computed QM/MM activation enthalpies (kcal/mol)a
• Method HF B3LYP LMP2 LCCSD LCCSD(T0)
  Experiment
• CM 28.3 10.2 9.5 18.7 13.1 12.7
• PHBHb 36.7 8.4 10.7 20.2 13.3
  12.0
• Average of 16 (CM) or 10 (PHBH) single-point
  calculations at B3LYP/MM optimized geometries,
  zero-point energy and 300 K thermal corrections
  from QM calculations on cluster models,
  aug-cc-VTZ basis on oxygen and cc-pVTZ basis on
  all other atoms, MMCHARMM for CM and MMGROMOS
  for PHBH.
• Average AM1/GROMOS values for PHBH 22.8 kcal/mol
• Accurate electronic structure methods and
  transition state theory describe enzymatic
• reactions quantitatively.

F. Claeyssens, J. N. Harvey, F. R. Manby, R. A.
Mata, A. J. Mulholland, K. E. Ranaghan, M.
Schütz, S. Thiel, W. Thiel, and H.-J. Werner,
Angew. Chem. 118, 7010 (2006).
32
How to introduce classical explicit polarisation
33
COS Model Overview

• Charge-on-spring (COS) model
• Virtual site with qv attached to polarizable
  center adapts position to electric field Ei ?
  induced dipole ?ind,i
• Positions of charges-on-spring and electric field
  components depend on each other ? iterative
  scheme employed
• 2-3 iterations per step ? MD 3-4 times more
  expensive

34
Solvent effects on an SN2 reaction
D. P. Geerke, S. Thiel, W. Thiel and W. F. van
Gunsteren, JCTC 3, 1499 (2007).
35
Comparison of (free) energy profiles
PMF in DMEpol
PMF in DMEnonpol
PMF in vacuo
PES in vacuo
36
Acknowledgement

• Richard Catlow
• Shimrit Cohen
• Karl-Erich Jaeger
• Christian Lennartz
• Frank Neese
• David OHagan
• Manfred Reetz
• Ansgar Schäfer
• Sason Shaik
• Paul Sherwood
• Wilfred van Gunsteren
• Hans-Joachim Werner

• Ahmet Altun
• Iris Antes
• Dirk Bakowies
• Salomon Billeter
• Marco Bocola
• Johannes Kästner
• Hai Lin
• Nikolaj Otte
• Jan Schöneboom
• Hans Martin Senn
• Frank Terstegen
• Stephan Thiel
• Alexander Turner
• Tell Tuttle
• Dongqi Wang
• Jingjing Zheng

Support from Schweizerischer Nationalfonds Europe
an Commission (ESPRIT/QUASI) German-Israeli
Foundation for Scientific Research Volkswagenstift
ung Deutsche Forschungsgemeinschaft (SFB 663)
37
(No Transcript)
38
PHBH QM/MM approach
39
PHBH Comparison of QM/MM and full QM results
Geometries of reactant, transition state, and
product taken from optimization of the complete
system (without water solvent) at the AM1/GROMOS
level. Charges (e) taken from single-point AM1
calculations(A) QM region (102 atoms)(B)
Complete system (7004 atoms)
a) Reactant and TS p-OHB, product 3,4-DOHBb)
Reactant and TS FADHOOH, product FADHO
Full QM calculations with our linear scaling
implementation of the conjugate gradient density
matrix search AM1 barrier of 15.1 kcal/mol (B)
compared with an AM1/GROMOS value of 21.3
kcal/mol.
40
PHBH References to QM/MM studies
1 L. Ridder, A. J. Mulholland, J. Vervoort and
I. M. C. M. Rietjens, J. Am. Chem. Soc. 120,
7641-7642 (1998). 2 L. Ridder, A. J.
Mulholland, I. M. C. M. Rietjens and J. Vervoort,
J. Mol. Graphics Modell. 17, 163-175
(1999). 3 L. Ridder, B. A. PalfeyI, M. C. M.
Rietjens, J. Vervoort and A. J. Mulholland, FEBS
Lett. 478, 197-201 (2000). 4 L. Ridder, J. N.
Harvey, I. M. C. M. Rietjens, J. Vervoort and A.
J. Mulholland, J. Phys. Chem. B 107, 2118-2126
(2003). 5 S. R. Billeter, C. F. W. Hanser, T.
Z. Mordasini, M. Scholten, W. Thiel and W. F. van
Gunsteren, Phys. Chem. Chem. Phys. 3, 688-695
(2001). 6 H. M. Senn, S. Thiel and W. Thiel, J.
Chem. Theory Comp. 1, 494-505 (2005).
41
PHBH Optimized B3LYP/GROMOS structures (TZVP
basis)
a) Snapshots labeled TI taken from MD of
thermodynamic integration.
42
Free energy changes from simulations overview

• Fixed constraint
• Thermodynamic integration sampling the mean
  force
• Continuously changing constraint
• Slow growth sampling the mean force
• Fast growth fast changing constraint,
  exponential average of the energy change
• Restraint (bias)
• Umbrella sampling sampling the distribution of
  the reaction coordinate
• Free-energy perturbation instantaneous changes,
  exponential average of the energy change
• And other methods

43
Thermodynamic integration

• The reaction is split into windows with different
  ?
• The force on the constrained ? is sampled
• Mean force force of constraint Fc
• Numerical integration along ? yields ?A
• Metric tensor correction accounts for constraint
  on the momentum canonically conjugated to ?

? constrained to this value
potential energy surface
J. G. Kirkwood, J. Chem. Phys. 3, 300 (1935).
44
Error analysis strategy

• Determine error bars for the mean and the
  variance (sib)2 in each window.
• Apply error propagation in each step of
  umbrella integration (data collection,
  combination of windows, integration) to calculate
  the sampling error.
• Use the insight gained to choose the
  parameters of umbrella simulations in an
  optimum manner.
• Test approximate expressions for the
  statistical error against exact results
  available for an analytical example potential.

45
Tests for equilibration

• MD trajectories are correlated
• De-correlation through coarse-graining
• Tests for
• Lack of trend in the mean
• Lack of trend in the variance
• Normality
• Lack of correlation
• Tests provide well-defined error bars for the
  mean and the variance

S. K. Schiferl and D. C. Wallace, J. Chem. Phys.
83, 5203 (1985).
46
Error analysis data collection

• Statistical tests provide variances for
  and (sib)2 in each window
• Error propagation leads to

• Analytic test potential Error bar of the
  mean force for the window centered at .
  Black Exact curve (sampled) Red Curves
  obtained from the given formula (10
  simulations).

47
Error analysis combining the windows

• Variations in the weights pi can be neglected.

Analytic test potentialVariance var (?A/??)
over the whole range of ?. Black Exact
curve.Red Curves obtained from the given
formula.
48
Error analysis integration

• var(?A/??) is defined on the bin-grid (width h)
• Integration from ?a to ?b according to Simpsons
  rule
• Taking into account the correlation between the
  bins
• Bins are correlated if influenced by the same
  window. An approximation of the covariance leads
  to
• sb average of sib (width of window) over the
  integration range

49
Choice of umbrella potential

• bias
• K 2? recommended (? is the maximum curvature of
  A(?))
• Global histogram enough, but not too much overlap

50
Number and range of windows
Analytic test exampleError in the free-energy
barrier decreases with increasing number of
windows.

• Overlap between the windows not required in UI,
  but advantageous to reduce the sampling error.
• Distance between the window centers should be
  .
• Stronger bias (larger K) requires more windows.

51
FEP formalism

• States A and B are part of the reaction profile
• Perturbation
• Energy of state A (EA) is calculated
• QM-atoms are perturbed moved to their places in
  state B, MM-atoms remain
• Perturbed energy EB is calculated
• Sampling
• QM part frozen
• ?EEB-EA is sampled at state A

52
Hysteresis in the optimization PHBH

• The change between two hydrogen-bond patterns
  leads to a hysteresis in the energies of the
  optimized structures.
• Significant structural changes challenge FEP
  sampling intermediate states required!

53
FEP energy contributions (kJ/mol) in PHBH
J. Kästner, H.-M. Senn, S. Thiel, N. Otte, W.
Thiel, J. Chem. Theory Comput. 2, 452 (2006)
54
Sampled free-energy changes in PHBH
55
Double iterative scheme for QM/MM-pol