University of Southern California Dissertations and Theses

Two state free energy barrier evaluation technique for dehalogenation reaction

(USC Thesis Other)

Two state free energy barrier evaluation technique for dehalogenation reaction

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content TWO STATE FREE ENERGY BARRIER
EVALUATION TECHNIQUE FOR
DEHALOGENATION REACTION
A Master Thesis
Presented to the Faculty of the Graduate School
of the University of Southern California
in Partial Fulfillment of the Requirements for the Degree of
Master of Science (CHEMISTRY)
by
Ilya Kupchenko
August 2016
c
2016 Ilya Kupchenko
ALL RIGHTS RESERVED
ABSTRACT
The enzymatic reactions play a key role in broad range of life processes. At
the same time theoretical description of enzyme catalysis processes represents
a challenging task due to the complexity of biomolecular systems of interest. In
the present work we validate and enhance a newly developed method: Parady-
namics(PD), targeted at evaluation of the free energy barriers of enzymatic reac-
tions. We present the simulation of the sample molecular reaction - haloalkene
dehalogenation to show, how our approach of using two state framework, in-
volving reference and target potentials, allows us to reduce the complexity of
task and brings an insight into the details of mechanisms of unexplored micro-
scale phenomena.
BIOGRAPHICALSKETCH
Ilya Kupchenko is a student of Master of Science in Chemistry program. His
interest to simulation of molecular systems emerged during his Bachelor of Sci-
ences degree program in Chemistry at Moscow State University(MSU) in 2003-
2008. After the graduation he continued to work as a research assistant for a few
years in MSU. During this time he also tried himself as a high-school teacher of
computer programming. He started his program at the University of Southern
California back in 2011. Today he’s career interests span not only statistical me-
chanics and physical chemistry, but also such fields like data science and prob-
ability theory, algorithms, parallel computing and data structures. Motivated
by possible applications in the field of physical chemistry, he additionally pur-
sued a Master of Science in Computer Science degree in Spring 2016. The results
of his work has been published in The Journal of Physical Chemistry B, in the
his paper “Enhancing Paradynamics for QM/MM Sampling of Enzymatic Re-
actions” of 2016 and in the his presentation at American Chemical Society 2015
Conference in Boston. Ilya is among of the winners of a prestigious Fulbright
Award. Today he is mainly focused on statistical mechanics.
iii
This thesis is dedicated to my parents.
iv
ACKNOWLEDGEMENTS
Acknowledgments: I would like to thank all USC community members who
helped me on my challenging path of becoming an authentic scientific re-
searcher. The base for my academic capacity was set during my graduate classes
program. The first year at USC was hard, but overcoming those hardships
boosted my analytical thinking skills, understanding the micro-scale natural
phenomena to the entirely different level. And the people who helped me were
the USC professors. They were guiding me through the strenuous “boot-camp”
consisting of hardcore masters-level courses. They were suggesting a wise tac-
tics of time-management and helping me with figuring out the concepts of vi-
tal importance for surviving during the course. And, finally they adopted my
needs, and balanced the load of the classes when the situation was dangerous
for me. I remember their kindness, readiness to help and availability of those
people at all times.
Another crucial component o my success at USC was my family and friends.
My family members were thousands of miles away from me but were always
keeping in touch with me. Whenever I needed some house-keeping advice,
difficult situation, or was just feeling lonely, my mother and father were ready to
talk to me. Additionally, I met a few people that were giving me a lot of support
and shared the hardest moments of the studies. They were my classmates, who
came from various countries, and just casual people whom I have met in USC
and who impressed me by their kind and warm attitude. Finally, I feel deeply
thankful to to my fiance who was near me while I was writing my thesis and
broght a lot of light and an entirely new sense in my life.
v
TABLEOFCONTENTS
Biographical Sketch . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iii
Dedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iv
Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . v
Table of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vi
List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vii
List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viii
1 Introduction 1
1.1 QM/MM simulations of condensed phases . . . . . . . . . . . . . 1
1.2 Paradynamic approach . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.2.1 General concepts . . . . . . . . . . . . . . . . . . . . . . . . 2
1.2.2 Previous studies . . . . . . . . . . . . . . . . . . . . . . . . 4
1.2.3 Inverted potential . . . . . . . . . . . . . . . . . . . . . . . . 5
1.2.4 Alternative approaches . . . . . . . . . . . . . . . . . . . . 7
1.3 The benchmark system: haloalkene dehalogenase . . . . . . . . . 9
2 TheoryandBackground 10
2.1 Studying the free energy surfaces . . . . . . . . . . . . . . . . . . . 10
2.1.1 Free energy computing techniques . . . . . . . . . . . . . . 10
3 Methods 14
3.1 Paradynamics: details . . . . . . . . . . . . . . . . . . . . . . . . . 14
4 Results 20
4.1 Modeling: Inverted EVB as Reference Potential . . . . . . . . . . . 20
4.1.1 Selecting the reaction coordiate . . . . . . . . . . . . . . . . 21
4.1.2 Validating the cycle using PMF . . . . . . . . . . . . . . . . 21
4.2 Comparing LRA to FEP . . . . . . . . . . . . . . . . . . . . . . . . 24
4.2.1 Using the LRA for the barrier calculation . . . . . . . . . . 25
4.2.2 Moving from EVB to DFT Potentials . . . . . . . . . . . . . 26
5 Conclusion 27
A Inputdata 29
A.1 Dicloroethane and acetic acid geometric parameters . . . . . . . . 29
A.2 EVB parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
B Simulationresults 35
Bibliography 40
vi
LISTOFTABLES
4.1 Calculating the reaction barrier using the two-state PD scheme
in water and in the active site of the protein. All values are given
in kcal/mole, and the LRA results are shown in the parenthesis.
Taken from [24]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
4.2 Components of the PD cycle computed using LRA approach for
the DFT B3LYP and 6-31g* basis set. Taken from [24]. . . . . . . . 26
A.1 Bonding EVB parameters: U
bond
= D

1expfa(bb
0
)g

2
.
Taken from [24]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
A.2 Angular EVB parameters: U =
1
2
K
0
(
0
)
2
. Taken from [24]. . . . 32
A.3 Non-bonded EVB parameters: U
nb
= A
i
A
j
r
12
B
i
B
j
r
6
. Taken
from [24]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
A.4 Atomic types and charges for EVB potential. Taken from [24]. . . 34
vii
LISTOFFIGURES
1.1 The thermodynamic cycle by knowing the three components of
which we can evaluate activation free energies. To compute the
barrier we first refined the EVB parameters, then we measured
the perturbations from Reference Potential to Target Potential in
Reactant and Product states and, finally we perturbed from RS
to TS using FEP on cheap classical RP . Hence, we obtained the
same amount of free energy as if we would ingrate over upper
curve of the cycle (from RS to TS on TP). Taken from [24] . . . . . 3
1.2 The comparison of the two methods: early version of Parady-
namics (A) and Metadynamics (MTD), is shown. The Meta-
dynamics is consequently filling the minimum of the initial
state(RS) to flatten the free energy surface, which requires it to re-
sample the region near current state after each Gaussian function
was added. The Paradynamics is using renormalization based
on geometry optimization and sets of single point calculations
in implicit solvating models. This allows it to cut-down a huge
amount of sampling and the initial RP , which is already close to
TP is quickly optimized wit several iterations. Taken from [37]. . 4
1.3 A thermodynamic cycle that provides a pictorial way of calcu-
lating the QM/MM activation free energy, where the QM/MM
barrier obtained by moving from the EVB/MM to the QM/MM
potential is shown by a red arrow. EVB2 represents a reference
potential whose minimum is placed at the position of the TS of
the TP . Taken from [24] . . . . . . . . . . . . . . . . . . . . . . . . 6
1.4 (A) Schematic representation of the two-step dehalogenation re-
action mechanism with the charge being transferred from car-
boxyl group to one of the chlorine atoms. (B) Cartesian geome-
tries of RS, TS and PS complexes of the same tye that we used for
our simulations. Taken from [12]. . . . . . . . . . . . . . . . . . . 9
4.1 Dehalogenation reaction RS and TS geometries with marked dis-
tances in intital and SN-2 complex states. Taken from [24]. . . . . 20
4.2 Thermodynamic cycles for PM6 QM and inverted EVB poten-
tials for the water system (a) and for the protein system (b). Re-
covery of the free energy barrier by means of PD(c). Taken from
[24]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
B.1 PMF free energy reaction profile derived using PM6 semi-
empirical method in the enzyme. Taken from [24]. . . . . . . . . . 35
B.2 LRA and FEP of RP to TP transition plotted relatively to
(weight of the state) for EVB RS in protein. Taken from [24]. . . . 36
B.3 LRA and FEP of RP to TP transition plotted relatively to
(weight of the state) for EVB TS in protein. Taken from [24]. . . . 36
viii
B.4 LRA and FEP of RP to TP transition plotted relatively to
(weight of the state) for EVB RS in water. Taken from [24]. . . . . 37
B.5 LRA and FEP of RP to TP transition plotted relatively to
(weight of the state) for EVB TS in water. Taken from [24]. . . . . 37
B.6 Schematic representation of thermodynamic cycle for degalo-
genation molecular process, EVB1 and EVB2 correspond to RS
and TS, and used as a starting points for FEP perturbation to the
target quantum mechanical potential. Taken from [24]. . . . . . . 38
B.7 Transition from reactiant state to SN-2 complex, the negative
charge of the oxygen is equally split on chlorine and oxygen.
1
and
2
are the eigen functions of a mixed two-state Hamiltonian.
Taken from [24]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
ix
CHAPTER 1
INTRODUCTION
1.1 QM/MMsimulationsofcondensedphases
Mixing the Quantum Mechanical and Molecular Mechanical Hamiltonian’s
(QM/MM)[2] is a popular in theoretical chemistry technique for studying var-
ious processes in biological molecules[8]. Since the complexity of the quantum
calculations grows quickly with the size of the system, in order to reduce the
computational time, the environment of the molecule of interest can be rep-
resented by classical Newtonian mechanics. This results in the drastic reduc-
tion of computational time. Recent development of electronic structure calcula-
tion techniques alongside with the high-performance and distributed comput-
ing methods brought the QM/MM approach on the new level of performance
and reliability as a new scientific tool for researching a wide variety of systems:
biomolecules, solutions and surface processes[22, 20, 16, 26].
A popular belief, sometimes called a “Moore’s law”, tells us that the perfor-
mance of computer processors grow exponentially with time. And, therefore,
for the last decade the growth of computer power allowed researchers to sam-
ple the phase spaces of the molecular systems with increasing efficiency. In
other words, the brute force solution of boosting the computational hardware
can be a relieve from a lot of things worth of concern, however, the quantitative
investigations of free energy landscapes of molecular processes by means of ab
initio methods and mixed QM/MM methods[2] remain a huge challenge due
to enormous complexity of electronic part of molecular Hamiltonian. In this
context development of approaches that would allow to reduce the amount of
1
calculations is a matter of crucial importance.
In the this thesis the author might use the direct quotations of his work: [24].
1.2 Paradynamicapproach
1.2.1 Generalconcepts
Following the economy of computational time by reducing the sampled phase
space areas a new approach called Paradynamics was developed recently by
Warshel et. al.[37] The key idea of the method is in using a computationally
cheap classical mechanical potential as a reference potential and a high-level
quantum mechanical method as a target potential [29, 36, 18]. The authors of
the method used the Empirical-Valence Bond(EVB) [4] method as a reference
potential. Such scheme, with using two potentials and thermodynamic pertur-
bations in the key points of the trajectory, that are set to the geometry param-
eters(corresponding to reactant, product, intermediate and transition states),
which is shown on the Figure 1.1, is called Paradynemics [37, 12, 7, 39].
The EVB method was developed in earlier studies of Warshel group [5] and
was extensively used to sample the reaction paths of various molecular pro-
cesses in water and protein phases. The refinement of EVB RP and TP to mini-
mize the difference between simplified classical and high-level quantum poten-
tial surfaces was used before. For the perturbation from reference potential to
target potential we used Free Energy Perturbation (FEP) techniques as well as
Linear Response Approximation (LRA).[33, 13] In order to be able to converge
2
Figure 1.1: The thermodynamic cycle by knowing the three components of
which we can evaluate activation free energies. To compute the
barrier we first refined the EVB parameters, then we measured
the perturbations from Reference Potential to Target Potential
in Reactant and Product states and, finally we perturbed from
RS to TS using FEP on cheap classical RP . Hence, we obtained
the same amount of free energy as if we would ingrate over
upper curve of the cycle (from RS to TS on TP). Taken from [24]
LRA we would need to accurately refine the RP EVB parameters to achieve cor-
respondence in RP and TP . FEP , on the other side is relatively more tolerant to
deviations between two surfaces due to the multiple intermediate steps. How-
ever, FEP requires a good correspondence between the initial and final frames
as well, since this can significantly reduce the error of exponential averaging.
The Paradynamic method allows us to close the thermodynamic cycle with-
out conducting the full integration from RS to TS and save a huge amount of
computational time. To compute the activation barrier, the overall amount of
calls to the high-level expensive target potential can be reduced in almost 200
times, which gives around two orders of performance increase, especially com-
3
Figure 1.2: The comparison of the two methods: early version of Parady-
namics (A) and Metadynamics (MTD), is shown. The Meta-
dynamics is consequently filling the minimum of the initial
state(RS) to flatten the free energy surface, which requires it
to re-sample the region near current state after each Gaussian
function was added. The Paradynamics is using renormaliza-
tion based on geometry optimization and sets of single point
calculations in implicit solvating models. This allows it to cut-
down a huge amount of sampling and the initial RP , which is
already close to TP is quickly optimized wit several iterations.
Taken from [37].
pared to the alternative methods like Methadynamics (MTD), see Figure 1.2,
where the TP potential surface is sampled directly [37, 23].
1.2.2 Previousstudies
As a target potential the authors of the method used to involve semi-empirical
PM6 and PM3 methods as well as density functional theory (DFT) method. The
target potential energy profile along the reaction coordinate (RC) has a mini-
mum at the reactant state (RS) and the product state (PS). The two states are
separated by the reaction barrier. The potential energy profile of the reference
4
potential has two minimum. The first one is located at the target potential reac-
tant state profile minimum along the reaction coordinate. Another minimum of
RP free energy is located in the transition state of the target potential. Hence, the
reference potential is inverted relatively to the target potential at the transition
state.
During our previous studies the system was tested on dichloromethane in
vacuum and in water and dycloroethane with acetic acid in DhlA protein and
in water.[12, 39] The results of the study of dichloromethane revealed that the
the LRA gives a very good approximation of FEP when the potential surface
topologies of RP and TP are close. Additionally, the best results are achieved
by the averaging of points obtained on both TP and RP potentials. [37] Apart
from that, an increase in the convergence was gained by applying Gaussian po-
tential functions to correct the free energy profile in the direction of the reaction
coordinate. The impact of gamma function is related to the reduction of the ab-
solute values of the energies, flattering of the potential and allowed to sample
the degrees of freedom orthogonal to the reaction coordinate more thoroughly.
1.2.3 Invertedpotential
Although, invoking PD framework enables us to have a huge savings in com-
putational cost, to sample the free energy surface even for a simple reaction we
need at least several thousands points per frame, which results in hundreds of
thousands of costly TP points for complete sampling of thermodynamic cycle
perturbations. Regarding the high-level ab initio methods, like DFT or Hartree-
Fock, it would cost us thousands or even millions of hours of single-processor
5
Figure 1.3: A thermodynamic cycle that provides a pictorial way of
calculating the QM/MM activation free energy, where the
QM/MM barrier obtained by moving from the EVB/MM to
the QM/MM potential is shown by a red arrow. EVB2 rep-
resents a reference potential whose minimum is placed at the
position of the TS of the TP . Taken from [24]
time. One of the cut-downs that might be found straightly is normalization and
parameter calibration of the EVB potential in solution. After deriving the op-
timal parameters we can move the simulation in the protein and apply extra
minor calibration for complete convergence. We used this same tactics used in
one of our recent works. [32]
However, since the amino acids of the active site environment of the protein
could polarize the electronic wavefunction in a specific way, one might need
to perform the QM/MM simulation right in the protein. To be able to do this
arduous task we introduced a special modification of the PD framework. Earlier
in our work we added the Gaussian functions to the EVB potential to match the
6
TP[37], which was a very complex task, since the amount of dof’s for a model
of molecule in Cartesian coordinates is 3N-6. The iterative optimization of PES
topology for even 20 atom molecule might take a lot of time itself and may
become a tricky process that can converge to many possible solutions.
Our new technique is based on placing the minimum of the second EVB state
“under “ the TS region . For those purposes we create an artificial state, that has
a minimum on the profile of the free energy along the reaction coordinate, and
which EVB force field parameters are optimized in a way to generate during
the MD trajectory the same geometry as generated by running MD on TP (see
Figure 1.3). In this manner, while being on the RP surface, such type of poten-
tial “holds” the molecular geometry in the TS region. This allows us to save
time on RP optimization by just setting the EVB parameters minimum values
(angles, bonds and dihedral angles) equal to optimized solute in some rough
approximation of solute. After this we perform FEP perturbation from RS and
TS regions of RP to the corresponding regions of the TP [32].
1.2.4 Alternativeapproaches
A few recent studies can be considered as relevant in terms of using thermody-
namic perturbation from classical RP to quantum TP in order to bring down the
simulation time to a reasonable amount. In of the most well-known papers by
Houk et al. the authors use the artificial minimum at the TS state[28]. Despite of
the seeming similarity, our approach has numerous differences. Not going into
details, we can highlight the key differences: one is that the TS minimum of RP
in our case exactly corresponds to the peak on the TP (in terms of position of TS
7
on free energy surface projection along the reaction coordinate), and another is
that unlike the [28] the RS minimum is implemented only in order to improve
the convergence of FEP thermodynamic integration.
Hiemdal and coworkers have used another approach that by its concep-
tual idea can be closer to what we developed [28]. In their approach our idea
is evolving into a separate framework, which is called “Q2MM” (quantum to
molecular mechanics). They use an RP based on the same reference classical
potential that was used by Houk and their colleagues. The fixed solute approx-
imation that they use is different from Molaris PD/LD scheme that we used for
our simulations. This does not allow them to adequately capture the entropic
contribution in the kinetic barrier of the reactions. Despite of some attempts to
run trajectories with moving the solute the FEP integration was done only us-
ing the RP sampling, which contributed a lot to the error of the FEP calculation,
as our further research proved and in some cases of the work even prevented
the convergence completely. One of the biggest faults of the team was that the
thermodynamic cycle values were not verified by direct PMF calculation on the
pure target potential.
In this thesis we focus on the benchmarking the new PD approach on the first
step of dehydrogenation reaction in the DhlA enzyme and, correspondingly,
the non-enzymatic esterification of dycloroethane, which is shown on Figure
1.4. This reaction has been widely studied in the previous works related to PD
approach. [38, 37]
8
Figure 1.4: (A) Schematic representation of the two-step dehalogenation
reaction mechanism with the charge being transferred from
carboxyl group to one of the chlorine atoms. (B) Cartesian ge-
ometries of RS, TS and PS complexes of the same tye that we
used for our simulations. Taken from [12].
1.3 Thebenchmarksystem: haloalkenedehalogenase
To validate the Paradynamics approach we conducted the benchmark simula-
tion using one of the widely studied reactions. This was the first SN-2 step of
dehalogenation reaction catalyzed by haloalkene dehalogenase (DhlA). This re-
action was extensively experimentally studied and simulated by Rosta et al.[3].
The 2DHD [40] structure was taken from Protein Data Bank [9] as the starting
point for geometry parameters relaxation.
The experimental results had shown that the acceptable free energy barrier
values can be obtained via perturbation from EVB potential with finely tuned
force field parameters to DFT potential energy surface.
9
CHAPTER 2
THEORYANDBACKGROUND
2.1 Studyingthefreeenergysurfaces
2.1.1 Freeenergycomputingtechniques
LinearResponseTheory
The Linear Response Approximation (LRA) is based on assumption that the
ambiance of a molecule or a system of interacting molecules(like protein-ligand
system) acts as a coupled set of harmonic oscillator states [10]. Those har-
monic oscillator states correspond to equilibrium states of the solvent of the
protein. The free energy profile along the chosen coordinate for such system
has parabolic shape. The susceptibility of the system to the perturbing potential
can be extracted from the dynamics of the unperturbed state. The distribution of
the potential energies of the equilibrium and perturbed state should be of Gaus-
sian shape. Such approach is exact in case of classical mechanical ensembles
and implies that the response f the system to perturbation can be represented as
linear function of the perturbation potential (for instance electric field) and ex-
pansion of the free energy. One of the important peculiarities of this technique
is the ability to determine a systematic error and deduct it, extrapolating the
knowledge to different systems and states.
10
Freeenergyperturbation(FEP)
Free energy perturbation method is another way to compute the difference in
free energy between two states. It is based on interpreting the configuration
integrals ratio, which corresponds to the free energy function difference, as a
thermodynamic integral of states energy difference. The first implementation
of the method was presented in pioneering work dedicated to the simulation of
high-temperature argon and nitrogen mixture by Zwanzig [42, 10]. For given
two states the free energy is evaluated according to the following expression,
which is called “Zwanzig formula”:
F(A! B) = F
B
F
A
=k
B
T ln
*
exp

E
B
E
A
k
B
T
!+
A
(2.1)
T is the temperature andk
B
is Boltzmann’s constant. In practice, to compute
the difference between two states we run the simulation using the potential of
one state to sample the phase space. At the same time we compute the energy
of the other state. The states can represent different conformations of molecules,
charge transfer, reactions of various types or even phase transitions and alchem-
ical transformation of one atom into another. Having two sets of energies for the
same sampling points we make a histogram of energy values of state A. We then
calculate exponential differences weighted by then
E
A
exp(
E
A
k
B
T
), wheren
E
A
is the
number of occurrences of phase space points with energy equalE
A
. This is what
is called the “thermodynamic averaging”: it allows us not only account for the
energetic effects, but also include entropic contributions through the exponen-
tial weights. Practically FEP is frequently done using the splitting of the whole
perturbation into several windows, called the “intermediate states”, since the
quality of the calculation depends significantly on the overlap between the po-
11
tential energy histograms of the initial and final states.
FEP approach is widely used to study the kinetics of a vast variety of pro-
cesses, including proton transfer, ligand-receptor interactions, adsorption pro-
cesses and etc. The method has a great role in computer aided drug design and
protein engineering [10].
WHAM
The Weighted Histogram Analysis Method (WHAM) is an advanced technique
of free energy calculations suggested by Robert H. Swendsen [19, 10] as an ex-
tension to the system of Ferrenberg-Swendsen equations that he has previously
derived. This approach which allows to collect and use information from sev-
eral windows in a most efficient way. WHAM has a number advantages over
standard FEP of Potential of Mean Force methods. It can optimize the links be-
tween windows, which is a major problem for standard PMF approach, usually
resolved by linear extrapolation.
WHAM works as a natural generalization of the single histogram method.
Assume, that we have conducted R molecular dinamical simulations at tem-
peratures T
i
where i 2 f1;:::;Rg and with the set of constraints
˜
V
j
, where the
constraints are weighted by the set of coefficients
j;m
, where j2f1;:::;Lg and
m 2 f1;:::;Rg. Then the WHAM formula, for iterative reweighting of the his-
tograms of thoseR calculations will be:
12
P
fg;
(fVg;) =
P
R
k=1
N
k
(fVg;)exp

P
L
j=0

j
V
j

P
R
m=1
n
m
exp

f
m

m
P
L
j=0

j;m
V
j

exp(f
j
) =
X
fVg;
P
fg
j
(fVg;)
(2.2)
where N
k
(fVg;) is a number of sampling points taken by histogram bin al-
located for constraining potential configurationfVg and the value of the coordi-
nate of the histogram during the run number i, that same is for n, but using a
different counter.[19]
This allows us to calculate the mutual reweighing of the histograms. To
start calculation we can set the inital f values to zero. The method is relatively
fast, but, like the previous techniques, it requires significant overlapping of the
intermediate states’ histograms.
13
CHAPTER 3
METHODS
3.1 Paradynamics: details
Compared to calculating a common thermodynamical cycle, in Paradynamics
(PD) [29, 36, 18, 13]we are trying to use the cheap molecular mechanical refer-
ence potential as well as the perturbations from it to the high-level QM/MM
quality target potential to evaluate the reaction activation free energy. To do
this we are plotting the thermodynamic cycle based on the four states: Reac-
tant State (RS) on Reference Potential (RP), Reactant State on Target Potential
(TP), Transition State (TS) on Reference Potential, Transition State (TS) on Ref-
erence Potential (RP). On the lower transition (from RS to TS on RP) we are
using our EVB[4] classical force field that is well-parametrized to reproduce the
key qualities of the system and we calculate this transition using multi-window
(11-41 windows) FEP[14, 27] . On the vertical arrows of the cycle we are us-
ing multi-window FEP as well, but here we use alchemical transformation from
EVB quantum region to QM quantum region. The EVB is designed to repro-
duce the bond lengths and angles corresponding to the QM average values and
uses Mulliken or electro-static potential fitted (ESP) charges for solvent-solute
interactions [37]:
g
E
RP
!E
TP

,
E
TP

1
2

D
E
TP

,
E
TP

E
RP

,
E
TP
E
E
TP
+
D
E
TP

,
E
TP

E
RP

,
E
TP
E
E
RP
(3.1)
In the above formula
,
is the reaction coordinate (RC) of TS and E
RP
and
14
E
TP
is the potential energy of RP and TP correspondingly.
g
E
RP
E
TP

E
TP

1
2

E
TP

E
TP

E
RP

E
TP

E
TP
+

E
TP

E
TP

E
RP

E
TP

E
RP

(3.2)
Analogously to the previous expression this one describes the same value
but for reactant state. The corresponding potential energies for the case of har-
monic RC constraint are given by:
E
RP
= E
RP
+E
cons
E
TP
= E
TP
+E
cons
E
cons
= K(
0
)
2
(3.3)
The value for the is chosen in correspondence to the positions of the re-
actant and transition states, so that it matches both EVB and QM free energy
surface (FES) minimum. The above equations for the free energies of the reac-
tant and transition states imply the linear reply to the constraint perturbation
and good overlap between the RS and TS. However, in most cases [7, 39], such
assumption is far from being correct: the EVB can not adequately represent a lot
of degrees of freedom of complex QM Hamiltonian. To solve this problem, and
find the activation barrier:
g
E
RP
E
TP
() = g
TP
() g
RP
() (3.4)
we are using multiple intermediate states, which are the mixture of QM and
MM states with certain proportion. By gradually increasing the weight of one
15
potential and decreasing the weight of another, we achieve a good overlapping
of each of the pairs of adjacent simulation windows:
E
m
= (1
m
)E
RP
+
m
E
TP
+E
cons
(3.5)
The m counter moves from zero to the number of windows N, and moves
from
m=0
= 0:0 to
m=N
= 1:0. After this we sum up all the windows and deriving
the total free energy of the perturbation from RP to TP . Therefore, the first part
of the above equation will be:
g
RP
() = G
m
kT ln

(x)exp

E
RP
E
m
kT

E
m
(3.6)
and the second part:
g
TP
() = G
m
kT ln

(x)exp

E
TP
E
m
kT

E
m
(3.7)
The G
m
is absolute free energy of certain window, and if the constraint
contribution is very small, we can neglect this part of potential to simplify our
calculations:
E
m
= (1
m
)E
RP
+
m
E
TP
(3.8)
The procedure for the calculating the barrier of the RS to TS transition on the
RP potential is a standard FEP that was extensively performed in the previous
scientific research. With the above equations we can compute the three values
16
out of four in our cycle. Assuming that the closed cycle for isolated thermo-
dynamic system results in zero free energy change, we can easily calculate the
fourth value by simple arithmetics, which will be the transition from RS to TS
on the high-level QM potential:
g
,
QM=MM
() =g
EVB1!QM(TS)
+ g
EVB2!QM(TS)
+ g
EVB1!EVB2
(3.9)
Hence, we computed the reaction barrier using in high-level QM/MM ap-
proach without strenuous direct sampling of the entire QM surface. This shows,
how the chemical knowledge contained in empirical EVB potential can be used
for saving a lot of computational resources.
In regular EVB we usually mix two states: EVB1 and EVB2 corresponding to
RS and product state(PS). This give an adequate TS for the most cases, which is
located somewhere near = 0:5, i.e. at at the state which is 50% of RS and 50%
PS. However, if we want to perturb from RP to TP at transition state we need
somewhat reliable and precise approximation of the transition state geometry.
For those purposes the standard EVB mapping scheme was modified and we
created an artificial equilibrium state corresponding to geometry parameters of
optimized at QM saddle point (see the Figure 1.3) We have also scaled the abso-
lute value of TS EVB gas phase free energy to reduce thermodynamic averaging
error.
Following the previous work related to PD method, to benchmark our in-
novative computational approach we considered a well studied the first SN-2
step of two-step process of dehalogenation reaction involving acetic acid and
dycloroethane in water and the same reaction in protein. For protein simulation
17
we decided to choose 2DHC geometry structure from protein data bank (PDB).
[15]
Like in our previous work we were using MOLARIS [12, 34] program pack-
age for our research. This package contains both MD integration algorithm
and multi-window FEP algorithm. The protein amino acid chains were mod-
eled by means of ENZYMIX force field which is implemented in the MOLARIS
package. To validate the results of our thermodynamic cycle calculations we
were using the multi-window PMF scheme based on Umbrella Sampling(US)
method. In other words we were conducting multiple runs, gradually changing
the minimum position of the constraint potential, and, therefore, equally sam-
pling the FES profile from RS to TS. [5, 6]. Surface constrained all atom solvent
(SCAAS)[34] model proved itself to be a robust tool for solute-solvent interac-
tions modeling, so we used it this time as well. The molecule was surrounded
by 18 sphere of explicit all-atom water molecule models, and 2 outer sphere of
point dipoles. On further distances the solvent was modeled by dielectric con-
stant. Local reaction field (LRF)[34] method was used to simulate long range
effects, like amino acid to amino acid or ligand to amino acid interactions.
For sampling the phase space we were using 120 ps trajectories. We were
mainly using 11 windows for both RS to TS on RP FEP mapping and we used
21-windows frames only for validation purposes. We were relaxing system for
100 ps and collected the energies and geometric coordinates only for the last
20 ps. We used Berendsen thermostat and the temperature was set to 300 K.
For the RS(EVB1) to TS(EVB2) on RP FEP mapping we have also tried the 200
ps windows, with 100 ps of relaxation and 100 ps of effective trajectory. Those
frameworks were used for both water phase and for protein. We have also var-
18
ied the time step an discovered that the optimal time step for EVB1 and EVB2 to
PM3 as well as to PM6 FEP perturbations the best integration step is 0.5 fs. For
RS to TS on RP we used 1.0 fs step. For PMF, on RS to TS FES we used again 0.5
fs. The lower step-size is because of the sensitivity of semi-empirical methods to
the input geometry: given higher MD integration step size the electronic wave
function tend to have convergence problems.
We tried to ionize the residues around the active site. In DhlA the ioniz-
able residues in the nearest proximity to the active site are: ASP260, ARG229,
LYS176, ASP178, GLU293, ASP266, LYS221, LYS224, LYS261. The effect of ion-
ization was not more than 0.5-1.0 kcal/mole. The residues were ionized by MO-
LARIS program in an automatic way, the corresponding proton-donating oxy-
gen, nitrogen and sulfur atoms were assigned to have a reasonable point charge.
Further, it was decided that we will proceed our simulations in the neutral state,
since it did not influence on the FES a lot, and our goal was to compare the ther-
modynamic cycle with the PMF rather than to reach the perfect precision of free
energy barriers.
We used 250 kcal/mole harmonic constraints and additional small bond
length constraints on SN-2 [Cl–C–O] group of 2.0 kcal/mole to keep the SN-
2 complex together during the PMF calculations. To relax the system we also
imposed angle constraints on SN-2 complex angle. The main RC constraint was
imposed gradually, starting 50 kcal/mole and going up to 250 kcal/mole. The
initial geometries for water phase simulations were generated using MOPAC
constrained optimization, were we varied constraints to create several inputs
for different RC values. For QM part of our QM/MM we used B3LYP funcional
[35, 11] with 6-31g* basis set. The electronic Hamiltonian was calculated us-
19
ing Qchem 4.0 [41] for semi empirical energy and force calculations we used
MOPAC2009 and its PM3 and PM6 methods.
20
CHAPTER 4
RESULTS
4.1 Modeling: InvertedEVBasReferencePotential
Our work is dedicated to enhancement and development of PD approach. From
Figure 4.1 you can see that the system of interest is relatively small (13 heavy
atoms) but, at the same time, big enough to provide solute entropic contribu-
tions, unlike smaller systems and processes involving less atoms, like water de-
protonation or SN-2 methyl chloride and iodine reaction, that were used for the
previous works dedicated to the development of Paradynamics. [7, 39, 21] The
overall performance and accuracy of PD depends on convergence of the FEP
calculations, which, in turn depend on the likelihood of RP and TP in terms of
FES topology.
Figure 4.1: Dehalogenation reaction RS and TS geometries with marked
distances in intital and SN-2 complex states. Taken from [24].
We used PM3 bondlength and angles to parametrize the EVB states for all
the cycles: PM3, PM6 and B3LYP . EVB has a very well refined parameters set so
we did not care much about the parameters for atoms and bonds that are not
involved in SN-2 complex. [17] As you can see on Figure 4.2 we pulled the TS
21
EVB state down by 45 kcal/mole to make the RP to TP perturbations approxi-
mately equal for both RS and TS. The off diagonal element for RP FEP mapping
was set to 15 kcal/mole. After tuning all the parameters and establishing a reli-
able workframe we tried to see if our modified PD scheme provides us the same
results as direct PMF on TP .
4.1.1 Selectingthereactioncoordiate
Our fist simulation was simulation of DhlA and water-phase dehalogenation
using PM6 as a QM part[30] and EVB as RP and PM6 as TP . The two EVB states
corresponding to the intial and transition state are shown on Fig. 4.1. Following
our previous research the reaction coordinate was set to be equal to the differ-
ence: = d(ClC)d(CO). Using the data from previous work [39] and several test
runs, as well as quasi-static multi-step optimization in MOPAC we determined
that the RC corresponding to RS is equal to 1.3 A, and for TS it is 0.2 A. The
angle constraint we used was equal to 160.0
o
.
4.1.2 ValidatingthecycleusingPMF
After running all the trajectories we calculated the free energies for the cycle
parts. For the evidence of convergence of PMF and FEP methods we did not
need to reproduce the exact results of the experimental studies. However, it
was preferable for us that the barrier would be high enough to keep the pro-
cess in the limits of transition state theory, i.e. higher than 10 kcal/mole. The
overall algorithm f the PMF integration was described in the previous papers
22
[39]. However, during our calculus we found a few significant problems in our
old PMF algorithm. The new algorithm with linear regression used to connect
piecewise FES PMF profile was implemented. We derived the PMF profiles for
both water and protein environment. For the PM6/MM approach they were
26.0 and 14.0 kcal/mole in water and in protein, correspondingly. The resulted
curves can be seen on Figure 4.2.
Figure 4.2: Thermodynamic cycles for PM6 QM and inverted EVB poten-
tials for the water system (a) and for the protein system (b).
Recovery of the free energy barrier by means of PD(c). Taken
from [24].
23
Table 4.1: Calculating the reaction barrier using the two-state PD scheme
in water and in the active site of the protein. All values are given
in kcal/mole, and the LRA results are shown in the parenthesis.
Taken from [24].
System g
EVB1!PM6
g
EVB2!PM6
g
EVB1!EVB2

g
,
RS!TS

PD
c

g
,
RS!TS

PMF
c
Water 77.0 (78.0) 59.0 (59.0) 8.0 26.0(27.0) 19.0
b
DhlA 89.0 (89.0) 75.0 (74.0) 0.0 14.0 (15.0) 10.0
b
9.4
d
a
RS! EVB1 and TS! EVB2
b
Experimental data is from [39]
c
The statistical errors of our method are discussed elsewhere.
d
Taken from the current project.
The obtained values correspond well with the presented in the experimental
work[25], they are 27.0 kcal/mole for the water and 15.3 kcal/mole for the en-
zymatic phase. However, the previous results for the analogous system derived
in one of our papers were a bit lower: 19.0 kcal/mole for water and 10.0 kcal/-
mole for protein. See the summary of the resulting free energies in the Table 4.1.
The reasons for such difference were unclear at this stage of modeling.
To understand the nature of the deviation from our previous results we tried
to vary the minimum positions of the RS and TS constraints for “vertical” RP to
TP perturbations. We varied them within [-0.2;+0.2] gap around the RC equilib-
rium value for each state (EVB1(RS) or artificial EVB2(TS)). This did not result in
more than 1.0 kcal/mole error. Our next attempt was to try the longer trajecto-
ries. We tried to take 50 ps out of 100 ps of trajectory time as meaningful points.
Again, as in the previous case this did not give us more than 1.0 kcal/mole
change in the energy. Our last step was to visualize the Gaussian distributions
24
of the potential energy. The energy distributions had a good overlap for all the
windows of FEP: at least 25-50% Gaussian curve square were overlapping. For
the future calculus we might even reduce the number of windows, since some
of adjacent FEP frames overlapped on almost 90%.
4.2 ComparingLRAtoFEP
As a part of our research project we evaluated the PMF and FEP-based PD meth-
ods in terms of accuracy and tried to set a reasonable error bars for our values.
This step was important for us in order to analyze the possibility of minimizing
the number of windows for our FEP RP to TP perturbations and, as an ultimate
goal, deducing it to two windows - one on RP and one on TP . This would cor-
respond to one-step LRA, assuming that we are using the simplified formula
for LRA[31]: E
TP
(
TP
) E
RP
(
RP
). Our first goal was to calculate the error of
the multi-window FEP and PMF error. To do this we used two different ap-
proaches, and conducted our validation for water phase only. In the first one
we run 100 ps trajectories on EVB potential constraining system at both RS and
TS. After that we took a snapshots of every 10 ps of run allowing the solvent
structure to renew. From each of those runs we calculated the free energy val-
ues and estimated the standard deviation(SD). For the FEP and the PMF the SD
was equal to 1.8 kcal/mole and 1.6 kcal/mole at TS. However, to make sure
that everything is stable and we can safely try the LRA we have also computed
several trajectories which were run from different optimized initial structures
both in water and in protein. This set of runs might be presented as some ana-
log of a piece of replica-exchange run[1], but being done manually and using a
small amount of exchanges. The small amount of trajectories is done because
25
of computational limitations. For the case of different starting geometries the
resulting errors were 2.1 kcal/mole and 1.3 kcal/mole. However, in practice
the divergence between the PMF and the PD cycle was around 6 kcal/mole for
PM6 QM part. It is still unclear, what is the source of this error. Most likely it is
a combination of different factors, one of them is a limited sampling due to the
computationally expensive QM potential.
4.2.1 UsingtheLRAforthebarriercalculation
Being confident in a reasonable error boundaries for FEP and PMF results we
decided to compare the FEP with LRA. As a two-state approach LRA gives enor-
mous boost in the speed of calculation. However, one should take care of build-
ing the appropriate RP , where the probability distributions of RS and especially
TS phase spaces correspond to each other and the energy Gaussian overlap well.
The same PD cycle with LRA instead of the FEP gave us the activation bar-
riers of 15.0 kcal/mole for the water and 27.0 kcal/mole for the protein. The
difference in the barriers between FEP and LRA turned out to be around 1.0
kcal/mole. This gave us a solid ground for using the LRA technique for further
simulations, where QM part wavefunction was be computed by means of DFT
when computing the TP energy. The results are shown in Table [?]
26
Table 4.2: Components of the PD cycle computed using LRA approach for
the DFT B3LYP and 6-31g* basis set. Taken from [24].
System g
EVB1!DFT
g
EVB2!DFT
g
EVB1!EVB2

g
,
RS!TS

PD

g
,
RS!TS

PMF
Water 42.0 28.0 8.0 22.0 -
DhlA 44.0 33.0 0.0 11.0 8.8
4.2.2 MovingfromEVBtoDFTPotentials
After our experience with the PM6 semi-empirical potential we understood that
even relatively basic QM technique takes a lot of processing time, when it comes
to thousands of steps. More advanced methods, like B3LYP take even more time
and using our 11-window FEP RP to TP perturbations becomes totally imprac-
tical. At the last step of our theoretical research we can make use of LRA advan-
tages. Reducing the number of windows to two, instead of 11 drastically drops
down the cost of PD cycle computing.
Following the LRA method we calculated the PD cycle for the case where RP
is an EVB classical potential and TP is DFT with B3LYP functional and 6-31g*
basis set. The overall results are presented in Tab. 4.2. The values of the acti-
vation free energy for the water environment was 22.0, while for the enzimatic
environment it constituted 10.7 kcal/mol. The overall catalytic effect of trans-
fering from water to the protein turned out to be 11.0 kcal/mole which is close
to the experimentally determined value. [39]
27
CHAPTER 5
CONCLUSION
The outcome of this work is a renovated Paradynamics protocol of molec-
ular dynamic simulation that can serve for studying the reactions in the water
and protein environments. The major contribution to the new version of PD
our new method is that instead of representing the transition state as a mixture
of reactant state potential (EVB1) and product state potential (EVB2), we create
a new system of two EVB states, where reaction state minimum remains on its
place, and at the phase space area of the transition state we create a new artificial
minimum. Both reaction state and transition state EVB potential parameters are
fitted in a way that they will reproduce the probability distribution of degrees of
freedom of the target high-level quantum mechanical or semi-empirical poten-
tial. In this was we can compute the two vertical perturbations: FEP value from
RP to TP begin constrained near the RS phase space region, FEP value from RP
to TP begin constrained near the TS phase space region, and, finally FEP value
of EVB1(RS) to EVB2(TS) on RP . Thus, we can close the cycle and obtain the bar-
rier avoiding calculation of upper PMF RS to TS curve, which needs about 40-60
windows of Umbrella Sampling on expensive QM potential. Additionally, we
proved that having well-fitted EVB parameters for RS and TS one can avoid us-
ing multi-window FEP for vertical transitions from RP to TP , and perform only
two-state LRA calculus. In this way the amount of time to close PD cycle is
several times less than the time needed to perform the full PMF.
Such framework enables us to calculate the free energy barriers using the
high-level ab initio or DFT methods for the ligand simulation in QM/MM mod-
els and finalizing the sampling within a reasonable time. However, to achieve a
28
close correspondence between the probability distributions of the solute molec-
ular coordinates we still need to do a preliminary optimization in the water, be-
fore moving to protein. One of the key competitive advantages of our method,
except its convergence speed, is that it allows to capture the solute probability
distribution over the reaction coordinate, and, therefore, solute entropic effect,
compared to a bunch of previously developed techniques that describes solvent
entropy impact only. Concluding the aforementioned facts, our renewed PD
framework based on the two-state thermodynamic cycle is proved to be robust,
fast and perspective tool for investigation of QM/MM FES topologies [16] of
bimolecular systems of various complexity.
29
APPENDIX A
INPUTDATA
A.1 Dicloroethaneandaceticacidgeometricparameters
Listing A.1: PDB Cartesian geometry of the RS after 100 ps relaxation in PM6
potential. Taken from [24].
ATOM 1 CL1 DCE 1 2.020 -0.398 1.628
ATOM 2 C1 DCE 1 1.412 -0.001 -0.003
ATOM 3 C2 DCE 1 2.230 -0.521 -1.155
ATOM 4 CL2 DCE 1 2.272 0.542 -2.702
ATOM 5 H11 DCE 1 0.321 -0.349 -0.104
ATOM 6 H12 DCE 1 1.371 1.170 -0.124
ATOM 7 H21 DCE 1 1.779 -1.412 -1.559
ATOM 8 H22 DCE 1 3.268 -0.767 -0.797
ATOM 9 C001ACC 2 -1.978 2.223 -1.228
ATOM 10 H001ACC 2 -2.962 2.641 -1.441
ATOM 11 H002ACC 2 -1.156 2.917 -1.028
ATOM 12 H003ACC 2 -2.054 1.741 -0.230
ATOM 13 C002ACC 2 -1.509 1.166 -2.218
ATOM 14 O001ACC 2 -2.262 0.205 -2.625
ATOM 15 O002ACC 2 -0.288 1.188 -2.610
30
Listing A.2: PDB Cartesian geometry of the TS after 100 ps relaxation in PM6
potential. Taken from [24].
ATOM 1 CL1 DCE 1 2.363 2.143 0.471
ATOM 2 C1 DCE 1 1.551 0.116 -0.091
ATOM 3 C2 DCE 1 2.439 -0.361 -1.146
ATOM 4 CL2 DCE 1 1.870 -0.239 -2.855
ATOM 5 H11 DCE 1 0.561 0.550 -0.391
ATOM 6 H12 DCE 1 1.738 -0.232 0.922
ATOM 7 H21 DCE 1 2.647 -1.453 -1.065
ATOM 8 H22 DCE 1 3.385 0.270 -1.330
ATOM 9 C001ACC 2 -0.119 -4.019 0.077
ATOM 10 H001ACC 2 0.648 -4.412 -0.656
ATOM 11 H002ACC 2 -1.096 -3.900 -0.424
ATOM 12 H003ACC 2 -0.171 -4.786 0.890
ATOM 13 C002ACC 2 0.306 -2.721 0.608
ATOM 14 O001ACC 2 0.757 -2.647 1.806
ATOM 15 O002ACC 2 0.284 -1.699 -0.160
31
A.2 EVBparameters
Table A.1: Bonding EVB parameters: U
bond
= D

1expfa(bb
0
)g

2
.
Taken from [24].
Bond type D b
0
a
C0–H0 100.40 0.99 2.0
C0–L0 79.00 1.82 0.8
C1–L- 80.00 2.10 0.8
O1–C1 80.00 2.00 0.8
C0–C0 96.00 1 1.54 0.8
C+–O+ 92.50 1.40 2.0
P0–O- 120.0 1.50 2.0
P0–O+ 120.0 1.50 2.0
P0–O0 95.0 1.58 2.0
C0–N0 93.0 1.44 2.0
CT–NT 93.0 2.00 2.0
CT–ST 93.0 2.20 2.0
C1–O1 80.00 2.00 0.8
32
Table A.2: Angular EVB parameters: U =
1
2
K
0
(
0
)
2
. Taken from [24].
Angle type
0
K

X-C+-Y 120.0 50.00
X-C0-Y 109.5 50.00
X-C1-Y 109.0 50.00
X-O1-Y 109.0 50.00
X-O-Y 120.0 50.00
X-P0-Y 109.0 100.00
X-CT-Y 120.0 50.00
X-S+-Y 103.5 100.00
X-ST-Y 107.0 100.00
X-O+-Y 120.0 50.0
33
Table A.3: Non-bonded EVB parameters: U
nb
= A
i
A
j
r
12
B
i
B
j
r
6
. Taken
from [24].
Atom type A B
H0 7.00 0.00
C0 632.00 24.00
C1 632.00 24.00
O+ 774.00 24.00
O1 774.00 24.00
L0 774.00 24.00
L- 20500.00 24.00
P0 1250.00 24.00
S+ 1022.00 24.00
ST 1022.00 24.00
O0 774.00 24.00
34
Table A.4: Atomic types and charges for EVB potential. Taken from [24].
EVB1 EVB2
Atom number Atom Type Charge Atom Type Charge
1 C0 -0.5001 C0 -0.5424
2 H0 0.1564 H0 0.2085
3 H0 0.1564 H0 0.2085
4 C+ 0.8442 C+ 0.7885
5 O+ -0.7861 O+ -0.6795
6 O- -0.7861 O1 -0.6006
7 L0 -0.2230 L- -0.5838
8 C0 -0.1217 C1 0.2074
9 C0 -0.1217 C0 -0.3096
10 L0 -0.2230 L0 -0.2254
11 H0 0.1511 H0 0.0698
12 H0 0.1511 H0 0.0698
13 H0 0.1511 H0 0.1943
14 H0 0.1511 H0 0.1943
35
APPENDIX B
SIMULATIONRESULTS
Figure B.1: PMF free energy reaction profile derived using PM6 semi-
empirical method in the enzyme. Taken from [24].
36
Figure B.2: LRA and FEP of RP to TP transition plotted relatively to
(weight of the state) for EVB RS in protein. Taken from [24].
Figure B.3: LRA and FEP of RP to TP transition plotted relatively to
(weight of the state) for EVB TS in protein. Taken from [24].
37
Figure B.4: LRA and FEP of RP to TP transition plotted relatively to
(weight of the state) for EVB RS in water. Taken from [24].
Figure B.5: LRA and FEP of RP to TP transition plotted relatively to
(weight of the state) for EVB TS in water. Taken from [24].
38
Figure B.6: Schematic representation of thermodynamic cycle for degalo-
genation molecular process, EVB1 and EVB2 correspond to RS
and TS, and used as a starting points for FEP perturbation to
the target quantum mechanical potential. Taken from [24].
39
Figure B.7: Transition from reactiant state to SN-2 complex, the negative
charge of the oxygen is equally split on chlorine and oxygen.

1
and
2
are the eigen functions of a mixed two-state Hamil-
tonian. Taken from [24].
40
BIBLIOGRAPHY
[1] Replica-exchange molecular dynamics method for protein folding. Chemi-
calPhysicsLetters, 314(12):141 – 151, 1999.
[2] Warshel A. Theoretical studies of enzymic reactions: Dielectric, elec-
trostatic and steric stabilization of the carbonium ion in the reaction of
lysozyme. J.Mol.Biol., 103:227, 1976.
[3] Warshel A. Theoretical studies of enzymic reactions: Dielectric, elec-
trostatic and steric stabilization of the carbonium ion in the reaction of
lysozyme. J.Mol.Biol., 103:227, 1976.
[4] Warshel A. An empirical valence bond approach for comparing reactions
in solutions and in enzymes. J.Am.Chem.Soc., 102:6218, 1980.
[5] Warshel A. Computer Modeling of Chemical Reactions in Enzymes and Solu-
tions. John Wiley & Sons, New York, 1991.
[6] Warshel A. Electrostatic basis for enzyme catalysis. Chem. Rev., 106:3210,
2006.
[7] Warshel A. Molaris-xg, v 9.11. University of Southern California, Los An-
geles, CA, 2012.
[8] Warshel A. Multiscale modeling of biological functions: From enzymes to
molecular machines (nobel lecture). Angew.Chem.,Int.Ed., 53:10020, 2014.
[9] Helen M. Berman, John Westbrook, Zukang Feng, Gary Gilliland, T. N.
Bhat, Helge Weissig, Ilya N. Shindyalov, and Philip E. Bourne. The protein
data bank. NucleicAcidsResearch, 28(1):235–242, 2000.
[10] Christophe Chipot.FreeEnergyCalculations. Springer Series in CHEMICAL
PHYSICS, Berlin, 2007.
[11] Becke A. D. Density-functional thermochemistry 0.3. the role of exact ex-
change. J.Chem.Phys., 98:5648, 1993.
[12] Rosta E. Towards accurate ab initio qm/mm calculations of free-energy
profiles of enzymatic reactions. J.Phys.Chem.B, 110:2934, 2006.
41
[13] Rosta E. Towards accurate ab initio qm/mm calculations of free-energy
profiles of enzymatic reactions. J.Phys.Chem.B, 110:2934, 2006.
[14] King G. Investigation of the free energy functions for electron transfer re-
actions. J.Chem.Phys., 93:8682, 1990.
[15] Verschueren K. H. G. Crystallographic analysis of the catalytic mechanism
of haloalkane dehalogenase. Nature, 363:693, 1993.
[16] Hu H. Free energies of chemical reactions in solution and in enzymes with
ab initio quantum mechanics/molecular mechanics methods. Annu. Rev.
Phys.Chem., 59:573, 2008.
[17] Nakayama H. Structure of a hyperthermophilic archaeal homing endonu-
clease, i-tsp061i: Contribution of cross-domain polar networks to ther-
mostability. J.Mol.Biol., 365:362, 2007.
[18] Bentzien J. Hybrid ab initio quantum mechanics/molecular mechanics cal-
culations of free energy surfaces for enzymatic reactions: The nucleophilic
attack in subtilisin. J.Phys.Chem.B, 102:2293, 1998.
[19] Shankar Kumar, John M. Rosenberg, Djamal Bouzida, Robert H. Swendsen,
and Peter A. Kollman. The weighted histogram analysis method for free-
energy calculations on biomolecules. i. the method. JournalofComputational
Chemistry, 13(8):1011–1021, 1992.
[20] Kamerlin S. C. L. Progress in ab initio qm/mm free-energy simulations
of electrostatic energies in proteins: Accelerated qm/mm studies of pka,
redox reactions and solvation free energies. J.Phys.Chem.B, 113:1253, 2009.
[21] Kamerlin S. C. L. The empirical valence bond model: Theory and applica-
tions. WileyInterdiscip.Rev.-Comput.Mol.Sci., 1:30, 2011.
[22] Kamerlin S. C. L. Multiscale modeling of biological functions. Phys.Chem.
Chem.Phys., 13:10401, 2011.
[23] Mones L. The energy gap as a universal reaction coordinate for the simu-
lation of chemical reactions. J.Phys.Chem.B, 113:7867, 2009.
[24] Jernimo Lameira. Enhancing paradynamics for qm/mm sampling of enzy-
matic reactions. TheJournalofPhysicalChemistryB, 120(9):2155–2164, 2016.
42
[25] Olsson M. H. M. Solute solvent dynamics and energetics in enzyme catal-
ysis: The s(n)2 reaction of dehalogenase as a general benchmark. J. Am.
Chem.Soc., 126:15167, 2004.
[26] Senn H. M. Qm/mm methods for biomolecular systems. Angew. Chem.,
Int.Ed., 48:1198, 2009.
[27] Torrie G. M. Nonphysical sampling distributions in monte carlo free-
energy estimation: Umbrella sampling. J.Comput.Phys., 23:187, 1977.
[28] Norrby P . O. Rationalizing the stereoselectivity of osmium tetroxide asym-
metric dihydroxylations with transition state modeling using quantum
mechanics-guided molecular mechanics. J.Am.Chem.Soc., 121:10186, 1999.
[29] Muller R. P . Ab initio calculations of free energy barriers for chemical reac-
tions in solution. J.Phys.Chem., 99:17516, 1995.
[30] Stewart J. J. P . Optimization of parameters for semiempirical methods v:
Modification of nddo approximations and application to 70 elements. J.
Mol.Model., 13:1173, 2007.
[31] Straatsma T. P . Estimation of statistical errors in molecular simulation cal-
culations. Mol.Phys., 57:89, 1986.
[32] Prasad B. R. Quantitative exploration of the molecular origin of the activa-
tion of gtpase. Proc.Natl.Acad.Sci.U.S.A., 110:20509, 2013.
[33] Lee F. S. Calculations of antibody-antigen interactions: Microscopic and
semi-microscopic evaluation of the free energies of binding of phosphoryl-
choline analogs to mcpc603. ProteinEng.,Des.Sel., 5:215, 1992.
[34] Lee F. S. Microscopic and semimicroscopic calculations of electrostatic en-
ergies in proteins by the polaris and enzymix programs. J. Comput. Chem.,
14:161, 1993.
[35] Lee C. T. Development of the colle-salvetti correlation-energy formula into
a functional of the electron-density. Phys. Rev. B: Condens. Matter Mater.
Phys., 37:785, 1988.
[36] Luzhkov V . Microscopic models for quantum mechanical calculations of
chemical processes in solutions: Ld/ampac and scaas/ampac calculations
of solvation energies. J.Comput.Chem., 13:199, 1992.
43
[37] Plotnikov N. V . Paradynamics: An effective and reliable model for ab ini-
tio qm/mm free-energy calculations and related tasks. J. Phys. Chem. B,
115:7950, 2011.
[38] Plotnikov N. V . Exploring, refining, and validating the paradynamics
qm/mm sampling. J.Phys.Chem.B, 116:10342, 2012.
[39] Plotnikov N. V . Computing the free energy barriers for less by sampling
with a coarse reference potential while retaining accuracy of the target fine
model. J.Chem.TheoryComput., 10:2987, 2014.
[40] Koen H. G. Verschueren, Frank Seljee, Henriette J. Rozeboom, Kor H. Kalk,
and Bauke W. Dijkstra. Crystallographic analysis of the catalytic mecha-
nism of haloalkane dehalogenase. Nature, 363(6431):693–698, June 1993.
[41] Shao Y. Advances in methods and algorithms in a modern quantum chem-
istry program package. Phys.Chem.Chem.Phys., 8:3172, 2006.
[42] Robert W. Zwanzig. High temperature equation of state by a perturbation
method. i. nonpolar gases. J.Chem.Phys., 22:1420, 1954.
44

Abstract (if available)

Abstract The enzymatic reactions play a key role in broad range of life processes. At the same time theoretical description of enzyme catalysis processes represents a challenging task due to the complexity of biomolecular systems of interest. In the present work we validate and enhance a newly developed method: Paradynamics (PD), targeted at evaluation of the free energy barriers of enzymatic reactions. We present the simulation of the sample molecular reaction—haloalkene dehalogenation, to show how our approach of using two state framework, involving reference and target potentials, allows us to reduce the complexity of task and brings an insight into the details of mechanisms of unexplored micro-scale phenomena.

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

Ab initio methodologies in studying enzymatic reactions

PDF

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

PDF

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

PDF

Development and application of free energy based methods to understand important biological phenomena

PDF

Understanding electrostatic effects in the function of biological systems

PDF

Quantitative computer-aided studies of artificial and enantioselective enzymes

PDF

Understanding organic and inorganic semiconductor surfaces by coherent nonlinear spectroscopy

PDF

Computational materials design for organic optoelectronics

PDF

Monochromatic exact solution to Maxwell's equations with finite energy in the vacuum

PDF

Photoinduced redox reactions in biologically relevant systems

PDF

Capturing the complexity of ion channels: simulations of long-time dynamics, conformational changes and the effect of the membrane potential

PDF

Studies on iron-chloride redox flow battery for large scale energy storage

PDF

Vorticity in superfluid helium nanodroplets

PDF

Electrocatalytic thiolate- and selenolate-based coordination polymers for solar energy conversion

PDF

Synthesis, characterization and reaction chemistry of polyazides and cyanometallates

PDF

Molecular orientation and structure of hydrogen bonds at water interfaces

PDF

Iridium and ruthenium complexes for catalytic hydrogen transfer reactions

PDF

Optical properties of lead and bismuth halide perovskites for photovoltaic applications

PDF

Computational analysis of the spatial and temporal organization of the cellular environment

PDF

High-frequency electron-electron double resonance techniques and applications

Asset Metadata

Creator Kupchenko, Ilya (author)

Core Title Two state free energy barrier evaluation technique for dehalogenation reaction

School College of Letters, Arts and Sciences

Degree Master of Science

Degree Program Chemistry

Publication Date 08/04/2016

Defense Date 06/01/2016

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

Tag free energy barrier reaction,OAI-PMH Harvest

Format application/pdf (imt)

Language English

Contributor Electronically uploaded by the author (provenance)

Advisor Warshel, Arieh (committee chair), Nakano, Aiichiro (committee member), Vilesov, Andrey (committee member)

Creator Email ikupchenko@gmail.com,kupchenk@usc.edu

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

Unique identifier UC11281504

Identifier etd-KupchenkoI-4737.pdf (filename),usctheses-c40-297489 (legacy record id)

Legacy Identifier etd-KupchenkoI-4737-0.pdf

Dmrecord 297489

Document Type Thesis

Format application/pdf (imt)

Rights Kupchenko, Ilya

Type texts

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

Access Conditions The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the a...

Repository Name University of Southern California Digital Library

Repository Location USC Digital Library, University of Southern California, University Park Campus MC 2810, 3434 South Grand Avenue, 2nd Floor, Los Angeles, California 90089-2810, USA