[PPT] - QM/MM methods (Slides from Edina Rosta, Todd Martinez, Marc PowerPoint Presentation

SLIDE 1

QM/MM methods

(Slides from Edina Rosta, Todd Martinez, Marc Vanderkamp)

SLIDE 2

Enzymes lower the activation energy for reaction

Energy Reactants Products Reaction without enzyme Reaction with enzyme: lower energy barrier

SLIDE 3

Enzyme catalytic cycles

Chemical conversion

Michaelis complex is formed (induced fit)

Difficult to dissect experimentally: use computer simulation
To simulate chemical change, we need to go beyond molecular

mechanics

SLIDE 4

Enzyme Catalysis

Rate acceleration of 1010 is

typical e.g. protease hydrolysis of peptide bonds

kcat/kuncat can reach 1016

e.g. uroporphyrinogen III decarboxylase

Lewis & Wolfenden PNAS 105, 17328 (2008)

SLIDE 5

Enzyme engineering

Finding more efficient enzymes for producing biofuels.

SLIDE 6

Enzyme engineering

Kendall Houk, UCLA David Baker, University of Washington

Nature vol 453, pages190–195 (08 May 2008)

SLIDE 7

Enzyme engineering

Kendall Houk, UCLA David Baker, University of Washington

Science 07 Mar 2008:

Vol. 319, Issue 5868, pp. 1387-1391

DOI: 10.1126/science.1152692

SLIDE 8

Enzyme engineering

Kendall Houk, UCLA David Baker, University of Washington

Science 16 Jul 2010:

Vol. 329, Issue 5989, pp. 309-313

DOI: 10.1126/science.1190239

SLIDE 9

Nobel Prize in Chemistry 1998

SLIDE 10

The recipients

Walter Kohn

Developed the Density Functional Theory (DFT), which became the most widely used quantum chemistry method due to its efficiency and accuracy.

John Pople

Developed numerous algorithms for quantum chemistry methods. Main founder behind the currently most widely used quantum chemistry software, Gaussian.

SLIDE 11

Their work enables us to

Obtain the solution of the

Schrödinger equation using approximate methods

Find energies and wave

functions of small to medium sized molecules

Provide accurate models of

chemical structure and reactivity

H E  = 

SLIDE 12

Nobel Prize in Chemistry 2013

SLIDE 13

Their work enables us to

Find the structures of

complex biomolecules by calculating their Newtonian dynamics

Find reaction mechanisms of

enzymes

Model and predict structure

and function of biological systems

SLIDE 14

A Hybrid QM/MM Approach

The development of hybrid QM/MM approaches is guided by the general idea that large chemical systems may be partitioned into an electronically important region which requires a quantum chemical treatment and a remainder which only acts in a perturbative fashion and thus admits a classical description.

SLIDE 15

Quantum Mechanics

In theory, a very accurate treatment
f the system
Largely ab initio, i.e. parameter-free
Very expensive — typically scales as

O(N4) or worse

Limited to very small systems at high

accuracy (e.g. DFT)

Can be used for larger systems at

lower accuracy (e.g. semi-empirical)

Entire proteins cannot be simulated

without enormous supercomputer power

SLIDE 16

Software Packages

QM software packages:

– ORCA – Q-Chem – Gaussian – Turbomole – MOPAC

MM software packages with

QM/MM:

– NAMD (NEW!) – CHARMM – ChemShell – AMBER – GROMACS – Q Site - Schrodinger

SLIDE 17

QM/MM implementations

SLIDE 18

The Simplest Hybrid QM/MM Model

Hamiltonian for the molecular system in the Born-Oppenheimer approximation:

              

es Ch External

f

Effect nuclei i es ch k ik k i electrons i es ch k ik k nuclei i j ij j i nuclei i nuclei j electrons i j ij electrons i ij j electrons i electrons i nuclei i j ij j i nuclei i nuclei j electrons i j ij electrons i ij j electrons i electrons i

R Q Z R Q R Z Z r R Z H R Z Z r R Z H

arg arg arg 2 2

1 2 1 1 2 1

                 

+ − + + −  − = + + −  − =

   

The main drawbacks of this simple QM/MM model are:  it is impossible to optimize the position of the QM part relative to the external charges because QM nuclei will collapse on the negatively charged external charges.  some MM atoms possess no charge and so would be invisible to the QM atoms  the van der Waals terms on the MM atoms often provide the only difference in the interactions of one atom type versus another, i.e. chloride and bromide ions both have unit negative charge and only differ in their van der Waals terms. “Standard” QM hamiltonian

SLIDE 19

A Hybrid QM/MM Model

So, it is quite reasonable to attribute the van der Waals parameters (as it is in the MM method) to every QM atom and the Hamiltonian describing the interaction between the QM and MM atoms can have a form: The van der Waals term models also electronic repulsion and dispersion interactions, which do not exist between QM and MM atoms because MM atoms possess no explicit electrons.

A. Warshel, M. Levitt // Theoretical Studies of Enzymic Reactions: Dielectric,

Electrostatic and steric stabilization of the carbonium ion in the reaction of

lysozyme. // J.Mol.Biol. 103(1976), 227-49

     

          − + + − =

nuclei i atoms MM j ij ij ij ij nuclei i atoms MM j ij j i electrons i atoms MM j ij j MM QM

R B R A R Q Z r Q H

6 12 /

ˆ

SLIDE 20

QM/MM implementations

SLIDE 21

Simple QM/QM implementation

SLIDE 22

QM/MM implementations

SLIDE 23

QM/MM implementations

SLIDE 24

QM/MM implementations

SLIDE 25

QM/MM implementations

SLIDE 26

The Hybrid QM/MM Model

Now we can construct a “real” hybrid QM/MM Hamiltonian: A “standard” MM force field can be used to determine the MM energy. For example, CHARMM force field has a form:

/ 12 6 /

ˆ ˆ

MM atoms MM atoms MM atoms electrons nuclei nuclei j i j ij ij QM MM i j i j i j ij ij ij ij MM atoms MM atoms electrons nuclei j i j el QM MM i j i j ij ij

Q Z Q A B H r R R R Q Z Q H r R     = − + + −       = − +

         

MM MM QM QM

H H H H ˆ ˆ ˆ ˆ

/

+ + =

/ /

ˆ ˆel

vdw QM QM MM QM MM MM

E H H E E =  +  + +

( )

 

      = − + + − + − + + +      

   

2 2 12 6

( ) ( ) 1 cos( ) 2 2

ij ij i j b MM nonbonded bonds angles dihedrals ij ij ij

A B q q K K E R R V n R R R

SLIDE 27

The Hybrid QM/MM Model

A “standard” MM force field can be used to determine the MM energy. For example, CHARMM force field has a form:

( )

 

      = − + + − + − + + +      

   

2 2 12 6

( ) ( ) 1 cos( ) 2 2

ij ij i j b MM nonbonded bonds angles dihedrals ij ij ij

A B q q K K E R R V n R R R

SLIDE 28

What does the QM program provide?

We need the QM/MM gradients corresponding to the total energy. To evaluate the gradient (and energy) of the QM/MM potential at given nuclei positions, the following practical information is required from the QM calculations:

Gradient of the wave function in the position of the QM atoms
Derivative of the electrostatic potential arising from the QM wave function at the

position of the MM atoms (Field)

(Energy of polarized QM system)

/ /

ˆ ˆel

vdw QM QM MM QM MM MM

E H H E E =  +  + +

Classical gradient Forces on QM nuclei (gradient of the energy) Gradient of the potential (field) at the position of the MM atoms (x qi) For derivatives On MM atoms:

SLIDE 29

Dividing Covalent Bonds across the QM and MM Regions

In many simulations it is necessary to have the QM/MM boundary cut covalent bonds, and a number of additional approximations have to be made.

SLIDE 30

QM/MM implementations

SLIDE 31

QM/MM implementations

SLIDE 32

Implementation of “link” atom approach

The link atom is placed along the bond vector joining the QM and MM atom The default link atom type is hydrogen It interacts with MM region only electrostatically (no VDW term). VdW interaction between QM and MM atoms which form 1-2 and 1-3 “bonded” pairs is not calculated. Bond stretching, angle bending, and torsion interactions between QM and MM regions are calculated as those in MM if 1-2, 1-2-3, or 1-2-3-4 terms contain at least one MM atom

SLIDE 33

Hints for running QM/MM calculations Choosing the QM region

SLIDE 34

QM/MM implementations

SLIDE 35

QM/MM implementations

SLIDE 36

QM/MM implementations

SLIDE 37

Classical alternative: EVB

Hrp is the approximately the

same in gas phase and in solution

Calculated free energies: εr

εp with FDFT and Eg with DFT

) )( (

g p g r rp

E E H − − =  

Hong, Rosta, Warshel; J. Chem. Phys., 2006

3 3

Cl CH Cl ClCH Cl

− −

+ → +

SLIDE 38

QM/MM implementation in NAMD

Electrostatic &

Mechanical embedding

Possible switching

function to cut off long range electrostatics

SLIDE 39

QM/MM implementation in NAMD

NAMD uses partial charges for the

QM atoms to calculate the electrostatic interaction with classical point charges.

There are two possibilities for the
rigin of the QM partial charges:

– the original partial charges found in the force field parameter files can be used – updated partial charges can be gathered at each step form the QM software output (controllable through the "qmChargeMode" keyword). The continuous update in charge distribution allows for a partial re- parameterization of the targeted molecule at each time step.

SLIDE 40

QM/MM implementation in NAMD

Link atom treatment
keywords

"qmBondDist" to scale the bond length

"qmLinkElement“

to change the terminating atom type

SLIDE 41

QM/MM implementation in NAMD

Several charge

redistribution scheme to remove partial charges from nearby QM atoms

"qmBondScheme"

keyword

No charge or

Redistributed Charge and Dipole (RCD) methods

SLIDE 42

QM/MM implementation in NAMD

ORCA and MOPAC

natively supported

python wrapper scripts

for GAUSSIAN, TeraChem and Q-Chem

SLIDE 43

QM/MM implementation in NAMD

Multiple QM regions

are supported

More on keywords,

etc: http://www.ks.uiuc.edu /Research/qmmm/

SLIDE 44

QM/MM reaction modelling

To overcome activation energy, need to ‘force’ a reaction

happening: apply ‘bias’ along a ‘reaction coordinate’

Potential Energy Surface QM/MM minimization Free Energy Surface QM/MM MD + statistics Use enhanced sampling e.g. ‘umbrella sampling’ ‘adiabatic mapping’

SLIDE 45

Dynamics on QM/MM potential

Minimizations corresponding to 0K is often problematic for larger systems:

Missing entropy effects are important in e.g., protein folding, ligand binding, etc.
Even is enthalpic contributions dominate, only local minima are found

There are several methods to model the thermostatted movements of the atoms

n a given potential. These should provide a proper thermodynamic ensemble at

long time scales. A popular one is Langevin dynamics to obtain canonical ensemble: R(t) is a delta-correlated stationary Gaussian process with zero-mean and a constant variance. g: damping constant Others: Andersen, Nose-Hoover, etc…

SLIDE 46

1. Umbrella Sampling

Run parallel

simulations with harmonic constraints moving along the reaction coordinate

Recover the unbiased

free energy surface from combined data using e.g., WHAM

( )

( ) ( )

2

1 2

pot A i A i

i A

U q k

E q

  = + −

Ferrenberg, Swendsen; Phys. Rev. Lett. 1989 - Cited: ~ 2800 times Kumar, Rosenberg, Bouzida, Swendsen, Kollman; J. Comput. Chem. 1992 (WHAM) - Cited: ~ 3600 times Torrie and Valleau; J. Comp. Phys. 1977 Cited: ~ 3200 times

SLIDE 47

1. WHAM

( )

( ) ( ) 1

Pr

=



k i

Nbin n k k i i

p

( )

( ) 1 1

ln

k i

NSim Nbin n k i k i

L p

= =

= 

( ) ( ) ( )

( ) ( ) 1

k i k i bin k j

U U k k i i i N U j j

e p p f e p e p

   − − − =

= =

 ( )

( ) ( ) ( )

( ) ( ) ( ) 1 1 1 1

ln 1

k i k k i i

NSim Nbin NSim Nbin n U U k k k i i k i k i

L f e p f e p

 



− − = = = =

  = + −    

  

Probability of observing a trajectory in the k-th simulation.

( ) k i

p

Equilibrium probability for bin i

( ) k i

n

Histogram count in bin i

( ) k i B

U k T

e

−

Boltzmann factor for biasing energy

( )

( ) 1

1

bin k j

k N U j j

f e p

 − =

=



Normalizing factor for the equilibrium probability

SLIDE 48

1. WHAM

( )

( ) ( ) 1

Pr

=



k i

Nbin n k k i i

p

( )

( ) 1 1

ln

k i

NSim Nbin n k i k i

L p

= =

= 

( ) ( ) ( )

( ) ( ) 1

k i k i bin k j

U U k k i i i N U j j

e p p f e p e p

   − − − =

= =

 ( )

( ) ( ) ( )

( ) ( ) ( ) 1 1 1 1

ln 1

k i k k i i

NSim Nbin NSim Nbin n U U k k k i i k i k i

L f e p f e p

 



− − = = = =

  = + −    

  

Probability of observing a trajectory in the k-th simulation.

( ) k i

p

Equilibrium probability for bin i

( ) k i

n

Histogram count in bin i

( ) k

L f  = 

i

L p  = 

(l)

( ) 1 (l) (l) 1

i

Msim k i k i Msim U l

n p N f e 

= − =

=

 

SLIDE 49

1. Free Energy Simulations Using the

String Method

R1 R2 RS TS PS

String-type methods have become

highly successful in recent years

Optimized a 1D string in the

multidimensional space of the internal reaction coordinates to

btain minimum free energy path
Hamiltonian replica exchange

between string images

E, Ren, Vanden-Eijnden, Phys. Rev. B, 2002

Henkelman, Jonsson, W. Yang, Brooks, … Rosta, Nowotny, Yang, Hummer, J. Am. Chem. Soc., 2011

SLIDE 50

1. Free Energy Simulations Using the

String Method

Start with a guess for the string
Run Umbrella Sampling

simulations

Determine forces for the images

along the string

Fit new string

R1 R2 RS TS PS

SLIDE 51

Start with a guess for the string
Run Umbrella Sampling

simulations

Determine forces acting on the

images along the string

Fit new string
Redistribute images
Run next iteration

R1 R2 RS TS PS

1. Free Energy Simulations Using the

String Method

SLIDE 52

Converged string:

– Forces are parallel to string

We use all data from all

string simulations with

Histogram Free implementation of WHAM: works with very high dimensionality

R1 R2 RS TS PS R1 R2 RS TS PS

( ) ( )

( ) ( ) (l) (l) (l) 1 1

1  

=  =

=  



k k NSim A NSim k k k A NSim A l

c f N f c

1. Free Energy Simulations Using the

String Method

Rosta, Nowotny, Yang, Hummer, J. Am. Chem. Soc., 2011

SLIDE 53

1. Free Energy Simulations Using the

String Method

SLIDE 54

Hamiltonian Replica Exchange

K. Hukushima and K. Nemotto, J. Phys. Soc. Japan, 1996

Fukunishi, Watanabe & Takada, J. Chem. Phys., 2002

Reaction Coordinate ξi

Replica 1 Replica 2



1



Running MD at different

temperatures in parallel

Couple the runs in order to

speed up lowest temperature’s dynamics

Preserve Peq at each

temperature

Detailed balance condition

has to be satisfied

Temperature

RS TS PS

SLIDE 55

1. RNase H: Free Energy Surface

Qe=r1-r2 Qp=r3-r4+r5-r6+r7-r8 Rosta, Nowotny, Yang, Hummer, J. Am. Chem. Soc., 2011

SLIDE 56

HDV ribozyme

Reaction mechanism

Ganguly, Thaplyal, Rosta, Bevilacqua, Hammes-Schiffer, JACS 2014

SLIDE 57

A: concerted

mechanism with phosphorane intermediate

B: sequential

pathway with a proton- transferred intermediate

C: sequential

pathway with a phosphorane intermediate

Ganguly, Thaplyal, Rosta, Bevilacqua, Hammes-Schiffer, JACS 2014

SLIDE 58

HDV ribozyme

Reaction mechanism

A B C PT

Ganguly, Thaplyal, Rosta, Bevilacqua, Hammes-Schiffer, JACS 2014

SLIDE 59

1. Umbrella Sampling Simulations & WHAM

Rosta, Hummer, JCTC, 2015

MC simulations
7 Umbrellas with 50 kcal/mol biasing force each

SLIDE 60

1. Systematic error when using WHAM in

conjunction with small biasing force

Rosta, Hummer, JCTC, 2015

6 Umbrellas with 50 kcal/mol biasing force each
1 Umbrella with 1 kcal/mol biasing force each

SLIDE 61

Master equation:

1. Markov-state models of conformational

transitions

1 1 ( ) ( )

( ) ( ) ( )

N N i j i j j j i j i j j i i

dP t P t k P t dt k

= =    

= −

 

( ) ( )

U F

k T k T

F U ⎯⎯⎯ → ⎯⎯ ⎯

Folding@home, V. Pande, Stanford

Voelz, Bowman, Beauchamp and Pande, J. Am. Chem. Soc., 2010

State U State F

SLIDE 62

DHAM: Dynamic Histogram Analysis Method

Our approach:

Discretized reaction coordinates
Short lagtimes

Most molecular dynamics integrators are Markovian!

SLIDE 63

1. DHAM: Dynamic Histogram Analysis Method

( )

( ) ( ) 1 1

Pr

k ji

Nbin Nbin T k k ji i j

M

= =



L = ln M ji

(k)

( )

Tji

(k )

j=1 Nbin

Õ

i=1 Nbin

Õ

k=1 NSim

Õ

( ) ( ) ( ) ( ) ( ) 1 =

= =



bin

k ji ji k k k ji i ji ji N k li li l

c M M f c M c M

Rosta, Hummer, JCTC, 2015

Biased Dynamical Trajectories Histogram

f Transitions

Markov Model Free Energy & Kinetics

( ) 1 ( ) ( ) ( ) 1 Msim k ji k ji Msim k k k i i ji k

T M n f c

= =

=

 

2 2

/ ln

lag

t   = −

( )

( ) ( ) ( ) ( ) ( )

exp / 2

k ji k k k k i ji j i B ji

M f c u u k T M =  − −

SLIDE 64

Applications using Umbrella Sampling

DHAM DHAM

SLIDE 65

1. Applications for “downhill” unbiased non-

equilibrium trajectories

DHAM

MC trajectories initiated from around x=0.63

SLIDE 66

1. Applications for “downhill” unbiased non-

equilibrium trajectories

( ) 1 ( ) 1 Msim k ji k ji Msim k i k

T M n

= =

= 



1 Nbin ij j i j

M p p

=



DHAM

SLIDE 67

1. Umbrella sampling QM/MM simulations of

catalytic reactions

R. Suardiaz, P. G. Jambrina, L. Masgrau, A. Gonzalez-Lafont, E. Rosta, J.M. Lluch, JCTC, 2016

Hydrogen abstraction from arachidonic acid in the catalytic reaction of human 15-LOX-2 lipoxygenase.

SLIDE 68

1. Umbrella sampling QM/MM simulations of

catalytic reactions

R. Suardiaz, et al, JCTC, 2016
P. Saura, R. Suardiaz, et al, ACS Catalysis, 2017

k=0.22 s-1 @T=300K Calculated prefactor using a barrier of 18.7 kcal/mol: A = 8.6x1012 s-1 kBT/h = 6.3x1012 s-1 Experimental kcat = 0.7 s-1

Lag time Rate 1 fs 0.24 s-1 10 fs 0.17 s-1 50 fs 0.29 s-1