Next-generation DFT-based quantum models for simulations of biocatalysis Darrin M. York U niversity of M innesota Minneapolis, Minnesota USA http:// theory.chem.umn.edu

Outline

•

AM1/d-PhoT model for RNA catalysis

•

Efficient treatment of long-range electrostatics in semiempirical calculations

•

Improved charge-dependant response properties

•

Selected applications

• • • • •

…in words

Study phosphate reactivity comprehensively (using small models) with high-level quantum models ( ab initio and DFT) Construct accurate semiempirical quantum models capable of being used in linear-scaling electronic structure and QM/MM simulations Develop improved (accurate, fast and general) models for electrostatics, solvation and generalized solvent boundary potentials.

Investigate how to improve next-generation semiempirical quantum models to account for charge-dependent response properties without significant sacrifice of efficiency.

Validate methods with respect to known reactions in solution, then apply them to the important problem of RNA catalysis in a realistic system consisting of many thousands of particles, and simulated for many tens of nanoseconds.

Phosphates and phosphoranes

Mechanisms for phosphoryl transfer

Dissociative D N

O RO lg P O O

Concerted A N D N

O RO lg P O O

Associative A N +D N

O RO lg P O O O nuc R O nuc R R O lg O P O O RO lg  O P O O  O nuc R RO lg O P O O O nuc R

QCRNA

– Online!

http://theory.chem.umn.edu/QCRNA Molecule (2000+) Reaction Mechanism (300+)

Giese et al., J. Mol. Graph. Model. 25, 423 (2006).

QCRNA

– Online!

http://theory.chem.umn.edu/QCRNA Potential Energy Surface

Reaction Tables Graphical Interface

Giese et al., J. Mol. Graph. Model. 25, 423 (2006).

Phosphate isomerization (Migration)

movie Liu et al., J. Phys. Chem. B, .109 , 19987 (2005);

Chem. Commun

.

,

31

, 3909 (2005). Silva-Lopez

et al., Chem. Eur. J.

, 11 , 2081 (2005); Mayaan

et al., J. Biol. Inorg. Chem

.

,

9

, 807 (2004). Range

et al., J. Am. Chem. Soc.

, 126 , 1654 (2004).

Parameter Optimization: AM1/d Methods

  2 (

λ

)    Mol Prop 

i w i

  

w iα

( 

i

 

Y i Semi

 )  2 (

λ

) 2 (

λ

) 

γ

T

 (

C



λ

 

Y i DFT



b

)  

0

2  Training set included a wide variety of biological phosphates and phosphoranes, hydrogen bonded complexes, proton affinities and reaction paths of associative and dissociative mechanisms in different charge states.

Nam

et al., J. Chem. Theory Comput., submitted

.

Why use a semiempirical model?

It is important to note that for the ribozyme systems of interest, the details of the mechanisms remain topics of considerable debate. Hence the goal is to test multiple mechanisms with a model that is sufficiently predictive to discern the most probable path.

A consensus has emerged that, in certain ribozymes such as HHR and HDV, a large scale conformational change either precedes or is concomitant with the chemical step of the reaction.

This necessitates the use of a quantum model that is able to be used with extensive conformational sampling (i.e., simulation) while providing an accurate description, in terms of energy, structure and charge distribution, along multiple mechanistic paths (i.e., not a single pre determined 1-D reaction coordinate) in order to be predictive.

Modification for AM1/d-PhoT Model

Want a

d

-orbital method for hypervalent species, but one that also describes reasonably hydrogen bonding interactions. Combine MNDO/d framework with a modified core-core term similar to AM1 (and retaining some AM1 parameters unmodified) to build a semiempirical model for phosphoryl transfer reactions: AM1/d-PhoT

Core-Core Repulsion

E MNDO AB



Z A Z B s A s A s B s B



1



e

 

A R AB



e

 

B R AB



E AB



E MNDO AB



Z A Z B R AB i 4

 

1



a i A e



b i A ( R AB



c i A ) 2



a i B e



b i B ( R AB



c i B ) 2



MNDO AM1 and PM3

Modified Core-Core Repulsion

E AB



E MNDO AB



Z A Z B R AB G A G B i 4

 

1



a i A e



b i A ( R AB



c i A ) 2



a i B e



b i B ( R AB



c i B ) 2



If G

A



and G

B

= 0, MNDO Hamiltonian If G

A



and G

B

= 1, AM1 and PM3

AM1/d-PhoT Model for Phosphoryl Transfer

AM1/d-PhoT Model for Phosphoryl Transfer

AM1/d-PhoT Model for Phosphoryl Transfer

AM1/d-PhoT Model for Phosphoryl Transfer

AM1/d-PhoT Model for Phosphoryl Transfer

Reaction Energies and Barrier Heights Error*

Neutral Rxn AM1/d AM1 PM3

Reaction Energy

No. Rxn MSE MUE

5 2.07

2.86

-7.32 -10.78

7.39 10.78

Monoanionic Rxn AM1/d AM1 PM3

4 0.84

1.96

-2.48

9.79

-4.94

8.80

Dianionic Rxn AM1/d AM1 PM3

2 -1.44

2.28

-9.00

9.00

-2.96

5.65

Dissociative Rxn AM1/d AM1 PM3

3 5.25

5.25

-23.24 -12.35

23.24 12.35

Activation Energy

No. TS MSE MUE

13 0.76

3.61

3.48 -18.76

6.62 18.76

11 -2.91

3.57

-0.36 -12.74

12.23 16.23

4 -3.33

3.33

-22.58 -31.77

22.58 31.77

3 3.35

6.60

10.08 -10.38

10.08 10.38

Relative Intermediate Energy

No. Int MSE MUE

8 -1.06

2.36

-42.29 -26.61

42.29 26.61

7 -6.59

6.59

-42.34 -34.10

42.34 34.10

*Errors are computed against

“B3LYP/6-311++G(3df,2p) adiabatic energies”

Linear Free Energy Relations Transphosphorylation of a cyclic phosphate with enhanced leaving groups.

Slope of plot is the Brønsted correlation parameter β reactions.

lg often used to characterize phosphoryl transfer The logk values were calculated from DFT and are contained in

QCRNA.

Gas Phase Proton Affinity I Molecule

H 3 O + HOH CH 3 OH CH 3 CH 2 OH C 6 H 5 OH CH 3 CO 2 H P(O)(OH)(OH)(OH) P(O)(O)(OH) P(O)(O)(OH)(OH) P(O)(O)(O)(OH) 2 P(O)(OH)(OH)(OCH 3 ) P(O)(O)(OH)(OCH 3 ) P(O)(OH)(OCH 3 )(OCH 3 ) P(O)(OH)(OCH 2 CH 2 O)

MSE MUE Ref.

165.0

390.3

381.5

378.2

350.1

347.2

330.5

310.6

458.9

581.1

329.3

454.9

329.4

329.5

B3LYP

-1.1

0.1

-2.2

-2.2

-2.4

-0.8

-2.4

-0.1

-1.1

-1.7

0.4

-1.4

0.7

-0.1

-1.0

1.1

AM1/d

3.8

5.4

2.0

2.9

-3.4

-2.7

-3.4

1.5

-1.9

10.4

0.3

0.7

1.8

-0.2

0.9

2.4

Error AM1

-2.0

20.5

2.7

4.7

-3.1

6.1

7.6

20.6

16.8

33.7

7.2

16.5

7.3

7.6

9.4

9.8

Range et al., Phys. Chem. Chem. Phys. 7,

B3LYP :

3070 (2005).

B3LYP/6-311++G(3df,2p)//B3LYP/6-31++G(d,p) PM3

-11.8

11.3

-1.9

-0.4

-6.9

0.9

15.0

35.1

24.7

36.4

14.9

22.8

12.3

11.8

8.5

11.0

MNDO/d

5.6

30.6

1.8

5.2

0.0

9.6

-12.2

-3.6

-2.8

16.3

-12.0

-7.6

-14.1

-17.1

-5.1

11.4

Gas Phase Proton Affinity II: Phosphorane Compounds

Error Molecule Ref.

P(OH)(OH)(OH)(OH)(OH) P(OH)(OH)(OH)(OH)(OH) P(OH)(OH)(OCH 2 CH 2 O)(OH) P(OH)(OH)(OCH 2 CH 2 O)(OH) P(OH)(OCH 2 )(OCH 2 CH 2 O)(OH) P(OH)(OCH 2 )(OCH 2 CH 2 O)(OH) P(OH)(OH)(OCH 2 CH 2 O)(OCH 2 )

MSE MUE

351.0

341.0

351.9

343.2

345.2

352.0

343.5

B3LYP

-0.4

-1.8

-0.9

-1.1

-0.7

-0.8

-1.1

-1.0

1.0

AM1/d

3.0

1.8

1.2

-2.5

-3.5

2.3

-0.7

0.2

2.1

AM1

9.0

13.6

5.9

8.0

3.6

5.4

6.2

7.4

7.4

PM3

8.3

9.0

1.7

-0.5

-2.3

-0.4

-0.9

2.1

3.3

MNDO/d

-1.3

-8.7

-11.8

-17.4

-20.2

-27.0

-19.5

-15.2

15.2

Range et al., Phys. Chem. Chem. Phys. 7,

B3LYP :

3070 (2005).

B3LYP/6-311++G(3df,2p)//B3LYP/6-31++G(d,p)

Example: QM/MM of Di-anionic Reactions in Water

-6 -5 -4

Comparison with DFT and Expt.

in kcal/mol

-3 35 30 25 20 15 10 5 -2 -1 0 -5 0 1 2 -10

q = R(P-O l ) - R(O n -P)

EP(-)….OH(-) DMP(-)…OH(-) TMP(-)...OH(-) 3 4 5 6

TS 1

*DFT:

B3LYP/6-311++G(3df,2p)

Dejaegere and Karplus,

JACS 1993

Cox and Ramsay,

Chem. Rev. 1964

Problems

•

Dispersion interactions

•

Relative conformational energies: sugar puckering and pseudorotation transition states

•

Proper treatment of polarizability and multiple charge states

The Problem of Charge-dependent Response Properties with Semiempirical Methods

Atoms are of course an

extreme

case: but typically polarizabilities of

neutral molecules are typically off by 25%, and anions by significantly more…

Giese et al., J. Chem. Phys., 123 , 164108 (2005).

Goal:

Improve charge-dependent response properties of semiempirical methods without significantly increasing computational cost.

Possible solutions: •

Reparameterize models

•

Increase minimal basis-set representation

•

Make basis set exponents charge dependent

DFT-based model…



E

[  ] 

E

[   ]

F

[    ]       (

r

)

v

(

r

)

d

(

r

)

d

3

r



N

3

r

   0 Giese et al., J. Chem. Phys. 123 , 164108 (2005).

E

[  ]    

E

[ 

ref

]  

E

( 1 , 

ref

2 ,  ) [  ]  

ref



E

( 1 , 

ref

2 ,  ) [  ] 

E

[  ] 

E

[ 

ref

]     

E

[   (

r

) ]   

ref

 (

r

)

d

3

r

 1 2   (

r

)    2

E

[  ]  (

r

)  (

r

 )   

ref

 (

r

 )

d

3

rd

3

r

  

E

[  ] 

E

[  ]  

E

( 1 , 

ref

2 ) [  ]    (

r

)  

k c k



k

(

r

) 

E

( 1 , 

ref

2 ) [  ] 

c

T



m

 1 2

c

T



η



c

m i

     

E

[   (

r

) ]    

ref



i

(

r

)

d

3

r



ij

  

i

(

r

)      (

r

2 )

E

[   ] (

r

 )    

ref



j

(

r

 )

d

3

rd

3

r



m i



ij

   

i D ij

(

r

)   

v

(

r

) (

r ij

)   

j D ij

(

r ij

)  

i

| (

r r

)  

r



j

| (

r

 )

d j

(

r

) 3

rd

  

d

3

r

 3

r D

(

r ij

) 

D

(

r ij

;

C i

,

W i

) 

D

(

r ij

;

C j

,

W j

) 

i

(

r

)   2 

i

(

q i

) 2   

i

( 

q i

) 2   (

u



U i

)

e

 

i

(

q i

) 2 |

r



R

i

| 2 

i

(

q i

)    

i

3 (

q i

)   1 / 3 2  

i

(

q i

)   

i

( 0 ) 

e

 3

B i q i

A Variational Electrostatic Projection (VEP) Method for QM/MM Calculations

Goal:

Model large-scale electrostatic effects of solvent-shielded macromolecular environment -

it’s linear response and

– in hybrid QM/MM calculations for a fraction of computational cost of explicit simulation

Method:

Green’s function approach that involves variational projection and reduced dimensional mapping of surrounding solvent-shielded macromolecular environment onto the dynamical reaction zone Gregersen and York, J. Phys. Chem. B, 109, 536-556 (2005).

Gregersen and York, J. Comput. Chem., 27, 103 (2006).

Multi-scale Quantum Models

External potential of solute and solvent Stochastic boundary Reaction Region

QM active site + MM surrounding

(Newtonian dynamics) Buffer Region (Langevin dynamics)

Linear-scaling QM/MM-Ewald Method

Nam et al., J. Chem. Theory Comput., 1, 2 (2005).

Applications to enzymes and ribozymes

• Hammerhead ribozyme Best characterized ribozyme – but complicated: role of metals, chemical/conformational steps, non-inline native structure • Hairpin ribozyme No metal cofactor, in-line configuration

General acid/base mechanism

Tai-Sung Lee

et al., submitted.

Mg 2+ ion is observed to coordinate the O2’ of G8 increasing it’s acidity in the early TS and then migrate closer to the leaving group O5’ position of the scissile phosphate in the late TS.

Simulations help to explain the long-standing disconnect between available structures and biochemical data (in particular, thio effect studies).

Early TS Late TS

Other Projects…

•

Parameters for RNA reactive intermediates

•

DNA bending

•

Polarization-exchange coupling

•

Linear-scaling electronic structure

• • • • • • • George Giambasu Dr. Tim Giese Yun Liu Dr. Evelyn Mayaan Adam Moser Dr. Kwangho Nam Dr. Kevin Range

Acknowledgements

• • • • • Dr. Olalla Nieto Faza Dr. Francesca Guerra Dr. Carlos Silva Lopez Prof. Xabier Lopez Dr. Anguang Hu • • • • • Prof Bill Scott Prof. Qiang Cui Dhd Marcus Elstner Prof. Jiali Gao Prof. Walter Thiel •

Funding/Resources:

• University of Minnesota • NIH • ACS-PRF Army High-Performance Computing Research Center • Minnesota Supercomputing Institute