University of Southern California Dissertations and Theses

Ab initio methodologies in studying enzymatic reactions

(USC Thesis Other)

Ab initio methodologies in studying enzymatic reactions

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content i

AB INITIO METHODOLOGIES IN STUDYING ENZYMATIC REACTIONS

by
Edina Rosta

A Dissertation Presented to the
FACULTY OF THE GRADUATE SCHOOL
UNIVERSITY OF SOUTHERN CALIFORNIA
In Partial Fulfillment of the
Requirements for the Degree
DOCTOR OF PHILOSOPHY
(CHEMISTRY)

August 2007

Copyright 2007 Edina Rosta
ii

Acknowledgments

I would like to thank Prof. Arieh Warshel his guidance and help through my
PhD. I am very thankful for the opportunity that I could learn from Prof. Warshel
about biological systems, scientific applications and methods. His insight for the
value of scientific work is outstanding in the academic world and I am grateful for
any feedback or judgment he commented my research during my PhD work.
I would like to thank all (current and former) members and visitors of Prof.
Warshel’s group: Dr. Pankaz Sharma, Dr. Hanbin Liu, Prof. Sonja Braun-Sand,
Dr. Marek Strajbl, Dr. Anton Burykin, Dr. Mitsunori Kato, Dr. Mats Olsson, Dr.
Gongyi Hong, Dr. Zhen Tao Chu, Prof. Jan Florian, Prof. Avital Shurki, Prof.
Janez Mavri, Dr. Claudia Nili Schutz, Dr. Marco Klahn and Dr. Yun Xiang.
I would also like to thank Prof. Michael Arbib for all the interesting hours I
could spend with his group. I am very grateful to Prof. Leonard Adleman for
showing me how beautiful science can be.

iii
Table of Contents
Acknowledgments ii
List of Tables v
List of Figures vi
List of Abbreviations xi
Abstract xii
Chapter 1. Introduction to the Simulations of Enzymatic Reactions 15
1.1 Motivation 15
1.2 Overview of the dissertation 20
Chapter 2. Simulations of Chemical Reactions in Solutions: Implicit
Solvation Models 22
2.1 Introduction 22
2.2 Continuum Solvation 24
2.3 Langevin Dipoles Solvation Model 27
Chapter 3. All atom simulations 42
3.1 Introduction 42
3.2 QM/MM 43
3.3 Frozen DFT Approach 47
3.4 Empirical Valence Bond (EVB) method 50
Chapter 4. Free energy calculations 58
4.1 The free energy perturbation (FEP) and the linear response
approximation (LRA) methods 58
4.2 QM/MM FEP Solvation Model 63
Chapter 5. QM/MM free energies using reference potentials 68
5.1 Introduction 68
5.2 Methods and systems 70
5.3 Results and Discussion 79
5.4 Conclusions 88
iv
Chapter 6. Linear Free Energy Relationships 91
6.1 Introduction 91
6.2 Transferability of the off diagonal elements between different
phases 95
6.3 LFERs in Phosphate Hydrolysis: an Introduction 108
6.4 Free energy surfaces for the hydrolysis of phosphate-monoesters in
solution 120
6.5 The Free Energy Surface for the Reaction of Ras 137
6.6 Summary 140
Chapter 7. Future Directions 151
Bibliography 153
Appendix A: The CPHF equations 172

v
List of Tables

Table 1. Comparison of the calculated pK
a
values for the QM/MM-FEP and
the COSMO methods. The experimental values are significantly
closer to the all atom model than to the continuum results........................... 66

Table 2. The EVB energetics of the S
N
2 reaction in water and in haloalkane
dehalogenase for simulations with five different initial conditions.
a
........... 80

Table 3. The LRA estimate of the free energy of transfer from the EVB to
the QM/MM surface in the RS and PS for the reaction of DhlA.
a
................ 82

Table 4. Results of the LRA simulation for the transition state of the
solution reaction and the overall QM/MM activation barrier.
a
..................... 84

Table 5. Results of LRA simulations at the transition state of the enzyme
reaction and the overall QM/MM activation barrier.
a
................................... 86

Table 6. Gas-phase ab initio activation energies for the attack of a
carboxylate on DCE.
a
.................................................................................... 88

Table 7. The correlation between the rate constants (and activation barriers)
and the observed pK
a
values for the hydrolysis reaction of
phosphomonoesters with different leaving groups.
a
................................... 129

Table 8. The charges of the leaving group in different states of the
hydrolysis of some selected systems.
a
........................................................ 131

Table 9. Difference of zero point energies (relative to the reactant state) for
different TSs of the magnesium mono methyl pyrophosphate trianion.
a
... 133

Table 10. Calculated isotope effects for the systems considered in Table 9.
a
.... 133

Table 11. Comparing the energies obtained by different solvent models and
different ab initio methods for characteristic points on the PES of the
hydrolysis of the magnesium mono methyl pyrophosphate trianion.
a
....... 141

vi
List of Figures

Figure 1. The figure shows the solvent excluded volume in the case of the
NaF. Each atom type has a van der Waals radius that is defined based
on experimental solvation energies............................................................... 25

Figure 2. Illustration of the Langevin dipoles solvation model: the water is
represented as a coarser and a denser grid of dipoles around the
solute. ............................................................................................................ 27

Figure 3. The Langevin Function. The figure shows the proportion of the
actual dipole compared to the total dipole moment in the function of
the interaction energy. This energy depends on the external electric
field and the temperature............................................................................... 28

Figure 4. Illustration of the smoothing function used to obtain continuous
derivatives. In the figure on the left side one can see the two potential
surfaces with and without the smoothing function. The right side
figure shows the smoothing of the deletion of a dipole, thus changing
the step function to a sigmoid curve. ............................................................ 30

Figure 5. Illustration of the use of potential energy scans to find transition
state barriers. The disadvantage here can easily seen if one tries to
imagine this procedure in higher than 2 dimensions..................................... 31

Figure 6. Test case of the [F…H…F]
-
system. The green solid line shows
the calculated potential energy (note that it is smooth) and the red
dots are the transition state points that have been found, starting from
different starting points. The F atoms were kept in one place while
the H was allowed to move. The van der Waals radius of the atoms is
illustrated on the bottom of the figure........................................................... 33

Figure 7. Transition state corresponding to the S
N
2 reaction of ........................... 34

Figure 8. Transition state search for the reaction of methylchloride and
ammonia. We have simultaneously mapped the potential energy
surface on a grid of 0.01 Å resolution and used gradient search to find
the transition state. ........................................................................................ 36

vii

Figure 9. Illustration of the Umbrella Sampling technique. A good
sampling can be achieved even at the transition state region due to the
fact that some of the mapping potentials have minima, rather than
maxima at the transition state along the reaction coordinate. In
practice we first evaluate the free energy changes corresponding to
the change between the mapping potentials and then we can use the
m
G Δ we obtained to find free energy surfaces corresponding to other
potential energies such as E
g
......................................................................... 54

Figure 10. Illustration of the EVB method: the free energy of the reaction
is built from the force field defining the reactants and another force
field defining the products. The ground state energy is obtained as the
eigenvalue of an empirical 2 by 2 Hamiltonian using the diabatic
states' energies and an empirical off-diagonal element................................. 56

Figure 11. Showing the free energy functions of the reactant and product
states, ( )
a
g
X
and ( )
b
g
X
, as two Marcus' type parabolas of equal
curvature. The reorganization energies at the minima of
a
g and
b
g
are given by
b b a a b
a
U U G λ
→
= − − Δ and
a a b a b
b
U U G λ
→
= − − Δ .
By assuming
b a
λ λ = one can obtain the LRA estimate of
b a
G
→
Δ
(see (Sham, Chu et al. 2000))........................................................................ 62

Figure 12. The thermodynamic cycle used to evaluate the QM/MM
activation barrier. .......................................................................................... 70

Figure 13. A schematic description of the intersection of
EVB
E and
/ QM MM
E at the TS region. The figure can be used as a qualitative
illustration (see text) for the validity of the LRA treatment of Eq.
(5.9) and (5.12).............................................................................................. 74

Figure 14. Showing the overall structure of the DhlA with the active site
Asp124 and a bound DCE............................................................................. 76

Figure 15. A schematic description of the reaction steps in DhlA (4A) and
a stick diagram describing the
N
S 2 step for the reaction in water with
an illustration of the link atom (4B).............................................................. 77

viii
Figure 16. Showing the SCAAS simulation system for the reference
reaction in water............................................................................................ 79

Figure 17. The EVB free energy profile of the
N
S 2 reaction in water and
in the haloalkane dehalogenase enzyme. The profiles were obtained
using FEP/US for the EVB ground state energy........................................... 81

Figure 18. Illustrating the fluctuations of
/ QM MM EVB
E E − for runs on the
transition state of the EVB surface at
‡
EVB
x . The figure depicts the
results of 5 runs, each of 10 ps...................................................................... 81

Figure 19. Showing the distributions of
/ QM MM EVB
E E − for trajectories on
the EVB and QM/MM reactant state. ........................................................... 83

Figure 20. Showing the active site of the DhlA for the TS of the
N
S 2
reaction step. ................................................................................................. 83

Figure 21. Obtaining the QM/MM activation free energy by moving from
the EVB to the QM/MM surfaces. ................................................................ 85

Figure 22. Diabatic and adiabatic FDFT energy profiles for the reaction,
3 3
Cl CH Cl ClCH Cl
− −
+ → + , in the gas phase and in solution. The
reaction coordinate is defined as the energy difference between the
diabatic surfaces,
1 2
ε ε ε Δ = − . ................................................................... 105

Figure 23. Plot of H
ij
of the above reaction in gas phase and solution. The
data are obtained from the diabatic and the adiabatic curves of Figure
22, using ( ) ( )( )
12 1 2 g g
H E E ε ε ε Δ = − − . Note that the adiabatic
barrier is underestimated by about 9 kcal/mol due to the use of the
DFT-BP method(Yang, Fleurat-Lessard et al. 2004).................................. 106

Figure 24. A schematic description of the potential surface for the
hydrolysis of phosphomonoesters. The figure provides a clear
definition of the associative and dissociative paths and also defines
the three reaction coordinates R
1
and R
2
and X. ....................................... 114

ix

Figure 25. A VB description of the energy surfaces for phosphate
hydrolysis in three extreme cases, associative (A), dissociative (B)
and concerted (C). The figure generates the surfaces for the
hydrolysis reaction by mixing four states using an EVB formulation
(the effect of the off-diagonal term is in included for simplicity) .
When the surface follows one of the extreme cases it is relatively
easy to describe the corresponding LFER by shifting parabolas in the
lower part of the surface as was done in ref. (Aqvist, Kolmodin et al.
1999). RS, IS and PS represent, respectively, reactant, intermediate
and product states. The indices as and ds stand for associative and
dissociative, respectively. The dots on the surfaces designate the
minima of the corresponding diabatic surfaces........................................... 116

Figure 26. The calculated surface for the hydrolysis of the mono methyl
phosphate dianion in the space defined by R
1
and R
2
. The figure
depicts the reaction path that connects the RS with the PS and the
corresponding TS. Energies are given in kcal/mol and distances in Å....... 120

Figure 27. The calculated TS structure for the hydrolysis of the mono
methyl phosphate dianion. .......................................................................... 121

Figure 28. The calculated surface for the hydrolysis of the phenyl
phosphate dianion in the space defined by R
1
and R
2
. The figure
shows the reaction path that connects the RS with the PS and the
corresponding TS. Energies are given in kcal/mol and distances in Å....... 122

Figure 29. The calculated surface for the hydrolysis of the mono methyl
pyrophosphate trianion in the space defined by R
1
and R
2
. The
figure shows the reaction path that connects the RS with the PS and
the corresponding TS. An alternative associative path with higher
energy barrier is shown (dashed line). Energies are given in kcal/mol
and distances in Å. ...................................................................................... 123

Figure 30. The calculated TS structure for the hydrolysis of the mono
methyl pyrophosphate trianion.................................................................... 124

Figure 31. The calculated surface for the hydrolysis of the mono methyl
pyrophosphate dianion in the space defined by R
1
and R
2
. The figure
shows the reaction path that connects the RS with the PS and the
x
corresponding TS. An alternative associative path with higher energy
barrier is shown (dashed line). Energies are given in kcal/mol and
distances in Å. ............................................................................................. 125

Figure 32. The calculated surface for the hydrolysis of the magnesium
mono methyl pyrophosphate trianion in the space defined by R
1
and
R
2
. The figure shows three possible reaction paths with similar
energy barriers that connect the RS with the PS and the
corresponding TSs. Energies are given in kcal/mol and distances in
Å.................................................................................................................. 126

Figure 33. Three possible calculated TSs structures for the hydrolysis of
the magnesium mono methyl pyrophosphate trianion. ............................... 127

Figure 34. Comparison of calculated (square boxes) with observed (filled
circles) LFER i.e. the correlation between the log k of the different
reactions and the corresponding pK
a
s of the leaving groups for
phosphate monoester dianions with different leaving groups. The rate
constants are given at 39 C° in units of s
1 −
. The corresponding data is
given in Table 7. The linear curve is taken from ref. (Lad, Williams et
al. 2003) where it was obtained as a least-square fit log k=1.0-1.26
pK
a
to measured hydrolysis rate constants of several
phosphomonoesters..................................................................................... 130

Figure 35. Defining the isotopic labeling of the pyrophosphate oxygens
used in Table 9and Table 10. The oxygen atoms marked with an
asterisk are substituted with an
18
O isotope whereas all other oxygen
atoms are left with
16
O................................................................................ 135

Figure 36. The calculated surface for the hydrolysis of GTP in RasGAP in
the space defined by R
1
and R
2
. The figure shows a part of the
reaction path (dashed line). The PS region is shown in Figure 37.
Energies are given in kcal/mol and distances in Å. .................................... 137

Figure 37. The calculated surface for the hydrolysis of GTP in RasGAP in
the space defined by the proton transfer coordinate, X and R
2
, fixing
1
R at the TS distance at about 2.1 Å. The figure shows the reaction
path that connects the RS with the PS and the corresponding TS.
Energies are given in kcal/mol and distances in Å. (See method
section for more details about our mapping procedure.)............................. 139
xi

List of Abbreviations

ai – ab initio;
CPHF – Coupled Perturbed Hartree-Fock;
DFT – Density Functional Theory;
ESP – Electrostatic Potential;
EVB – Empirical Valence Bond;
FDFT – Frozen Density Functional Theory;
FEP – Free Energy Perturbation;
LFER – Linear Free Energy Relationship;
LRA – Linear Response Approximation;
MM – Molecular Mechanics;
MO – Molecular Orbital;
PES – Potential Energy Surface;
QM/MM – Quantum Mechanics/Molecular Mechanics;
QSAR – Quantitative Structure Activity Relationships;
US – Umbrella Sampling;

xii
Abstract

Modeling chemical reactions in condensed phases, especially in biological
environment, such as in enzymes, is a major task of computational chemistry. At
the current computer resources one can usually not afford to solve the
Schrödinger equation for the molecules, and not even for all the electrons of the
system. Simplified classical physical models are used for most of the processes
such as protein folding, binding calculations or self assembly of molecules. For
chemical reactions, however these classical models are insufficient and quantum
treatment is necessary for the atoms participating in the reactions.
Our goal is to provide computational methods that can quantitatively
reproduce experimentally observed reaction rates and free energy changes, such
as activation free energies. For this purpose we developed methods that describe
the reactive part of the system with ab initio calculations, providing accurate
quantum mechanical treatment, and at the same time account for the
environmental effects in a consistent way. My contribution by this PhD work is to
offer several alternatives and examples how reactions in complex systems can be
studied.
The concept behind solvation models can be looked at as the first version of
"coarse graining": simplifying the fast degrees of freedom and keeping the slow
ones that contribute to the high energy transition states.
xiii
Main results of the dissertation include: (1) major developments for the
Langevin dipoles solvation method: new implementation in the Gaussian program
package, derivation and implementation of analytical gradients, new formulation
of the model that allows a variational treatment of the wave function; (2)
development and application of a new all atom solvation model and
implementation in the Q-Chem program package; (3) theoretical investigation of
the concept of Linear Free Energy Relationships (LFERs) in the context of the
transferability of the off-diagonal elements between different phases for the
Empirical Valence Bond (EVB) method and in applications for phosphate
hydrolysis reactions in water and in proteins; (4) development and applications of
a method using reference potential energy surfaces to allow for better sampling
and accuracy using Quantum Mechanics/Molecular Mechanics (QM/MM).
It is found in our work that ab initio methodologies are very powerful in
accurately modeling potential energy surfaces. The usual total energy
minimization is essential and it works well in most cases, however it is often
difficult, especially in larger systems, to include all the degrees of freedom to
follow a reaction path. In these cases one can easily end up in local minima. In
other words, not having adequate sampling of all degrees of freedom to find better
minima leads to artificially too high reaction barriers. We could overcome this
problem by reducing the dimensionality of the problem: introducing simplified
solvation models is shown to reproduce experimental activation energies for a
xiv
series of reacting molecules. Alternatively, we also used reference potentials,
where much better sampling could be achieved and we obtained converging
results for reactions in protein and water, reproducing experimental findings.

15

Chapter 1.
Introduction to the Simulations of
Enzymatic Reactions

1.1 Motivation
All atom simulations in biological systems are often restricted in size and
accuracy due to the high demand in computer time to evaluate thermodynamic
properties with proper sampling and chemical accuracy. To save computer time
empirical potential energy functionals, molecular mechanics (MM) potentials can
be constructed and used. However these potential surfaces are unable to describe
chemical reactions. Possible solutions are, to either apply accurate quantum
mechanical description using ab initio methods, in which case the calculations are
more limited in terms of system size and simulation length or one can use various
approaches with mixed level of accuracy, such as implicit description of the
environment, QM/MM or Frozen DFT (FDFT) methods.
It is of basic interest to examine what is the level of the separability in a
complex system, when one can still reproduce experimental measurements of a
reaction quantitatively, but treat some parts of the system using simplified
models. In solution these models can provide an effective potential, the solvation
energy, in addition to the accurate ab initio description of the solute. However, in
16
an enzyme, such a simplification is often not possible, so one has to include all
the atoms, but assume that their degrees of freedom are separable from the
quantum system where the reaction takes place. We try to separate the sampling
problem of the molecular mechanics part, which requires long trajectories and
many different initial conditions, from the evaluation of the wave function in the
active space. By doing this, we make sure that the two systems are properly
coupled, which means that the wave function is polarized by the environment and
also, the energy of the environment includes the electrostatic and van der Waals
interaction energies with the solute. We have followed this approach in all the
solvation models used here (including the QM/MM FEP approach) and in our
EVB as a reference model calculations.
To measure the predictive power of our approaches, there is always a great
deal of effort spent on validating methods and results of our simulations. In most
cases we validated our methods on experimental activation free energies, which
we reproduced both in protein environment and in the reference reaction in water.
One very important feature of our calculations is that they provide in depth,
detailed description of the mechanism, which is most of the times impossible to
obtain from experimental results, where one often only has a few indirect clues
about the reaction mechanism.
One of the main questions that we asked was how certain mechanisms
proceed. For example we examined the case of the phosphate hydrolysis, where a
17
phosphate ester bond is broken by the assist of a water molecule. Here one of the
debates was whether there is an associative or dissociative mechanism. We also
asked the questions how do the protonation states of the participating species
affect the reactions? What is the protonation state of the ionizable residues in the
protein? In the majority of proteins that catalyze phosphate hydrolysis, there are
metal ions in the active site that are essential to obtain catalysis. What is their
role? How do they help the catalysis?
The goal of this dissertation is to provide a methodology with which reliable
energetics can be obtained for enzymatic reactions, and based on these calculated
free energy profiles of the reactions, one could finally be able to give answers to
similar questions as the ones asked above. Understanding chemical reactions of
biological systems at the molecular level is crucial for further progress in
biochemistry and it is a task that best can be achieved by the use of computer
simulations and their validations on experimental evidence.
In more details the following section will guide the reader through the basic
structure of this dissertation. We will provide there more specific explanation on
what steps we took to accomplish our goals.

In terms of software development, the following computer programs have
been developed or modified:
• MOLARIS molecular modeling package (Chu, Villa et al. 2004): New
modules for QM/MM calculations were added;
18
• Q-CHEM ab initio program package (Shao, Molnar et al. 2006): The
MOLARIS program was coupled together with Q-Chem and QM/MM module
was created;
• GAUSSIAN (Frisch, Trucks et al. 1998): New modules of the Langevin
Dipoles solvation model and its analytic gradient evaluation task was created;
The main results of this dissertation are discussed in the following papers and
abstracts:
1. Rosta, E.; Klahn, M.; Warshel, A., Towards accurate ab initio QM/MM
calculations of free-energy profiles of enzymatic reactions. Journal of Physical
Chemistry B 2006, 110, (6), 2934-2941.
2. Klahn, M.; Rosta, E.; Warshel, A., On the mechanism of hydrolysis of
phosphate monoesters dianions in solutions and proteins. Journal of the American
Chemical Society 2006, 128, (47), 15310-15323.
3. Shao, Y.; Molnar, L. F.; Jung, Y.; Kussmann, J.; Ochsenfeld, C.; Brown,
S. T.; Gilbert, A. T. B.; Slipchenko, L. V.; Levchenko, S. V.; O'Neill, D. P.;
DiStasio, R. A.; Lochan, R. C.; Wang, T.; Beran, G. J. O.; Besley, N. A.; Herbert,
J. M.; Lin, C. Y.; Van Voorhis, T.; Chien, S. H.; Sodt, A.; Steele, R. P.; Rassolov,
V. A.; Maslen, P. E.; Korambath, P. P.; Adamson, R. D.; Austin, B.; Baker, J.;
Byrd, E. F. C.; Dachsel, H.; Doerksen, R. J.; Dreuw, A.; Dunietz, B. D.; Dutoi, A.
D.; Furlani, T. R.; Gwaltney, S. R.; Heyden, A.; Hirata, S.; Hsu, C. P.; Kedziora,
G.; Khalliulin, R. Z.; Klunzinger, P.; Lee, A. M.; Lee, M. S.; Liang, W.; Lotan, I.;
19
Nair, N.; Peters, B.; Proynov, E. I.; Pieniazek, P. A.; Rhee, Y. M.; Ritchie, J.;
Rosta, E.; Sherrill, C. D.; Simmonett, A. C.; Subotnik, J. E.; Woodcock, H. L.;
Zhang, W.; Bell, A. T.; Chakraborty, A. K.; Chipman, D. M.; Keil, F. J.; Warshel,
A.; Hehre, W. J.; Schaefer, H. F.; Kong, J.; Krylov, A. I.; Gill, P. M. W.; Head-
Gordon, M., Advances in methods and algorithms in a modern quantum chemistry
program package. Physical Chemistry Chemical Physics 2006, 8, (27), 3172-
3191.
4. Hong, G. Y.; Rosta, E.; Warshel, A., Using the constrained DFT approach
in generating diabatic surfaces and off diagonal empirical valence bond terms for
modeling reactions in condensed phases. Journal of Physical Chemistry B 2006,
110, (39), 19570-19574.
5. Klahn, M.; Braun-Sand, S.; Rosta, E.; Warshel, A., On possible pitfalls in
ab initio quantum mechanics/molecular mechanics minimization approaches for
studies of enzymatic reactions. Journal of Physical Chemistry B 2005, 109, (32),
15645-15650.
6. Rosta, E.; Warshel, A., Analytical derivatives for langevin dipoles
solvation model. Abstracts of Papers of the American Chemical Society 2003,
225, U771-U771.
For more details please visit the author’s web site:
http://www-rcf.usc.edu/~rosta/cveng0.htm

20

1.2 Overview of the dissertation

My dissertation is built mostly on the published and presented results of my
works (see (Rosta and Warshel 2003; Rosta and Warshel 2004; Klahn, Braun-
Sand et al. 2005; Rosta and Warshel 2005; Hong, Rosta et al. 2006; Klahn, Rosta
et al. 2006; Rosta, Klahn et al. 2006; Shao, Molnar et al. 2006)) and for a
complete understanding, details and explanations are included of the methods that
I used, and most of which was developed in our Laboratory (Villa and Warshel
2002).
The material in this Dissertation is organized as follows. Chapter 2 describes
implivit solvation models. We start with the continuum dielectric theory and its
applications. In Section 2.3 the Langevin Dipoles model is reviewed and new
developments are also discussed.
In Chapter 3 we would like to give a general introduction about available all
atom simulations for biological systems. Here we review the MM force field, the
QM/MM methods in general, a special QM/QM method: the Frozen DFT method
and finally the EVB. We will discuss how the potential energy is calculated and
how the free energy can be constructed from molecular dynamics simulations
using the above mentioned approaches.
21
Chapter 4 reviews the free energy perturbation (FEP) approach that was used
extensively in all EVB and FDFT calculations. We will also give a derivation here
for the Linear Response Approximation (LRA) and we can see how it can be
obtained from the FEP by introducing Taylor expansion. In Section 4.2 the
QM/MM - FEP method will be introduced, giving some preliminary results of its
performance as a solvation model.
In Chapter 5 we will guide through (Rosta, Klahn et al. 2006) how one can
obtain ab initio QM/MM free energies using model potentials. We chose the
haloalkane dehalogenase enzyme and its reference solution reaction for our
studies. The main advantage of the work discussed in this section is the use of a
reference potential which allows a fast evaluation of the free energy surface.
In the next Chapter we introduce the Linear Free Energy Relationships and
their use in computational chemistry. A brief derivation is also shown how one
can obtain such relationships. Related to the Marcus picture, we will also discuss
here the transferability of the diabatic coupling matrix elements that are
represented in the EVB method as the off-diagonal elements of the EVB
Hamiltonian. In Sections 6.3-6.6 we will show a very interesting application
where we were able to reproduce experimental LFER at several points. We
concluded that different molecules had different mechanism although they all fit
the same experimental LFER curve.

22

Chapter 2.
Simulations of Chemical Reactions
in Solutions: Implicit Solvation
Models

2.1 Introduction

Solvent molecules represent an enormous amount of extra degrees of freedom.
Treating them quantum mechanically is usually beyond the scope with the current
computer technologies. One might think that a possible solution is to include only
a few solvent molecules from the first solvation shells. This approach, however,
still represents a gas phase cluster, which behaves significantly differently from
the real solution. Although the interaction strength decreases with distance, the
number of molecules, thus the number of interactions increases cubically, so their
effect must be considered.
It was only recently realized by many researchers (Tomasi 2004) that some
treatment of the environment is unavoidable using ab initio calculations and it can
be best achieved by simplified models that account for a proper electrostatic
picture.
23
There are several options on what level one wishes to include the solvent. The
least approximation is done by QM/QM methodology, such as the ONIOM
(Maseras and Morokuma 1995) or frozen DFT type of approaches (Wesolowski
and Warshel 1993). It is usually still very expensive to calculate these systems
and the number of configurations is still very large, which is required to obtain
meaningful, convergent free energies. The sampling problem is also present, and
it is very difficult to overcome for other all atom systems, such as QM/MM
methods, however there is much progress in this direction in the recent years
(Riccardi, Schaefer et al. 2005; Kastner, Senn et al. 2006; Rosta, Klahn et al.
2006). These all atom models will be discussed in the next Chapter and
applications will be shown in Chapters 5 and 6.
Since the solvent is normally only constituted of one type of molecule, one
can easily introduce further simplifications based on the isotropy of the system.
This is done for implicit solvent models. We are going to discuss two types of
these models in this chapter: the continuum solvation models and the Langevin
Dipoles model.
In case of the Langevin Dipoles, the solvent is represented as a grid of dipoles
that can change their orientation and can be polarized due to the external electric
field. Here the solvents' translational and rotational degrees of freedom are not
explicitly handled. During my PhD work, I took this model and implemented it
into the Gaussian ab initio program package(Frisch, Trucks et al. 1998). We also
24
modified the original algorithm(Florian and Warshel 1997) and added a term,
which provides smooth surfaces and a differentiable formulation. This will be
reviewed in Section 2.3.
Continuum solvation models treat the solvent as a continuum dielectric and do
not take into account their discreteness. The solvent molecules' electrostatic
effects are still included, which is the most important and does not cancel even at
long distances. There are many variations of this model, we will review some of
them in the next section.

2.2 Continuum Solvation

According to the theory of electricity (Purcell 1965), the energy density of the
electrostatic field is expressed by:

π
ε
ρ
8
2
E
dV
dU
= =
(2.1)

where E is the magnitude of the field and ε is the dielectric constant of the
medium. For the case of the point charge the magnitude of the electrostatic field is
equal to
2
r
q
E
ε
= , where q is the charge of the ion and r is the radial distance
from its center. The volume element (in spherical coordinates) can be written as:
25
dr r dV
2
4π = . The energy of electromagnetic field outside of the ion of radius a
surrounded by an infinite medium of dielectric ε therefore is,
a
q
r
dr q
dr r
r
q
dV r U
a a a
ε ε
π
ε π
ε
ρ
2 2
4
8
) (
2
2
2
2
2
2 ∫ ∫ ∫
∞ ∞ ∞
= =








= = (2.2)








− = − = Δ
1 2
2
1 2
1 1
2
) ( ) (
ε ε
ε ε
a
q
G G G (2.3)
Equation (2.3) is the famous formula for solvation free energy of an ion
proposed first by M. Born (Born 1920).
The solvation free energy of an ion in water, which is the difference between
the free energy of an ion in water 80 = ε and in vacuum 1 = ε is given by
a
q
G G G
w
sol
2
1
80
1
2
) 1 ( ) 80 ( 





− = = − = = Δ ε ε (2.4)
The parameter a in Eq. (2.4) is the so-called solvation (Born) radius of the
ion. The Born radius is in general different from the Pauling’s crystal radius of the
ion. It can be calculated from (2.4) if the experimental value of the solvation
energy of an ion in water is known. In this sense Eq. (2.4) can be considered as a
definition of the solvation radius.

Figure 1. The figure shows the solvent excluded volume in the case of the NaF. Each atom
type has a van der Waals radius that is defined based on experimental solvation energies.
26

For molecules, one can build a collection of possibly overlapping spheres,
where the center is defined by the position of the atoms and the radius is the
above defined solvation radius. This space that is not accessible for the solvent is
called cavity. Once the shape and size of the cavity are defined, one may pass to
the solution of the electrostatic problem. Here we assume that the dielectric
constant is 1 inside the cavity and
sol
ε outside of the cavity.
The next step is solving the Poisson equation for the electrostatic potential, φ :

ρ
φ
ε
Δ =−

(2.5)
using the above defined position dependent ε and the charge density ρ . Then
the electron density is polarized according to the electrostatic potential obtained in
Eq. (2.5), and the free energy of solvation can be accounted for as half of the
interaction energy:
( ) ( )
3
1
2
cavity
solv
V
G r r dr φ ρ Δ =
∫
. (2.6)
There are several implementations of basically the above described idea in
different quantum chemistry codes. We used the COSMO and PCM models
(Mierts, Scrocco et al. 1981) that are some of the most well known approaches
that are readily available in the Gaussian program package. We have found in our
studies of phosphate hydrolysis (see Chapter 6.3) that unfortunately, none of these
methods have correctly implemented gradients in Gaussian, so geometry
27
optimizations are problematic. We also used PQS program package, which
allowed us to optimize geometries in solution using their continuum solvation
model.

2.3 Langevin Dipoles Solvation Model
The Langevin dipoles solvation model was developed in our group in 1997
(Florian and Warshel 1997). It belongs to the family of implicit solvation models,
however unlike the widely used continuum methods, this model represents the
solvent in a discrete manner.

Figure 2. Illustration of the Langevin dipoles solvation model: the water is represented as a
coarser and a denser grid of dipoles around the solute.

28
The model has the following simplifications: it assumes that the water
molecules are restricted to a grid, so their translational degrees of freedom are not
taken into account. There is a coarser grid outside and a denser grid inside. The
inner structure of the molecules is simplified to represent the water as a dipole.
This dipole is allowed to orient towards the electric field and the length of the
dipole is determined according to the Langevin function. This function can be
derived by assuming a Boltzmann average for the orientation of fixed length
dipoles at the temperature of the simulation,

Figure 3. The Langevin Function. The figure shows the proportion of the actual dipole
compared to the total dipole moment in the function of the interaction energy. This energy
depends on the external electric field and the temperature.

The dipole expression has the following form:

0 0
ˆ
( , ) ( )
total total
s s s scaling
L f r

β
μ ξ μ μ ξ = (2.7)
29
Here L is the Langevin function,
0
μ is the original dipole moment of the
solvent.
scaling
f in is a smoothing function that takes the distance to the closest
solute atom as an argument. These quantities will be discussed below in more
details.
total
s
ξ is the overall electric field at the position of the dipole:
' '
' ' 3 5
' 1 ' '
'
only for the iterative model
3 ( )
NDip
total QM s ss
s s ss s
s ss ss
s s
r
r
r r

μ
ξ ξ μ
=
≠
 
= + − +
 
 
∑

(2.8)
and
ˆ
total
s
ξ is the unit vector of the direction of this dipole:
ˆ
total
total s
s total
s

ξ
ξ
ξ
= .
We have two variations of the Langevin dipoles model. In the simpler case,
the interaction among the solvent dipoles is neglected. In this case the total
electrostatic energy only consists of the electrostatic energy coming from the QM
or solute field. The iterative model uses the total field including the electrostatic
potential from all the dipoles, this can be determined in a self consistent approach
using iterations until the dipoles converge.
The Langevin function depends on the electric field in the following way:
0
0
1
( , ) ,
total x x
total s
s x x
e e
L x
e e x k T
ξ μ
ξ μ
−
−
+
= − =
−
(2.9)
Here the temperature, T and the Boltzmann constant k appears in the
equation. A derivation of Langevin's formula can be obtained by integrating the
rotational angle of the dipole using the Boltzmann weighting of the interaction
energy with the total electric field.
30

Figure 4. Illustration of the smoothing function used to obtain continuous derivatives. In the
figure on the left side one can see the two potential surfaces with and without the smoothing
function. The right side figure shows the smoothing of the deletion of a dipole, thus
changing the step function to a sigmoid curve.

During my PhD work, I implemented a smoothing function illustrated in
Figure 4 to allow for a continuous energy functional with respect of the solute
atomic coordinates. In the original model, while the solute cavity changed, the
Langevin dipoles were removed once their distance fell below a threshold. This
caused a discontinuous behavior. Our aim was to introduce a scaling function to
allow for the continuous deletion of the dipoles with the change of the solute
cavity. We used the following function, where the
s
r
β
is the distance of the dipole
s to the closest solute atom β and
vdw
r
β
is the van der Waals radius corresponding
to the same solute atom.

( ) 1 tanh( (log( ) log( )))
vdw
scaling s
f r M r r
β β β
= + −
(2.10)
31

Figure 5. Illustration of the use of potential energy scans to find transition state barriers.
The disadvantage here can easily seen if one tries to imagine this procedure in higher than 2
dimensions.
Having well defined differentiable potential energy is the first step towards
analytical derivatives. Much effort was then spent during this work to derive
formulation for gradients and implement these analytical derivatives. We give
here a brief derivation of the formulation.
The free energy change related to solvation can be obtained by taking half of
the interaction energy between the solute and the solvent electrostatic energy.
From now on, we refer to the solute as the molecule, where the electrons are
treated by quantum mechanics; however the nuclei are treated classically using
the Born-Oppenheimer approximation. We can divide the interaction energy
between the solute and solvent into two parts: one that comes from the interaction
with the nuclei (
nuc
V ) and one that comes from the interaction with the electrons
(
elec
V ):
32
3
1 1 1
2
NDip NDip NAt
QM eff s s
s s nuc elec
s s
s
r
E q V V
r

α
α
α
α
μ
μ ξ
= = =
= = = +
∑ ∑ ∑
(2.11)
here NDip represents the number of dipoles and NAt represents the number of
solute atoms.
nuc
V can simply be calculated by the above sum, using the nuclear
charge,
nuc
q
α
of the α
th
atom.
3
1 1
NDip NAt
nuc s s
nuc
s
s
r
V q
r

α
α
α
α
μ
= =
=
∑ ∑
(2.12)
The electronic contribution to the interaction energy is not so straightforward
to obtain; one has to evaluate a one-electron integral given in Eq. (2.13) to find
the interaction with the dipoles.
1
3
1 , 1 1 , 1
1
NDip NDip nbasis nbasis
s s
elec s ij i j s ij ij
s i j s i j
s
r
V P P h
r

μ ϕ ϕ μ
= = = =
= − = −
∑ ∑ ∑ ∑
(2.13)
This integral is typically evaluated by replacing the solvent dipole by two
opposite charges that are in a very close distance to each other.
Using the chain rule, the derivative of the electronic energy part of the
interaction with respect to a nuclear coordinate can be written in the following
form:
1 1 1
2 2 2
s
ij s s s kl s
s ij ij ij ij ij
i j s i j s k l i j s
kl
h
P dE
P h P h P
dr r P r r

α α α α
μ μ
μ
∂
∂ ∂ ∂
= + +
∂ ∂ ∂ ∂
∑∑ ∑∑ ∑ ∑∑
(2.14)
Out of the three terms of Eq. (2.14), the first and third term can be obtained by
simple algebra, however the second term, requires more thoughts. Since the
33
density is implicitly dependent on the nuclear coordinate, its derivative,
ij
P
r
α
∂
∂

requires the solution of the Coupled-Perturbed Hartree-Fock (CPHF) equations.
The evaluation of these equations is quite expensive and a better solution is to
redefine the Fock matrix of the wave function, to contain the interaction with the
solvent in a variational manner. Before showing a small derivation for the
expression of the gradients in case of a variational Fock matrix, we would like to
present some results that we obtained by approximating
ij
P
r
α
∂
∂
from numerical
differences.

Figure 6. Test case of the [F…H…F]
-
system. The green solid line shows the calculated
potential energy (note that it is smooth) and the red dots are the transition state points that
have been found, starting from different starting points. The F atoms were kept in one place
while the H was allowed to move. The van der Waals radius of the atoms is illustrated on the
bottom of the figure.
F
Fixed at 3.5 Å
34
Our aim is to be able to find transition states via analytical gradients and
numerical second derivatives. Transition states can be characterized by all zero
first derivatives and all but one positive eigenvalues of the second derivatives
matrix. The one negative second derivative defines the coordinate along the
reaction path. Our first test case was the [F…H…F]
-
system, where the two
fluorides had different van der Waals radius for the purpose of obtaining a non-
symmetrical shape of the potential energy surface. In this way all the effects of
the geometry change (com-pared to the symmetrical structure) come from the
solvation. The positions of the F atoms are fixed and the H atom is allowed to
move to find the transition state point. We have started two optimizations to see
how stable the results are numerically. The converged structures are shown by the
red dots in Figure 6. It can be seen that the structures are right at the transition
state point: at the top of the hill of the potential surface.

Figure 7. Transition state corresponding to the S
N
2 reaction of
NH
3
+ CH
3
Cl [NH
3
CH
3
Cl] NH
3
CH
3
+
+ Cl
-
.
35

We also examined a more realistic case, the S
N
2 reaction of methylchloride
attacked by ammonia, resulting a methylamine cation and a chloride anion. This
reaction is an excellent prototype not only of the very important reaction types of
the S
N
2 reactions, but also for studying solvent effects on the geometrical
structure of transition states in chemistry. This is due to the fact that the reactants
are neutral molecules, however the products are ions which have very strong
interaction with the solvent and thus solvation causes a major change in the
energetics and structures in this reaction compared to the gas phase equivalent.
The same system was studied earlier by the PCM method(Amovilli, Mennucci et
al. 1998) and the transition state was not reported to be found using explicitly the
gradients, but rather a map of the energy profile was calculated as shown on
Figure 5. The disadvantage of this approach is that it takes significantly more time
to compute such maps compared to just finding transition sItates by using gradient
search.
36

Figure 8. Transition state search for the reaction of methylchloride and ammonia. We have
simultaneously mapped the potential energy surface on a grid of 0.01 Å resolution and used
gradient search to find the transition state.

This point is illustrated by Figure 8. We have mapped the potential energy
surface at the transition state on a grid defined by the attacking N and C distance
and by the leaving Cl and C distance. The grid resolution was 0.01 Å, thus we
calculated about 150 points to obtain the surface shown on the figure. At the same
time the gradient search required less than 10 points to locate the TS. The
accuracy of the transition state structure was about 0.02 Å and 0.2 kcal/mol,
which is considered within the limits of chemical accuracy.
The above examples used numerical gradients for the derivative of the density
matrix, which is very expensive to evaluate and our aim is to have fully analytical
37
gradients. For this purpose we need to modify the Fock matrix of the wave
function by adding the terms from the solvent dipoles in a way that the formal
derivative of the energy with respect to the density matrix element would be
identical to the corresponding element of the Fock matrix:
( ) E f P = (2.15)
and
( ) f
F
P
P
μν
μν
∂
=
∂
. (2.16)
In the following section we give a short derivation to find the gradients in case
when the above conditions are satisfied. The total energy derivative can be written
in the following form:
( ) ( ) ( ) ( )
F
P
E f f f f
r r r r P r
P P P P P
P

μν
μν
μν
α α α α μν α
∂
∂ ∂ ∂ ∂ ∂ ∂
= + = +
∂ ∂ ∂ ∂ ∂ ∂ ∂
∑
(2.17)
In Eq. (2.17) the density matrix, P is defined by the occupied one electron
orbitals, we use here the standard notation, the greek letters symbolize atomic
basis functions while arabic letters symbolize the molecular orbitals (MOs):
occ.
i* i
c c
i
P
μν μ ν
=
∑
(2.18)
The derivative of the molecular orbital coefficients can be expressed with a
unitary transformation:
MO
m
mi
U c
i
m
c
r
μ α
μ
α
∂
=
∂
∑
(2.19)
38
Applying this expansion on the derivative of the density matrix, one gets:
( )
occ. occ. MO
i m i i m
mi mi
c U c c U c c
i
i
i
i i m
P c
c
c
r r r
μν μ α α ν
μ ν μ ν μ ν
α α α
  ∂ ∂
∂
= + = +
 
 
∂ ∂ ∂
 
∑ ∑∑
(2.20)
Substituting Eq. (2.20) into Eq. (2.17) we get:
( )
occ. MO
m i i m
mi
( )
U c c c c
i m
E f
F
r r
P
α
μν μ ν μ ν
μν
α α
∂ ∂
= + +
∂ ∂
∑ ∑∑
(2.21)
The last term of this equation can be rearranged by reordering the
summations:
( ) ( )
occ. MO occ. MO
m i i m m i i m
mi mi
U c c c c U c c c c
i m i m
F F
α α
μν μ ν μ ν μν μ ν μ ν
μν μν
+ = +
∑ ∑∑ ∑∑ ∑
(2.22)
Here we can recognize the Fock matrix in the basis of the MOs:
( ) ( )
occ. MO occ. MO occ. MO
m i i m
mi mi mi
U c c c c U 2 U
mi im im
i m i m i m
F F F F
α α α
μν μ ν μ ν
μν
+ = + =
∑∑ ∑ ∑∑ ∑∑
(2.23)
We can use that the Fock matrix is diagonal in the basis of the MOs, so the
total energy derivative finally becomes:
occ. occ. occ.
ii ii
ii
S S ( ) ( ) ( )
2 U
ii ii ii
i i i
E f f f
F F
r r r r r r
P P P
α
α α α α α α
ε
∂ ∂ ∂ ∂ ∂ ∂
= + = − = −
∂ ∂ ∂ ∂ ∂ ∂
∑ ∑ ∑

(2.24)

These results allow us to only use the explicit derivative of the energy
expression with respect to the derivative of the atomic coordinates, the occupied
orbital energies and the derivative of the diagonal elements of the overlap matrix
of the occupied orbitals. These quantities are readily available in all ab initio
programs and thus enforcing the above form of the Fock matrix offers immense
advantages compared to the original formalism.
39
To actually realize this advantage, we need to derive the corresponding Fock
matrix element. This can be done by simply evaluating the
elec
ij
E
P
∂
∂
derivative:
1 , 1
1 1
1
2 1
2
NDip nbasis
s
s kl kl
NDip NDip
s k l s elec elec s
ij s ij s
s s s ij ij ij
P h
E
F h
P P P

μ
μ
μ ξ
= =
= =
∂ −
 
∂ ∂
= = = − +
 
 
∂ ∂ ∂
 
∑ ∑
∑ ∑ ∑

(2.25)
For the practical purposes we would like to give the detailed formulas that can
be directly implemented into any code. Given the Langevin function,
0 0
ˆ
( , ) ( )
total total
s s s scaling
L f r

β
μ ξ μ μ ξ = , we can use the chain rule to find the derivative
of this function as follows:
total total
s s s s s
total
total
ij s ij ij
s
P P P

μ μ ξ μ ξ
ξ
ξ
∂ ∂ ∂ ∂ ∂
= +
∂ ∂ ∂ ∂
∂
(2.26)
and further applying the chain rule one gets:
0
0
( , )
( , )
total total total total
s s s s s s s s
total total total total
ij s ij s s ij s ij
L
P P L P P

μ ξ μ ξ μ μ ξ μ ξ
ξ ξ μ ξ ξ
∂ ∂ ∂ ∂ ∂
= − +
∂ ∂ ∂ ∂ ∂
(2.27)
which is the same as the following:
0
0
( , )
( , )
total total total
s s s s s s s
total total total total
ij s s s ij s ij
L
P L P P

μ μ ξ μ μ ξ μ ξ
ξ μ ξ ξ ξ
  ∂ ∂ ∂ ∂
= − +
 
∂ ∂ ∂ ∂
 
. (2.28)
40
What is left is to find the derivative of the electric field with respect of the
density matrix element,
total
s
ij
P

ξ ∂
∂
. In case of the noniterative version of the
Langevin dipoles model, we can easily present the working formulas:
, 1
nbasis
s
kl kl total QM
k l s s s
ij
ij ij ij
P h
h
P P P

ξ ξ
=
∂
∂ ∂
= = =
∂ ∂ ∂
∑
(2.29)
and
ˆ ˆ
total total
total total s s s
s s ij
ij ij
h
P P

ξ ξ
ξ ξ
∂ ∂
= =
∂ ∂
(2.30)
Resulting in the following final expression for
s
ij
P

μ ∂
∂
:
( )
0
0
( , )
ˆ
( , )
total
total s s s s s s s
s ij ij total total total total
ij s s s s
L
h h
P L

μ μ ξ μ μ μ
ξ
ξ μ ξ ξ ξ
  ∂ ∂
= − +
 
∂ ∂
 
. (2.31)

Similarly, one can obtain an expression for
total
s
ij
P

ξ ∂
∂
in case of the iterative
version, however there one has to apply the same method in writing the derivative
of the dipole expression as it is given in the Appendix A. The Z-vector formalism
can be used in a straightforward manner, however in this case its use does not add
in any significant expense to the calculations, because it only concerns classical
terms and one electron integrals for the electrostatic potential.
41
We can start the case of the iterative dipoles by writing the expression for the
electric field used in this case:
' '
' ' 3 5
' 1 ' '
'
3 ( )
NDip
total QM s ss
s s ss s
s ss ss
s s
r
r
r r

μ
ξ ξ μ
=
≠
 
= + − +
 
 
∑
(2.32)
We are looking for the gradient of Eq. (2.32) with respect to the density
matrix elements:
'
' 1
'
total QM total NDip
s s s s
s
ij ij s ij
P P P

ξ ξ ξ μ
μ
=
∂ ∂ ∂ ∂
= +
∂ ∂ ∂ ∂
∑
(2.33)
To obtain the derivative of the dipole, we can start by writing the equation for
the dipole noting that it depends on the electric field coming from the solute and
the other dipoles, all other parameters can be taken as constants with respect to
the density matrix elements:
( ) '
,
QM
s s s s

μ μ ξ μ =
(2.34)
Writing the derivative then becomes a set of linear equations that has the
following form:
( )

'
'
'
'
' 1
'
,
ss
s s
QM
QM NDip
s s s
s s s s
QM
s
ij ij s ij
s
A
t t
d
dP P P

μ ξ μ
μ ξ μ μ
μ
ξ =
∂ ∂ ∂ ∂
= +
∂ ∂ ∂
∂
∑
(2.35)
After solving Eq. (2.35), the results can be substituted into Eq. (2.33) and a
expression for the Fock matrix element can be obtained that will satisfy Eq.(2.25).
42

Chapter 3.
All atom simulations

3.1 Introduction
According to the previous chapter, we were able to show how successful
implicit solvation models are. For biological systems, however one does not deal
with a pure solvent, but rather a complex chemical environment, like proteins,
DNA or lipids. These systems are no longer isotropic: evolution used the fact that
certain charge distribution, polarities, molecular cavities and flexibilities can
enhance the chemical and physical behavior of biological systems. This means
that modeling these systems is much more difficult and a simplified description
often does not reflect reality. In this chapter we would like to review the all atoms
models that are generally used for studying enzymatic reactions.
We start this by introducing the QM/MM (Warshel and Levitt) method that is
now widely used in the broad field of molecular modeling and just very recently
ab initio codes also implement and use this method.
We would also like to review a less well known approach: the FDFT or
Constrained DFT method. This QM/QM approach offers very accurate
description of the system under study, but it does not cost as much as the full DFT
43
treatment. It also has a key feature: it can describe diabatic states as well, not only
adiabatic ones, like the usual QM wave function.
Another approach that we would like to describe here, is the EVB method
(Warshel and Weiss 1980). This method has been developed and used extensively
in our group and recently it has many new implementations by others as well.
This approach is not formally a quantum mechanical approach, it is rather a
molecular mechanical treatment, however it can be used to describe chemical
reactions unlike the regular MM. Its main advantage over the ab initio description
is that it costs only as much as the classical force field calculations do. The
method has been successfully applied for a number of problems that otherwise
could not have been tackled by ab initio approaches (Warshel 1991; Shurki and
Warshel 2003).

3.2 QM/MM
First principle approaches are at present too expensive for effective modeling
of very large molecules like enzymes. One can utilize the fact that
macromolecules are assembled by the same type of bonds that connect the atoms
in small molecules. Quantum mechanical description can be successfully replaced
by a molecular force field to describe the environment, assuming that it is
properly coupled to the reactive quantum system. The force field potential surface
44
is expressed as a sum of contributions from bonded atoms and interactions
between nonbonded atoms:
( ) ( ) ( )
,
,
b nb
U U b U U r
φ
φ
Θ
= Θ + + (3.1)
Here b,Θ ,φ and r are, respectively, the vectors of bond lengths, bond angles,
torsional angles, and the vector of Cartesian coordinates which are used to
evaluate the nonbonded distances. The first term defines a very deep potential
well, and since the molecule stays in most cases inside this well (except in
extreme cases of bond dissociation), it is reasonable to approximate this part of
the potential surface by its quadratic expansion, which is given by
( ) ( ) ( )
2 2
, , 0, , 0,
1 1
, Cross terms
2 2
b b i i i i i i
i i
U b K b b K
Θ Θ
Θ = − + Θ −Θ +
∑ ∑
(3.2)
where U is usually given in kcal/mol, b in Å and Θ in radians. The torsional
potential U
φ
is a periodic function, which can be described by the leading terms
in the Fourier expansion of the potential.
( ) ( )
,
1
1 cos
2
i i
i
U K n
φ φ
φ φ = −
∑
(3.3)
The nonbonded potential,
nb
U can be described by an atom-atom interaction
potential of the form
( )
12 6 i j
nb ij ij ij ij ind
ij ij
q q
U A r B r C U r
r
− −
= − + +
∑
(3.4)
45
where
ij
r is the distance between the indicated atoms, the
i
q 's are the residual
atomic charges, and
ind
U is the many body inductive effect of the electronic
polarizabilities.
Evaluating the potential energy surface at some extrema usually does not
allow one to obtain macroscopic observables, although in some cases energy
minimization can lead to very good estimate of the enthalpies. To be able to
quantitatively predict experimentally measurable macroscopic quantities, one has
to evaluate thermodynamic functions, such as the free energy, enthalpy and
entropy. This can only be done by the sampling of the configurational space. To
find relevant configurations for evaluating average properties, one can use several
methods. These methods include for example Monte Carlo, Molecular Dynamics
(MD), accelerated MD, Replica Exchange approaches. Throughout our
calculations, we used the MD simulation technique.
Molecular Dynamics (Hill, Tokmakoff et al.) simulations evaluate the motion
of the atoms in a given system and provide the positions or trajectory of these
atoms as a function of time. The trajectories are calculated by solving the classical
equation of motion (Newton's law) for the molecule under consideration.
i i i
i
U
m r F
r
∂
= =−
∂
(3.5)
The strength of MD approaches is associated with the fact that they have the
ability to simulate, at least in principle, the true microscopic behaviors of
46
macromolecules. The weakness comes from the fact that some properties reflect
extremely long time processes that are not accessible by current computers.
The time needed to reach accurate results for average properties depends on
the model used and number of local minima. Eventually, with the increase of
computer power it has become feasible to reach simulation times of nanoseconds
and to start to obtain meaningful average properties of macromolecules.
For example, models that trim the protein to a sphere with the proper spherical
boundary conditions (Aqvist 1996) converge much faster than models that involve
periodic boundary conditions (van Gunsteren and Mark 1992) since the latter
involves more degrees of freedom and more minima. As much as accuracy is
concerned, the proper treatment of long-range effects is crucial and improved
convergence is usually associated with more proper treatment of long-range
forces (Lee, Chu et al. 1993).
QM/MM approaches have provided a general scheme for studies of chemical
processes in proteins e.g. ref. (Warshel and Levitt 1976; Field, Bash et al. 1990;
Bakowies 1996; Gao 1996; Friesner and Beachy 1998; Monard and Merz 1999;
Garcia-Viloca, Gonzalez-Lafont et al. 2001; Field 2002; Shurki and Warshel
2003). A significant progress has been made with calibrated semiempirical
QM/MM approaches (e.g. ref. (Shurki and Warshel 2003)) that includes careful
evaluations of the relevant activation free energies by free energy perturbation
approaches that date back to the 80’s (e.g. ref. (Warshel, Sussman et al. 1988)). In
47
current works, an ab initio representation is often used within the framework of a
QM/MM treatment, since such QM(ai) representations have been shown to
provide “chemical accuracy” in studies of gas phase reactions of small molecules.
Unfortunately it represents at present a very challenging or even unpractical
technique to try to evaluate the potential of mean force (PMF) for enzymatic
reactions by QM(ai)/MM approaches due to the requirement of a very extensive
sampling and thus extremely expensive repeated evaluation of the QM energies.

3.3 Frozen DFT Approach
A very promising option for treating larger systems by first principle
approaches is provided by the frozen DFT (FDFT) and constraint DFT (CDFT)
approaches (Wesolowski and Warshel 1993; Hong, Strajbl et al. 2000). The basic
idea behind these approaches is to treat the environment quantum mechanically,
but freezing (or constraining) their electron density. In this way the entire system
is treated by an ab initio DFT approach, and a formally rigorous nonadditive
kinetic energy functional evaluates the coupling between regions of the active site
and the environment.

The FDFT approach presents a general way of coupling two subsystems by
means of an orbital-free and first-principle effective potential. This makes it
possible to cast the concept of "embedding potential" in DFT terms. The FDFT
48
approach is related in a formal way to the work of Cortona (Cortona 1991), who
did not deal, however, with the issue of embedding a subsystem in a larger
system, which is described by a more approximate method. Wesolowski &
Warshel (Wesolowski and Warshel 1993) realized that the coupling term in any
hybrid method could be obtained by a partial minimization of the total energy
functionals. The CDFT embedding idea has been adopted in related fields such as
studies of molecules on metal surfaces (Kluner 2001).
The total energy of an N-electron system is defined within the Kohn-Sham
formulation as:
1 ( ) ( ')
[ ] [ ] ( ) ( ) ' [ ( )]
2 | ' |
ext xc
E T V d d d E
ρ ρ
ρ ρ ρ ρ = + + +
−
∫ ∫∫
r r
r r r r r r
r r
(3.6)
where ( )
ext
V r is the external potential, [ ( )]
xc
E ρ r is a density functional of the
exchange-correlation energy and
*
1
( ) ( ) ( )
N
i i
i
ρ ϕ ϕ
=
=
∑
r r r (3.7)
* 2
1
1
[ ( )] ( )( ) ( )
2
N
i i
i
T d ρ ϕ ϕ
=
= − ∇
∑
∫
r r r r (3.8)
In the FDFT/CDFT approach, we divide the system into sub-systems, and
apply Kohn-Sham formulation to each subsystem separately, taking into account
the interaction with the other subsystems. For instance, if we divide the whole
system into two sub-systems with N
I
and N
I'
electrons respectively, the total
energy of the system is defined as:
49
'
' '
[ , ]
1 ( ) ( ')
[ ] [ ] [ , ] ' ( ) ( ) [ ( )]
2 | ' |
CDFT
I I
nadd
I I s I I ext xc
E
T T T d d V d E
ρ ρ
ρ ρ
ρ ρ ρ ρ ρ ρ
=
+ + + + +
−
∫∫ ∫
r r
r r r r r r
r r
(3.9)
( ) d
I I
r N ρ =
∫
r (3.10)
( )
' '
d
I I
r N ρ =
∫
r (3.11)
'
( ) ( ) ( )
I I
ρ ρ ρ = + r r r (3.12)
' ' '
[ , ] [ ] [ ] [ ]
nadd
s I I s I I s I s I
T T T T ρ ρ ρ ρ ρ ρ = + − − (3.13)
where [ ]
s
T ρ is defined as the orbital-independent kinetic energy functional of
the electron density of the system (Wesolowski 1997). Taking the electron
density when region I’ is isolated and minimizing the total energy with respect to
the orbitals of region I leads to the following modified Kohn-Sham equation:
( ) ( ) ( )
2 #
1
2 , eff I i i i
V ϕ ε ϕ   − ∇ + =
 
r r r , (3.14)
where the effective potential ( )
#
, eff I
V r is defined as:

( ) ( )
( ) ( )
( )
'
#
,
,
nadd
I I
KS
eff I eff
I
T
V V
δ ρ ρ
δρ
 
 
= +
r r
r r
r
, (3.15)
with the Kohn-Sham potential, ( )
KS
eff
V r given by:

( ) ( )
( ) ( )
( )
'
' '
d ' d '
' '
I I KS
eff ext xc
V V V
ρ ρ
ρ = + + +  
 
− −
∫ ∫
r r
r r r r r
r r r r
. (3.16)
Once ( )
I
ρ r is known a counterpart equation for region I’ could be derived.
These two coupled equations can be solved iteratively. In case a part of the
50
system is described classically we use a FDFT - QM/MM formulation where the
effect of the MM part is incorporated in the FDFT Hamiltonian as described in
ref. (Olsson, Hong et al. 2003).

3.4 Empirical Valence Bond (EVB) method
The EVB (Empirical Valence Bond) method (Warshel and Weiss 1980) is a
simple and effective way of introducing bond forming and bond breaking in the
formalism of force field, MM methods. This is very important since the modeling
of chemical reactions otherwise requires a quantum mechanical treatment, which
is usually not affordable at the level of MD.
The EVB method represents the potential energy surfaces of proteins by a
combination of a classical empirical force field and an empirical valence bond
force field. The classical force field is used to simulate the parts of the protein
removed from the actual chemical reaction being studied since there is no bond
breaking or making in this region and the classical force field is extremely well
suited for the representation of molecules that are not undergoing chemical
transformations. In the small region of the protein where there is a chemical
reaction taking place, an empirical valence bond formalism is used. The energy is
simply built by mixing two force field potentials, one that describes the reactant
states and one that describes the product states.
51
The ground state potential energy surface (PES) is constructed by creating a
two by two empirical Hamiltonian, where the diagonal elements consist of the
force field energies that describe the reactant or the product and the off diagonal
element is either an exponential decay function of the distance or similar
parameter or sometimes taken as a constant, which is fitted to experimental
results. In certain cases, one can have more than two "diabatic" states that will
make it possible to describe several reaction steps or intermediate
structures(Shurki 2006). In this case, of course, the empirical Hamiltonian will be
of higher dimension, according to how many states we will have.
More specifically, the diagonal elements are described by
) ( 6 12
2
166 ) ( 332 ) (
) 2)( ( )}] ( exp{ [1
2
i
Δ + + + − − + − +
− + − − − = =
∑ ∑ ∑
∑ ∑
′
−
′
′
−
′
′
−
′
′
ss
(i)
Ss
k' k,
4
k' k, k k,
r μ (i)
k
k k,
(i)
k k k,
k k,
k k,
(i)
m 0,
m
(i)
m θ
m
(i)
m 0,
(i)
m i ii
U U r α r / exp 1 Q Q r B r A
θ θ / K b b β D ε H
k k, s

(i)
ss
(i)
Ss
0
i
U U ε Δ + + + = (4.2)
( ) [ ]
ij
r Δ − = μ exp A H
ij ij

where the b
(i)
’s and θ
(i)
's are, respectively, the bond lengths and bond angles
in the quantum mechanical system composed of the n bases and the excess
protons. The r
k,k’
runs over all the nonbonded distances in the quantum system.
The r
-4
term represents an approximation for the inductive interaction between the
solute charges.
(i)
Ss
U describes the interaction between the quantum system (the
“solute”) in its i-th state and the surrounding classical system (the “solvent”),
52
which includes water molecules and/or protein atoms.
ss
U is the solvent-solvent
classical potential surface. Finally,
) (i
Δ is the so-called “gas phase shift” that
determines the relative energy of the diabatic states(Warshel 1991). Other than the
usual force field parameters of an MM force field, there are two more extra
parameters, the off diagonal elements (the H
ij
s) and the gas phase shift. These are
generally fitted to experimental information about the reaction barrier and the
reaction free energy. The ground state energy, E
g
, is obtained by diagonalizing the
resulting EVB Hamiltonian.
g g g
E C HC = (4.3)
We would like to emphasize one very important advantage of the EVB
potential. One can define the reaction coordinate in a Hilbert space fashion: there
are two well defined diabatic Hamiltonians, which can be used to drive the
reaction from one state to the other. This is especially important in condensed
phases, since it helps capturing the solvent reorganization energy and includes
important entropic effects of both the solute and the solvent.
The free energy profile is thus built by using a free energy
perturbation/umbrella sampling (FEP/US) method. We start the evaluation of the
free energy by running molecular dynamic simulation on a mapping potential of
the form
( )
11 22
1
m m m
U H H λ λ = − + (3.17)
53
Here λ
m
gradually increases from 0.0 to 1.0 with a small constant increment in
n+1 fixed increments ( 0 / , 1/ , , /
m
n n n n λ = … ) and the system is driven gradually
from reactant state to product state. The free energy contribution associated with
the change of U
m
to U
m’
can be regarded as a perturbation, and the corresponding
free-energy difference [ ]
m m'
G δ λ λ → can be obtained with free energy
perturbation (FEP) method:
[ ] { } ( ) { }
' '
exp exp
m m m m
m
G U U δ λ λ β β − → = − − (3.18)
where
m
<> is a arithmetic average of the quantity defined inside the bracket
over the configurations sampled by moving the entire system on potential surface
U
m
. The overall free energy profile with respect to λ is a simple summation:
( ) ( )
1
1
0
m n
n m m
m
G G λ δ λ λ
= −
+
=
Δ = →
∑
(3.19)
This approach will force the system to fluctuate near the TS at some of the
m
λ
values, while allowing all geometrical coordinates to move freely, including the
solulte and solvent degrees of freedom.
The free energy, Δ G
m
, associated with changing θ from 0 to m/n is evaluated
by the FEP procedure that will be described in chapter 4.1 (see, also chapter 3.3.2
in (Warshel 1991)). The free energy functional that corresponds to the adiabatic
ground state surface, E
g
, Eq. (4.3) is then obtained along the reaction coordinate,
1 2 22 11
x H H ε ε = − = − , which can be written as:
54
( ) ( ) ( ) ( ) ( )
1
' ln ' exp
m
m g m
g x G x x E x x
ε
β δ β ε
−
 
Δ = Δ − − − −
 
(3.20)
where x is the generalized reaction coordinate, taken in the EVB approach as
the energy gap between the reactant and product diabatic potentials x=ε
j
-ε
i
. ε
m

in Eq. (3.20) is the mapping potential that keeps x in the region of x’. If the
changes in ε
m
are sufficiently gradual, the free energy functional g(x’) obtained
with several values of m overlap over a range of x’, and patching together the full
set of Δ g(x’) gives the complete free energy curve for the reaction.
-300 -200 -100 0 100 200 300
0
50
100
150
200
Free energy (kcal/mol)
Δε
50% ε
react
+ 50% ε
prod
70% ε
react
+ 30% ε
prod

Figure 9. Illustration of the Umbrella Sampling technique. A good sampling can be achieved
even at the transition state region due to the fact that some of the mapping potentials have
minima, rather than maxima at the transition state along the reaction coordinate. In practice
we first evaluate the free energy changes corresponding to the change between the mapping
potentials and then we can use the
m
G Δ we obtained to find free energy surfaces
corresponding to other potential energies such as E
g
.
The FEP/US approach may also be used to obtain the free energy functional
of the isolated diabatic states. For example, the diabatic free energy, Δ g
1,
of the
reactant state can be calculated as
55
(3.21)

The diabatic free energy profiles of the reactant and product states represent
microscopic equivalents of the Marcus’ parabolas (Marcus 1964; Marcus 1993).
The reorganization energy, λ
ij
, of the i→j reaction step can be obtained directly
from the plots of the diabatic curves (as in Figure 11) or by using
( )
j
i j
i
i j j i
ε ε − − − =
→
ε ε 2 1 λ (3.22)
The use of the energy gap as a reaction coordinate appears to provide a very
powerful tool for many studies since the electrostatic contribution to the energy
gap corresponds to the generalized solvent coordinate, which is hard to explore by
other approaches. The EVB approach also provides convenient analytical first
derivatives that make it easy to explore the actual dynamics of the system (being
limited, however, by the available computer time).
The EVB method satisfies some of the main requirements for reliable studies
of enzymatic reactions. Among the obvious advantages of the EVB approach is
the facilitation of proper configurational sampling and converging free energy
calculations. This includes the inherent ability to evaluate non-equilibrium
solvation effects (Villa and Warshel 2001). Another important feature of the EVB
method is the ability to capture correctly the linear relationship between activation
free energies and reaction energies (LFER) observed in many important reactions
(for example, (Warshel 1991)). Furthermore, the EVB benefits from the
aforementioned ability to treat consistently and conveniently the solute-solvent
56
coupling. This feature is essential not only in allowing one to properly model
charge-separation reactions, but also in allowing for a reliable and convenient
calibration. Calibrating EVB surfaces using ab initio calculations was found to
provide quite reliable potential surfaces.

Figure 10. Illustration of the EVB method: the free energy of the reaction is built from the
force field defining the reactants and another force field defining the products. The ground
state energy is obtained as the eigenvalue of an empirical 2 by 2 Hamiltonian using the
diabatic states' energies and an empirical off-diagonal element.

The empirical valence bond (EVB) model thus provides an extremely
powerful way for modeling and analyzing chemical reactions in solutions and
proteins. This model, as we use it in practice, is based on the unverified
-300 -200 -100 0 100 200 300
0
50
100
150
200
Free energy (kcal/mol)
Δε
12
12 prod
react
H
H
H ε
ε  
=
 
 
57
assumption that the off diagonal elements do not change significantly upon
transfer of the reacting system from one phase to another. This ad hoc assumption
has been rationalized by its relationship to empirically observed linear free energy
relationships as well as by other qualitative considerations. Nevertheless, this
assumption has not been established in a more rigorous way. This present work
in Chapter 6.2 explores the validity of the above EVB key assumption by a
rigorous numerical approach. This can be done by exploiting the ability of the
frozen density functional theory (FDFT) and the constraint density functional
theory (CDFT) models to generate convenient diabatic states for QM/MM
treatments and thus to examine the relationship between the diabatic and adiabatic
surfaces, as well as the corresponding effective off diagonal elements. It has been
found that, at least for the test case of S
N
2 reactions, the off diagonal element does
not change significantly upon moving from gas phase to solutions and thus the
EVB assumption is valid and extremely useful.

58

Chapter 4. Free energy calculations

4.1 The free energy perturbation (FEP) and the
linear response approximation (LRA) methods
With a given set of potential functions one can simulate experimentally
observed macroscopic properties using microscopic models. According to the
theory of statistical mechanics, one should consider all the quantum mechanical
energy levels of the system in order to evaluate different average properties
(McQuarrie 1976). However, in the classical limit it is possible to approximate the
average of a given property, A (which is independent of the momentum of the
system), by (McQuarrie 1976):

(4.1)

(4.2)
where dr designates the volume element of the complete space spanned by the
3n vector r associated with the n atoms of the system. The evaluation of Eq. (4.1)
requires us to explore all points in the entire configuration space of the given
system. Such a study of solvated macromolecules is clearly impossible with any
of the available computers. However, one can hope that the average over a limited
59
number of configurations will give similar results to those obtained from an
average over the entire space. With this working hypothesis we can try to look for
an efficient way of spanning phase space.
Evaluation of free energies by statistical mechanical approaches is extremely
time-consuming due to sampling problems. Fortunately, it is possible in some
cases to obtain meaningful results using perturbation approaches. Such
approaches exploit the fact that many important properties depend on local
changes in the macromolecules so that the effect of the overall macromolecular
potential cancels out. Such calculations are usually done by the so-called free-
energy perturbation (FEP) method (Zwanzig 1954; Valleau and Torrie 1977) and
the related umbrella sampling method (Valleau and Torrie 1977). This method
evaluates the free energy associated with the change of the potential surface from
U
1
to U
2
by gradually changing the potential surface using the relationship

(4.3)
The free-energy increment associated with the change of U
m

to U
m’
can be obtained by (Valleau and Torrie 1977)

(4.4)

where
m
indicates that the given average is evaluated by propagating
trajectories over U
m
. The overall free energy change is now obtained by changing
the
m
λ in n equal increments and evaluating the sum of the corresponding G δ :
60

(4.5)
The FEP approach has been used extensively in studies of free energies of
biological systems (e.g., (Kollman 1993) ). It must be emphasized here that the
convergence of FEP approaches is quite slow and that obtaining meaningful
results requires proper treatment of long-range effects (see below).
FEP calculations of solvation free energies or related properties for large
cofactors or ligands are expected to converge extremely slowly. A promising
strategy may be provided by exploiting the fact that electrostatic effects in
solutions (and probably proteins) seem to follow the linear response
approximation (Kalra, Garde et al.). Our derivation of the LRA method is based
on studying the functions that describe the free energies of the reactant (0) and
product (1) states. The LRA (Kalra, Garde et al.) can be derived from the FEP
formulas by using the first order Taylor expansion of the following equation for a
single step reverse and forward FEP (see e.g. (Lu, Singh et al. 2003), eq. (6)):

( ) ( )
( )
( )
1 0
0 1
0 1
/ 2
0
/ 2
1
U U
G U U
U U
e
e
e
β
β
β
− −
− Δ →
− −
= (4.6)
from this equation we express the free energy change
1 0
A A − :
( ) ( )
( )
( )
( )
( )
1 0 0 1
/ 2 / 2
0 1
0 1
ln ln
U U U U
G U U e e
β β
β
− − − −  
− Δ → = −
 
 
(4.7)
61
Now we examine the expression
( )
( )
1 0
/ 2
0
ln
U U
e
β − −
and use the first order
Taylor expansion of the exponential function:
( )
( )
( )
( ) ( )
1 0
/ 2
1 0
0 0 0
0
ln ln 1 / 2 ln 1 / 2 / 2
U U
e U U U U
β
β β β
− −
≅ − − = − Δ ≅− Δ
(4.8)
Where we also used the first order Taylor expansion of the logarithm function,
to we obtain part of the final LRA expression:
( ) ( ) ( )
0 1
0 1
1/ 2 G U U U U β β β − Δ → ≅ − Δ − Δ (4.9)
So the working formula simplifies to:
( ) ( )
0 1
0 1
1/ 2 G U U U U Δ → = Δ + Δ (4.10)
It is important to note that the above derivation does not require that U Δ
would be very small, but rather that the fluctuations from the average value of
U Δ would be small, since the main component of U Δ can be evaluated
outside of the averaging and it results the same formula as (4.10).
62

Figure 11. Showing the free energy functions of the reactant and product states, ( )
a
g
X

and ( )
b
g
X
, as two Marcus' type parabolas of equal curvature. The reorganization
energies at the minima of
a
g and
b
g are given by
b b a a b
a
U U G λ
→
= − − Δ and
a a b a b
b
U U G λ
→
= − − Δ . By assuming
b a
λ λ = one can obtain the LRA estimate of
b a
G
→
Δ (see (Sham, Chu et al. 2000)).

The microscopic validity of the LRA was established in (Warshel 1982;
Hwang and Warshel 1987; Kuharski, Bader et al. 1988; King and Warshel 1990)
and subsequent studies, and the general result of Eq. (4.10) was introduced in
(Lee and Warshel 1992) (for an excellent discussion, see (Hummer and Szabo
63
1996)). This equation has been found to provide a powerful way of estimating the
results of FEP calculations of proteins.

4.2 QM/MM FEP Solvation Model
While the accuracy of the simplified solvent models mentioned above has
been found to be reasonable in earlier studies, it is important to examine the
possible effect of moving to a more atomistic solvent model. This is done here by
using a special combined ab initio quantum mechanics molecular mechanics free
energy perturbation treatment (this approach is referred to here as a QM(ai)/MM-
FEP approach). In general we advocate the evaluation of the free energy surface
of reactions in solution by using a QM(ai)/MM approach that uses the EVB as a
reference potential(Rosta, Klahn et al. 2006). This approach is particularly
effective when one allows the solute coordinates to change during the mapping
process. Here, however, we restrict ourselves to a simpler approach where the
solute coordinates are fixed at each point (the missing entropic effects will be
considered below). In this case we are interested primarily in the solvation of the
fixed solute and we can use the relationship:
( )
1 2 1 2 1 2
( , ) , ( , , ) ( ) (1)
g S g s
gas solv S pol S S
G R R E R R G R R G Δ = Δ + Δ + ΔΔ → Q Q Q
(4.11)
Where
g
S
Q and
s
S
Q are the vectors of the partial charges of the solute in the
gas phase and in solution (when the solvent potential is included in the ab-initio
Hamiltonian). These partial charges are taken as the charges derived from the
64
electrostatic potential (ESP) according to the Merz-Singh-Kollman scheme(Singh
and Kollman 1984) (Besler, Merz et al. 1990).
S
pol
G ΔΔ is the solute energy upon
change of its electron density from gas phase to its solution value. The solvation
free energy term ) , , (
2 1
s
S solv
R R G Q Δ is evaluated by using a standard adiabatic
charging (AC) free energy perturbation (FEP) procedure (e.g. refs (Warshel
1991)), changing the solute residual charges from zero to their value in the gas
phase in n+1 integration steps using:
m
s
S m
λ Q Q =
(4.12)
and representing the mapping potential by
) (
ss Ss m m
U U U + = Q
(4.13)
where
Ss
U and
ss
U are the solute-solvent and solvent-solvent potentials parts
of the potential surface as represented by the ENZYMIX force field. With this
mapping potential we use the standard FEP equation
1
1 1
1
1 2 1
1
( ) ln exp( ( ) / ) )
( , , ) ( ) (4)
m
solv m m m m U
n
g
solv S solv m m
m
G U U RT
G R R G
λ λ β
λ λ
−
+ +
+
+
=
ΔΔ → = − < − − >
Δ = ΔΔ →
∑
Q
(4.14)
where (1/ )
B
k T β = ,
B
k is the Boltzmann constant and T is the absolute
temperature. In this QM(ai)/MM FEP approach we fix the solute charges ,
g
S
Q .
Thus the effect of polarizing the solute (the last term in Eq. (4.11)) is only
65
evaluated in a mean field using the LD solvent model as discussed in ref. (Florian
and Warshel 1997). We also have an alternative approach where
)) ( , , ( ) , ( ) , (
2 1 2 1 2 1
s s
S solv gas
r R R G R R E R R G Q Δ + Δ = Δ
(4.15)
where now the potential
m
U fluctuates upon fluctuations of
s
S
Q as a result of
the solvent fluctuations. Our preliminary studies indicate that both approaches
give similar results due to the reasonable physical picture provided by the last
term of Eq. (4.11).
Similar methods were used by several groups (Ishida and Kato 2003; Kastner,
Senn et al. 2006); (Zhang, Liu et al. 2000). It was found (Kastner, Senn et al.
2006) that it is a very good approximation to replace the QM density by point
charges and only use these charges to evaluate the QM -- MM interaction energy.
It was also investigated what is the error that is caused by not polarizing the QM
wave function and this error was slightly greater than the error caused by
replacing the density with point charges, however in the case of p-
hydroxybenzoate hydrolase (PHBH) it was less than 3 percent of the activation
energy, also less than 3 kcal/mol. Their results were also compared to
thermodynamic integration and energy minimization of the total energy and all
the three results gave similar results.
The above mentioned two QM/MM FEP approaches evaluate the minimum
energy path by optimization of the QM/MM total energy and then use the FEP to
find the free energy from the MM part between the minimized QM structures
66
along the reaction path. Our method differs in the sense that we obtain absolute
solvation energy estimates, thus solvation free energies for each QM structure and
do not need, at least in principle to evaluate the FEP along the whole path.
We present here some preliminary results using our solvation method of the
calculation of the pK
a
s of small organic acids.
Table 1. Comparison of the calculated pK
a
values for the QM/MM-FEP and the COSMO
methods. The experimental values are significantly closer to the all atom model than to the
continuum results.
pK
a

Molecule QM/MM-FEP COSMO Experiment
HCOOH 5.8 12.4 3.75
CH
3
OH 19.3 34.9 15.5

The pKa values of the two different acids were calculated using the QM/MM-
FEP solvation model and the results were compared to the COSMO solvation
model implemented in the Gaussian 98 program (Frisch, Trucks et al. 1998). In
case of both methods, we applied 6-311++G** B3LYP/DFT level of theory to
evaluate the gas phase molecular geometries and then determined the solvated
polarized wave function using the COSMO model. The ESP charges were used
for determining the QM/MM-FEP solvation energy that were taken from the
COSMO polarized wave function by fitting the atomic charges to the electron
density of the system as described earlier in this Section. The total energy was
constructed by using the polarized
pol gas pol
H Ψ Ψ
energy of the wave function
and adding the solvation contribution.
67
The results are shown in Table 1. The COSMO model clearly overestimates
the experimental pK
a
values by 10-20 pK
a
units while the QM/MM-FEP only has
about 2-4 pK
a
unit error.

68

Chapter 5. QM/MM free energies
using reference potentials

5.1 Introduction
Recent realization of the importance of the proper sampling of QM(ai)/MM
sampling led to several advances(Muller and Warshel 1995; Sakane, Yezdimer et
al. 2000; Zhang, Liu et al. 2000; Strajbl, Hong et al. 2002; Iftimie, Salahub et al.
2003; Liu, Wood et al. 2003; Olsson, Hong et al. 2003; Crespo, Marti et al. 2005;
Hu and Yang 2005; Pradipta 2005; Rod and Ryde 2005). The main direction of
the recent advance has been based on different adaptations (Sakane, Yezdimer et
al. 2000; Iftimie, Salahub et al. 2003; Liu, Wood et al. 2003; Crespo, Marti et al.
2005; Pradipta 2005; Rod and Ryde 2005) of the idea presented in our group's
work (Muller and Warshel 1995; Strajbl, Hong et al. 2002) of using a classical
potential as a reference for the QM/MM calculations. However, a key point that
has been overlooked in some of the recent studies (e.g. ref. (Liu, Wood et al.
2003; Rod and Ryde 2005)), is the challenge of obtaining the free energy of the
solute and not just the surrounding solvent. That is, simulations with a fixed
solute miss the crucial entropic contribution of the solute configuration and this
contribution cannot be estimated by the harmonic approximation at least when
one deals with activation free energies of reactions in condensed phases.
69
Furthermore, calculations that estimate the QM/MM solute path along a fixed
energy minimized path may reflect artificial minima in the complex landscape of
the enzyme active site(Klahn, Braun-Sand et al. 2005). Now, as demonstrated in
our early works it is rather simple to obtain converging results using a classical
reference potential once the solute is fixed(Bentzien, Florian et al. 1998). The
problems start when the solute is allowed to fluctuate. In this case even with the
EVB as a reference potential there can be a significant difference between the
QM/MM potential and the reference potential.
A significant part of our effort in the last few years has been focused on the
challenge of obtaining converging results when the solute is allowed to fluctuate.
This led to a development of a LRA treatment that allowed us to obtain a
reasonable convergence even in the challenging case of the auto dissociation of
water in water(Strajbl, Hong et al. 2002). However we were less successful in
obtaining the activation barrier since the LRA formulation was not derived for
this case. Thus we were unable to provide a general prescription for evaluating
QM(ai)/MM activation free energies.
In the present work we report a modification of our LRA treatment that finally
allows us to obtain converging QM(ai)/MM activation free energies while
considering the effect of the substrate motion. This presents, in our view, a very
significant progress in QM/MM simulations, where it is finally possible to
70
evaluate the activation free energies of enzymatic reactions by ab initio QM/MM
approaches and doing this without enormous computational effort.
Our approach is demonstrated in a study of the S
N
2 reaction of haloalkane
dehalogenase and its corresponding reference solution reaction. It is shown that
we can now obtain the catalytic effect of enzymatic reactions by QM(ai)/MM
approaches.

5.2 Methods and systems
Our approach is based on the thermodynamic cycle of Figure 12. In this cycle
we have a reference free energy surface
EVB
g Δ , which is used to evaluate the
QM/MM activation free energy.

Figure 12. The thermodynamic cycle used to evaluate the QM/MM activation barrier.

In order to progress in a rigorous way, we start by defining the configurational
partition function of the system
71

( ) ( ) , E x g x
Q e dxd c e dx
α α
τ β
α τ
τ
∞
−Δ −Δ
−∞
= =
∫∫ ∫
(5.1)
where E
α
is the given potential surface (α is taken as a or b for the EVB or
QM/MM surfaces, respectively), x is the reaction coordinate and τ are all the
coordinates perpendicular to the reaction coordinate, g
α
Δ is the resulting free
energy profile along the reaction coordinate and c
τ
is a constant.
Next we define the partition function and free energy of the reactant (RS) (or
product (PS)) state by
( )
( )
( )
( )
‡ ‡
- g -
e e
x x
x G RS
x
x
q x
Q RS c dx c dx c c
c
α α
β β α
α τ τ τ
Δ Δ
−∞ −∞
= = =
∫ ∫
(5.2)
This expression can provide the free energy of the reaction by
( ) ( ) ( ) G RS PS G PS G RS
α α α
Δ → = Δ − Δ (5.3)
We can also use the ( ) g x Δ function to evaluate the rate constant in our
multidimensional surface by (for details see e.g. ref. (Warshel 1991))

( ) { }
( ) { }
{ }
‡
‡ ‡ ‡
‡ ‡
1
/ exp / exp
2
1
/ exp
2
x
TS
TS
k x x g x x g x dx
x x G

α α α α
α
κ β β
κ β
−∞
 
 
= Δ −Δ Δ −Δ  
 
 
 
 
 
= Δ −Δ
 
 
∫

(5.4)

72
where κ is the transmission factor, x is the velocity along the reaction
coordinate and
‡
x Δ is the width of transition state (TS) region.
‡
G
α
Δ is defined by
Eq. (5.4) and can be expressed as
( ) ( ) ( ) ( ) ( ) { }
‡
‡ ‡ 0 1 0
‡
1
ln exp -
x
RS RS
G g x g x g x g x dx
x
α α α α α α
β β
−
−∞
 
Δ = Δ − Δ + Δ − Δ  
Δ
 
 
∫
(5.5)
Where
0
RS
x is the lowest point at the RS. The correction by the last term in the
equation provides a convenient way of incorporating the entropic effects of the
ground state into the transition state theory rate constant (Warshel 1991). This
correction is in fact the origin of our
‡ ‡
w cage
g g Δ − Δ term (see e.g. ref. (Strajbl,
Florian et al. 2000; Shurki, Štrajbl et al. 2002)).
In order to evaluate the ( )
/ QM MM
G RS Δ and
( )
‡
/ QM MM QM
g x Δ terms of Eq.(5.5),
we exploit the availability of the EVB reference surface and write:
( ) ( ) ( )
b a a b
G RS G RS G RS
→
Δ = Δ + ΔΔ (5.6)
and

( ) ( ) ( )
‡ ‡ ‡
b b a b a b b
g x g x g x
→
Δ = Δ + ΔΔ (5.7)
The ( ) G RS ΔΔ term is evaluated by the LRA end point approach of ref.
(Strajbl, Hong et al. 2002) using

( ) ( ) ( ) ( )
( ) { } ( )
1
1
ln /
1
ln exp
2
a b b a
b a b a b a
a b
a
G RS Q RS Q RS
E E E E E E
β
β β
−
→
−
ΔΔ = − =
 
− − − ≅ − + −
 
(5.8)
73
where we approximated the free energy perturbation treatment by the LRA
expression. It is useful to note that the exponential average ( ) { }
exp
b a
a
E E β − −
in Eq. (5.8) is just a formal expression that is usually evaluated by a free energy
perturbation treatment.
The
( )
‡
g x ΔΔ is evaluated by the same LRA approach, but now considering
the partial partition function at
‡
x . This gives

( ) ( ) ( ) ( ) ( ) ( ) ( ) { }
( ) ( ) ( ) ( )
( )
‡ 1 ‡ ‡ 1 ‡ ‡
‡ ‡ ‡ ‡
ln / ln exp
1
2
a b b b b a b b b a b
a
b b a b b b a b
a b
g x q x q x E x E x
E x E x E x E x
β β β
− −
→
 
ΔΔ = − = − − −
 
 
≅ − + −
(5.9)
Up to this point in the derivation the main advance is the idea of fixing the
system at
‡
x . However even with this seemingly simple idea we have a
significant problem, since the transition state for the EVB ( a ) and QM/MM (b )
free energy surfaces can be different. In order to resolve thus problem we start by
modifying Eq. (5.9) by using

( ) ( ) ( ) ( )
‡ ‡ ‡ ‡ ‡
b b a a a b a b a b
g x g x g x g x x
→
Δ = Δ + ΔΔ + ΔΔ → (5.10)
The last term can be evaluated by calculating the Potential of Mean Force
(PMF) for moving on the QM/MM surface from
‡
EVB
x to
‡
/ QM MM
x . The PMF
calculations can be performed by defining two sets of constraints (one for
‡
/ QM MM
x
74
and the other for
‡
EVB
x ) and pulling the system from one constraint to another on
the QM/MM surface. Furthermore, we are able to use the LRA formulation and to
derive an expression that evaluates the last two terms of Eq. (5.10) in one step, in
which case the reaction coordinates,
EVB
x and
/ QM MM
x can also be different. In
this case we write:

( ) ( ) ( )
‡ ‡ ‡ ‡
b b a a a b a b
g x g x g x x
→
Δ = Δ + ΔΔ → (5.11)
This can be done by including the coordinate
‡
x in the LRA treatment and
writing:
( ) ( ) ( ) ( ) ( ) ( ) ( )
‡ ‡ ‡ ‡ ‡ ‡ ‡ ‡
1
2
b b a a a b a b b a a a b b a b
a b
g x g x g x x E x E x E x E x
→
 
Δ − Δ = ΔΔ → ≅ − + −
 
(5.12)
This new expression will be used in the present work in the evaluation of
( )
‡
/ / QM MM QM MM
g x Δ .

Figure 13. A schematic description of the intersection of
EVB
E and
/ QM MM
E at the TS
region. The figure can be used as a qualitative illustration (see text) for the validity of the
LRA treatment of Eq. (5.9) and (5.12)

75
In order to clarify some of the elements of our new LRA concept, we illustrate
in Figure 13 a typical behavior of
EVB
E and
QM/MM
E at the transition state (TS)
region. The figure considers these two surfaces in the simplifying case when the
reaction coordinate x is similar for both surfaces (this is not a formal requirement
in our treatment) and the same set of coordinates, τ, which span the space
perpendicular to x. Our task is to find the free energy change associated with the
move from
‡ ‡
EVB EVB EVB
( , ) E x τ to
‡ ‡
QM/MM QM/MM QM/MM
( , ) E x τ . Now, the motion
along the τ axis can clearly be described by the LRA approximation since it
corresponds to two parabolas leading to Eg. (5.9). Similarly the motion along the
x axis can be described by the LRA since it corresponds to the intersection of two
inverted parabolas. The overall application of Eq. (5.12) corresponds to the
motion from
‡ ‡
EVB EVB
( , ) x τ to
‡ ‡
QM/MM QM/MM
( , ) x τ and is also described by the
LRA approach since it corresponds to transitions between two surfaces that can be
approximated by
2
EVB
( ) ( )
i i i
i
E Q Q δ λ δ ≈
∑
and
' 2
QM/MM
( ) ( )
i i i
i
E Q Q δ λ δ ≈ + Δ
∑

respectively, where only one of the λ−coefficients is negative in each case (
i
Q δ
are the displacements from the TS after principal axis transformation).
It is useful to point out at this point that the calculations of
( )
‡ ‡
/ EVB QM/MM EVB QM MM
g x x
→
ΔΔ → allow the solvent to relax at the TS of the solute
coordinate. Such an approach misses the so called non equilibrium solvation
76
(NES) effect (see discussion in ref. (Villa and Warshel 2001)). Fortunately the
( )
‡
EVB
EVB
g x Δ is evaluated by the EVB mapping procedure(Warshel 1991) that
captures the NES free energy. Assuming that the NES contribution is similar on
the QM/MM and EVB surfaces, we can treat
( )
‡
/
/
QM MM
QM MM
g x Δ as the correct
activation free energy that contains the NES effects.

Figure 14. Showing the overall structure of the DhlA with the active site Asp124 and a bound
DCE.

In order to validate our approach we took the first step in the reaction of
haloalkane dehalogenase (DhlA) as a benchmark (see references in ref. (Olsson
and Warshel 2004)) using the crystal structure of the enzyme (Verschueren,
Franken et al. 1993) as a starting point (Protein Data Bank entry code 2DHD).
This system catalyzes the cleavage of the halogen from the dichloro-ethane
77
(DCE) and similar alkyl-halides by an S
N
2 nucleophilic attack of the oxygen of
Asp 124 on the carbon of the DCE (see Figure 14). In the second step of the
enzymatic reaction the formed ester is hydrolized by a water molecule resulting
an alcohol as a final product. These reaction steps are illustrated in Figure 15.

Figure 15. A schematic description of the reaction steps in DhlA (4A) and a stick diagram
describing the
N
S 2 step for the reaction in water with an illustration of the link atom (4B).

The S
N
2 step of the reaction of the DhlA has been subject to many recent
theoretical and experimental studies (Pries, Kingma et al. 1995; Lau, Kahn et al.
2000; Kmunicek, Luengo et al. 2001; Lewandowicz, Rudzinski et al. 2001;
Bohac, Nagata et al. 2002; Bruice 2002; Shurki, Štrajbl et al. 2002; Bosma,
Pikkemaat et al. 2003; Devi-Kesavan and Gao 2003; Devi-Kesavan, Garcia-
Viloca et al. 2003; Paneth 2003; Sun, Kalju et al. 2003; Nam, Prat-Resina et al.
2004; Olsson and Warshel 2004; Kmunicek, Hynkova et al. 2005; Marti, Moliner
et al. 2005; Soriano, Silla et al. 2005) and of course the S
N
2 reaction has also been
78
a central subject in physical organic chemistry (e.g. ref. (Shaik, Schlegel et al.
1992)). The EVB potential surface has been described in details in previous works
(Olsson and Warshel 2004). The QM(ai)/MM surface was evaluated as described
in details in (Klahn, Braun-Sand et al. 2005). Our treatment used a standard link
atom technique rather than the more reliable hybrid orbital type treatments
because we found it very convenient to link our MM program MOLARIS (Lee,
Chu et al. 1993; Chu, Villa et al. 2004) to standard QM programs (in this case
GAUSSIAN 98 (Frisch, Trucks et al. 1998)) through an automated external script.
The QM system consisted of the dichloro-ethane and the formate-group, treating a
total of 12 atoms quantum mechanically. The hydrogen atom of the formate ion is
treated as a link atom to the carbons of the Asp 124 residue and the methyl group
for the reactions in the protein and in the reference solution reaction, respectively.
The quantum system was treated by the B3LYP/DFT method with the 6-31G*
basis set.
The specific EVB parameters used in this work constitute a small
modification of the parameters that were used in ref (Olsson and Warshel 2004).
The EVB-MD simulations were performed using the MOLARIS program with a
simulation sphere of 22 Å described by the ENZYMIX force field and subject to
the Surface Constraint All Atom Solvent (SCAAS)(Warshel, Sussman et al. 1988)
spherical boundary conditions and the local reaction field (LRF) long range
treatment(Lee and Warshel 1992). The EVB free energy profile of the reacting
79
systems (water and protein) was evaluated by performing the EVB FEP/umbrella
sampling (US) procedure with 31 windows of 10 ps length each and a time step of
1 fs. The simulations were performed at 300K.

5.3 Results and Discussion

Figure 16. Showing the SCAAS simulation system for the reference reaction in water.
As stated in the previous section we took the first step in the reaction of DhlA
as our benchmark system. The calculations started with the evaluation of the EVB
free energy surface for the reaction in DhlA and in the reference solution system.

80
Table 2. The EVB energetics of the S
N
2 reaction in water and in haloalkane dehalogenase for
simulations with five different initial conditions.
a

Water Protein
EVB
mapping
( ) G RS TS ΔΔ →

( ) G RS PS ΔΔ →

( ) G RS TS ΔΔ →

( ) G RS PS ΔΔ →

Sim. 1. 26.97 -1.11 19.82 4.47
Sim. 2. 28.32 -6.19 19.53 3.67
Sim. 3. 29.48 -5.43 19.98 4.52
Sim. 4. 28.86 -5.27 20.37 5.87
Sim. 5. 28.09 -0.47 19.39 3.75
Average 28.34 -3.69 19.82 4.46
a
The energies are given in kcal/mol.
The final free energy profiles were then obtained by averaging the results of
five mapping calculations from different initial conditions, each obtained after
100 ps relaxation at the reactant state. The simulation system for the reference
reaction in water is depicted in Figure 16. The results of the simulations are
summarized in Table 2. The corresponding EVB surfaces are described in Figure
17.
81

Figure 17. The EVB free energy profile of the
N
S 2 reaction in water and in the haloalkane
dehalogenase enzyme. The profiles were obtained using FEP/US for the EVB ground state
energy.
In the next step we perform the LRA calculations of the
b a
E E − terms of
Eq. (5.8) (where a and b designate the EVB and QM/MM surfaces,
respectively).

Figure 18. Illustrating the fluctuations of
/ QM MM EVB
E E − for runs on the transition state of
the EVB surface at
‡
EVB
x . The figure depicts the results of 5 runs, each of 10 ps.

82
The fluctuations involved in these simulations are illustrated in Figure 18, the
corresponding distributions are shown in Figure 19 and the resulting averages are
summarized in Table 3. These results were then used to convert the EVB reaction
energy to the QM/MM free energy using Eq. (5.6). This treatment gave
/
( )
QM MM
G RS PS Δ → of -10.42 and -2.25 kcal/mol for the reaction in the protein
and water, respectively.
Table 3. The LRA estimate of the free energy of transfer from the EVB to the QM/MM
surface in the RS and PS for the reaction of DhlA.
a

Water Protein
/
( )
QM MM EVB
EVB RS
E E −
5.49 6.81
/
( )
QM MM EVB
EVB PS
E E −
5.30 2.36
/
/ ( )
QM MM EVB
QM MM RS
E E −
-16.60 -19.27
/
/ ( )
QM MM EVB
QM MM PS
E E −
-29.88 -28.25
( )
/ EVB QM MM
G PS
→
ΔΔ
-6.73 -6.71
a
The energies are given in kcal/mol. G ΔΔ (PS) is given
relative to G ΔΔ (RS) and the results reflect the average of two
simulation runs.
83

Figure 19. Showing the distributions of
/ QM MM EVB
E E − for trajectories on the EVB and
QM/MM reactant state.

Finally we turn to the QM/MM TS for the reaction in water and in the protein
active site (the structure of the TS in the protein active site is described in Figure
20).

Figure 20. Showing the active site of the DhlA for the TS of the
N
S 2 reaction step.
84
Table 4. Results of the LRA simulation for the transition state of the solution reaction and
the overall QM/MM activation barrier.
a

Water TS
simulations
( ) ( )
‡ ‡
b a a a
b
E x E x −
( ) ( )
‡ ‡
b a a a
a
E x E x −
( ) ( )
‡ ‡
b b a b
b
E x E x −
Sim. 1. -15.70 3.52 -17.64
Sim. 2. -18.24 4.01 -17.83
Sim. 3. -16.98 5.77 -18.33
Sim. 4. -11.58 6.42 -17.32
Sim. 5. -13.34 6.17 -17.43
Average -15.17 5.18 -17.71
Average with
RS as a
reference
1.43 -0.31 -1.11

( )
‡
a b a
g x
→
ΔΔ =
(1.43+ -0.31) / 2 = 0.56
( )
‡ ‡
a b a b
g x x
→
ΔΔ → =
(-0.31 + -1.11) / 2 = -0.71
Overall
‡
Δ G Δ G Δ G Δ G
28.34+0.56 = 28.90 (obs.:
27.0)
28.34 -0.71 = 27.63 (obs.:
27.0)
a
The energies are given in kcal/mol. The subscripts a and b designate, respectively the EVB
and QM/MM potential surfaces. The lower part of the table provides two estimates for
‡
G Δ ;
the first involves the use of the first two terms of Eq. (5.10) assuming that
‡ ‡
/ EVB QM MM
x x ≅ .
The second one involves the use of Eq. (5.11) and (5.12), while taking
‡
/ QM MM
x from the
corresponding gas phase estimates. The observed activation barrier is taken from the estimate
of ref. (Olsson and Warshel 2004).
Although our approach allows us to evaluate the ab initio barrier by using the
EVB barrier, we face the challenge of finding configurations that constitute
‡
QM/MM
x without actually calculating the entire
/
( )
QM MM
g x Δ . Our way to overcome
the challenge is to exploit the fact that in many cases (e.g. S
N
2 and related
reactions)
‡ ‡ ‡
p w gas
x x x ≈ ≈ (see e.g. (Maulitz, Lightstone et al. 1997)). Taking
‡
EVB
x
to be as our approximation for the TS of both the QM and the EVB potentials we
85
obtained the results that are illustrated in Figure 21. Obtaining the QM/MM
activation free energy by moving from the EVB to the QM/MM surfaces. and
summarized in the first two columns of Table 4 and Table 5, respectively, for the
solution and protein reactions.

Figure 21. Obtaining the QM/MM activation free energy by moving from the EVB to the
QM/MM surfaces.
Now, if our assessment that the TSs are similar for the QM/MM and EVB
potentials is correct, there is no need to evaluate the last two terms in Eq. (5.10).
With this in mind we obtained a first rough estimate of the effect of moving to
‡
/ QM MM
x , by evaluating the average force along the reaction coordinate at
‡
EVB
x for
trajectories on the EVB and on the QM/MM surface. In both cases we found the
average forces similar indicating that the above approximation is reasonable in the
present case. The second and more rigorous estimate of the QM/MM estimate
involves the use of Eqs. (5.11) and (5.12) for a direct evaluation of the QM/MM
free energy at the gas phase estimate of
‡
/ QM MM
x . This gave
‡
/ , . QM MM enz
G Δ =18.6
86
kcal/mol and
‡
/ , . QM MM wat
G Δ =27.6 kcal/mol. The corresponding calculated catalytic
effect is
‡ ‡ ‡
. . wat enz
G G G ΔΔ = Δ − Δ ≅ 9 kcal/mol, whereas the estimate based on the
first two terms of Eq. (5.10) gives
‡
/
11
QM MM
G ΔΔ ≅ kcal/mol. Thus our results are
in good agreement with the corresponding observed value (
‡
. obs
G ΔΔ =11.7
kcal/mol).
Table 5. Results of LRA simulations at the transition state of the enzyme reaction and the
overall QM/MM activation barrier.
a

Protein TS
simulations
( ) ( )
‡ ‡
b a a a
b
E x E x −
( ) ( )
‡ ‡
b a a a
a
E x E x −
( ) ( )
‡ ‡
b b a b
b
E x E x −
Sim. 1. -16.51 -2.18 -15.09
Sim. 2. -16.59 -0.85 -14.00
Sim. 3. -17.17 -2.74 -15.08
Sim. 4. -16.45 0.43 -15.19
Sim. 5. -16.14 5.97 -15.10
Average -16.61 0.13 -14.89
Average
with RS as a
ref. point
2.70 -6.68 4.38
Apply Eq.
(5.9) or Eq.
(5.12)
( )
‡
a b a
g x
→
ΔΔ =
(2.70+-6.68) / 2 = -1.99
( )
‡ ‡
a b a b
g x x
→
ΔΔ → =
(-6.68 + 4.38) / 2 = -1.15
Overall
‡
Δ G Δ G Δ G Δ G
19.82-1.99 = 17.83 (obs.: 15.3
a
) 19.82-1.15 = 18.67 (obs.: 15.3
a
)
a
The energies are given in kcal/mol. The subscripts a and b designate, respectively the EVB and
QM/MM potential surfaces. The lower part of the table provides two estimates for
‡
G Δ ; the first
involves the use of the first two terms of Eq. (5.10) assuming that
‡ ‡
/ EVB QM MM
x x ≅ . The second
one involves the use of Eq. (5.11) and (5.12), while taking
‡
/ QM MM
x from the corresponding gas
phase estimates. The observed activation barrier is taken from the estimate of ref. (Olsson and
Warshel 2004).

87
The present study focuses on the evaluation of the free energy of a given ab
initio QM/MM approach without addressing the reliability of the specific
QM/MM potential surface. Since we used here the DFT B3LYP method with a 6-
31G* basis set, it is interesting to ask how accurate the results of this method are
assuming that we obtain its correct free energy surface. In order to explore this
issue we examine the gas phase ab initio activation free energies of the S
N
2
reaction with formate and acetate as nucleophiles and with different ab initio
models (see Table 6) .We found that using larger basis sets, including diffuse
functions does not change the results significantly, so the choice of 6-31G* is
very reasonable. The DFT/B3LYP method seems to underestimate the transition
state barrier by about 8 kcal/mol in the gas phase for the formate system, relative
to the MP2/6-311G** results. However, the QM/MM calculations correspond to
solution calculations and the difference between the B3LYP/6-31G* and MP2/6-
311G** calculations of formate in a Langevin Dipoles (LD) solvent model
(Florian and Warshel 1997) was found to be about 4 kcal/mol. Furthermore, the
gas phase estimate obtained by CCSD(T)/6-31G* is only about 23 kcal/mol for
the formate system and the G2 results for acetate are only 21.3 kcal/mol. Thus we
estimate that the error associated with the quantum method used and with the use
of formate instead of acetate is about 2-3 kcal/mol. We do not expect that a
smaller error range can be obtained without addressing such issues as the effect of
88
charge transfer to the solvent and the effect of using polarizable force field (that
was not used in this study).
Table 6. Gas-phase ab initio activation energies for the attack of a carboxylate on DCE.
a

B3LYP
6-31G*
B3LYP
6-611++G**
MP2
6-31G*
MP2
6-311++G**
G2 (Devi-
Kesavan,
Garcia-
Viloca et al.
2003)
Formate
+DCE
b

17.04
(+1.86)
17.22
(+2.43)
24.70
(-1.65)
25.36
(-2.26)
-
Acetate
+DCE
18.13 18.35 25.88 26.18 21.3
a
The energies are given in kcal/mol.
b
The difference between the solvation free energies of the transition state and the reactant state are
given in parenthesis. These solvation free energies were evaluated using the LD model(Florian and
Warshel 1997).
With this in mind we believe that it is still important to calibrate or at least to
validate the QM/MM results using experimental information about the reference
solution reaction. However, the difference between the barrier in the enzyme and
in the reference solution reaction should be sufficiently reliable to let one reach
clear conclusions about the catalytic effect of the enzyme.

5.4 Conclusions
The evaluation of the free energy changes in proteins by classical force fields
requires usually a very extensive averaging over the configurational space of the
protein. Thus it is quite likely that the same requirement will hold for QM/MM
89
calculations. This means that simple minimization approaches of the type used in
gas phase QM calculations might not be effective in evaluations of activation
energies of chemical reactions in proteins (e.g. see ref (Klahn, Braun-Sand et al.
2005)). Unfortunately, evaluating the activation free energies for QM(ai)/MM
surfaces is extremely challenging due to the need of extensive evaluation of the
QM(ai) energies.
Here we present a refinement of our previous treatments(Strajbl, Hong et al.
2002), that allows one to evaluate the QM(ai)/MM activation free energies of
reactions in condensed phases. Our treatment uses, as before, the EVB as a
reference potential and then uses the LRA approach to evaluate the free energies
of transfer from the EVB to the QM/MM surfaces in the reactant and product
state. However, at the transition state we constraint the system to be at a fixed
value of the reaction coordinate, using the EVB and QM/MM partial partition
functions (the q ‘s of Eqs. (5.1) and (5.9)) to evaluate the QM/MM free energy at
the transition state.
The new version allows us to evaluate the QM(ai)/MM activation free
energies for reactions in proteins and solutions in an effective and economical
way .This advance is particularly significant since it evaluates the free energy
associated with the both the substrate and the solvent motions. The evaluation of
the free energy associate with the solvent (or protein) configurational space
papered to be a relatively simple task once one uses a classical reference potential.
90
The main problem has been to use the reference potential in evaluation of the free
energy contributions associated with the solute motions where the difference
between the reference EVB potential and the QM/MM potential can be large. The
present refinement finally allows us to overcome the problems with the solute
fluctuations and thus to present a general way for evaluation QM(ai)/MM
activation free energies for reactions in proteins and solutions.
It might be useful to estimate at this stage the time saving offered by the
present approach, relative to a full QM/MM PMF treatment. Here we start by the
observation that the EVB umbrella-sampling procedure, which is comparable in
terms of the time steps needed to standard PMF approaches, requires at least 30
windows with around 10000 configurations in each (assuming 10ps long
simulations with the usual 1fs time step). Since each QM/MM step costs (e.g. in
our case with B3LYP/6-31G*) around 2 min on one 64-bit Intel Xeon 3.2 GHz
processor, we would need
4
10 hours for an ab initio PMF calculation. This can be
compared to a total of about 100 hours (involving 5000 single point calculations)
that are required in our LRA approach that samples only the reactant, product and
the transition state (in contrast to a sampling along the entire reaction path). Of
course, our approach can also benefit greatly from the similarity between the EVB
and ab initio surface. Thus, by requiring that the EVB potential will be similar to
the corresponding ab initial we can reduce the simulation length required for
convergence.
91
Despite of the advance of the present QM(ai)/MM treatment, it still involves
an extensive and complicated procedure. Thus it is encouraging to see that the
catalytic effect obtained by the EVB method is quite reliable. Thus we believe
that many issues that require extensive configurational averaging (e.g. the effect
of mutations and the relationship between folding energy and catalysis) should
still be explored using the EVB method. In our view the QM/MM-FEP treatment
should be particularly important in cases where the exact mechanism in the
enzyme presents an open question, and also as a tool for refining and calibrating
EVB studies of mutational effects.

Chapter 6. Linear Free Energy
Relationships

6.1 Introduction
LFERs are one of the most important empirical relationships in organic
chemistry; they are on some level a unifying theory of organic reactivity. LFERs
are used to predict equilibrium or rate constants for specific reactions by
substituent effect. Most typically the Hammett equation is used for substituent
effects in the following form:

log
X
X
H
k
k
ρσ =
(6.1)

92
where the parameter
X
σ is the substituent constant and rho is a reaction
specific constant.
Although these concepts originate from chemistry, they are now widely used
as important empirical descriptors in biological systems as well, in the so called
quantitative structure activity relationships (QSARs). The most important
application of the QSARs is currently in pharmaceutical drug design: the whole
cheminformatics research area evolved from the simple concepts of LFERs to the
complex use of QSARs.
Despite the importance of the LFERs, there is no rigorous derivation of these
equations, only few attempts were made based on quantum mechanical
ideas(Mezey 1987; Cioslowski and Fleischmann 1991; Girones, Carbo-Dorca et
al. 2003).
A very promising possibility towards understanding the basic nature of the
LFERs comes from the transition state theory of electron transfer (ET) reactions.
Here the rate constant can be related to the reaction free energy,
0
G Δ , the so
called driving force by the Arrhenius relationship:

( )
2
0
exp
4
B
G
k A
k T
λ
λ
 
− Δ +
 
=
 
 
 
(6.2)

Here
B
k is the Boltzmann constant, k is the rate constant, λ is the
reorganization energy and A is a prefactor that depends for example on the
93
frequency of crossing the barrier top. Considering λ and A constant, we obtain a
relationship between the driving force and the rate:
( )
2
0
log const. k G λ =− Δ + +
(6.3)

which implies a quadratic dependence on the reaction free energy that is
equivalent to for example pK
a
of the reaction. One can then argue that for a small
change in the driving force the parabola can be very well estimated by a constant
slope linear equation when
0
G λ Δ . Similar derivation with the same
conclusion can be found with more details in (Warshel 1991).
The validity of LFERs for chemical processes in solutions has been the
subject of many studies (e.g. (Leffler and Grundwald 1963; Marcus 1964; Albery
1980; Kreevoy and Kotchevar 1990; Warshel 1991; Schweins, Geyer et al.
1996)). Furthermore experimental(Fersht, Leatherbarrow et al. 1986; Toney and
Kirsch 1989; Mihai, Kravchuk et al. 2003) and theoretical (Warshel 1984;
Warshel, Schweins et al. 1994) studies have indicated that such relationships are
also valid in enzymes and more importantly, for the purpose of this work, with a
similar correlation coefficients in the enzyme and in the reference solution
reactions (refs. (Warshel 1984; Warshel, Schweins et al. 1994).
Related to the validity of the LFERs, i.e. that the LFERs can be used with the
same correlation coefficient in different environments; we would like to examine
the validity of the assumption that the H
ij
s are transferable. Based on the
experimental evidences on the LFERs, there are good qualitative reasons to
94
assume that the off diagonal elements are transferable between different phases,
but quantitative examinations of this issue are clearly long overdue.
Experimental LFER studies have been used to try to predict reaction
mechanism in solutions and in enzymes(Jencks 1987; Warshel 1991; Warshel,
Hwang et al. 1992; Warshel, Schweins et al. 1994; Aqvist, Kolmodin et al. 1999).
Some of the most important groups of reactions in these studies were the
phosphate ester hydrolysis reactions. The nature of the hydrolysis of phosphate
monoester dianions in solutions and in proteins is still a problem of a significant
current interest. Section 6.3, based on the publication (Wagner, DeMarco et al.)
explores this problem by systematic calculations of the landscape of the potential
surfaces of the reactions of a series of phosphate monoesters with different
leaving groups. These calculations involve computational studies ranging from ab
initio calculations with implicit solvent models to ab initio QM/MM free energy
calculations. The calculations reproduce the observed linear free energy
relationship (LFER) for the solution reaction and thus are consistent with the
overall experimental trend and can be used to explore the nature of the transition
state (TS) region, which is not accessible to direct experimental studies. It is
found that the potential surface for the associative and dissociative paths is very
flat and that the relative height of the associative and dissociative TS is different
in different systems. In general, the character of the TS changes from associative
to dissociative upon decrease in the pK
a
of the leaving group. It is also
95
demonstrated that traditional experimental markers such as isotope effects and the
LFER slope cannot be used in a conclusive way to distinguish between the two
classes of transition states. In addition it is found that the effective charges of the
TS do not follow the previously assumed simple rule. Armed with that experience
we explore the free energy surface for the GTPase reaction of the RasGAP
system. In this case it is found that the surface is flat but that the lowest TS is
associative. The present study indicates that the nature of the potential surfaces
for the phosphoryl transfer reactions in solution and proteins is quite complicated
and cannot be determined in a conclusive way without the use of careful
theoretical studies that should of course, reproduce the available experimental
information.

6.2 Transferability of the off diagonal elements
between different phases
The present Section describes an attempt to quantify the above hypothesis of
the transferability of the off diagonal coupling elements between different phases.
We face the problem here that the transferability of the off diagonal terms cannot
probably be proven in a rigorous way, for which one of the main reasons is that
the diabatic states are not uniquely defined. Furthermore, we also realize that
there is no current way of rigorously formulating the solvent dependence of the
96
off diagonal term, despite some attempts to progress in this direction (Schmitt and
Voth 1998). Fortunately, however, it is clearly possible to try to validate the EVB
assumptions by using a numerical approach, and such a strategy will be adopted
here.
In trying to establish the validity of the EVB assumptions numerically, we
will have to choose a specific zero order diabatic representation. This can be done
by using different quantum mechanical approaches, but doing so by the EVB and
the frozen DFT (FDFT) or the constraint DFT (CDFT) methods (Wesolowski and
Warshel 1993; Hong, Strajbl et al. 2000) is probably the most natural way (and
currently perhaps the most rigorous (Wesolowski 2002)) of representing
interacting fragments with a full flexibility of allowing or restricting charge
transfer between these fragments. Thus these approaches provide an effective way
for examining the validity of the key EVB assumption (although alternative ways
(Mo and Gao 2000) are also useful). That is, there are of course different ways for
obtaining H
12
theoretically (Shaik and Shurki 1999) (Mo and Gao 2000) and
estimating the effective H
12
experimentally (Hush 1961) (Hush 1968) (Gwaltney,
Rosokha et al. 2003). However, the EVB intentionally avoids the rigorous VB
formulation, since such a treatment suffers from major complications due to the
overlap between the diabatic wave functions. Instead the EVB chose arbitrary
diabatic states and then force them (by using the proper H
12
) to reproduce
"reality", namely the adiabatic surface and charges. Of course, the selection of the
97
diabatic states is not unique and rather arbitrary (except the requirement that the
separate fragments at infinity will be reproduced by this diabatic states).
However, taking the energies and charges of these states as those of the
corresponding FDFT states provides a well defined diabatic selection, without the
need to specify a wave function and the resulting formal conflicts. Basically we
have a formulation that can be examined rigorously within its framework. Using
this formulation we will try to show that with this selection H
12
is transferable and
the basic EVB assumption is, in fact valid.

Our starting point is the generic reaction:
A B C A B C + − → − + (6.4)
Where we can describe the reactant (r) and product (p) states by diabatic
valence bond like wave-functions:
[ ] [ ]
[ ] [ ]
r
p
A B C
A B C
Ψ = Ψ Ψ −
Ψ = Ψ − Ψ
(6.5)
We can also consider the adiabatic ground state wave-function, which
describes the entire system. Comparing the diabatic and adiabatic energies of our
system will allow us to examine the validity of the EVB approach. This point will
be clarified below by considering briefly the EVB and CDFT/QM/MM methods.
98
The EVB approach describes reactions by using diabatic surfaces and mixing
them by effective off diagonal terms. In the simple 2-state case, for e.g. i = 1 or 2,
we can write:
gas intra Ss
i i i
ii i
H U U ε α = = + + (6.6)
0
exp{ ( )}
ij
H A R R μ = − − (6.7)
where
intra
i
U is the intramolecular potential of the solute system, which is
represented by a molecular mechanics-type force field that describes the current
bonding structure of the state i.
Ss
i
U represents the interaction energies (van der
Waals and electrostatic interactions) between the solute (S) and the environment
(s). Finally
gas
i
α is a constant that reflects the difference between the gas phase
energy of the diabatic states when all the corresponding fragments are at infinite
separation(Warshel 1991). The off diagonal element H
12
depends typically on the
distance R between the A and C atoms. The environment is sometimes also
modeled by a molecular mechanics force field, but the parameters of this potential
energy are common for all states.
The ground state energy,
g
E of the system is obtained as the lowest
eigenvalue of the above described by a 2x2 EVB Hamiltonian.
1 12
12 2
0
g
g
E H
H E
ε
ε
−
=
−
(6.8)
99
Now the EVB approach assumes that H
12
does not change significantly upon
moving the solute from the gas phase to solution and our task is to examine this
approximation in a rigorous way. This can be done by using the CDFT approach
for generating a diabatic representation (see below) and then extracting the
relevant H
ij
by comparing the diabatic surface to the corresponding adiabatic
surface obtained by the full DFT method. Performing such calculations in the gas
phase and in solution should allow us to extract the solvent effect on H
ij
and thus
to validate our assumption. More specifically, since the selection of the diabatic
states is always a mathematical convenience rather than a reflection of any
absolute reality, we may choose
1
ε and
2
ε that provide the most practical utility.
However, here we will focus on diabatic states that can be obtained from quantum
mechanical calculations (see below). In addition we will also evaluate the actual
adiabatic ground state energy, E
g
, by the same Hamiltonian used to determine the
diabatic states. The 2-state EVB model leads through Eq. (6.8) to:
( )( )
12 1 2 g g
H E E ε ε = − − . (6.9)
Thus, having quantum mechanical values for
1
ε and
2
ε as well as E
g
provides
an operational definition of H
12
.

With the operational definition of Eq. (6.9), we can look for an ab initio
treatment that can produce wave-functions both in the diabatic and adiabatic
states where the relationship between these states is defined the same way in the
100
gas phase and in solution. One of the most unique ways of doing so is provided by
the FDFT/CFDF method (Hong, Strajbl et al. 2000), and since more details on the
frozen density functional method can be found elsewhere, we will provide below
only a brief description.

The FDFT approach can be used to define the diabatic energies of the
reactant state (H
11
) and product state (H
22
) of our system. This can be done by
decomposing the density as either ( ) ( ) ( )
1
A BC
N N
r r r ρ ρ ρ = + (for one diabatic
state) or ( ) ( ) ( )
2
AB C
N N
r r r ρ ρ ρ = + (for the other diabatic state).
Since we are interested here in the validity of the EVB in studies of reactions
in solutions, we will have to evaluate the relevant activation barriers by both the
EVB and the FDFT approach. The general procedure was described in Chapter
3.3.
The actual adiabatic ground state free energy profile along the reaction
coordinate,
1 2 22 11
x H H ε ε = − = − , is evaluated by the umbrella sampling / FEP
method:
( ) ( ) ( ) ( ) ( )
1
' ln ' exp
m
m g m
g x G x x E x x
ε
β δ β ε
−
 
Δ = Δ − − − −
 
(6.10)
Similar formulas can be used to evaluate the diabatic free energies associated
with H
11
and H
22
.
101
In the present case we consider, as a general benchmark for charge transfer
reactions, the S
N
2 reaction:

3 3 A B A B
Cl CH Cl Cl CH Cl
− −
+ → + . (6.11)
The VB resonance structures for the diabatic reactant (r) and product (p) states
are, respectively:
3
3
[ ]
[ ]
r A B
p A B
Cl H C Cl
Cl CH Cl
−
−
Ψ = −
Ψ = −
(6.12)
In order to represent the VB diabatic reactant state by the FDFT approach, we
divide the whole system into the
A
Cl
−
ion and the
3 B
CH Cl subsystems, and
express the total energy in the gas phase as :
3
3 3
3 3
3 3
11 1
[ , ]
[ ] [ ] [ , ]
( ( ) ( ))( ( ') ( '))
1
'
2 | ' |
( )( ( ) ( )) [ ( ) ( ))]
B
A
B B
A A
B B
A A
B B
A A
FDFT
CH Cl
Cl
nadd
CH Cl s CH Cl
Cl Cl
CH Cl CH Cl
Cl Cl
ext CH Cl xc CH Cl
Cl Cl
H E
T T T
d d
V d E
ε ρ ρ
ρ ρ ρ ρ
ρ ρ ρ ρ
ρ ρ ρ ρ
−
− −
− −
− −
= =
= + +
+ +
+ +
−
+ + +
∫∫
∫
r r r r
r r
r r
r r r r r r
(6.13)

The modified Kohn-Sham equation was applied to obtain the densities of the
subsystems. Similarly, we described the product state by:
102

3
3 3
3 3
3 3
22 2
[ , ]
[ ] [ ] [ , ]
( ( ) ( ))( ( ') ( '))
1
'
2 | ' |
( )( ( ) ( )) [ ( ) ( ))]
A
B
A A
B B
A A
B B
A A
B B
FDFT
CH Cl
Cl
nadd
CH Cl s CH Cl
Cl Cl
CH Cl CH Cl
Cl Cl
ext CH Cl xc CH Cl
Cl Cl
H E
T T T
d d
V d E
ε ρ ρ
ρ ρ ρ ρ
ρ ρ ρ ρ
ρ ρ ρ ρ
−
− −
− −
− −
= =
= + +
+ +
+ +
−
+ + +
∫∫
∫
r r r r
r r
r r
r r r r r r
(6.14)
Finally we obtain the adiabatic ground state energy
full dft
g
E
−
, by treating the
whole system at the same level with a regular DFT calculation.
The same treatment is also implemented in solution, where the solvent effect
is simulated by a QM/MM approach (e.g. ref. (Warshel and Levitt 1976)), using
the formulation of ref. (Olsson, Hong et al. 2003). The actual QM/MM treatment
has been implemented in the program MOLARIS (Chu, Villa et al. 2004) (Lee,
Chu et al. 1993) by using a very convenient "weak coupling" philosophy where
any QM package can be called by MOLARIS. In this work, the QM region is
treated by the DFT and CDFT/FDFT method that were implemented in the
deMon program (St-Amant and Salahub 1990). The Becke88 exchange potential
(Becke 1988), the Perdew86 correlation potential (Perdew 1986) and the non-
additive kinetic energy functional were used in all the calculations. Standard 6-
31+G* basis sets were used for expanding molecular orbitals.
The MM part (the water molecules in the classical region) was represented by
the regular ENZYMIX force field. The simulation model includes an 18 Å
explicit simulation sphere, surrounded by a surface region whose average
polarization and radial distribution are determined by the surface constrained all-
103
atom solvent (SCAAS) model(Vaidehi, Wesolowski et al.). The surface region is
embedded in a bulk continuum region with a dielectric constant of 80. Long-range
electrostatic effects are treated by the local reaction field model.
The diabatic and adiabatic QM(FDFT)/MM free energy profiles were
obtained by using 21 windows of 2 ps each with 0.5 fs stepsize at 300
o
K, in each
window the runs were preequilibrated. The calculated profiles are summarized in
Figure 22 and show the expected trend in S
N
2 reactions (Shaik, Schlegel et al.
1992), where the adiabatic barrier is much larger in solution due to the loss of
solvation upon formation of the
1/ 2 1/ 2
3
Cl CH Cl
− −

delocalized transition state.
We also examined longer dynamics and runs that started with long MM
relaxations, and obtained similar results, thus found that the free energy
calculations converged. The adiabatic barrier for this reaction in gas phase in
known to be overestimated by about 9 kcal/mol due to the use of the DFT-BP
method (Yang, Fleurat-Lessard et al. 2004). Assuming similar effect in solution,
we expect that the barriers would be around 16+9=25 kcal/mol with a better ab
initio method. Now, since the calculated profile starts when the reactant
molecules are in the same solvent cage and thus the calculated barrier should be
increased by about 2.3 kcal/mol (e.g. see ref. (Strajbl, Sham et al. 2000)) and be
around 27.3 kcal/mol. This result is in a good agreement with the experimental
barrier of about 26.6 kcal/mol (Albery 1978).
104
It is important to note that our calculations represent a complete free energy
profile for an ab initio QM/MM simulation of a reaction in solution with a
complete flexibility of the solute and the solvent. Such studies are extremely
challenging and quite rare (see discussion in ref. (Rosta, Klahn et al. 2006)), even
gas phase free energy calculations on this relatively simple system have
significant difficulties (Yang, Fleurat-Lessard et al. 2004). Also note that our
model includes the non equilibrium solvation effects (which are automatically
obtained by Eq.(6.10)) and the solute entropic effects, which are typically missing
in PMF type calculations with fixed solute coordinates for each mapping window
(see discussion in ref. (Shurki and Warshel 2003)).
It is important to note that the gas phase barrier is evaluated for the same
structures as those used in the solution study. This means that the barrier is
evaluated more or less from the minimum of the gas-phase profile (rather than the
configuration of the separated fragments), where the corresponding gas-phase
barrier is expected to be around 4 kcal/mol in the DFT-BP method (Yang, Fleurat-
Lessard et al. 2004), as compared to our calculated barrier of about 6 kcal/mol.
This 2 kcal/mol deviation reflects the fact that the gas-phase maximum is
calculated using the structures of the TS in solution.
105

-100 -50 0 50 100 150
0
10
20
30
Energy (kcal/mol)
Δε (kcal/mol)
Ε
gas
g
ε
gas
1
ε
gas
2
Ε
sol
g
ε
sol
1
ε
sol
2

Figure 22. Diabatic and adiabatic FDFT energy profiles for the reaction,
3 3
Cl CH Cl ClCH Cl
− −
+ → + , in the gas phase and in solution. The reaction coordinate is
defined as the energy difference between the diabatic surfaces,
1 2
ε ε ε Δ = − .

Inserting the free energy curves of Figure 22 into Eq. (6.9) allows us to obtain
the effective H
12
for both the gas phase and solution reaction for any point along
the reaction coordinate between 60 ε Δ = − kcal/mol to 60 ε Δ = kcal/mol. At
larger absolute values of ε Δ the calculations do not seem to be sufficiently
reliable in view of the asymmetry of the results for positive and negative ε Δ , thus
we did not evaluate H
12
in this range. Obtaining fully symmetric results would
require significantly longer runs and sampling with different initial conditions.
The results obtained from the above analysis are summarized in Figure 23.
The most striking conclusion from Figure 23 is that H
12
is very similar in the gas
106
phase and the solution cases, despite the very large change in the corresponding
diagonal energies. This suggests that the EVB off diagonal elements are indeed
transferable between different phases.
-60 -40 -20 0 20 40 60 80
0
2
4
6
8
10
12
14
16
18
H
ij
(kcal/mol)
Δε (kcal/mol)
H
ij
gas phase
H
ij
solution

Figure 23. Plot of H
ij
of the above reaction in gas phase and solution. The data are obtained
from the diabatic and the adiabatic curves of Figure 22, using
( ) ( )( )
12 1 2 g g
H E E ε ε ε Δ = − − . Note that the adiabatic barrier is underestimated by
about 9 kcal/mol due to the use of the DFT-BP method(Yang, Fleurat-Lessard et al. 2004).

This work explores the long-standing question of the validity of the EVB's
key assumption that the off diagonal mixing terms are transferable between
different environments. This examination exploited the ability of the FDFT
approach to generate convenient and physically consistent diabatic states for
QM/MM treatments and this to be used in numerical studies of the relationship
between the diabatic and adiabatic surfaces in different environments. Obtaining
107
these surfaces allows us to evaluate the corresponding off diagonal elements in
different environments and thus to examine the effect of the environment on these
crucial parameters. It has been found that at least for the test case of S
N
2 reactions
the off diagonal element does not change significantly upon moving from gas
phase to solution. This is quite significant since the S
N
2 reaction is a prototype to
many charge transfer reactions and since the present study found that the off
diagonal element does not depend significantly on the environment for the entire
reaction coordinate and thus for a different degree of charge transfer. Of course
more studies will be required in the future, examining for example charge
separation reactions, but the present study is the first one that actually examines
the crucial problem addressed in this work. Now from the indications obtained so
far, it seems that the EVB assumption is valid and extremely useful.
The transferability of the off diagonal elements between different phases is
important not only for the quantitative EVB studies but also as a justification of
the assumptions made for a long time by the physical; organic chemistry
community in LEFR studies, where it has been assumed implicitly that the
correlation coefficients do not change significantly in different environments.
Finally, it is important to note that even if in some cases, the off diagonal term
changes upon moving from the gas phase to solution are larger than the current
estimate, the changes upon moving from solutions to enzymes are expected to be
significantly smaller since in both cases we deal with polar environments.
108

6.3 LFERs in Phosphate Hydrolysis: an
Introduction
Many biological molecules contain phosphate and the reactions that break,
connect or exchange the bonding between such groups play a major role in key
processes ranging from the replication of the genetic material to energy and signal
transduction (for reviews see for example (Cox and Ramsay 1964; Kirby and
Varvoglis 1967; Westheimer 1987; Vetter and Wittinghofer 1999; Cleland and
Hengge 2006)). One of the most important reactions in such systems is the
hydrolysis of phosphate monoesters such as GTP and ATP that play a major role
in energy and signal transduction processes. These systems have been explored
extensively in the past by different experimental and theoretical approaches that
examined factors that can change or modulate the rate relevant constant (e.g.
(Langen, Schweins et al. 1992; Ahmadian, Stege et al. 1997; Aqvist, Kolmodin et
al. 1999; Glennon, Villa et al. 2000; Strajbl, Shurki et al. 2003; Shurki, Strajbl et
al. 2004) (Schweins, Geyer et al. 1995; Schweins, Geyer et al. 1996; Vetter and
Wittinghofer 1999; Allin, Ahmadian et al. 2001; Allin and Gerwert 2001;
Dittrich, Hayashi et al. 2003; Du, Black et al. 2004; Grigorenko, Nemukhin et al.
2005; Klahn, Braun-Sand et al. 2005) (Bourne and Williams 1984; Cui, Elstner et
al. 2001; Yang, Gao et al. 2003; Li and Cui 2004)). Furthermore, significant
109
insight has been obtained from studies of the mechanism of the corresponding
hydrolysis reactions in solution (e.g. (Cox and Ramsay 1964; Kirby and Jenks
1965; Guthrie 1977; Pecoraro, Hermes et al. 1984; Barnes, Wilkie et al. 1994;
Bourne 1997; Hoff and Hengge 1998; Wang, Topol et al. 2003; Cleland and
Hengge 2006; Grigorenko, Rogov et al. 2006)). Additional insight has been
provided by studies of related systems (e.g. (Florián, Åqvist et al. 1998; Zhou and
Taira 1998; Hu and Brinck 1999; Iche-Tarrat, Barthelat et al. 2005; Liu, Lopez et
al. 2005)). Nevertheless, at present we still have a limited understanding of the
detailed reaction paths in solution and in enzyme active sites. The problem stems
to a large extent from the fact that the energies of the transition states of
phosphate hydrolysis reactions are relatively high, and it is very hard to
characterize such transition states experimentally via key intermediates. Thus,
despite frequent implications that many of the mechanistic issues have been
resolved experimentally (e.g. refs (Admiraal and Herschlag 1995; Admiraal and
Herschlag 2000)), it is extremely hard to deduce the nature of the transition state
(TS) in a unique way from the available experiments. For example, a careful
analysis (Florián, Åqvist et al. 1998; Aqvist, Kolmodin et al. 1999) has indicated
that the difference between associative and dissociative paths cannot be deduced
from the observed LFER or the observed isotope effect, as well as from other
traditional markers. In view of the difficulty of obtaining a unique mechanistic
insight by direct experimental studies, it is important to try to exploit
110
computational approaches in studies of the hydrolysis of phosphate monoesters. It
is crucial, however, to demand that the computational analysis will reproduce all
key experimental findings.
Since the initial work in our group about the systematic evaluation of the
surfaces for associative and dissociative phosphate hydrolysis (Florián, Åqvist et
al. 1998) in solution, there were several attempts to explore this crucial problem.
We will not consider here gas phase studies that were reviewed in ref. (Zhou and
Taira 1998), since these studies cannot be used to resolve mechanistic problems in
condensed phases. Here we would like to clarify that gas phase calculations that
focused on the existence of a pentacoordinated intermediate with a clear
minimum on the potential surface are not particularly useful. That is, this issue
might be relevant in some solution reactions but is of a little significance when
one considers gas phase calculations with say 80 kcal/mol barriers where the
possible difference between having a minimum of a few kcal/mol and not having
it is questionable.

Some studies yielded interesting insights in the reaction of monoanions as
was done for instance in ref. (Hu and Brinck 1999), which overlooked however
the entropy cost associated with water bridge dissociative states (see discussion in
ref. (Glennon, Villa et al. 2000)). While this issue is of a significant interest our
focus here is on the studies of the hydrolysis of phosphate dianions. In this respect
111
it is important to clarify that the controversy about the nature of the TS
(associative or dissociative) is arguably the most important open mechanistic
problem in the field. Thus attempts to imply that such controversy does not exists
at least in the case of diesters (Lopez, Dejaegere et al. 2006) and that all workers
support the associative mechanism are not justified. In fact even the associative
nature of the hydrolysis of phosphodiesters has not been established or widely
accepted (e.g. for a recent review see (Hengge 2005)).
Recently reported studies vary in their sophistication and the relevance to the
key issue of the energetics of the associative and dissociative path. For example
the systematic study of Wang et al. (Wang, Topol et al. 2003) examined the
associative path of the dianion of mono methyl phosphate but focused on the
monoanion dissociative path that was then discussed as a support for a general
dissociative mechanism (although the reaction in Ras involves probably a dianion
system). A recent Car-Parrinello DFT study of Akola and Jones (Akola and Jones
2003) found similar energy barriers of about 40 kcal/mol for the dissociative and
associative paths of ATP with magnesium ion in solution but argued, without
sufficient justifications, that the calculations for the dissociative barrier must
represent a major overestimate. These workers also provide estimates for the
reaction free energy but the corresponding PMF seems to reflect a far too short
sampling time. A systematic QM/MM study (Grigorenko, Rogov et al. 2006) of
the hydrolysis of triphosphate in solution was reported recently. This study was
112
instrumental in identifying an associative mechanism with a late proton transfer.
However the calculated barrier was about 20 kcal/mol rather than the 27 kcal/mol
barrier that is observed for such a system.
Although the above studies were quite instructive we feel that it is essential to
generate a more quantitative picture that will reproduce the overall experimental
trend in a series with different leaving groups(Klahn, Rosta et al. 2006). It is also
crucial to evaluate the complete free energy surface rather than only isolated
transition states in order to understand the overall nature of the associative or
dissociative reaction paths. The present Chapter will focus on the analysis of the
reaction surface of the hydrolysis of phosphate monoesters dianions and related
systems, using state of the art computational approaches. This study can be
considered as a follow up to the earlier attempt to map theoretically the
associative and dissociative potential energy surface (PES) in solution (Florián
and Warshel 1998), with a focus on the availability of more computer power and
several alternative approaches. We will also exploit the recent availability of more
specific experimental information about the magnitude of the activation barriers
of the hydrolysis in key systems (e.g. (Williams and Wyman 2001; Lad, Williams
et al. 2003; Wolfenden 2006)).
It will be shown that our free energy surfaces reproduce the observed
experimental trend over a wide range of rate constants and basically reproduce the
experimental LFER. The experimental validation of the calculated surfaces will
113
be taken as an indication that the calculated surfaces and reaction paths are
reasonable. It will be shown that the lowest TS changes gradually from
associative to dissociative upon increase in the acidity of the leaving group. We
will also explore the effect of a magnesium ion on the solution surface and then
move to the Mg
+ 2
catalyzed RasGAP reaction. Here we will argue that the very
flat TS region of the solution reaction can be easily changed by the protein active
site and that the TS in the RasGAP reaction involves an associative TS.

The transition states of phosphate hydrolysis reactions have been traditionally
classified (e.g. ref (Basaif, Davis et al. 1989; Barnes, Wilkie et al. 1994;
Schroeder, Lad et al. 2006)) as associative or dissociative transition states
according to the distance between the reacting phosphate and the leaving
group,
1
R , and the distance to the attacking nucleophile,
2
R (see Figure 24). This
definition is related to the More O’Ferrall-Jencks (MOFJ) diagrams (More
O'Ferrall 1970; Jencks 1985) but it is drawn in terms of bond lengths rather than
the bond order description used e.g. by Williams (Basaif, Davis et al. 1989).
114

Figure 24. A schematic description of the potential surface for the hydrolysis of
phosphomonoesters. The figure provides a clear definition of the associative and dissociative
paths and also defines the three reaction coordinates R
1
and R
2
and X.

The definition of Figure 24 is sometimes extended to “loose” TS (Basaif,
Davis et al. 1989), but the main features are uniquely defined in terms of the 2
dimensional surface outlined in Figure 24. Here it is crucial to state that the
overall surface and reaction path provide a more unique definition than alternative
attempts to clarify the hydrolysis mechanism in terms of the sum of distances (e.g.
ref (Mildvan 1997)) or other structural features. It is also important to note that
the conventional definition in terms of bond order does not provide any clear
insight about the energetics on the reaction surface, nor information about the
charge distribution ( e.g. see ref. (Barnes, Wilkie et al. 1994)). More important is
(as will be clarified below) that definitions in terms of the observed LFER and the
corresponding effective charges are quite problematic.
It is also useful to clarify that the proton of the attacking water is always
allowed to find its lowest energy at any point on the map. Thus the position of the
115
proton is not a part of the definition of the associative or dissociative path
(although it is of course an interesting parameter). With this in mind, it is not so
productive to include the proton position in the definitions of the path (e.g. ref
(Du, Black et al. 2004)). Similarly, defining the substrate assisted mechanism as a
distinguishable path (e.g. ref (Du, Black et al. 2004)) is problematic, since even
when the proton transfer from the attacking water to the phosphate oxygen atoms
is concerted with the nucleophile attack (e.g. ref (Florián and Warshel 1998;
Glennon, Villa et al. 2000)) we still have a substrate as a base mechanism. At any
rate, in view of the above points we believe that the role of theoretical analysis of
phosphate hydrolysis is to provide a clear description of the landscape of the
reacting system and to locate key reaction paths and key transition states. It is also
important to examine whether various experiments can actually help in
discriminating between different assumed mechanistic options or, more likely,
that the experimental information is less unique than what has been commonly
assumed. The present work will address the above problems by evaluating free
energy surfaces of representative systems and exploring the balance between
associative and dissociative pathways.
116

Figure 25. A VB description of the energy surfaces for phosphate hydrolysis in three extreme
cases, associative (A), dissociative (B) and concerted (C). The figure generates the surfaces
for the hydrolysis reaction by mixing four states using an EVB formulation (the effect of the
off-diagonal term is in included for simplicity) . When the surface follows one of the extreme
cases it is relatively easy to describe the corresponding LFER by shifting parabolas in the
lower part of the surface as was done in ref. (Aqvist, Kolmodin et al. 1999). RS, IS and PS
represent, respectively, reactant, intermediate and product states. The indices as and ds
stand for associative and dissociative, respectively. The dots on the surfaces designate the
minima of the corresponding diabatic surfaces.
The conceptual part of the analysis will be guided by the empirical valance
bond (EVB) description (Warshel 1991; Aqvist and Warshel 1993; Warshel and
Florian 2004). This description basically considers the four corners of Figure 24
as the centers of parabolic surfaces (see Figure 25), that correspond to different
zero order diabatic states. The mixing of these surfaces by the EVB off-diagonal
terms leads to the actual adiabatic reaction surface. If the surface follows more or
less one of the extreme paths (associative or dissociative) it is easy to quantify the
117
nature of the observed LFER by shifting the parabolas in the lower part of Figure
25, as was done in ref. (Aqvist, Kolmodin et al. 1999); the inverse is not true and
the observed LFER cannot tell us about the nature of the parabolas, since the
shifts of the parabolas in the different cases can lead to a similar LFER, (see ref.
(Aqvist, Kolmodin et al. 1999) and the present work).
This work tries to generate “consensus“ surfaces for phosphate hydrolysis in
solution by combining the results of several methods in order to establish the error
range of the calculated results, and by insisting that the calculated barriers will be
consistent with the corresponding observed rate constants.

To generate the effective free energy surfaces (PES) for the hydrolysis
reaction we started by constraining the values of the two chosen reaction
coordinates,
1
R and
2
R (see Figure 24) and optimizing the structures of the
constrained system in the gas phase (note that we later map the surface in
solution) using the density functional theory (DFT) based method B3LYP/6-
31+G*. The resulting structures are subsequently solvated using the COSMO
continuum model (Klamt and Schuurmann 1993; Barone and Cossi 1998) and
verified in some cases using the LD model in its standard implementation. An
additional local relaxation of the energy of the given solvent model for points
{
1
R ,
2
R } on the reaction surface of the solvated system is also applied in some
cases. The error range associated with the choice of the quantum method and the
118
basis set is also established. This is done by evaluating selective points using the
methods MP2/6-311++G** and B3LYP/Lanl2DZ. The evaluation of the
corresponding PESs required us to address two problems. The first problem is the
dependence of the optimized structures on their initialization. This problem is
handled by starting the optimization procedure from several initial structures for a
given point {
1
R ,
2
R }, and taking the lowest value as the corresponding effective
energy. The second problem arises from the neglect of additional possible
reaction coordinates, in our case the coordinate, X, that describes the proton
transfer from the attacking water to the phosphate oxygen atoms, (Figure 24).
This problem is handled by examining the structures along the reaction paths
carefully, after finding the initial reaction paths on our PES, to ensure that all
reaction paths develop continuously, i.e. that all points on the PES are actually
connected. In cases where we find discontinuities, especially for the proton
transfer reaction coordinate, we estimate the energy barrier for the proton transfer
with additional energy scans in the direction of this coordinate. At any rate, at the
end of our mapping procedure we end up with a carefully scanned solution
surface.
In principle we should include in Eq. (4.11) zero point energy effects and
activation entropies. However, the zero point contributions to the activation
barriers were found to be smaller than 1 kcal/mol. The activation entropy effect
was found to be much smaller than their gas phase barrier estimates and was
119
neglected (see discussion in (Strajbl, Florian et al. 2001)). Furthermore the
experimental studies of Wolfenden (Wolfenden 2006) established that for the
reactions of the type studied here (hydrolysis by water) the contributions of the
activation entropies are small. A more careful study that will use a restraint
release approach is left to subsequent studies, see Chapter 7.

In our QM/MM treatment used standard link atom technique rather than the
more reliable hybrid orbital type treatments (e.g. (Assfeld and Rivail 1996))
because we found it very convenient to link our MM program MOLARIS (Chu,
Villa et al. 2004) (Lee, Chu et al. 1993) to standard QM programs (in this case
GAUSSIAN 98 (Frisch, Trucks et al. 1998)) through an automated external script.
The calculations were performed with the B3LYP/DFT method with 6-31+G*
basis set for the QM/MM calculations.
The specific FEP parameters used in this work are the standard ENZYMIX
parameters and the simulations were performed using the MOLARIS program
with a surface constraint all atom solvent (SCAAS) (King and Warshel 1989)
spherical boundary conditions of 22 Å radius and the LRF long range treatment
(Lee and Warshel 1992). The solvation QM(ai)/MM free energies were evaluated
by performing the AC-FEP calculations (Warshel, Sussman et al. 1986) with 31
windows of 10ps length each and a time step of 1fs. The simulations were
performed at 300K.
120

6.4 Free energy surfaces for the hydrolysis of
phosphate-monoesters in solution

(i) Monomethyl phosphate dianion
We started our systematic study by evaluating the free energy surfaces for the
monomethyl phosphate dianion. The results of the calculation are summarized in
Figure 26 and Figure 27.

Figure 26. The calculated surface for the hydrolysis of the mono methyl phosphate dianion
in the space defined by R
1
and R
2
. The figure depicts the reaction path that connects the
RS with the PS and the corresponding TS. Energies are given in kcal/mol and distances in Å.

This system was explored systematically in our previous study(Florián, Åqvist
et al. 1998) and the present results are similar to those obtained before, including
the results considered briefly in the study of ref. (Wang, Topol et al. 2003), that
focused, however, on the monoanion. Ref. (Florián, Åqvist et al. 1998) gave two
associative TSs of 40 kcal/mol and 43 kcal/mol and a dissociative TS of 56
121
kcal/mol, while ref. (Wang, Topol et al. 2003) gave an associative barrier of 42
kcal/mol for the dianion.

Figure 27. The calculated TS structure for the hydrolysis of the mono methyl phosphate
dianion.

As seen from the figures the hydrolysis reaction is associative with a relatively
high barrier of 47 kcal/mol. This result is in a good agreement with the
experimental results (Lad, Williams et al. 2003) for alkyl phosphates (about 44
kcal/mol for methylphosphate at 25 C°). Interestingly the calculated activation
barrier for hydrolysis of methylphosphate by OH
−
is about 40 kcal/mol. Thus it
would be interesting to assess the contribution from this OH
−
pathway to the
hydrolysis of monomethyl phosphate in water and to examine whether the
experimental estimate of the contribution to the barrier from the hydrolysis by
water should be slightly modified. The nature of the overall surface can be
described by case A in Figure 25.

122
(ii) Phenylphosphate dianion
Next we evaluate the surface for a case that involves a leaving group with a
lower pK
a
, namely the phenyl phosphate dianion. The results of these calculations
are summarized in Figure 28. As seen from the figure we still have as associative
TS but the activation barrier is much smaller than in the case of the monomethyl
phosphate. The calculated barrier of 35 kcal/mol is again in very good agreement
with the observed value of about 35 kcal/mol (Lad, Williams et al. 2003).

Figure 28. The calculated surface for the hydrolysis of the phenyl phosphate dianion in the
space defined by R
1
and R
2
. The figure shows the reaction path that connects the RS with the
PS and the corresponding TS. Energies are given in kcal/mol and distances in Å.

(iii) Mono methyl pyrophosphate trianion
After establishing the strong dependence of the barrier on the pK
a
of the
leaving group we returned to the systems that are more relevant to the hydrolysis
of GTP and ATP although in these systems the electrostatic repulsion is screened
by a Mg
+ 2
ion. As a representative system we considered cases with methyl
phosphates as a leaving group, since the effect of the a phosphate group on the
123
hydrolysis of the g-b bond is relatively small (Aqvist, Kolmodin et al. 1999; Allin
and Gerwert 2001). We start the study of these systems by evaluating the surface
for the hydrolysis of the mono methyl pyrophosphate trianion. The results of the
calculations are summarized in Figure 29 and Figure 30. As seen from Figure 29
the dissociative TS is now at a lower energy than the associative TS.

Figure 29. The calculated surface for the hydrolysis of the mono methyl pyrophosphate
trianion in the space defined by R
1
and R
2
. The figure shows the reaction path that
connects the RS with the PS and the corresponding TS. An alternative associative path with
higher energy barrier is shown (dashed line). Energies are given in kcal/mol and distances in
Å.

Thus we conclude that the decrease in the pK
a
of the leaving group changes
the character of the TS from associative to dissociative. The EVB behaviour of
the surface in the path through the dissociative TS can be described by case B of
Figure 25. The calculated barrier is about 34 kcal/mol which is in a very good
agreement with the observed barrier for the pyrophosphate trianion (about 33
kcal/mol) (Wolfenden 2006). The success of our calculations of this highly
charged system indicates that the strong electrostatic repulsion is properly
124
screened by the solvent model used. It is also apparent that the real solvent (i.e.
water) screens theses interaction quite effectively since as will be shown below
the trianion system falls nicely of the observed LFER.

Figure 30. The calculated TS structure for the hydrolysis of the mono methyl pyrophosphate
trianion.

(iv) Mono methyl pyrophosphate dianion
In view of the relatively high barrier of the trianion system it is important to
consider the effect of protonation of the pyrophosphate system. This is done in the
calculations summarized in Figure 31. Now the difference in energy between the
associative and the dissociative TS becomes quite small. Although the protonated
system has a lower barrier than the deprotonated system it is very likely that the
hydrolysis of GTP and ATP in proteins involves deprotonated systems. It is
interesting to note that in the reactant state the proton is always transferred to the
phosphate that is hydrolyzed even when the position of the proton is initialized at
the leaving group. This proton is eventually transferred back to the leaving group
as expected on the reaction path toward the product state.
125

Figure 31. The calculated surface for the hydrolysis of the mono methyl pyrophosphate
dianion in the space defined by R
1
and R
2
. The figure shows the reaction path that connects
the RS with the PS and the corresponding TS. An alternative associative path with higher
energy barrier is shown (dashed line). Energies are given in kcal/mol and distances in Å.

(v) Magnesium mono methyl pyrophosphate trianion
Despite our success in reproducing the barrier for the hydrolysis of mono
methyl pyroposphate trianion, we felt that a good reference system for the
GTPase and ATPase reactions in proteins can be provided by the hydrolysis of
mono methyl pyrophosphate in the presence of an Mg
+ 2
ion. The position of the
magnesium ion was initialized between the two phosphates to form a bidentate
complex (see Figure 33), resembling the usual coordination of the magnesium ion
with GTP or ATP in enzymes like Ras p21 or ATPases. Furthermore we included
into the QM system four additional water molecules around the magnesium ion to
complete its first solvation shell. The results of these calculations are summarized
in Figure 32, Figure 33.
126

Figure 32. The calculated surface for the hydrolysis of the magnesium mono methyl
pyrophosphate trianion in the space defined by R
1
and R
2
. The figure shows three possible
reaction paths with similar energy barriers that connect the RS with the PS and the
corresponding TSs. Energies are given in kcal/mol and distances in Å.

This solution system has basically three TSs with similar energy but the
dissociative TS is likely to be with higher energy in active sites of proteins. That
is, although it is absolutely essential to examine the situation inside the protein by
QM/MM calculations that include the protein it seem from the surface of Figure
32 that moving to the dissociative TS, which involves a large increase in the
distance between the neuclophile and leaving group, costs only a few kcal/mol.
Thus even small steric forces in an active site that is originally complementary to
the reactant state are likely confine the reaction surface to the associative region
where the distance between the nucleophile and the leaving group is small.
127

Figure 33. Three possible calculated TSs structures for the hydrolysis of the magnesium
mono methyl pyrophosphate trianion.
At any rate, the calculations indicate that the potential surface involves a flat
ridge at the TS region and that small perturbations can easily shift the TS position.
The surface of Figure 32 is quite complicated and can be considered as a
combination of case A, B and C in Figure 25. This can be represented with a
single EVB surface with special mixing terms. Interestingly the lowest activation
barrier for the associative path is around 32 kcal/mol as compared to the observed
barrier of 28 kcal/mol (Kotting and Gerwert 2004).

128
(vi) Other systems
In addition to the systems considered above we also examined additional
systems to gain more insight about the overall trend in phosphate hydrolysis
reactions. The additional systems include the hydrolysis of the acetyl phosphate
dianion, with a very good acetate leaving group. The corresponding results are
included in our LFER considerations that will be described below.

The present study allows us to explore the validity of the long standing
assumption that the observed LFER can be used to determine the nature of the TS
in phosphate hydrolysis reactions (e.g see ref (Admiraal and Herschlag 1995)).
Our previous studies (Florián, Åqvist et al. 1998; Aqvist, Kolmodin et al. 1999)
already clarified that the traditional assumptions about the relationship between
the observed LFER and the nature of the TS are very problematic. The main
problem is the LFER should be based on at least three VB states (e.g. Figure 25)
and thus does not follow the traditional assumptions in physical organic
chemistry. This point has also been demonstrated in a pioneering ab initio study
of phosphate monoesther dianions in solution(Florián, Åqvist et al. 1998), that
involved methanol as a nucleophile and substituted phenyl leaving groups. Here
we continue along the previous line of reasoning by comparing the calculated and
observed LFER for the hydrolysis of different phosphomonoesters. The results of
129
the present study are summarized in Figure 34 and show that we could reproduce
the observed trend over a very wide range of rate constants.
Table 7. The correlation between the rate constants (and activation barriers) and the
observed pK
a
values for the hydrolysis reaction of phosphomonoesters with different
leaving groups.
a

Leaving group pK
a

calc
G
≠
Δ
exp
G
≠
Δ
log k
calc
log k
exp
CH
3
O
-
15.5
g
47 44
b
-20.1 -18.2
b

PhO
-
10.0
g
35 35
b
-11.9 -12.0
b

CH
3
PO
4
2-
6.3
g
34 33
c
-11.0 -10.3
c

HPO
4
2-
7.2
g
32 29
f
-9.8 -7.7
f

(Mg
2+
)CH
3
PO
4
2- i
6.5
h
31 28
e
-8.9 -6.8
e

CH
3
COO
-
4.8
g
28 28
d
-6.8 -6.6
d

PhCOO
-
4.2
g
25 24
d
-4.4 -4.3
d

a
The rate constants are given in s
-1
and calculated for 39C°. The energies are given in kcal/mol
and the observed values are obtained from the observed rate constants using transition state theory.
The reacting system includes an attacking water molecule, a phosphate dianion group and the
indicated leaving group. The sources of the observed rate constants are refs. (Lad, Williams et al.
2003), (Admiraal and Herschlag 1995), (Di Sabato and Jencks 1961), (Kotting and Gerwert 2004)
and (Wolfenden 2006) for b, c, d, e and f, respectively. The observed pK
a
values are taken from
ref. (Jencks and Regenstein 1976) and (Sigel, Bianchi et al. 2001) for g and h.
i
This system involves a Mg
2+
ion and is compared to measurements of solvated GTP (Kotting
and Gerwert 2004).

Apparently as pointed out above, the rate determining TSs are very different
for the different molecules considered in the LFER study. Thus it seems clear that
the fact that a given molecule that follows this LFER cannot be used to deduce the
nature of the corresponding TS.
130

Figure 34. Comparison of calculated (square boxes) with observed (filled circles) LFER i.e.
the correlation between the log k of the different reactions and the corresponding pK
a
s of
the leaving groups for phosphate monoester dianions with different leaving groups. The rate
constants are given at 39 C° in units of s
1 −
. The corresponding data is given in Table 7. The
linear curve is taken from ref. (Lad, Williams et al. 2003) where it was obtained as a least-
square fit log k=1.0-1.26 pK
a
to measured hydrolysis rate constants of several
phosphomonoesters.
The LFER observed in phosphate hydrolysis reactions has been traditionally
used to define effective TS charges (e.g. refs (Leffler and Grundwald 1963;
Williams 1984)). This definition relates the change in the charge of the leaving
group to the Brønsted a or to the ratio between b
lg
and b
eq
(Williams 1984).
Such an analysis is very reasonable and would be physically consistent in the case
of intersection of two VB states. Otherwise it may turn out to be just a way for
classifying the given reaction. Here we perform a preliminary analysis of the
nature of the effective charges by calculating the change in the charges of the
leaving group from the reactant to the product state for the associative and
131
dissociative paths of the hydrolysis of the mono methyl pyrophosphate trianion.
The corresponding results are summarized in Table 8.
Table 8. The charges of the leaving group in different states of the hydrolysis of some
selected systems.
a

mono
methyl
phos-
phate
dianion
mono
methyl
pyrophos-
phate
trianion
mono
methyl
pyrophos-
phate
trianion
b

GTP in
solution
tetraanion
GTP in
RasGAP
tetraanion

leaving groups
state
CH
3
O PO
4
CH
3
PO
4
CH
3
PO
4
PO
3
CH
3
PO
4
PO
3
CH
3
RS
-0.29
(-1.71)
-1.44
(-1.56)
-1.42
(-1.58)
-2.48
(-1.52)
-2.51
(-1.49)
associative TS
-0.51
(-1.49)
-1.74
(-1.26)
-1.86
(-1.14)
-2.73
(-1.27)
-2.79
(-1.21)
dissociative
TS
-0.97
(-1.03)
-1.89
(-1.11)
-1.72
(-1.28)
-2.87
(-1.13)
-2.83
(-1.17)
PS
-1.88
(-1.12)
-1.69
(-1.31)
-2.88
(-1.12)
-2.89
(-1.11)
a
The charges are given in the unit of the elementary charge. The reported charges of the leaving
group are taken as the sum of the partial charges of all the atoms of the leaving group, for different
points on the PES for the hydrolysis reaction. The charge of the remaining part of the system (i.e.
the phosphate and hydrolyzing water molecule) is given in parenthesis. The charges are evaluated
according to the Merz-Singh-Kollman scheme.
b
This system contained a Mg
2+
ion with its hydration shell (i.e. four additional water
molecules). The charge transferred to the magnesium ion and its hydration shell was added in
equal halves to the leaving group and the remaining system consisting of the phosphate and the
hydrolyzing water molecule.
The first point that emerges from the table is the fact that the change in the TS
charge upon moving from the RS to the TS does not follow a simple rule. For
example when we have CH
3
O
−
as a leaving group the changes are -0.22 and -
0.68 for the associative and dissociative TSs, respectively, while for the case with
PO
4
CH
3
−
as a leaving group the changes are found to be -0.30 and -0.45 for the
associative and dissociative TSs respectively. This means that the effective TS
132
charges cannot be used to describe the different mechanisms. Similarly and quite
significantly it indicates (in contrast to the implications of ref. (Admiraal and
Herschlag 1995)) that the electrostatic interaction of the environment in Ras and
related systems can stabilize the associative and dissociative TSs in a similar way
(Glennon, Villa et al. 2000). This point will be quantified below. The second
point is that the change in ab initio charges between the RS and TS are not similar
to the changes in charge deduced by traditional physical organic chemistry (e.g.
refs (Leffler and Grundwald 1963; Williams 1984; Davis, Hall et al. 1988)). For
example in the case of ATP hydrolysis(Admiraal and Herschlag 1995) and GTP
hydrolysis(Maegley, Admiral et al. 1996) in solution it is assumed that the change
in charge is -0.83. On the other hand, our ab initio calculations give -0.30 for the
low energy TS and -0.45 for the dissociative TS. The situation becomes even
more striking in the case of the hydrolysis reaction in RasGAP when the changes
in charge are -0.28 and -0.32 for the associative and dissociative TSs,
respectively. Here one could argue that the effective charges should not be equal
to the ab initio charges but the problem is, in fact, more fundamental. That is,
traditional physical organic chemistry approaches provide an excellent estimate of
the effective charge in cases of two intersecting energy parabolas, which might be
a reasonable description in the concerted case (case C in Figure 25). However, in
the more general case the most unique definition of group charge can probably be
obtained by the VB approach that provide the natural definition of the degree of
133
charge transfer (see e.g. ref. (Warshel 1991), chapter 2) and can quantify the
effect of the environment and various substantiates on such charges (see also ref.
(Shaik, Schlegel et al. 1992)). Now the charge changes obtained by well
calibrated EVB Hamiltonians should be equal to the corresponding ab initio ESP
charge changes, since this is the way we calibrate the EVB surface. This means
that the ESP charges provide what is probably the best description of the effective
charges and that our conclusions about the complex nature of these charges are
valid.
Table 9. Difference of zero point energies (relative to the reactant state) for different TSs of
the magnesium mono methyl pyrophosphate trianion.
a

associative TS

dissociative TS
isotopic
composition GA SA GD SD
γ
16
O
4
-0.130 -0.609 -1.037 -1.973
γ
18
O
3
-0.136 -0.625 -1.040 -1.978
γ
18
O
4
-0.137 -0.628 -1.039 -1.980
a
The zero point energies are given in kcal/mol. γ
16
O
4
denotes the phosphate with the natural
16
O
isotope, γ
18
O
3
and γ
18
O
4
are the
18
O labeled atoms according to the notation of Figure 35. The
calculations considered structures that differ from each other only slightly for both, the associative
and the dissociative mechanism (this allows us to examine the sensitivity of our calculated values
to small structural changes). GA and GD are gas phase optimized and subsequently solvated TSs,
while SA and SD are solution optimized TSs (here G and S designates gas and solution
respectively) . The corresponding structures are depicted in the supplementary material.
Table 10. Calculated isotope effects for the systems considered in Table 9.
a

associative structures dissociative structures isotopic
composition GA SA GD SD
γ
18
O
3
1.0092 1.0256 1.0061 1.0081
γ
18
O
4
1.0122 1.0308 1.0041 1.0112
a
γ
18
O
3
and γ
18
O
4
are the
18
O labeled ions according to the scheme of Fig. 12.
The different structures are defined in the caption of Table 4.

134
The O16/O18 isotope effect has been used extensively in attempts to
determine the nature of the transition state in phosphoryl transfer reactions (e.g.
refs (Du, Black et al. 2004; Cleland and Hengge 2006)). Although it is tempting
to accept the assignments obtained from isotope effect studies it has been argued
that the corresponding analysis is far less unique than what has been commonly
assumed (Aqvist, Kolmodin et al. 1999). Basically the assumption of the
relationship between the isotope effect and the reaction path are in some respects
circular since the nature of the TSs have never been established by a direct
experimental approach and the basic arguments about the expected isotope effect
are problematic (Aqvist, Kolmodin et al. 1999). The general difficulty of
obtaining unique mechanistic interpretation of isotope effects experiments has
also been noted by the experimental community (e.g. (Anderson, Cassano et al.
2006)). Interestingly, a recent theoretical study (Liu, Gregersen et al. 2006) of
transferification found what is basically an associative TS with isotope effects
larger than one that could have been probably interpreted traditionally as an
evidence for a dissociative mechanism.
135
O P
O*
O*
O P
OH
O
O*
O* P
O*
O*
O P
OH
O
O*
-
-
γ
18
O
3
-
-
γ
18
O
4

Figure 35. Defining the isotopic labeling of the pyrophosphate oxygens used in Table 9and
Table 10. The oxygen atoms marked with an asterisk are substituted with an
18
O isotope
whereas all other oxygen atoms are left with
16
O.

In view of the above problems we calculated the isotope effects for the
hydrolysis of the magnesium mono methyl pyrophosphate trianion (shown in
Figure 32), considering different TSs as the rate limiting states. The calculations
were done using the B3LYP/6-311++G method for the gas phase TSs (described
in the supplementary information), with the isotopic composition described in
Figure 35. The calculated zero point energies are summarized in Table 9. These
zero point energies were used to evaluate the corresponding isotope effects, which
are summarized in Table 10. The calculations indicate that the isotope effect
expected from associative and disociative TSs is quite similar and thus cannot be
used to establish the nature of the TS.
Interestingly, the calculated isotope effects are in the range found
experimentally for the GTPase reaction of Ras (Du, Black et al. 2004) which are
1.0013 for g
18
O
3
/ g
16
O
3
and 1.0128 for g
18
O
4
/ g
16
O
4
, respectively. More
accurate results could have been obtained by evaluating the isotope effects in
136
solutions and proteins using the centroid path integral methods (Gillan 1987; Voth
1996) and in particular the robust quantized classical path (QCP) simulation
approach (Olsson, W. W. Parson et al. 2006). However, the required evaluation of
realistic vibrational spectra in solution or proteins is quite demanding and is left to
subsequent works since our main aim is to establish the fact that the very
qualitative reasoning that has been used in analyzing phosphate hydrolysis
reactions is not supported by actual quantum mechanical calculations. Of course,
it would be interesting to use more careful calculations to try to see which TSs in,
for example, the reaction of Ras give a better agreement with the observed isotope
effect, but such a study cannot be done by the current use of circular arguments
and should be based on extremely careful simulation studies.

137
6.5 The Free Energy Surface for the Reaction of
Ras

Figure 36. The calculated surface for the hydrolysis of GTP in RasGAP in the space defined
by R
1
and R
2
. The figure shows a part of the reaction path (dashed line). The PS region is
shown in Figure 37. Energies are given in kcal/mol and distances in Å.

Armed with the insight from the above solution studies we can explore related
enzymatic reactions. In particular we are interested in understanding the GTPase
reaction of the RasGAP system. This GTP hydrolysis involves a large catalytic
effect (about 11 kcal/mol) that has been the central subject of biochemical and
structural studies of signal transduction (for reviews see (Sprang 1997; Vetter and
Wittinghofer 1999)). Although significant progress has been made in the
theoretical description of this reaction it has been hard to asses whether the
reaction involves an associative or dissociative mechanism (Maegley, Admiral et
al. 1996; Schweins, Geyer et al. 1996; Sprang 1997; Vetter and Wittinghofer
1999; Glennon, Villa et al. 2000). The main difficulties have been associated with
138
the fact that both pathways can be catalyzed by the protein in a similar way
(Glennon, Villa et al. 2000). Here we perform a careful and systematic QM/MM
study of the surface for GTP hydrolysis in Ras/GAP. The calculations started with
an extensive relaxation using EVB potential surfaces similar to those used in ref
(Shurki and Warshel 2004). The relaxed structures were used as a starting point
for the QM(ai)/MM FEP calculations that were done for the
1
R and
2
R space. Note
that in each point of this surface we perform extensive FEP simulations that allow
for a full relaxation of the protein environment. The results of this simulation
study are summarized in Figure 36 and Figure 37. Apparently the shape of the
surface is similar to that obtained for magnesium mono methyl pyrophosphate
trianion in water, but it seems to have more associative character. The calculated
barrier height is found to be 14 kcal/mol, which is quite similar to the 16 kcal/mol
observed barrier. The very good agreement may be somewhat coincidental since
the calculations did not involve an actual evaluation of the potential of mean force
(PMF). However, the mechanistic conclusions are probably valid. A completely
conclusive study might require additional exploration of the effect of different
mutations and the evaluation of the free energy surface by an approach that also
consider the solute entropy (Rosta, Klahn et al. 2006). It would also be interesting
to reproduce the observed LFER in Ras (Schweins, Geyer et al. 1996). However,
the mechanistic picture that emerges from the present study is probably
significantly more conclusive than the assumed mechanistic information that
139
emerged from interpretation of experimental results, such as isotope effects (Du,
Black et al. 2004) or even from structural studies of transition state analogues
(whose direct relevance to the exact transition state is hard to establish).

Figure 37. The calculated surface for the hydrolysis of GTP in RasGAP in the space defined
by the proton transfer coordinate, X and R
2
, fixing
1
R at the TS distance at about 2.1 Å.
The figure shows the reaction path that connects the RS with the PS and the corresponding
TS. Energies are given in kcal/mol and distances in Å. (See method section for more details
about our mapping procedure.)

Finally it should be kept in mind, that the premise of several experimental
attempts to explore the nature of the TS in Ras is very problematic. That is, it has
been assumed in several recent works (e.g. refs (Maegley, Admiral et al. 1996; Li
and Cui 2004)) that the transition state of the reference reaction in solution is
dissociative and thus presumably the same is true for the TS in Ras. Apparently as
shown in this work this assumption is not justified. That is, ref (Admiraal and
Herschlag 1995) envisioned a surface with an enormous barrier for the
140
dissociative path (Fig 1 of ref (Admiraal and Herschlag 1995)),yet the careful
study reported here indicates that all the experimental problems that were used to
support the proposal of a dissociative TS are equally valid for supporting the
proposal of an associative TS.

6.6 Summary
At this point it might be useful to consider the corresponding error range of
our calculations and it s impact on our overall conclusions. We start by comparing
in Table 11 the energetics of different key points on the surfaces of the
magnesium mono methyl pyrophosphate trianion obtained by both the LD and the
QM/MM –FEP approach. As seen from the table the very different models give
similar results with deviations between 1 to 5 kcal/mol. The effect of using
different quantum mechanical Hamiltonians was estimated to be relatively small
since the B3LYP/6-311++G** method used here is quite reliable. Finally, we
examined the magnitude of the zero point energy corrections (see Table 9) and
found them to be smaller than 2 kcal/mol for the energy barrier in the most
extreme case and smaller than 1 kcal/mol in most cases. We did not include the
gas phase entropic contributions since our earlier studies (e.g. ref (Strajbl, Florian
et al. 2001)) indicated that these contributions are small and cannot be deduced
from the gas phase calculations. Furthermore the experimental studies of
Wolfenden (Wolfenden 2006) established that for the reactions of the type studied
141
here (hydrolysis by water molecule) the contributions of the activation entropies
are small. A more careful study that will use a restraint release approach is out of
the scope of this work, see Chapter 7.
Table 11. Comparing the energies obtained by different solvent models and different ab
initio methods for characteristic points on the PES of the hydrolysis of the magnesium mono
methyl pyrophosphate trianion.
a

COSMO LD FEP

structure
B3LYP/
6-311++G**
MP2/
6-311++G**
B3LYP/
Lan2DZ

RS 0 0 0 0 0
associative TS 32 36 33 33 35
concerted TS 31 36 35 32 34
dissociative TS 31 33 36 31 35
completely
dissoc. TS
31 33 34 31 34
PS 2 -3 3 -1 7
a
The energies are given in kcal/mol. The different solvent models (COSMO,
LD and FEP) are discussed in the method section. The B3LYP/6-311++G** has
been the method of for the evaluation of all the PESs.

It is important to reemphasis that our calculations are not base on taking gas
phase transition states and solvating them but on a very systematic construction of
the free energy surface for the solution reaction by a careful interpolation of
values obtained by consistent calculations of the energies of the solvated quantum
system (see for example discussion in the method section). Having the complete
surface is also much more informative than having only some specific TSs since it
provides crucial information about the nature of the overall landscape. This
142
should allow one to better appreciate the nature of the reaction paths. It is also
important to point out that our conclusions on the rather late proton transfer step
were determined by a careful scanning and not by gas phase calculations and
these results were contained by the different approaches used here.
Another point of concern is the ability of the present study to reproduce the
proper screening for the strongly interacting negative charges in say the
hydrolysis of the mono methyl pyrophosphate. Here we would like to clarify that
we are not dealing with gas phase calculations but with proper calculations in
solutions. In particular we are encouraged by obtaining similar results for the
simplified solvent models: the COSMO and LD models and the microscopic
QM/MM calculations. Furthermore it is important to point out that we have
reproduced the actual experimental trend in such systems.
At this stage it is important to examine the relations between the calculations
of the reference reactions and the corresponding observed results. As established
in Table 7 and Figure 34 we obtained a very encouraging agreement between the
calculated and observed barriers where the calculated rate constant follows nicely
the observed LFER. This lends significant credibility to the calculated results.
One might still argue that with an estimated error range of 2-4 kcal mole there
is no reason to accept our conclusions. However, our main point is not the exact
height of the associative and dissociative barriers but the fact that they are
143
frequently similar and that their relative height changes with the pK
a
of the
leaving group. In our view this point has been established by the present study.

This work presented what is perhaps the most extensive study of the potential
surfaces for the hydrolysis of phosphate monoester dianions in solution, providing
complete 2 dimensional surfaces for key representative systems. The resulting
TSs reproduced the observed LFER in a semiquantitative way, thus lending
credibility to the subsequent analysis. The calculated surfaces were found to be
quite flat at the TS region, changing between associative and dissociative
character depending on the nature of the leaving group. This established that it is
unjustified to use the observed LFER to deduce the actual mechanism in a unique
way.

Evaluation of the transition state charges establish that the calculated charges
and their correlation with the nature of the TS are quite different than what has
been assumed in the early insightful studies of physical organic chemists (Leffler
and Grundwald 1963; Williams 1984) (see below). It is also demonstrated that
isotope effects cannot be used in a unique mechanistic analysis of the hydrolysis
of phosphate monoester dianions.
Apparently, although early experimental studies were extremely insightful, the
resulting insight cannot be used for a detailed analysis and must be augmented by
144
ab initio QM/MM calculations. This problem becomes even more serious when
one tries to evaluate the effect of enzyme active sites on the reaction surface for
phosphate hydrolysis and /or the effect of different mutations.
This work explored the uniqueness of the effective charge extracted from
LFER analysis (Leffler and Grundwald 1963; Williams 1984; Davis, Hall et al.
1988). As mentioned earlier, the effective charges provide a useful way for
describing the change in charge of the leaving group in cases which can be
represented uniquely by two intersecting parabolas. However, in more complex
cases they can only be considered as a book keeping description of the observed
LFER. This point has been established by evaluating the change in the ab initio
charge of the leaving group for associative and dissociative mechanisms. It was
found that the effective charges are very similar in both cases. Further studies
should involve an EVB analysis of the proper effective charges.
After validating our approach by studies of the solution reaction and in
particular the reaction in the presence of the Mg
+ 2
ion, we evaluated the surface
for the GTPase reaction in the Ras/GAP system. Here it is found that the surface
is associative in nature and that again it is impossible to establish by simple
experiments (e.g. isotope effect or LFER) whether it is associative or dissociative.
Nevertheless, it would be useful to try to reproduce the observed LFER in Ras
and to attempt to reproduce the effect of different mutations on the vibrational
frequencies of the bound GTP in the associative and dissociative TSs.
145

The nature of the charge distribution in the TSs of the GTP hydrolysis
reaction has been an issue of significant interest as much as the catalytic reaction
of Ras is concerned. EVB studies (e.g. refs (Glennon, Villa et al. 2000)) that
reproduced the experimental catalytic effect found that there is a significant shift
of the positive electrostatic potential of the protein towards the leaving group
upon transfer of the substrate to the protein active site in both the associative and
dissociative paths (Glennon, Villa et al. 2000; Shurki and Warshel 2004). In fact,
the idea and more importantly the demonstration of a shift of the electrostatic
potential was introduced first in ref (Muegge, Schweins et al. 1996), where it was
correlated with the change of the protein between its GDP-bound to GTP-bound
structures. At any rate, the shift of the electrostatic potential should lead to a
charge shift, and indeed the existence of such a shift has been supported, at least
in the reactant state, by careful FTIR studies (Allin and Gerwert 2001) and
quantified in a theoretical work (Klahn, Schlitter et al. 2005). The present work
focused on the implication that the shift of the electrostatic potential and the
corresponding catalytic effect indicates that we have a dissociative TS (e.g. ref
(Maegley, Admiral et al. 1996)). First, as it seems to be clear from Table 8, the
change in charge shift upon moving from the RS to the TS in the RasGAP system
is virtually equal for the associative and the dissociative path. Second, true
assessment of the catalytic effect in both paths cannot be accomplished without
146
actual calculations of the corresponding energetics for both paths. To the best of
our knowledge the only study that addressed this issue has been reported in ref.
(Glennon, Villa et al. 2000). This study found a similar electrostatic catalysis for
the associative and dissociative mechanisms.
Recent studies have examined the energetics of the GTP hydrolysis reaction in
solution (Grigorenko, Rogov et al. 2006) and in the RasGAP system (Grigorenko,
Nemukhin et al. 2005) by a QM/MM approach which involved energy
minimization and TS search. These studies were very instructive in pointing out
towards an associative TS where the proton is still attacked to the nucleophilic
water molecule. However, these studies have not provided the overall free energy
surface and underestimated significantly the barrier for the solution reaction.
Furthermore, the rate determining steps was found to be very different for the
solution and the protein reactions (these rate limiting steps involved a formation
of an associative TS in the solution reaction and a proton transfer to the leaving
group in the protein reaction). At present it is hard to pinpoint the source of the
difference between the results of these studies and the results of the present study.
We can only point out that in view of the complexity of the system studied, we
placed so much emphasis on calibrating and validating our calculations using the
solution LFER and examining the results using several approaches. No such
validation effort was reported by other studies.
147
As should be clear from the above discussion we believe that the present study
can be considered as the most systematic available computational analysis of the
hydrolysis of phosphate monoesters dianions in water. The reason for our view is
not only the fact that we reproduce nicely the observed rate constants for very
different systems, but also our estimate of the error range by using quite different
approaches. We also like to reemphasis that our study is not base on any gas
phase transition state search but on a very careful mapping of the surface of the
solvated system. Having a full surface rather than just isolated TSs is in our view
a major advantage of the present treatment. It is also important to mote that that
we do not base our claim for reasonable accuracy on the accuracy of the DFT
calculations which is clearly not perfect but on the compounded experience from
modeling the overall experimental trend. This point is particularly important when
we address the effect of the protein on the calculated trend, In this respect it is
useful to mention that we were able to reproduce the every large solvent screening
of the repulsion between the negative charges of the reacting system since we
reproduce the experimental barrier in the highly charged mono methyl
pyrophosphate trianion system.
There has been a significant interest in the validity of the substrate as a base
proposal (e.g. refs (Schweins, Langen et al. 1994; Schweins, Geyer et al. 1995))
for Ras and related systems. In some studies (e.g. ref (Admiraal and Herschlag
2000)) it was argued that there are experimental evidences against such a
148
mechanism in solution. However, it has been clarified (e.g. ref (Glennon, Villa et
al. 2000)) that there is no direct experiment that can establish this issue. There is,
nevertheless, a point of clarification. That is, the substrate as a base idea addresses
the fact that the proton transfer (PT) step in Ras does not involve Gln61 as a base
(Langen, Schweins et al. 1992) and that the proton is transferred in one way or
another to the g-phosphate. As discussed in refs (Florián, Åqvist et al. 1998) and
(Glennon, Villa et al. 2000), this transfer can be done in a concerted way with the
nucleophilic attack and still be considered as a phosphate as a base mechanism (as
long as the mechanism is an associative mechanism). The phosphate as a base
mechanism does not require any preequilibrium PT step as implied in some works
(e.g. ref (Du, Black et al. 2004)) and only requires that the PT will occur on the
way or near the associative TS. It is also important to clarify that, in contrast to
early assumptions (see discussion in ref. (Florián and Warshel 1997)), it is
impossible to determine experimentally whether the water attack occurs in a
stepwise mechanism or a concerted way (since it is impossible to examine the
OH
-
attack on a protonated phosphate due to conflicting pH requirements). At any
rate, as already clarified in ref (Glennon, Villa et al. 2000) the main issue is
whether we have an associative or dissociative TS in the Ras and RasGAP
reactions, and whether Gln61 participates in the reaction in a direct way.
The dissociative mechanism has been frequently associated with a transition
state with a protonated TS (Maegley, Admiral et al. 1996). It has also been
149
suggested that this is a well established experimental fact for the solution reaction
(e.g. refs (Maegley, Admiral et al. 1996) (Li and Cui 2004) ) and the it consistent
with the observed isotope effect (Pecoraro, Hermes et al. 1984). While in the case
of a leaving group with a very large pK
a
we may have a fully dissociative
mechanism and the proton of an attacking water may be transferred to the
negatively charged leaving group, the situation is quite different in cases of S
N
2
and/or associative reactions, where the first TS involves a concerted PT to the
phosphate or a transition state where the proton is not yet transferred.
Obviously the analysis presented here is not the last word in the exploration of
the nature of the TS in the hydrolysis of phosphomonoesters in solutions and
proteins. Further QM/MM calculations should explore the effect of the solute
entropy by FEP approaches that do not fix the solute (Rosta, Klahn et al. 2006).
However, we believe that the present systematic study has helped to advance our
understanding of the nature of the reaction path and that the analysis of the type
presented here is crucial for further progress in studies of G-proteins and related
systems.
It is also useful to note that the present study provided compelling reminders
about the major difference between experimental information and the
interpretation of such experiments. As shown here in several key examples (e.g.
the analysis of LFERs and isotope effects) it is not justified to accept traditional
interpretations of experimental result as experimental facts. In the case of high
150
energy TSs it is crucial to add computational studies as a key tool in the
interpretation of experimental results.
151

Chapter 7. Future Directions

We would like to continue our studies of phosphate hydrolysis to extend the
reaction mechanism to the biologically also very important phosphate diesters,
which participate in DNA and RNA hydrolysis. Here we wish to examine the
reaction of an attacking OH
-
, compared to the attacking H
2
O molecule. We
already have very promising preliminary results for the LFER in this case. Also,
we plan to investigate an enzymatic system where such LFER was measured for a
series of compounds (Zalatan and Herschlag 2006). This future study on one hand
would test the performance of the QM/MM-FEP method and on the other hand it
is an ideal system for studying entropic effects, since in this case large entropy
contributions were measured (Schroeder, Lad et al. 2006), while the same reaction
with water as an attacking group had very small entropy change. The entropy is
calculated by defining a classical EVB system that mimics the ab initio structures,
but allows one for a much better sampling due to the lower cost of the
calculations. We use a quasiharmonic approximation for the calculation of the
entropy for high energies, while applying a relatively large force constant to keep
the molecule in place. Then in a subsequent simulation we release this constraint
that kept the molecule in place and use a FEP procedure to account for the free
energy change corresponding to this restraint release (Xiang and Warshel 2006).
152
We also plan to conduct a systematic study using the FDFT approach to
generate diabatic and adiabatic potential surfaces for a series of molecules to
reproduce the LFER type of relationships. In this study the change in the off-
diagonal terms would be of significant interest. The FDFT method will allow us
to study the properties of Marcus parabolas of the different systems that exhibit
the LFERs. We can also examine parameters other than the off-diagonal elements
of these Marcus parabolas, such as the reorganization energy or the reaction free
energy and study the effects of the solvent or protein system on these parameters.
We also plan to implement the gradients of the variational Langevin Dipoles
model that have been derived during this PhD work. This work is already in
progress using the Gaussian 98 program package.
Further implementations of the QM/MM and the QM/MM-FEP protocols are
also underway in the Q-Chem program (Shao, Molnar et al. 2006).

153

Bibliography

Admiraal, S. J. and D. Herschlag (1995). "Mapping of the Transition State for
ATP Hydrolysis: Implications for Enzymatic Catalysis." Chem. Biol. 2:
729.

Admiraal, S. J. and D. Herschlag (2000). "The Substrate-Assisted General Base
Catalysis Model for Phosphate Monoester Hydrolysis: Evaluation Using
Reactivity Comparisons." J. Am. Chem. Soc. 122(10): 2145-2148.

Ahmadian, M. R., P. Stege, et al. (1997). "Confirmation of the Arginine-Finger
Hypothesis for the GAP-Stimulated GTP-Hydrolysis Reaction of Ras."
Nat. Struct. Biol. 4(9): 686-689.

Akola, J. and R. O. Jones (2003). "ATP hydrolysis in water - A density functional
study." J. Phys. Chem. B 107(42): 11774-11783.

Albery, J. K., M. M. (1978). "Methyl transfer reactions." Adv. Phys. Org. Chem.
16: 87-157.

Albery, W. J. (1980). "The Application of the Marcus Relation to Reactions in
Solution." Annu. Rev. Phys. Chem. 31: 227-263.

Allin, C., M. R. Ahmadian, et al. (2001). "Monitoring the GAP catalyzed H-Ras
GTPase reaction at atomic resolution in real time." Proc. Nat. Acad. Sci.
USA 98: 7754-7759.

Allin, C. and K. Gerwert (2001). "Ras Catalyzes GTP Hydrolysis by Shifting
Negative Charges from - to -Phosphate As Revealed by Time-Resolved
FTIR Difference Spectroscopy." Biochemistry 40: 3037-3046.

Amovilli, C., B. Mennucci, et al. (1998). "MCSCF Study of the SN2 Menshutkin
Reaction in Aqueous Solution within the Polarizable Continuum Model."
J. Phys. Chem. B 102(16): 3023-3028.

Anderson, V. E., A. G. Cassano, et al. (2006). Isotope Effects in Chemistry and
Biology. Boca Raton, Taylor&Francis.
154
Aqvist, J. (1996). "Calculation of Absolute Binding Free Energies for Charged
Ligands and Effects of Long-Range Electrostatic Interactions." J. Comp.
Chem. 17(14): 1587-1597.

Aqvist, J., K. Kolmodin, et al. (1999). "Mechanistic Alternatives in Phosphate
Monoester Hydrolysis: What Conclusions Can Be Drawn From Available
Experimental Data?" Chem. & Biol. 6(3): R71-R80.

Aqvist, J. and A. Warshel (1993). "Simulation of Enzyme Reactions Using
Valence Bond Force Fields and Other Hybrid Quantum/Classical
Approaches." Chem. Rev. 93: 2523-2544.

Assfeld, X. and J. L. Rivail (1996). "Quantum chemical computations on parts of
large molecules: The ab initio local self consistent field method."
Chemical Physics Letters 263(1-2): 100-106.

Bakowies, D., Thiel, W. (1996). "Hybrid Models for Combined Quantum
Mechanical and Molecular Approaches." J. Phys. Chem 100: 10580-
10594.

Barnes, J. A., J. Wilkie, et al. (1994). "Transition-state structural variation and
mechanistic change." J. Chem. Soc. Faraday Trans. 90: 1709 - 1714.

Barone, V. and M. Cossi (1998). "Quantum calculation of molecular energies and
energy gradients in solution by a conductor solvent model." J. Phys.
Chem. A 102(11): 1995-2001.

Basaif, S. A., A. M. Davis, et al. (1989). "Effective Charge-Distribution for
Attack of Phenoxide Ion on Aryl Methyl Phosphate Monoanion - Studies
Related to the Action of Ribonuclease." Journal of Organic Chemistry
54(23): 5483-5486.

Becke, A. D. (1988). "Density-Functional Exchange-Energy Approximation With
Correct Asymptotic-Behavior." Physical Review A 38(6): 3098-3100.

Bentzien, J., J. Florian, et al. (1998). "Quantum mechanical-molecular mechanical
approaches for studying chemical reactions in proteins and solution." ACS
Symposium Series 712(Combined Quantum Mechanical and Molecular
Mechanical Methods): 16-34.
Besler, B. H., K. M. J. Merz, et al. (1990). "Atomic charges derived from
semiempirical methods." Journal of Computational Chemistry 11(4): 431-
439.
155

Bohac, M., Y. Nagata, et al. (2002). "Halide-Stabilizing Residues of Haloalkane
Dehalogenases Studied by Quantum Mechanic Calculations and Site-
Directed Mutagenesis." Biochemistry 41(48): 14272-14280.

Born, M. (1920). "Volumen and Hydrationswärme der Ionen." Z. Phys. 1: 45.

Bosma, T., M. G. Pikkemaat, et al. (2003). "Steady-State and Pre-Steady-State
Kinetic Analysis of Halopropane Conversion by a Rhodococcus
Haloalkane Dehalogenase." Biochemistry 42(26): 8047-8053.

Bourne, H. R. (1997). "G proteins - The arginine finger strikes again." Nature
389: 673-674.

Bourne, N. and A. Williams (1984). "Effective Charge on Oxygen in Phosphoryl
(-Po3-2-) Group Transfer from an Oxygen Donor." Journal of Organic
Chemistry 49(7): 1200-1204.

Bruice, T. C. (2002). "A View at the Millennium: the Efficiency of Enzymatic
Catalysis." Acc. Chem. Res. 35(3): 139-148.

Chu, Z. T., J. Villa, et al. (2004). "MOLARIS version beta9.05". University of
Southern California, Los-Angeles.

Cioslowski, J. and E. D. Fleischmann (1991). "Assessing molecular similarity
from results of ab initio electronic structure calculations." J. Am. Chem.
Soc. 113(1): 64-67.

Cleland, W. W. and A. C. Hengge (2006). "Enzymatic Mechanisms of Phosphate
and Sulfate Transfer." Chem. Rev. 106(8): 3252-3278.

Cortona, P. (1991). "Self-Consistently Determined Properties of Solids Without
Band-Structure Calculations." Phys. Rev. B 44: 8454-8458.

Cox, J. R., Jr. and O. B. Ramsay (1964). "Mechanisms of Nucleophilic
Substitution in Phosphate Esters." Chem. Rev. 64: 317-352.

Crespo, A., M. A. Marti, et al. (2005). "Multiple-Steering QM-MM Calculation of
the Free Energy Profile in Chorismate Mutase." J. Am. Chem. Soc.
127(19): 6940-6941.

156
Cui, Q., M. Elstner, et al. (2001). "A QM/MM Implementation of the Self-
Consistent Charge Density Functional Tight Binding (SCC-DFTB)
Method." J. Phys. Chem. B 105: 569-585.

Davis, A. M., A. D. Hall, et al. (1988). "Charge Description of Base-Catalyzed
Alcoholysis of Aryl Phosphodiesters - a Ribonuclease Model." J. Am.
Chem. Soc. 110(15): 5105-5108.

Devi-Kesavan, L. S. and J. Gao (2003). "Combined QM/MM Study of the
Mechanism and Kinetic Isotope Effect of the Nucleophilic Substitution
Reaction in Haloalkane Dehalogenase." J. Am. Chem. Soc. 125(6): 1532-
1540.

Devi-Kesavan, L. S., M. Garcia-Viloca, et al. (2003). "Semiempirical QM/MM
potential with simple valence bond (SVB) for enzyme reactions.
Application to the nucleophilic addition reaction in haloalkane
dehalogenase." Theoretical Chemistry Accounts: Theory, Computation,
and Modeling (Theoretica Chimica Acta) 109(3): 133.

Di Sabato, G. and W. P. Jencks (1961). "Mechanism and catalysis of reactions of
acyl phosphates. II. Hydrolysis." J. Am. Chem. Soc. 83: 4400 - 4405.

Dittrich, M., S. Hayashi, et al. (2003). "On the mechanism of ATP hydrolysis in
F-1-ATPase." Biophysical Journal 85(4): 2253-2266.

Du, X., G. E. Black, et al. (2004). "Kinetic isotope effects in Ras-catalyzed GTP
hydrolysis: Evidence for a loose transition state." PNAS 101(24): 8858-
8863.

Fersht, A. R., R. J. Leatherbarrow, et al. (1986). "Quantitative Analysis of
Structure-Activity Relationships in Engineered Proteins by Linear Free-
Energy Relationships." Nature 322: 284.

Field, M. (2002). "Stimulating Enzyme Reactions: Challenges and Perspectives."
J. Comp. Chem. 23: 48-58.

Field, M. J., P. A. Bash, et al. (1990). "A Combined Quantum Mechanical and
Molecular Mechanical Potential for Molecular Dynamics Simulations." J.
Comp. Chem. 11: 700-733.

Florián, J., J. Åqvist, et al. (1998). "On the Reactivity of Phosphate Monoester
Dianions In Aqueous Solution: Brønsted Linear Free-Energy
157
Relationships Do Not Have an Unique Mechanistic Interpretation." J. Am.
Chem. Soc. 120: 11524 -11525.

Florian, J. and A. Warshel (1997). "Langevin Dipoles Model for Ab Initio
Calculations of Chemical Processes in Solution: Parameterization and
Application to Hydration Free Energies of Neutral and Ionic Solutes, and
Conformational Analysis in Aqueous Solution." J. Phys. Chem B 101(28):
5583-5595.

Florián, J. and A. Warshel (1997). "A fundamental assumption about OH
-
attack
in phosphate ester hydrolysis is not fully justified." J. Am. Chem. Soc.
119: 5473 - 5474.

Florián, J. and A. Warshel (1998). "Phosphate Ester Hydrolysis in Aqueous
Solution: Associative Versus Dissociative Mechanisms." J. Phys. Chem. B
102: 719-734.

Friesner, R. and M. D. Beachy (1998). "Quantum Mechanical Calculations on
Biological Systems." Curr. Op. Struct. Biol. 8: 257-262.

Frisch, M. J., G. W. Trucks, et al. (1998). Gaussian 98. Pittsburgh, PA, Gaussian,
Inc.

Gao, J. (1996). "Hybrid quantum and molecular mechanical simulations: an
alternative avenue to solvent effects in organic chemistry." Acc. Chem.
Res. 29: 298-305.

Garcia-Viloca, M., A. Gonzalez-Lafont, et al. (2001). "A QM/MM study of the
Racemization of Vinylglycolate Catalysis by Mandelate Racemase
Enzyme." J. Am. Chem. Soc 123: 709-721.

Gillan, M. J. (1987). "Quantum-Classical Crossover of the Transition Rate in the
Damped Double Well." J. Phys. C. Solid State Phys. 20: 3621-3641.

Girones, X., R. Carbo-Dorca, et al. (2003). "Molecular Basis of LFER. Modeling
of the Electronic Substituent Effect Using Fragment Quantum Self-
Similarity Measures." J. Chem. Inf. Comput. Sci. 43: 2033-2038.

Glennon, T. M., J. Villa, et al. (2000). "How does GAP catalyze the GTPase
reaction of Ras? A computer simulation study." Biochemistry 39: 9641-
9651.

158
Grigorenko, B. L., A. V. Nemukhin, et al. (2005). "QM/MM modeling the Ras-
GAP catalyzed hydrolysis of guanosine triphosphate." Proteins-Structure
Function and Bioinformatics 60(3): 495-503.

Grigorenko, B. L., A. V. Rogov, et al. (2006). "Mechanism of triphosphate
hydrolysis in aqueous solution: QM/MM simulations in water clusters."
Journal of Physical Chemistry B 110(9): 4407-4412.

Guthrie, J. P. (1977). "Hydration and Dehydration of Phosphoric Acid
Derivatives: Free Energies of Formation of the Pentacoordinate
Intermediates of Phosphate Ester Hydrolysis and of Monomeric
Metaphosphate." J. Am. Chem. Soc. 99: 3991-4001.

Gwaltney, S. R., S. V. Rosokha, et al. (2003). "Charge-Transfer Mechanism for
Electrophilic Aromatic Nitration and Nitrosation via the Convergence of
(ab Initio) Molecular-Orbital and Marcus-Hush Theories with
Experiments." J. Am. Chem. Soc. 125(11): 3273-3283.

Hengge, A. (2005). "Mechanistic studies on enzyme-catalyzed phosphoryl
transfer." Advances in Physical Organic Chemistry 40: 49.

Hill, J. R., A. Tokmakoff, et al. (1994). "Vibrational dynamics of carbon
monoxide at the active site of myoglobin: picosecond infrared free-
electron laser pump-probe experiments." J. Phys. Chem. 98: 11213-11219.

Hoff, R. H. and A. C. Hengge (1998). "Entropy and enthalpy contributions to
solvent effects on phosphate monoester solvolysis. The importance of
entropy effects in the dissociative transition state." Journal of Organic
Chemistry 63(19): 6680-6688.

Hong, G., M. Strajbl, et al. (2000). "Constraining the electron densities in DFT
method as an effective way for ab initio studies of metal catalyzed
reactions." J. Comp. Chem. 21: 1554-1561.
Hong, G. Y., E. Rosta, et al. (2006). "Using the constrained DFT approach in
generating diabatic surfaces and off diagonal empirical valence bond
terms for modeling reactions in condensed phases." Journal of Physical
Chemistry B 110(39): 19570-19574.

Hu, C.-H. and T. Brinck (1999). "Theoretical studies of the Hydrolysis of the
Methyl Phosphate Anion." J. Phys. Chem. A 103: 5379-5386.

159
Hu, H. and W. Yang (2005). "Dual-topology/dual-coordinate free-energy
simulation using QM/MM force field." The Journal of Chemical Physics
123(4): 041102.

Hummer, G. and A. Szabo (1996). "Calculation of free energy differences from
computer simulations of initial and final states." J. Chem. Phys. 105(5):
2004-2010.

Hush, N. S. (1961). "Adiabatic theory of outer sphere electron-transfer reactions
in solution." Trans. Faraday Soc. 57: 557.

Hush, N. S. (1968). "Homogeneous and heterogeneous optical and thermal
electron transfer." Electrochimica Acta 13(5): 1005.

Hwang, J. K. and A. Warshel (1987). "Microscopic examination of free-energy
relationships for electron transfer in polar solvents." Journal of the
American Chemical Society 109(3): 715-20.

Iche-Tarrat, N., J. C. Barthelat, et al. (2005). "Theoretical studies of the
hydroxide-catalyzed P-O cleavage reactions of neutral phosphate triesters
and diesters in aqueous solution: Examination of the changes induced by
H/Me substitution." Journal of Physical Chemistry B 109(47): 22570-
22580.

Iftimie, R., D. Salahub, et al. (2003). "An efficient Monte Carlo method for
calculating ab initio transition state theory reaction rates in solution." J.
Chem. Phys. 119(21): 11285-11297.

Ishida, T. and S. Kato (2003). "Theoretical Perspectives on the Reaction
Mechanism of Serine Proteases: The Reaction Free Energy Profiles of the
Acylation Process." J. Am. Chem. Soc. 125(39): 12035-12048.

Jencks, W. P. (1985). Chem. Rev. 85: 511-527.

Jencks, W. P. (1987). Catalysis in Chemistry and Enzymology. New York, Dover.

Jencks, W. P. and J. Regenstein (1976). Handbook of Biochemistry and
Molecular Biology. Cleveland, CRC Press.

Kalra, A., S. Garde, et al. (2003). "Osmotic water transport through carbon
nanotube membranes." Proc. Nat. Acad. Sci. 100: 10175-10180.

160
Kastner, J., H. M. Senn, et al. (2006). "QM/MM Free-Energy Perturbation
Compared to Thermodynamic Integration and Umbrella Sampling:
Application to an Enzymatic Reaction." J. Chem. Theory Comput. 2(2):
452-461.

King, G. and A. Warshel (1989). "A Surface Constrained All-Atom Solvent
Model for Effective Simulations of Polar Solutions." J. Chem. Phys.
91(6): 3647-3661.

King, G. and A. Warshel (1990). "Investigation of the Free Energy Functions for
Electron Transfer Reactions." J. Chem. Phys. 93: 8682-8692.

Kirby, A. J. and W. P. Jenks (1965). "The Reactivity of Nucleophilic Reagents
toward the p-Nirophenyl Phosphate Dianions." J. Am. Chem. Soc. 87:
3209.

Kirby, A. J. and A. G. Varvoglis (1967). "The Reactivity of Phosphate Esters.
Monoester Hydrolysis." J. Am. Chem. Soc. 89: 415-423.

Klahn, M., S. Braun-Sand, et al. (2005). "On possible pitfalls in ab initio quantum
mechanics/molecular mechanics minimization approaches for studies of
enzymatic reactions." Journal of Physical Chemistry B 109(32): 15645-
15650.

Klahn, M., E. Rosta, et al. (2006). "On the Mechanism of Hydrolysis of
Phosphate Monoesters Dianions in Solutions and Proteins." J. Am. Chem.
Soc. 128(47): 15310-15323.
Klahn, M., J. Schlitter, et al. (2005). "Theoretical IR spectroscopy based on
QM/MM calculations provides changes in charge distribution, bond
lengths, and bond angles of the GTP ligand induced by the Ras-protein."
Biophysical Journal 88(6): 3829-3844.

Klamt, A. and G. Schuurmann (1993). "Cosmo - a New Approach to Dielectric
Screening in Solvents with Explicit Expressions for the Screening Energy
and Its Gradient." Journal of the Chemical Society-Perkin Transactions
2(5): 799-805.

Kluner, T., N. Govind, Y. A. Wang, E. A. Carter (2001). "Prediction of Electronic
Excited STates of Adsorbates on Metal Surfaces from First Principles."
Phys. Rev. Lett. 86(26): 5954.

161
Kmunicek, J., K. Hynkova, et al. (2005). "Quantitative Analysis of Substrate
Specificity of Haloalkane Dehalogenase LinB from Sphingomonas
paucimobilis UT26." Biochemistry 44(9): 3390-3401.

Kmunicek, J., S. Luengo, et al. (2001). "Comparative Binding Energy Analysis of
the Substrate Specificity of Haloalkane Dehalogenase from Xanthobacter
autotrophicus GJ10." Biochemistry 40(37): 11288-11288.

Kollman, P. (1993). "Free Energy Calculations: Applications to Chemical and
Biochemical Phenomena." Chem. Rev. 93: 2395-2417.

Kotting, C. and K. Gerwert (2004). "Time-resolved FTIR studies provide
activation free energy, activation enthalpy and activation entropy for
GTPase reactions." Chemical Physics 307(2-3): 227.

Kreevoy, M. M. and A. T. Kotchevar (1990). "Dynamics of Hydride Transfer
Between NAD
+
Analogues." J. Am. Chem. Soc. 112: 3579.

Kuharski, R. A., J. S. Bader, et al. (1988). "Molecular model for aqueous ferrous
ferric electron transfer." J. Chem. Phys. 89: 3248-3257.

Lad, C., N. H. Williams, et al. (2003). "The rate of hydrolysis of
phosphomonoester dianions and the exceptional catalytic proficiencies of
protein and inositol phosphatases." Proceedings of the National Academy
of Sciences of the United States of America 100(10): 5607-5610.

Langen, R., T. Schweins, et al. (1992). "On the Mechanism of Guanosine
Triphosphate Hydrolysis in ras p21 Proteins." Biochemistry 31: 8691-
8696.

Lau, E. Y., K. Kahn, et al. (2000). "The importance of reactant positioning in
enzyme catalysis: A hybrid quantum mechanics/molecular mechanics
study of a haloalkane dehalogenase." Proc. Natl. Acad. Sci. USA 97(18):
9937-9942.

Lee, F. S., Z. T. Chu, et al. (1993). "Microscopic and Semimicroscopic
Calculations of Electrostatic Energies in Proteins by the POLARIS and
ENZYMIX Programs." J. Comp. Chem. 14: 161-185.

162
Lee, F. S. and A. Warshel (1992). "A Local Reaction Field Method for Fast
Evaluation of Long-Range Electrostatic Interactions in Molecular
Simulations." J. Chem. Phys. 97: 3100-3107.

Leffler, J. E. and E. Grundwald (1963). Rates and Equilibria of Organic
Reactions. New York, Wiley: 156-161.

Lewandowicz, A., J. Rudzinski, et al. (2001). "Chlorine Kinetic Isotope Effects
on the Haloalkane Dehalogenase Reaction." J. Am. Chem. Soc. 123(19):
4550-4555.

Li, G. and Q. Cui (2004). "Mechanochemical Coupling in Myosin: A Theoretical
Analysis with Molecular Dynamics and Combined QM/MM Reaction
Path Calculations." J. Phys. Chem. B 108(10): 3342-3357.

Liu, W. B., R. H. Wood, et al. (2003). "Density and temperature dependences of
hydration free energy of Na+ and Cl- at supercritical conditions predicted
by ab initio/classical free energy perturbation." J. Phys. Chem. B 107(35):
9505-9513.

Liu, Y., B. A. Gregersen, et al. (2006). "Transesterification Thio Effects of
Phosphate Diesters: Free Energy Barriers and Kinetic and Equilibrium
Isotope Effects from Density-Functional Theory." Biochemistry 45(33):
10043-10053.
Liu, Y., X. Lopez, et al. (2005). "Kinetic isotope effects on thio-substituted
biological phosphoryl transfer reactions from density-functional theory."
Chemical Communications(31): 3909-3911.

Lopez, X., A. Dejaegere, et al. (2006). "Nucleophilic Attack on Phosphate
Diesters: A Density Functional Study of In-Line Reactivity in Dianionic,
Monoanionic, and Neutral Systems." J. Phys. Chem. B 110(23): 11525-
11539.

Lu, N., J. K. Singh, et al. (2003). "Appropriate methods to combine forward and
reverse free-energy perturbation averages." The Journal of Chemical
Physics 118(7): 2977-2984.

Maegley, K. A., S. J. Admiral, et al. (1996). "Ras-catalyzed hydrolysis of GTP: A
new perspective from model studies." Proc. Natl. Acad. Sci., USA 93:
8160-8166.

163
Marcus, R. A. (1964). "Chemical and electrochemical electron transfer theory."
Ann. Rev. Phys. Chem. 15: 155-196.

Marcus, R. A. (1993). "Electron-transfer reactions in chemistry - theory and
experiment (Nobel lecture)." Angew. Chem. Int. Ed. Engl. 32: 1111-1121.

Marti, S., V. Moliner, et al. (2005). "Improving the QM/MM Description of
Chemical Processes: A Dual Level Strategy To Explore the Potential
Energy Surface in Very Large Systems." J. Chem. Theory Comput. 1(5):
1008-1016.

Maseras, F. and K. Morokuma (1995). "IMOMM: A new integrated ab initio +
molecular mechanics geometry optimization scheme of equilibrium
structures and transition states." Journal of Computational Chemistry
16(9): 1170-1179.

Maulitz, A. H., F. C. Lightstone, et al. (1997). "Nonenzymatic and Enzymatic
Hydrolysis of Alkyl Halides: A Theoretical Study of the S
N
2 Reactions of
Acetate and Hydroxide Ions with Alkyl Chlorides." Proc. Natl. Acad. Sci.
USA 94: 6591-6595.

McQuarrie, D. A. (1976). Statistical Mechanics. New York, Harper and Row.
Mezey, P. G. (1987). "The shape of molecular charge distributions: Group theory
without symmetry." Journal of Computational Chemistry 8(4): 462-469.

Mierts, S., E. Scrocco, et al. (1981). "Electrostatic interaction of a solute with a
continuum. A direct utilizaion of AB initio molecular potentials for the
prevision of solvent effects." Chemical Physics 55(1): 117.

Mihai, C., A. V. Kravchuk, et al. (2003). "Application of Bronsted-Type LFER in
the Study of the Phospholipase C Mechanism." J. Am. Chem. Soc.
125(11): 3236-3242.

Mildvan, A. S. (1997). "Mechanisms of signaling and related enzymes." Proteins
29(4): 401-416.

Mo, Y. and J. Gao (2000). J. Chem. Phys. 104: 3012-3020.

Monard, G. and K. M. Merz (1999). "Combined Quantum Mechanical/Molecular
Mechanical Methodologies Applied to Biomolecular Systems." Acc.
Chem. Res. 32: 904-911.

164
More O'Ferrall, R. A. (1970). "Relationships between E2 and E1cB mechanisms
of β-elimination." J. Chem. Soc. B: 274 - 277.

Muegge, I., T. Schweins, et al. (1996). "Electrostatic Control of GTP and GDP
Binding in the Oncoprotein p21 ras." Structure 4: 475-489.

Muller, R. P. and A. Warshel (1995). "Ab Initio Calculations of Free Energy
Barriers for Chemical Reactions in Solution." J. Phys. Chem. 99: 17516-
17524.

Nam, K., X. Prat-Resina, et al. (2004). "Dynamics of an Enzymatic Substitution
Reaction in Haloalkane Dehalogenase." J. Am. Chem. Soc. 126(5): 1369-
1376.

Olsson, M. H. M., G. Hong, et al. (2003). "Frozen Density Functional Free
Energy Simulations of Redox Proteins: Computational Studies of the
Reduction Potential of Plastocyanin and Rusticyanin." J. Am. Chem. Soc.
125: 5025-5039.
Olsson, M. H. M., W. W. Parson, et al. (2006). "Dynamical Contributions to
Enzyme Catalysis: Critical Tests of A Popular Hypothesis." Chem. Rev.
106(5): 1737-1756.

Olsson, M. H. M. and A. Warshel (2004). "Solute Solvent Dynamics and
Energetics in Enzyme Catalysis: The S
N
2 Reaction of Dehalogenase as a
General Benchmark." J. Am. Chem. Soc. 126(46): 15167-15179.

Paneth, P. (2003). "Chlorine Kinetic Isotope Effects on Enzymatic
Dehalogenations." Acc. Chem. Res. 36(2): 120-126.

Pecoraro, V. L., J. D. Hermes, et al. (1984). "Stability-Constants of Mg-2+ and
Cd-2+ Complexes of Adenine-Nucleotides and Thionucleotides and Rate
Constants for Formation and Dissociation of Mgatp and Mgadp."
Biochemistry 23(22): 5262-5271.

Perdew, J. P. (1986). "Density-Functional Approximation For The Correlation-
Energy Of The Inhomogeneous Electron-Gas." Physical Review B 33(12):
8822-8824.

Pradipta, B. (2005). "Accelerating quantum mechanical/molecular mechanical
sampling using pure molecular mechanical potential as an importance
function: The case of effective fragment potential." The Journal of
Chemical Physics 122(9): 091102.
165

Pries, F., J. Kingma, et al. (1995). "Histidine 289 is essential for hydrolysis of the
alkyl-enzyme intermediate of haloalkane dehalogenase." J. Biol. Chem.
270(18): 10405-11.

Purcell, E. M. (1965). Electricity and Magnetism, McGraw-Hill.

Riccardi, D., P. Schaefer, et al. (2005). "pK
a
Calculations in Solution and Proteins
with QM/MM Free Energy Perturbation Simulations: A Quantitative Test
of QM/MM Protocols." J. Phys. Chem. B 109(37): 17715-17733.

Rod, T. H. and U. Ryde (2005). "Quantum mechanical free energy barrier for an
enzymatic reaction." Phys. Rev. Lett. 94(13): 138302.
Rosta, E., M. Klahn, et al. (2006). "Towards accurate ab initio QM/MM
calculations of free-energy profiles of enzymatic reactions." Journal of
Physical Chemistry B 110(6): 2934-2941.

Rosta, E., M. Klahn, et al. (2006). "Towards Accurate Ab Initio QM/MM
Calculations of Free-Energy Profiles of Enzymatic Reactions." J. Phys.
Chem. B 110(6): 2934-2941.

Rosta, E. and A. Warshel (2003). "Analytical derivatives for langevin dipoles
solvation model." Abstracts of Papers of the American Chemical Society
225: U771-U771.

Rosta, E. and A. Warshel (2004). "Ab initio QM/MM simulation with proper
sampling: A study of the reaction of orotidine-5 '-monophosphate
decarboxylase." Abstracts of Papers of the American Chemical Society
227: U320-U320.

Rosta, E. and A. Warshel (2005). "Ab initio free energies of enzymatic reactions
using model potentials." Abstracts of Papers of the American Chemical
Society 229: U743-U743.

Sakane, S., E. M. Yezdimer, et al. (2000). "Exploring the ab initio/classical free
energy perturbation method: The hydration free energy of water." The
Journal of Chemical Physics 113(7): 2583.

Schmitt, U. W. and G. A. Voth (1998). "Multistate Empirical Valence Bond
Model for Proton Transport in Water." J. Phys. Chem. B 102: 5547-5551.

166
Schroeder, G. K., C. Lad, et al. (2006). "The time required for water attack at the
phosphorus atom of simple phosphodiesters and of DNA." Proc. Nat.
Acad. Sci. USA 103(11): 4052-4055.

Schweins, T., M. Geyer, et al. (1996). "Linear Free Energy Relationships in the
Intrinsic and GTPase Activating Protein-Stimulated Guanosine 5'-
Triphosphate Hydrolysis of p21 ras." Biochemistry 35: 14225-14231.

Schweins, T., M. Geyer, et al. (1995). "Substrate-assisted Catalysis as a
Mechanism for GTP Hydrolysis of p21 ras and Other GTP-binding
Proteins." Nature Struct. Biol. 2(1): 36-44.
Schweins, T., R. Langen, et al. (1994). "Why Have Mutagenesis Studies Not
Located the General Base in ras p21." Nature Struct. Biol. 1(7): 476-484.

Shaik, S. and A. Shurki (1999). "Valence Bond Diagrams and Chemical
Reactivity." Angew. Chem.Int. Ed. Engl. 38: 586-625.

Shaik, S. S., H. B. Schlegel, et al. (1992). Theoretical Aspects of Physical Organic
Chemistry. Application to the SN2 Transition State. NY, Wiley
Interscience.

Sham, Y. Y., Z. T. Chu, et al. (2000). "Examining Methods for Calculations of
Binding Free Energies: LRA, LIE, PDLD-LRA, and PDLD/S-LRA
Calculations of Ligands Binding to an HIV Protease." Proteins: Struct.
Funct. Genet. 39: 393-407.

Shao, Y., L. F. Molnar, et al. (2006). "Advances in methods and algorithms in a
modern quantum chemistry program package." Physical Chemistry
Chemical Physics 8(27): 3172-3191.

Shurki, A. (2006). "Valence Bond – Rebirth of the Phoenix or Relic from the
Stone Age." Theoretical Chemistry Accounts: Theory, Computation, and
Modeling (Theoretica Chimica Acta) 116(1): 253.

Shurki, A., M. Strajbl, et al. (2004). Electrostatic basis for bioenergetics.
Energetics of Biological Macromolecules, Pt E. 380: 52-84.

Shurki, A., M. Štrajbl, et al. (2002). "How Much Do Enzymes Really Gain by
Restraining Their Reacting Fragments?" J. Am. Chem. Soc. 124: 4097-
4107.

167
Shurki, A. and A. Warshel (2003). "Structure/Function Correlations of Protreins
using MM, QM/MM and Related Approaches; Methods, Concepts, Pitfalls
and Current Progress." Advances in Protein Chemistry 66: 249-312.

Shurki, A. and A. Warshel (2004). "Why does the Ras switch "break" by
oncogenic mutations?" Proteins: Struct. Funct. Bioinf. 55(1): 1-10.

Sigel, H., E. M. Bianchi, et al. (2001). "Acid-base properties of the 5 '-
triphosphates of guanosine and inosine (GTP(4-) and ITP4-) and of several
related nucleobase derivatives." Journal of the Chemical Society-Perkin
Transactions 2(4): 507-511.

Singh, U. C. and P. A. Kollman (1984). "An approach to computing electrostatic
charges for molecules." Journal of Computational Chemistry 5(2): 129-
145.

Soriano, A., E. Silla, et al. (2005). "Dynamic and Electrostatic Effects in
Enzymatic Processes. An Analysis of the Nucleophilic Substitution
Reaction in Haloalkane Dehalogenase." J. Am. Chem. Soc. 127(6): 1946-
1957.

Sprang, S. R. (1997). "G PROTEIN MECHANISMS: Insights from Structural
Analysis." Annu. Rev. Biochem. 66: 639-678.

St-Amant, A. and D. R. Salahub (1990). "New algorithm for the optimization of
geometries in local density functional theory." Chemical Physics Letters
169(5): 387.

Strajbl, M., J. Florian, et al. (2000). "Ab initio/LD studies of chemical reactions in
solution: Reference free-energy surfaces for acylation reactions occuring
in serine and cysteine proteases." Int. J. Quantum. Chem. 77: 44 -53.

Strajbl, M., J. Florian, et al. (2001). "Ab Initio Evaluation of the Free Energy
Surfaces for the General Base/Acid Catalyzed Thiolysis of Formamide
and the Hydrolysis of Methyl Thioformate: A Reference Solution
Reaction for Studies of Cysteine Proteases." J. Phys. Chem. B. 105: 4471-
4484.

168
Strajbl, M., G. Hong, et al. (2002). "Ab-initio QM/MM Simulation with Proper
Sampling: "First Principle" Calculations of the Free Energy of the Auto-
dissociation of Water in Aqueous Solution." J. Phys. Chem. B 106:
13333-13343.

Strajbl, M., Y. Y. Sham, et al. (2000). "Calculations of Activation Entropies of
Chemical Reactions in Solution." J. Phys. Chem. B 104: 4578-4584.

Strajbl, M., A. Shurki, et al. (2003). "Converting conformational changes to
electrostatic energy in molecular motors: The energetics of ATP
synthase." Proc. Nat. Acad. Sci. USA 100(25): 14834-14839.
Sun, H., K. Kalju, et al. (2003). "Comparison of formation of reactive conformers
for the S
N
2 displacements by CH
3
CO in water and by Asp124-CO in a
haloalkane dehalogenase." Proc. Natl. Acad. Sci. USA 100(5): 2215-2219.

Tomasi, J. (2004). "Thirty years of continuum solvation chemistry: a review, and
prospects for the near future." Theoretical Chemistry Accounts: Theory,
Computation, and Modeling (Theoretica Chimica Acta) 112(4): 184.

Toney, M. D. and J. F. Kirsch (1989). "Direct Bronsted Analysis of the
Restoration of Activity to a Mutant Enzyme by Exogenous Amines."
Science 243: 1485-1488.

Vaidehi, N., T. A. Wesolowski, et al. (1992). "Quantum-Mechanical Calculations
of Solvation Free Energies. A Combined ab initio Pseudopotential Free-
Energy Perturbation Approach." J. Chem. Phys. 97: 4264-4271.

Valleau, J. P. and G. M. Torrie (1977). Modern Theoretical Chemistry. New
York, Plenum Press.

van Gunsteren, W. F. and A. E. Mark (1992). "On the Interpretation of
Biochemical Data by Molecular Dynamics Simulation." Eur. J. Biochem.
204: 947.

Verschueren, K. H. G., S. M. Franken, et al. (1993). "Refined X-ray Structures of
Haloalkane Dehalogenase at pH 6.2 and pH 8.2 and Implications for the
Reaction Mechanism." Journal of Molecular Biology 232(3): 856.

Vetter, I. R. and A. Wittinghofer (1999). "Nucleoside triphosphate-binding
proteins: different scaffolds to achieve phosphoryl transfer." Q Rev
Biophys 32(1): 1-56.

169
Villa, J. and A. Warshel (2001). "Energetics and Dynamics of Enzymatic
Reactions." J. Phys. Chem. B 105: 7887-7907.

Villa, J. and A. Warshel. (2002). "MOLARIS: Version alpha9.06, Theoretical
Background and Practical Examples." from
http://diana.imim.es/~jvilla/molaris/theory_molaris/theory.html.

Voth, G. A. (1996). "Path-integral centroid methods in quantum statistical
mechanics and dynamics." Adv. Chem. Phys. 93: 135-218.

Wagner, G., A. DeMarco, et al. (1976). "Dynamics of the aromatic amino acid
sidechains in the globular conformation of the basic pacreatic trypsin
inhibitor. 1.
1
H-NMR studies." Biochem. Struct. Mech. 2: 139-149.

Wang, Y. N., I. A. Topol, et al. (2003). "Theoretical studies on the hydrolysis of
mono-phosphate and tri-phosphate in gas phase and aqueous solution."
Journal of the American Chemical Society 125(43): 13265-13273.

Warshel, A. (1982). "Dynamics of Reactions in Polar Solvents. Semiclassical
Trajectory Studies of Electron-Transfer and Proton-Transfer Reactions." J.
Phys. Chem. 86: 2218-2224.

Warshel, A. (1984). "Simulating the energetics and dynamics of enzymic
reactions." Pontificiae Academiae Scientiarum Scripta Varia 55(Specif.
Biol. Interact.): 59-81.

Warshel, A. (1991). Computer Modeling of Chemical Reactions in Enzymes and
Solutions. New York, John Wiley & Sons.

Warshel, A. and J. Florian (2004). The Empirical Valence Bond. The
Encyclopedia of Computational Chemistry.

Warshel, A., J. K. Hwang, et al. (1992). "Computer Simulations of Enzymatic
Reactions: Examination of Linear Free-Energy Relationships and
Quantum-Mechanical Corrections in the Initial Proton-transfer Step of
Carbonic Anhydrase." Faraday Discuss. 93: 225.

Warshel, A. and M. Levitt (1976). "Theoretical studies of enzymic reactions:
dielectric, electrostatic and steric stabilization of the carbonium ion in the
reaction of lysozyme." J. Mol. Biol. 103: 227-249.

170
Warshel, A., T. Schweins, et al. (1994). "Linear Free Energy Relationships in
Enzymes. Theoretical Analysis of the Reaction of Tyrosyl-tRNA
Synthetase." J. Am. Chem. Soc. 116: 8437-8442.

Warshel, A., F. Sussman, et al. (1988). "Evaluation of Catalytic Free Energies in
Genetically Modified Proteins." J. Mol. Biol. 201: 139-159.
Warshel, A., F. Sussman, et al. (1986). "Free Energy of Charges in Solvated
Proteins: Microscopic Calculations Using a Reversible Charging Process."
Biochemistry 25: 8368-8372.

Warshel, A. and R. M. Weiss (1980). "An Empirical Valence Bond Approach for
Comparing Reactions in Solutions and in Enzymes." J. Am. Chem. Soc.
102(20): 6218-6226.

Wesolowski, T. A. (1997). "Density functional theory with approximate kinetic
energy functionals applied to hydrogen bonds." The Journal of Chemical
Physics 106(20): 8516-8526.

Wesolowski, T. A. (2002). "Development of Novel Computational Strategies to
Match the Challenges of Supramolecular Chemistry, Biochemistry, and
Materials Science." CHIMIA 56(12): 707.

Wesolowski, T. A. and A. Warshel (1993). "Frozen Density Functional Approach
for Ab Initio Calculations of Solvated Molecules." J. Phys. Chem. 97:
8050-8053.

Westheimer, F. H. (1987). "Why Nature Chose Phosphates." Science 235: 1173-
1178.

Williams, A. (1984). "Effective Charge and Leffler's Index as Mechanistic Tools
for Reactions in Solution." Acc. Chem. Res. 17: 425-430.

Williams, N. H. and P. Wyman (2001). "Base catalysed phosphate diester
hydrolysis." Chemical Communications(14): 1268-1269.

Wolfenden, R. (2006). "Degrees of Difficulty of Water-Consuming Reactions in
the Absence of Enzymes." Chem. Rev. 106(8): 3379-3396.

Xiang, Y. and A. Warshel (2006). Unpublished results.

171
Yang, S. Y., P. Fleurat-Lessard, et al. (2004). "Free Energy Profiles for the
Identity S
N
2 Reactions Cl
-
+ CH
3
Cl and NH
3
+ H
3
BNH
3
: A Constraint Ab
Initio Molecular Dynamics Study." J. Phys. Chem. A 108(43): 9461-9468.

Yang, W., Y. Q. Gao, et al. (2003). "The missing link between thermodynamics
and structure in F1-ATPase." Proc. Nat. Acad. Sci. USA 100(3): 874-879.

Zalatan, J. G. and D. Herschlag (2006). "Alkaline Phosphatase Mono- and
Diesterase Reactions: Comparative Transition State Analysis." J. Am.
Chem. Soc. 128(4): 1293-1303.

Zhang, Y., H. Liu, et al. (2000). "Free energy calculation on enzyme reactions
with an efficient iterative procedure to determine minimum energy paths
on a combined ab initio QM/MM potential energy surface." J. Chem.
Phys. 112(8): 3483-3492.

Zhou, D. M. and K. Taira (1998). "The hydrolysis of RNA: From theoretical
calculations to the hammerhead ribozyme-mediated cleavage of RNA."
Chemical Reviews 98(3): 991-1026.

Zwanzig, R. W. (1954). "High-Temperature Equation of State by a Perturbation
Method. I. Nonpolar Gases." J. Chem. Phys. 22: 1420.

172

Appendix A: The CPHF equations

We would like to find the total derivative of the energy with respect of a
coordinate, such as the atomic coordinates:
1 1 1
2 2 2
s
ij s s s kl s
s ij ij ij ij ij
i j s i j s k l i j s
kl
h
P dE
P h P h P
dr r P r r

α α α α
μ μ
μ
∂
∂ ∂ ∂
= + +
∂ ∂ ∂ ∂
∑∑ ∑∑ ∑ ∑∑
(8.1)
This energy is a functional of some parameters, such as in our case the
dipoles, which have a dependence on the density matrix and thus on the atomic
coordinates implicitly. Let us collect these nonvariational parameters,
, 1
m
t m N = … :
ˆ
( , , ) E E r t H
α
=
(8.2)
The equations that define the parameters can also be listed:
ˆ
( , , ) 0 1,...,
m
r t H m N
α
Ω = =
(8.3)
In our case these equations are the Langevin equations defining the dipoles
based on the electrostatic field. The total derivatives can be expressed using the
Hellman-Feynman theorem:
Hellman-Feynman
ˆ ˆ
ˆ ˆ
t m
m
m
t dE E H E E H
E t
dr r t r r H H

α
α α α α
∂ ∂ ∂ ∂ ∂ ∂
= + = +
∂ ∂ ∂ ∂ ∂ ∂
∑
(8.4)
173
and the
m
t
r
α
∂
∂
derivatives are unknown. These can be obtained by
differentiating their defining equations:

ˆ
0 1,...,
ˆ
mi
m i
m m m i
i
i
A
t
d t H
m N
dr r t r H

α α
α α α
Ω
Ω ∂Ω ∂Ω ∂ ∂
= + = =
∂ ∂ ∂ ∂
∑
(8.5)
We can write the above equation in a matrix formalism:
0 0
m mi i
i
A t At
α α α α
= Ω + ⇒ = Ω +
∑
(8.6)
Using the above obtained matrix, A, we can now complete the total derivative
of the energy:
1
ˆ
ˆ
t
dE E H
E A
dr r H
α
α α
−
∂ ∂
= − Ω
∂ ∂
(8.7)
Here we can introduce the Z-vector, which depends only on the method, but
not the r
α
.
1 t
Z E A
−
=
(8.8)
So our working formula becomes:
ˆ
ˆ
dE E H
Z
dr r H
α
α α
∂ ∂
= − Ω
∂ ∂
. (8.9)

Abstract (if available)

Abstract Modeling chemical reactions in condensed phases, especially in biological environment, such as in enzymes, is a major task of computational chemistry. At the current computer resources one can usually not afford to solve the Schrödinger equation for the molecules, and not even for all the electrons of the system. Simplified classical physical models are used for many of the processes. For chemical reactions, however these classical models are insufficient and quantum treatment is necessary for the atoms participating in the reactions.

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

High level ab initio computational studies on carbocations and related intermediates

PDF

Advancing ab initio QM/MM free energy calculations: refining, validating and quantifying the reference potential approach

PDF

Development of robust ab initio methods for description of excited states and autoionizing resonances

PDF

Electronic structure and spectroscopy in the gas and condensed phase: methodology and applications

PDF

Photoinduced redox reactions in biologically relevant systems

PDF

Modeling x-ray spectroscopy in condensed phase

PDF

Electronically excited and ionized states in condensed phase: theory and applications

PDF

Understanding electrostatic effects in the function of biological systems

PDF

Quantum mechanical description of electronic and vibrational degrees of freedom in laser-coolable molecules

PDF

Advancing computational methods for free energy calculations in proteins and the applications for studies of catalysis in nucleotide-associated proteins

PDF

Explorations in the use of quantum annealers for machine learning

PDF

The effect of microhydration on ionization energy and proton transfer in nucleobases: analysis and method development

PDF

Quantum molecular dynamics and machine learning for non-equilibrium processes in quantum materials

PDF

Simulations across scales: insights into biomolecular, mechanisms, magnetic materials, and optical processes

PDF

Adiabatic and non-adiabatic molecular dynamics in nanoscale systems: theory and applications

$Study of rotation and phase separation in ³He, ⁴He, and mixed ³He/⁴He droplets by X-ray diffraction$

PDF

Study of rotation and phase separation in ³He, ⁴He, and mixed ³He/⁴He droplets by X-ray diffraction

Asset Metadata

Creator Rosta, Edina (author)

Core Title Ab initio methodologies in studying enzymatic reactions

School College of Letters, Arts and Sciences

Degree Doctor of Philosophy

Degree Program Chemistry

Publication Date 07/02/2007

Defense Date 05/07/2007

Publisher University of Southern California (original), University of Southern California. Libraries (digital)

Tag ab initio,computer simulations,enzyme,free energy,OAI-PMH Harvest,QM/MM,quantum mechanics,solvation

Language English

Advisor Warshel, Arieh (committee chair), Krylov, Anna I. (committee member), Petruska, John A. (committee member)

Creator Email rosta@usc.edu

Permanent Link (DOI) https://doi.org/10.25549/usctheses-m567

Unique identifier UC1224000

Identifier etd-Rosta-20070702 (filename),usctheses-m40 (legacy collection record id),usctheses-c127-510167 (legacy record id),usctheses-m567 (legacy record id)

Legacy Identifier etd-Rosta-20070702.pdf

Dmrecord 510167

Document Type Thesis

Rights Rosta, Edina

Type texts

Source University of Southern California (contributing entity), University of Southern California Dissertations and Theses (collection)

Repository Name Libraries, University of Southern California

Repository Location Los Angeles, California

Repository Email cisadmin@lib.usc.edu