Friedrich-Alexander-Universität Erlangen-Nürnberg
Naturwissenschaftliche Fakultät II
Computer-Chemie-Centrum
QM/MM Docking and Simulations of FRET
DEN NATURWISSENSCHAFTLICHEN FAKULTÄTEN
DER FRIEDRICH-ALEXANDER-UNIVERSITÄT ERLANGEN-NÜRNBERG
ZUR
ERLANGUNG DES DOKTORGRADES
vorgelegt von
Frank Beierlein
aus Erlangen
Als Dissertation genehmigt von den Naturwissenschaftlichen Fakultäten der
Friedrich-Alexander-Universität Erlangen-Nürnberg
Tag der mündlichen Prüfung: 17.8.2005
Vorsitzender der Promotionskommission: Prof. Dr. D.-P. Häder
Erstberichterstatter: Prof. Dr. T. Clark
Zweitberichterstatter: Prof. Dr. Dr. h.c. R. van Eldik
Die vorliegende Arbeit entstand in der Zeit von September 2001 bis Juni 2005
am Computer-Chemie-Centrum der Friedrich-Alexander-Universität Erlangen-
Nürnberg.
Acknowledgements
“… I have prefaced these volumes with the names of my authorities. I have
done so because it is, in my opinion, a pleasant thing and one that shows an
honorable modesty, to own up to those who were the means of one’s
achievements …”
Adapted from Pliny the Elder, Natural History, Preface.
Prof. Dr. Tim Clark
Dr. Harald Lanig
Prof. Dr. S. Schneider
Dipl.-Chem. Olaf Othersen
Dipl.-Chem. Anselm Horn
Dr. Matthias Hennemann
Dr. Gudrun Schürer
Dr. Nico van Eikema Hommes
Dr. Jana Kurzawa
Dipl.-Phys. Dan Thomas Marian
Dipl.-Phys. Ciprian Roman
Dipl.-Biol. Heike Meiselbach
Dipl.-Chem. Jürgen Bulitta
Cand.-Chem. Marius Kroll
Isabelle Bundesmann
Clark group
Schneider group
HPC support group, RRZE
Eva Beierlein
My parents Karin and Rainer Beierlein
Deutsche Forschungsgemeinschaft (DFG): SFB 473 and SFB 583
Competence Network for Technical, Scientific High Performance Computing in
Bavaria (KONWIHR)
“Quia occupationibus tuis publico bono parcendum erat, quid singulis
contineretur libris, huic epistulae subiunxi summaque cura, ne legendos eos
haberes, operam dedi. Tu per hoc et aliis praestabis ne perlegant, sed, ut
quisque desiderabit aliquid, id tantum quaerat et sciat quo loco inveniat.”
C. Plinius Secundus, Naturalis Historia, Praefatio
Zusammenfassung
In dieser Arbeit werden zwei neue Ansätze aus der Computerchemie
vorgestellt, in denen gemischt quantenmechanische/molekülmechanische
(QM/MM) Verfahren verwendet werden, um tiefere Einblicke in die Struktur und
die Funktion von Protein-Ligand-Komplexen zu erhalten. Außerdem wird
gezeigt, wie umfangreiche Computersimulationen zur Erklärung komplexer
Befunde aus der Proteinspektroskopie verwendet werden können.
QM/MM Docking
Im ersten Teil der Arbeit wird ein kombinierter QM/MM-Docking-Ansatz
vorgestellt, mit dem Protein-Inhibitor-Komplexe untersucht werden. Mit Hilfe des
Docking-Programms AutoDock wird die Bindung kleiner organischer Moleküle
in der Bindetasche von Enzymen vorhergesagt. Hierbei wird der Ligand als
flexibles Molekül behandelt, das Protein ist jedoch starr und es kann keine
konformative Anpassung des Proteins (induced fit) als Reaktion auf die
Ligandbindung erfolgen. Die durch das Docking erhaltenen Strukturen werden
mit semiempirischen AM1 QM/MM-Optimierungen verfeinert. Hierbei wird eine
verbesserte Beschreibung des Bindungsmodus des Liganden erhalten.
Außerdem liefert die quantenmechanische Rechnung eine korrekte
Beschreibung der elektronischen Eigenschaften des Liganden. Durch die
Verwendung einer flexiblen Proteinumgebung bei den QM/MM-Rechnungen
können strukturelle Anpassungen des Enzyms als Reaktion auf die
Ligandbindung (induced fit) sogar im Rahmen einer einfachen
Geometrieoptimierung simuliert werden, und auf den Einsatz teuerer Methoden
(z.B. Moleküldynamik) kann verzichtet werden.
Die Methode wurde mit einem Satz durch Röntgenkristallographie
aufgeklärter Protein-Inhibitor-Komplexe validiert, deren Geometrien aus der
Proteindatenbank (PDB) entnommen wurden. Hierbei wurden auch
Metalloproteine, für die das Dockingprogramm keine Parameter hat, erfolgreich
gedockt. Hierzu wurden die Standard-Kraftfeldparameter um geeignete
Metallparameter ergänzt. Zusätzlich zu den Strukturen von zusammen mit dem
Inhibitor kristallisierten Proteinen wurden auch ligandfrei kristallisierte
Strukturen von HIV-1-Protease und Thrombin erfolgreich für Docking-
Experimente verwendet. Durch den Vergleich der erhaltenen Strukturen mit den
Ergebnissen, die unter Verwendung von mit einem Liganden kristallisierten
Proteinstrukturen erhalten wurden, konnte gezeigt werden, daß die Methode
begrenzte konformative Anpassungen eines Enzyms nach Bindung eines
Liganden auch durch einfache Geometrieoptimierungen simulieren kann und in
diesen Fällen auf teurere, rechenzeitintensivere Verfahren verzichtet werden
kann.
Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei
kristallisierte Enzym. Links: Ergebnis des Dockings mit AutoDock (starre Umgebung).
Rechts: Ergebnis der VAMP QM/MM-Optimierung.
Die quantenmechanische Beschreibung des Liganden liefert ein detailliertes
Bild seiner elektronischen Eigenschaften, z. B. seiner Polarisation durch die
Gruppen im aktiven Zentrum eines Enzyms.
Überlagerung der QM/MM-optimierten Struktur eines HIV-1-Protease-Inhibitors im
aktiven Zentrum des Enzyms mit der entsprechenden Struktur berechnet in der
Gasphase. Das durch die Bindetasche induzierte elektrostatische Potential wurde auf
die van der Waals Oberfläche der Moleküle projiziert. Nur die Extrema sind gezeigt, der
Bereich zwischen –1 und +7 kcal mol-1 wurde nicht abgebildet. Der resultierende
induzierte Dipolmomentvektor (1.64 D) wird durch den Pfeil repräsentiert.
Die Ergebnisse der vorgestellten Arbeit weisen darauf hin, daß ein QM/MM-
Moleküldynamik-Ansatz den Effekt des induced fit auch in solchen Fällen
simulieren kann, in denen größere konformative Anpassungen des Enzyms, wie
z. B. das Öffnen einer Bindetasche, erforderlich sind.
Simulationen zum Fluoreszenz-Resonanz-Energietransfer (FRET)
Im zweiten Teil der Arbeit wird eine Studie vorgestellt, mit Hilfe derer die
Ergebnisse zeitaufgelöster Fluoreszenzspektroskopie der Aminosäure
Tryptophan in Proteinen erklärt werden können. Bei derartigen Messungen wird
häufig mehr als eine Lebensdauer beobachtet, was gewöhnlich mit der Existenz
mehrerer Seitenkettenkonformationen von Tryptophan oder mehrerer
Proteinkonformationen erklärt wird („Rotamerenmodell“). In dieser Arbeit wird
das aus zwei Monomeren bestehende Tet-Repressor (TetR)-Protein untersucht,
welches ein interessantes Modellsystem für die Erforschung der Mechanismen
der transkriptionellen Regulation darstellt. In Komplexen des Tet-Repressors
mit Tetracyclin wird häufig Fluoreszenz-Resonanz-Energie-Transfer (Förster-
Energie-Transfer, FRET) von der Aminosäure Tryptophan zum Induktor
Tetracyclin beobachtet.
Fluoreszenz-Resonanz-Energie-Transfer (FRET) von Tryptophan 43 (blau) zum
Induktor Tetracyclin (grün). Nur eines der beiden TetR-Monomere ist gezeigt (PDB
Code 2trt).
Die Struktur und Dynamik sowie die spektroskopischen Eigenschaften des
TetR-Tetracyclin-Komplexes werden mit Hilfe eines gemischt
klassischen/quantenmechanischen Verfahrens untersucht. Eine klassische
Moleküldynamik (MD)-Simulation liefert die Eingabegeometrien für
semiempirische QM/MM CI-Rechnungen, die verwendet werden, um die zur
Simulation der spektroskopischen Eigenschaften der relevanten Chromophore
(Tryptophan und Tetracyclin) erforderlichen Größen zu berechnen.
Schematische Darstellung des kombinierten Moleküldynamik-QM/MM-Verfahrens.
Die Analyse der klassischen MD-Simulationen ergab für den Chromophor
Tryptophan 43, einer in den DNA-Bindeköpfen der beiden Monomere des TetR-
Proteins vorkommenden Aminosäure, zwei (Monomer 1) bzw. fünf (Monomer 2)
Rotamere. Zusätzlich zu den Änderungen der Seitenkettenwinkel, die die
Rotamere ineinander überführen, wurde eine schnelle Bewegung der
Tryptophan-Seitenkette beobachtet.
Seitenkettenwinkel χ1/ χ2 von Tryptophan 43 (Monomer 1) im Verlauf der MD-
Simulation.
Strukturen von Tryptophan 43, die die Clusterzentren der Clusteranalyse der χ1/ χ2-
Seitenkettenwinkel repräsentieren (links: Monomer 1, rechts: Monomer 2).
Für die weitere Analyse wurde das Monomer 1, das zwei Rotamere in der
MD aufweist, verwendet. Die vorgestellte Methode ist im Gegensatz zu den
meisten experimentellen Verfahren in der Lage, die Rotamere unabhängig
voneinander zu untersuchen.
Neben den stationären Absorptions- und Fluoreszenzspektren von
Tryptophan und Tetracyclin wurden auch zeitaufgelöste Fluoreszenzspektren
von Tryptophan simuliert. Hierbei wurden zum einen intrinsische (ohne
Löschung) Fluoreszenzabklingkurven von Tryptophan aus den Einstein-
Koeffizienten, die für jeden MD-Snapshot errechnet wurden, generiert. Die
Simulationen zeigen sowohl für die gesamte MD-Trajektorie (2 Rotamere) als
auch für die einzelnen Rotamere ein biexponentielles Abklingverhalten. Dieses
wird durch die schnelle Bewegung der Tryptophan-Seitenkette, die im Laufe der
Simulation zu einer Vielfalt von Tryptophan-Mikroumgebungen führt, verursacht.
Biexponentielle Fits der berechneten intrinsischen (ohne Löscher)
Fluoreszenzabklingkurven von Tryptophan 43 im Monomer 1. Rot: Gesamte MD-
Trajektorie (enthält beide Rotamere von Tryptophan 43 im Monomer 1), grün: Rotamer
1, blau: Rotamer 2.
Zusätzlich zur intrinsischen Fluoreszenz von Tryptophan 43 wurde auch
dessen Fluoreszenzlöschung durch Förster-Energietransfer (FRET) zum
Induktor Tetracyclin untersucht. Hierzu wurde für jeden MD-Snapshot die
Geschwindigkeitskonstante für den Fluoreszenz-Resonanz-Energie-Transfer
(FRET) berechnet. Es konnte gezeigt werden, daß sowohl die gesamte MD-
Trajektorie (2 Rotamere) als auch die einzelnen Rotamere ein biexponentielles
Abklingverhalten, verursacht durch FRET, aufweisen. Die schnellen
Bewegungen der Tryptophan-Seitenkettenwinkel führen auch hier zu einer
Vielfalt an relativen Donor-/Akzeptorgeometrien, was aufgrund der großen
Empfindlichkeit von FRET für derartige Geometrieänderungen das
biexponentielle Verhalten verursacht.
Die Analyse der relativen FRET-Geschwindigkeitskonstanten zeigt, daß die
beiden Rotamere eine unterschiedliche Wahrscheinlichkeit für eine Löschung
durch Energietransfer aufweisen. Entsprechend werden für die beiden
Rotamere unterschiedliche Lebensdauern gefunden.
Biexponentielle Fits der berechneten Fluoreszenzabklingkurven von Tryptophan 43
im Monomer 1 in Gegenwart des FRET-Akzeptors Tetracyclin. Rotamer 1 (grün) weist
längere, Rotamer 2 (blau) kürzere Lebensdauern auf als die Gesamttrajektorie (rot).
Diese Ergebnisse zeigen, daß das klassische „Rotamerenmodell“, welches
gewöhnlich zur Interpretation zeitaufgelöster Emissionsspektren von
Tryptophan dient, die tatsächlichen Verhältnisse zu stark vereinfacht, indem es
annimmt, daß jedes Rotamer einer diskreten Lebensdauer zuzuordnen ist. Die
vorgestellten Ergebnisse zeigen, daß das biexponentielle Abklingverhalten
zumindest im Fall des untersuchten TetR/Tetracyclin-Systems aus den
schnellen Bewegungen des Tryptophan-Chromophors und den resultierenden
mannigfaltigen Donor-/Akzeptorgeometrien folgt. Die unterschiedlichen FRET-
Geschwindigkeitskonstanten/Lebensdauern der Rotamere zeigen aber, daß das
Rotamerenmodell für Systeme mit FRET-Akzeptoren erweitert werden kann.
Abstract
The following presents two approaches that use combined quantum
mechanical/molecular mechanical (QM/MM) methods in order to gain deeper
insights into the structure and function of protein-ligand complexes and to
understand complex results of protein spectroscopy by extensive computer
simulations.
In the first part, a combined QM/MM docking approach for investigating
protein-inhibitor complexes is presented. Starting points for QM/MM
optimizations are generated with AutoDock. The subsequent semiempirical
AM1 QM/MM optimization of the complex obtained by the docking procedure
gives a more detailed description of the binding mode and the electronic
properties of the ligand. As a flexible protein environment is used in the QM/MM
optimizations, the method is able to simulate limited structural changes of the
enzyme upon binding a ligand, even within a simple geometry optimization. The
method was validated using a set of structurally known protein-inhibitor
complexes, whose crystallographic data were taken from the Protein Data
Bank. In addition to protein structures taken directly from complexes with the
inhibitors, structures of uncomplexed HIV-1-protease and thrombin were also
used successfully for QM/MM docking experiments. By comparing the resulting
structures with those obtained using protein structures from protein-inhibitor
complexes, it is shown that the method is able to simulate the effect of the
induced fit when a simple optimization is adequate to reproduce the protein
movement. Describing the ligand quantum mechanically gives a detailed view of
its electronic properties, e.g. its polarization within the active site of the enzyme.
This study suggests strongly that a QM/MM molecular dynamics approach will
be able to simulate the induced fit in general cases.
A computational model study designed to simulate the results of time-
resolved fluorescence spectra of tryptophan in proteins is presented in the
second part. In such measurements, the occurrence of more than one
fluorescence lifetime is generally attributed to the existence of several
tryptophan rotamers and/or structural conformations of the protein structure.
The system chosen for this initial study is the tetracycline repressor (TetR)
protein, an interesting model system for the investigation of the mechanisms of
transcriptional regulation. Fluorescence resonance energy transfer (FRET) from
tryptophan to tetracycline is frequently observed in complexes of the Tet-
repressor with the antibiotic tetracycline. A combined classical/quantum
mechanical approach is used to model the structure and the spectroscopic
properties of the TetR-tetracycline complex. A classical molecular dynamics
simulation provides input geometries for semiempirical QM/MM CI calculations,
which are used to calculate tryptophan absorption and fluorescence spectra as
well as fluorescence resonance energy transfer (FRET) rate constants. It is
shown how a biexponential tryptophan fluorescence decay can be obtained
from a molecular dynamics (MD) trajectory containing two tryptophan rotamers
by calculating the FRET rate constant for every MD snapshot. The results
indicate that the classical “rotamer model”, used to explain the
multiexponentiality of time-resolved tryptophan emission spectra, can be
extended to systems with FRET acceptors present in the protein matrix.
Table of Contents
1 Introduction .....................................................................................................1
1.1 Some General Remarks on Computational Chemistry ...............................1
1.2 Methods......................................................................................................4
1.2.1 Quantum Mechanics – Semiempirical Molecular Orbital Theory..........4
1.2.2 Classical Force Fields ........................................................................11
1.2.3 Molecular Dynamics...........................................................................13
1.2.4 QM/MM Methods................................................................................19
1.2.5 Molecular Docking..............................................................................21
2 QM/MM Docking: A New Protocol and its Evaluation for Known Test Systems .........................................................................................................23
2.1 Introduction ...............................................................................................23
2.1.1 Combining Docking Techniques with other Computational
Methods ...............................................................................................23
2.1.2 The Protein-Inhibitor Systems Investigated........................................24
2.1.3 Overview of the Method .....................................................................25
2.2 Computational Details...............................................................................26
2.3 Results and Discussion ............................................................................32
2.3.1 Validation of the QM/MM Docking Approach......................................32
2.3.2 Simulating the Induced Fit..................................................................34
2.3.3 Polarization of the Ligand by the Environment ...................................39
2.4 Conclusions and Outlook..........................................................................42
3 Simulating FRET: A Combined Molecular Dynamics-QM/MM Approach........................................................................................................44
3.1 Introduction ...............................................................................................44
3.1.1 The Tetracycline Class of Antibiotics..................................................44
3.1.2 The Tet-Repressor/Operator (TetR/tetO) System...............................47
3.1.3 Fluorescence Resonance Energy Transfer (FRET)............................50
3.1.4 Tryptophan Fluorescence ...................................................................51
3.1.5 Time-Resolved Tryptophan Emission: The Rotamer Model and
Alternative Approaches ........................................................................53
3.1.6 Spectroscopic and Computational Investigations of the Tet
Repressor/Tetracycline System............................................................57
3.1.7 Overview of the Method......................................................................58
3.1.8 Molecular Dynamics Simulations for the Simulation of Protein
Spectra and FRET................................................................................60
3.2 Computational Details ...............................................................................63
3.2.1 Validation: Glycyltryptophan Absorption and Emission Spectra .........63
3.2.2 Simulating FRET in TetR: Molecular Dynamics Simulations ..............66
3.2.3 Simulating FRET in TetR: Quantum Mechanical Simulations.............69
3.2.4 Calculation of FRET Data ...................................................................72
3.3 Results and Discussion.............................................................................78
3.3.1 Validation: Glycyltryptophan Absorption and Emission Spectra .........78
3.3.2 Simulating FRET in TetR: Molecular Dynamics Simulations ..............81
3.3.3 Simulating FRET in TetR: Quantum Mechanical Simulations.............86
3.3.4 Simulating FRET in TetR: Tryptophan Quenching by FRET...............96
3.4 Conclusions.............................................................................................102
3.5 Outlook....................................................................................................104
4 References ...................................................................................................105
1
1 Introduction
1.1 Some General Remarks on Computational Chemistry
“Computational chemistry has existed for half a century, growing from the
province of a small nucleus of theoretical work to a large, significant component
of scientific research. By virtue of the great flexibility and power of electronic
computers, basic principles of classical and quantum mechanics are now
implemented in a form which can handle the many-body problems associated
with the structure and behavior of complex molecular systems” (John A. Pople,
1997).[1]
In the last four decades a new field of chemistry has gained considerable
importance: computational chemistry (which includes applied theoretical
chemistry, molecular modeling and chemoinformatics). The foundations of
computational chemistry were laid with the development of theoretical models in
the early part of the 20th century. Although the field was restricted to a few
dedicated theoreticians in the early days, the rapid development of computer
technology and the availability of easy-to-use commercial and academic
software has now opened it to a large number of more application-oriented
scientists.[2-4]
The use of computers in chemistry helps save time – even computations of
large systems can be faster than experiments – and costs – as no chemicals
are required that otherwise would have to be disposed of, a “green” aspect of
this new field of chemistry. For experimental chemists, modeling studies can
give valuable guidance, especially as they allow compounds that exist only
virtually to be examined. Furthermore, even complex reaction mechanisms can
be investigated by calculations. Visualizing the results helps to understand
complex chemical issues. Therefore, computer simulations are essential in all
fields of modern chemistry: life sciences, materials sciences as well as the
classical disciplines such as inorganic and organic chemistry benefit from
computational studies.[5]
2
The increasing importance of computational chemistry results is visible from
the rising number of publications and citations. Figure 1 gives an overview of
the number of search results in Science Citation Index Expanded[6] for selected
keywords (AMBER, CHARMM, MM3, QM/MM, AM1, PM3, B3LYP and 6-31G).
The analysis of the keywords shows that more compute-expensive methods,
especially density functional theory (B3LYP), but also methods for treating large
protein systems with force fields (AMBER) or combined quantum
mechanical/molecular mechanical (QM/MM) methods become more important
today. The significance of computational chemistry is also emphasized by the
Nobel prizes awarded to theoretically oriented chemists (Robert S. Mulliken,
1966; John A. Pople, 1998; Walter Kohn, 1998).[7]
Figure 1. Number of search results in Science Citation Index Expanded[6] for the
keywords AMBER, CHARMM, MM3, QM/MM, AM1, PM3, B3LYP and 6-31G. The sum
is drawn in red.
Computational chemistry is also crucial for a field that plays a pivotal role in
the global economy: pharmaceutical research. Remarkable process has been
achieved in the different disciplines involved in this field since 1950. It
culminated in the elucidation of the sequence of the human genome. Today,
3
genomics, proteomics, bioinformatics, combinatorial chemistry and ultrahigh-
throughput screening provide an enormous number of targets and data. Hence,
the application of data-mining methods and virtual screening have a growing
impact on the validation of “druggable” targets, lead finding and the prediction of
ADMET (absorption, distribution, metabolism, excretion and toxicity) profiles.[8]
However, the development periods until a new drug comes to the market are
still long (12-15 years) and are coupled with high financial effort (up to $ 880
million).[8-10] Despite the introduction of combinatorial chemistry and the
establishment of high-throughput screening (HTS), the average number of new
chemical entities (NCEs) introduced annually to the world market was only
about 37 for the period 1991-1999.[8] In a report by the Boston Consulting
Group, it is concluded that genomics and other modern technologies, like
computational chemistry, could reduce the cost of developing a new drug from $
880 million by about $ 300 million and could also shorten the development time
by two years.[9, 10] In particular, in silico methods are expected to speed up the
drug discovery process, to provide a quicker and cheaper alternative to in vitro
tests, and to reduce the number of compounds with unfavorable
pharmacological properties in an early stage of drug development.[8, 9, 11]
The following thesis gives two examples of the use of molecular modeling
techniques to investigate complex biomolecular processes: the interaction of
proteins with potential drugs and the transfer of excitation energy from a donor
group in a protein to an acceptor. Especially the latter simulations are essential
in order to interpret the experimental results, which cannot be explained
satisfactorily by traditional means. This is a good example of how the synergy
between theoretical simulations and experiments has vastly accelerated
progress in many areas of interest, even the “wettest of wet chemists”[3] can
now benefit from a comparison of experimental results with those obtained from
simulation/theory.
4
1.2 Methods
1.2.1 Quantum Mechanics – Semiempirical Molecular Orbital Theory
Starting point for a discussion or a historical review of quantum mechanics[4,
12-17] is the Schrödinger equation.[18] The full, time-dependent form is
(1).
Provided that the external potential V is time-independent, Eq. (1) can be
written in the more familiar way
(2).
The most common way to write the Schrödinger equation is to abbreviate the
left side of the equation using the Hamiltonian operator H .
(3)
This results in the shortest possible notation of the Schrödinger equation
(4).
The Schrödinger equation cannot be solved exactly for more than two
particles because of the three-body-problem.[4] However, taking the Born-
Oppenheimer approximation[19] into account – separation of the electronic
motion from that of the nucleus – the Schrödinger equation itself can be
separated into two individual equations, where the one related to the nucleus is
often omitted. The resulting equation refers only to the electronic factors.
(5)
Another important approximation is given by the Hartree-Fock (HF)
approach.[20] It plays a significant role in modern quantum chemistry and is
usually a starting point for more accurate methods. It overcomes the main
problem, that Eq. (5) can only be solved exactly for the smallest molecule, +2H .
( ) ( )2 2 2 2
2 2 2
,,
2r t
V r t im x y z t
⎧ ⎫ ∂Ψ⎛ ⎞∂ ∂ ∂⎪ ⎪− + + + Ψ =⎨ ⎬⎜ ⎟∂ ∂ ∂ ∂⎪ ⎪⎝ ⎠⎩ ⎭
hh
( ) ( )2
2
2V r E r
m⎧ ⎫
− ∇ + Ψ = Ψ⎨ ⎬⎩ ⎭
h
22ˆ
2H V
m= − ∇ +h
H EΨ = Ψ
ˆe e e eH EΨ = Ψ
5
The complicated many-electron problem is simplified to a one-electron problem
by replacing the N-electron wavefunction by a set of N one-electron
wavefunctions. This means that one electron then interacts mathematically with
N-1 one-electron wavefunctions, i.e. the electron ‘sees’ an averaged field of the
remaining electrons. In polyelectronic systems, there are more interactions,
such as Coulomb- and exchange interactions, than in single electron systems.
These can be described by integrals. By concentrating on one electron, a
possible notation of the Hartree-Fock equation with the core Hamilton operator
Hcore, the Coulomb operator J, the exchange operator K and χi, representing the
spin orbital is
(6).
A simpler notation uses the Fock operator, an effective one-electron
Hamiltonian if for the electron in the polyelectronic system
(7)
or written in the standard eigenvalue form
(8).
Assuming that each electron moves in a ‘fixed’ field comprising the nuclei
and the other electrons has important implications in finding a way to solve the
equation for one electron, which naturally will affect the solutions for other
electrons. The general strategy is called the self-consistent field (SCF)
approach. First, a set of trial solutions χi to the Hartree-Fock eigenvalue
equations is obtained and used to calculate the Coulomb and exchange
operators. Then the Hartree-Fock equations are solved, which gives a second
set of solutions χi. These are used in the next iteration. The SCF method refines
the individual electronic solutions that correspond to lower and lower total
energies until a point is reached where the results for all electrons are
unchanged, which means that the field is now self-consistent.
( ) { } ( ) ( )1 1
ˆ ˆ ˆ1 (1) (1) 1 1N N
corej j i ij j
j jH J K χ ε χ
= =
⎡ ⎤+ − =⎢ ⎥
⎣ ⎦∑ ∑
i i i jf χ ε χ=
i i ij jj
f χ ε χ=∑
6
The Hartree-Fock equations are usually solved in different ways for atoms
and for molecules. Assuming a spherically symmetrical electron distribution, the
equations can be solved numerically, but are not particularly useful. Thus,
analytical approximations with the following form of the wavefunction ψ can be
obtained, where Y is a spherical harmonic and R is a radial function.
(9)
In 1930, Slater[21] suggested an analytical form for the radial functions, which
are now universally known as Slater-type orbitals (STOs).
(10)
Slater provided a series of empirical rules for choosing the orbital exponents
ζ, which are given by
(11).
Z is the atomic number, σ a derived shielding constant and n* is the effective
principal quantum number, which takes the same values as the true principal
quantum number for n = 1, 2, or 3, but has a value of 3.7, 4.0 and 4.2 for n = 4,
5 and 6, respectively.
As the direct solution of the Hartree-Fock equation is not a practical
proposition, an alternative approach is necessary. This is to write a spin orbital
as a linear combination of single-electron orbitals:
(12)
The one-electron orbitals φν are commonly called basis functions and often
correspond to the atomic orbitals.
The derivation of the Hartree-Fock equation for a closed-shell system was
first given by Roothaan[22] and Hall.[23] Finally written as a matrix equation, the
Roothaan-Hall equation is
( ) ( , )nl lmR r Yψ θ φ=
( ) ( ) 1/ 21/ 2 1( ) 2 2 !n n rnlR r n r e ζζ −+ − −= ⎡ ⎤⎣ ⎦
*Zn
σζ −=
1
K
i icν νν
ψ φ=
=∑
7
(13).
Unlike ab initio methods, semiempirical methods neglect or approximate
some of the integrals in the Fock matrix F to save time in the calculations. This
is done by reducing the number of electrons to the valence electrons only; the
remaining electrons are subsumed into the nuclear core. Furthermore, a
common feature is to set the overlap matrix S equal to the identity matrix I, which reduces the overlap matrix to 1 (diagonal elements) or 0 (off-diagonal
elements). This leads to the simplified Roothaan-Hall equation
(14).
This, however, does not mean that all overlap integrals are set to zero in the
calculation of the Fock matrix elements; even in the simplest semiempirical
models some overlaps are specifically included. Many semiempirical methods
are based on the ZDO (zero differential overlap) approximation, in which the
overlap between pairs of different orbitals is set to zero.
The Roothaan-Hall equations are not applicable to open-shell systems such
as radicals, which contain one or more unpaired electrons. Several approaches
have been devised to treat open-shell systems. The first of these is the spin-
restricted Hartree-Fock (RHF) theory, which uses combinations of singly and
doubly occupied molecular orbitals. The closed-shell approach as described
above is a special case in RHF theory. The doubly occupied orbitals use the
same spatial functions for electrons of both spins, α and β. The alternative
approach is the spin-unrestricted Hartree-Fock (UHF) theory of Pople and
Nesbet,[24] which uses two distinct sets of molecular orbitals: one for electrons
of α spin and one for electrons of β spin.
In the SCF method, electrons are assumed to move in an average potential
of the other electrons, which also means that the present position of an electron
is not influenced by the presence of a neighboring electron. Naturally, electrons
tend to avoid each other more than described by Hartree-Fock theory, which
gives rise to a lower energy. The difference between the exact energy and the
Hartree-Fock limit is called correlation energy. A way to incorporate the electron
= FC SCE
= FC CE
8
correlation is to perform configuration interaction (CI), where excited states are
included in the description of an electronic state. By taking the excited states
into account, the overall wavefunction can be described by a linear combination
of ground and excited state wavefunctions.
(15)
Ψ0 is the wavefunction obtained from solving the Hartree-Fock equations, Ψ1,
Ψ2, etc. are wavefunctions that represent configurations derived by replacing
one or more of the occupied spin orbitals by a virtual spin orbital. The energy of
the system is then minimized in order to determine the coefficients c0, c1, etc.,
using a linear variational approach, just as for the single-determinant
calculation. A CI calculation is therefore very complex and the total number of
ways to permute N electrons within K orbitals is (2K!)/[N!(2K-N)!]. Even at rather
small values of N and K, this leads to an enormous number of configurations.
Although it is the best way to describe a system, a full CI calculation is not
practicable for most systems, thus reduction to smaller CI calculations is
required. For example, in configuration interaction singles (CIS), only
wavefunctions that differ from the Hartree-Fock wavefunction by a single spin
orbital are included. Other examples are CID or CISD calculations, referring to
double substitutions or single and double substitutions, respectively. A final
word should be said about Brillouin’s theorem,[25] which states that single
excitations do not mix directly with single determinant, closed-shell ground-state
wavefunctions, only higher excitations, e.g. double, triple, do. However, since
single excitations mix with doubles, singles have an indirect effect on ground-
state wavefunctions.
In 1965, the first method to implement the ZDO approximation in a practical
way was published by Pople et al. This first three-dimensional self-consistent
field method was called CNDO[26, 27] (complete neglect of differential overlap).
Prior to this, almost only qualitative π-molecular orbital methods existed, e.g. the
PMO theory (perturbational molecular orbital), developed by Dewar in 1952.[28]
The only quantitative method so far was the Pariser-Parr-Pople method
(PPP).[29, 30] The differences to ab initio methods were obvious from the name of
the method, which already indicated approximation.
0 0 1 1 2 2 ...c c cΨ = Ψ + Ψ + Ψ +
9
Although these theories were initially successful, the lack of ability to
calculate e.g. differences in spin states led to the development of new methods.
In the successor to CNDO, INDO (intermediate neglect of differential
overlap),[31] the constraint of CNDO that monocentric two-electron integrals are
equal, was lifted. Such methods were unable to represent the repulsion
between two lone pairs correctly because the 4c2e integrals were calculated
using s-orbitals, even those involving p-orbitals. This gave need for another
model, which finally resulted in the development of the NDDO methods (neglect
of diatomic differential overlap),[26, 27] where all monoatomic overlaps were taken
into account.
The programs still had the major problem that geometry optimization was not
very practical because of the poor results. In this situation, Dewar and co-
workers took over and published the MINDO/3 method in 1975.[32] The most
significant improvement over earlier methods was the careful optimization of the
parameters used in the approximation. Although the basic equations in
MINDO/3 were similar to those in INDO, the origin of the parameters was
different. MINDO/3 contained several adjustable parameters instead of atomic
spectral data, which were optimized to give the best fit to experimental data of
molecules. Based on the INDO method, MINDO/3 was not able to represent
lone-pair lone-pair interactions adequately.
Not willing to accept this situation, Dewar and Thiel published a new method
– MNDO (modified neglect of differential overlap) – based on the NDDO method
in 1977.[33] The parameters in MNDO were obtained from experimental data and
optimized to reproduce several atomic and molecular properties, e.g. heats of
formation, dipole moments, ionization potentials and molecular geometries.
Instead of using diatomic parameters in the resonance integral and core-core
repulsion terms, which severely limited extending MINDO/3 to new elements,
MNDO used entirely monoatomic parameters. This new perspective led to the
development and derivation of parameters for most of the 2nd and 3rd row
elements over the next decade.
In the early 1980s, the first MOPAC versions were published in the Quantum
Chemistry Program Exchange.[34, 35] Both MINDO/3 and MNDO were
implemented and the program allowed various operations, such as geometry
10
optimization, transition state localization, gradient minimization and vibrational
frequency calculations. Over the next few years, the predictive power of
MOPAC rose as the methods improved their ability to predict molecular
properties. This program made semiempirical methods still more widely used,
but simultaneously various limitations became apparent. The most significant
error was the inability of MNDO to reproduce hydrogen bonds, a special interest
and major structure-determining feature in biological systems.
By modifying the core-core repulsion with additional Gaussian terms and new
parameters, the new method AM1 (Austin model 1) was developed by the group
of Dewar in 1985.[36] In another method – based on the same modifications as
AM1 and called PM3 (parametric method number 3)[37] – Stewart used only
adjustable parameters, which were optimized using a large set of molecular
reference data and an automated parameterization procedure.
Over the last couple of years new methods have included d-orbitals into
semiempirical methods in order to be able to treat transition metals or
hypervalent molecules more accurately. Some of these approaches are SAM1
(semi ab initio method version 1),[38, 39] MNDO/d,[40-46] AM1*,[47, 48] and
AM1(d).[49] Some of these approaches are available commercially, but have not
yet been published. Additionally, more applications have been made to a wider
range of ground- and excited-state atomic and molecular properties, to
polarizabilities and hyperpolarizabilities, to reaction mechanisms and also to
large systems such as proteins and enzymes.
The popularity and widespread use of AM1, PM3 and MNDO is mainly due to
their implementation in semiempirical program packages like VAMP,[50, 51]
MOPAC[52] and AMPAC.[53] There are also other methods in use, e.g. the
SINDO1 program developed by Jug and Nanda[54] or the ZINDO program by
Zerner and co-workers.[55-57] ZINDO for example, can perform a wide variety of
semiempirical calculations and has been particularly useful for calculations on
transition metal and lanthanide compounds and for predicting molecular
electronic spectra.
Quantum mechanical techniques are typically only used for moderate-sized
molecules (up to 100 atoms for ab initio or density functional theory (DFT) and
11
up to 500 for semiempirical MO calculations), because they scale badly, i.e. a
calculation for twice as large a molecule does not require twice as much
computer time and resources but rather 2n times as much (n varies between
about three for DFT calculations to four for Hartree-Fock and very large
numbers for ab initio techniques with explicit treatment of electron-correlation.[12,
17] Thus, large biological systems such as protein or DNA, which can consist of
tens of thousands of atoms, cannot be treated with conventional quantum
mechanics. Although some linear scaling methods for QM calculations have
been developed,[12, 17, 58-60] large systems are usually computed using classical
force fields or hybrid quantum mechanical/molecular mechanical approaches
(see following sections).
1.2.2 Classical Force Fields
In contrast to quantum mechanics, force field methods – also known as
molecular mechanics – ignore the electronic motions and calculate the energy
of a system as a function of the nuclear positions only.[4] As in the case of
semiempirical approaches, several assumptions must be made. One of them is
the well known Born-Oppenheimer approximation, without which it would be
impossible to contemplate writing the energy as a function of the nuclear
coordinates at all. Although force fields are based on rather simple functions,
e.g. Hooke’s law, and models of interaction, their accuracy sometimes can
compete with that of the highest-level quantum-mechanical calculations, in a
fraction of computer time. However, of course, molecular mechanics cannot
provide properties that depend upon the electronic distribution in a molecule.
Parameterization plays a key role in the development of force fields, but
parameters derived and tested on small molecules can be used for a wide
range of larger molecules such as polymers.[4]
Most of the force fields in use today can be interpreted in terms of a relatively
simple four-component picture of the intra- and intermolecular forces within the
system.[4] Penalty functions are applied to deviations of bonds and angles away
from their ‘reference’ or ‘equilibrium’ values. There is also a term that accounts
for barriers of rotation about chemical bonds (i.e. describing the energetics of
12
twisting the 1,4 atoms attached to the bonds formed by the atoms 2-3)[61] and
finally, terms to express the interaction between non-bonded parts of the
system are also implemented. These are the basic terms a force field is built of,
whereas more sophisticated force fields can contain additional terms. A
functional form to describe the potential energy as a function of the positions r
of N particles is[4]
(16)
or, to give a more general representation of the different terms used in a
force field:
(17)
As we can see in Eq. (16), the functions used in force fields are rather
simple. The first two terms are harmonic potential functions, which refer to the
bond lengths and angles, respectively. The third term is the torsional potential,
taking rotation barriers into account, whereas the fourth term describes the non-
bonded interactions between atoms separated by at least three bonds or
located in different molecules. In simple force fields the non-bonded term
usually uses Coulomb’s law to describe the electrostatic interactions and a
Lennard-Jones potential for the van der Waals interactions. The van der Waals
interaction described here consists of attracting induced dipole and dispersion
(London) interactions and repulsive exchange interactions.
Since it uses quite simple functions, a force field must be well parameterized
for the type of calculation intended in order to give good results. While it may be
tempting to try to predict other quantities that were not included in the
parameterization process, it is not necessarily a failing if a force field is unable
to do so. One important point worth thinking of is that force fields are empirical.
There is no ‘correct’ form of a force field, nevertheless most of the force fields in
use have a similar form.
2 2,0 ,0
12 6
1 1 0
( ) ( ) ( ) (1 cos( ))2 2 2
44
N i i ni i i i
bonds angles torsions
N Nij ij i j
iji j i ij ij ij
k k VV r l l n
q qr r r
θ θ ω γ
σ σε
πε= = +
= − + − + + −
⎛ ⎞⎡ ⎤⎛ ⎞ ⎛ ⎞⎜ ⎟⎢ ⎥+ − +⎜ ⎟ ⎜ ⎟⎜ ⎟ ⎜ ⎟⎜ ⎟⎢ ⎥⎝ ⎠ ⎝ ⎠⎣ ⎦⎝ ⎠
∑ ∑ ∑
∑∑
( ) bonds angles torsions non bondedV r V V V V −= + + +
13
A concept that is common in most force fields is that an atom type for each
atom must be defined in the input. While quantum mechanical methods only
need the atomic numbers of the nuclei, the nuclear coordinates, the overall
charge and the spin multiplicity, a molecular mechanics program requires the
definition of an atom type, which contains information of the hybridization state
or sometimes the local environment. This can lead to quite an extensive
parameterization for some atoms. For example, the MM2[62]/MM3[63] force fields
of Allinger and co-workers, which are widely used to calculate ‘small’ molecules,
distinguish between different types of carbon: sp3, sp2, sp, carbonyl,
cyclopropane, radical, cyclopropene and carbonium ion. In the AMBER force
fields of Kollman and co-workers (Weiner et al.,[64, 65] Cornell et al.[66]) the
carbon atom at the junction between a six and a five-membered ring (e.g. in the
amino acid tryptophan) is assigned an atom type that is different from the
carbon atom in an isolated five-membered ring, such as histidine, which in turn
is different from the atom type of a carbon atom in a benzene ring. Furthermore,
AMBER uses different types for a histidine amino acid, depending on its
protonation state.[4]
It is often found that force fields designed to model specific classes of
molecules (such as proteins and nucleic acids, in the case of AMBER) use
more specific atom types than force fields for more general purpose use.[4]
1.2.3 Molecular Dynamics
As shown in the previous section, the energy of a molecule can be calculated
easily using a simple mechanical model. In many cases, stationary points, e.g.
local minima or transition states, are obtained by performing a geometry or a
transition-state optimization. In the case of a geometry-optimization, the
optimizer searches for a local minimum near the starting point on the energy
surface. The real situation, however, is different: molecules do not stand still,
not even at 0 K (which is illustrated by the existence of the zero-point energy in
quantum mechanics). An approach for the simulation of the dynamic behavior of
molecules are molecular dynamics (MD) simulations, which sample the phase
space of a molecule (defined by the positions of the atoms and their velocities)
14
by integrating Newton’s equations of motion. Because MD simulations account
for thermal motion, the molecules simulated may possess enough thermal
energy to overcome potential barriers, which makes the technique in principle
suitable for conformational analyses, especially in the case of large
molecules.[67]
Starting from initial velocities (see below), the mechanical model is in
permanent movement and the movement of each atom influences the motion of
the other atoms, resulting in a complex mathematical problem of coupled
differential equations. This can only be solved by numerical integration and not
analytically. The protocol must allow sampling the fastest possible movement of
an atom within the simulation system. These are the vibrations of bonds
involving a hydrogen atom (e.g. a C-H bond in a methyl group), taking between
10 and 100 fs. Therefore, the integration steps must be at least one order of
magnitude smaller than the fastest motion, i.e. about 1 fs. Other motions in
proteins occur on a much longer time scale, e.g. 100 ps for rotations around
single bonds,[67] 4 ps for the reorientation of water, 10 ps - 100 ns for domain
motions, 100 µs - 1 s for aromatic ring-flips and allosteric shifts[68] and up to 1 s
for protein folding processes.[67] The latter process would require 1015
integration steps. Such simulations are of course not possible with the current
algorithms and computer power.
The starting point of an MD simulation is an initial set of coordinates (from X-
ray crystallography or NMR investigations) This structure is normally geometry
optimized prior to the MD simulation in order to remove bad contacts and initial
strain (usually van der Waals overlap of non-bonded atoms), which might
disturb the subsequent MD. Normally the same force-field function is used for
geometry optimization and the energy evaluations during the MD, e.g. Eq. (16).
After assigning velocities vi, which typically represent a low-temperature
Maxwell distribution, the simulation is started by calculating the acceleration ai
for each atom i according to Newton’s law Fi = mi ai, as written in Eq. (18).[67]
(18)
2
2
( )i ii i i i
i
V x xF m a mx t
∂ ∂− = = =∂ ∂
15
Fi is the force acting on the Cartesian coordinate xi for i=1, …, 3N for the N
atoms in the molecule or molecular system, mi is the atomic mass of atom i, V is
the potential energy function given in Eq. (16), and t is the time. The total
energy E of the system is the sum of all kinetic (1/2 mv2) and potential energy
V(x) contributions:
(19)
The position xi(t + ∆t) of an atom at the time (t + ∆t) can be calculated based
on the known position xi(t) at the time t, e.g. by the Leapfrog[69] algorithm (which
is a modification of the well-known Verlet[70] algorithm) given in Eq. (20).
(20)
By defining the atomic displacement and velocities at a distinct time, these
equations can be integrated.[67]
The temperature T of a system is related to the mean kinetic energy of all
atoms N via Eq. (21), where kB is the Boltzmann constant and <vi2> the average
of the squared velocities of atom i.
(21)
The computer time required for an MD simulation is proportional to the
square of the number of atoms in the system. The calculation of the non-
bonded interactions in Eq. (16) is the most time-consuming part of the energy
evaluations. One way to speed up calculations is the introduction of so-called
cutoffs to the van der Waals and Coulomb interactions. These cutoffs reduce
the number of non-bonded interactions by simply defining a maximum distance
at which two atoms are allowed to interact through space. Several cutoff
schemes have been introduced, from a simple sphere to switched or shifted
cutoffs.[67]
The fastest motions in a molecule are vibrations of bonds involving hydrogen.
If these are to be reproduced in an MD simulation, the integration time step
21 ( )2
xE m V xt
∂⎛ ⎞= +⎜ ⎟∂⎝ ⎠
2
2
( ) ( ) ( 2)
( )( 2) ( 2)
i i i
ii i
x t t x t v t t t
x tv t t v t t tt
+ ∆ = + + ∆ ∆
∂+ ∆ = − ∆ + ∆∂
21 32 2i i Bm v Nk T=∑
16
must be 1 fs or less. As it is often unnecessary to include hydrogen motions, the
time step can be increased by a method called SHAKE, which constrains the
motion of hydrogen atoms and therefore allows integration steps of up to 2 fs.[71]
In calculations with explicit solvent, water models like TIP3P[72] can be used,
which are intended to reproduce the experimentally known bulk water properties
with limited computational effort (the atoms of these water molecules have fixed
charges and a fixed relative orientation). Another way to speed up calculations
is the use of united-atom force fields which decrease the number of atoms by
contracting non-essential hydrogens (e.g. of methyl or methylene groups) to
pseudo-carbon (united) atoms with increased van der Waals radii and modified
charges.[67]
Considering the influence of the solvent is necessary as enzymes normally
work in an aqueous environment. The most straightforward approach is to use
explicit solvent molecules in the simulation system using a supermolecule
approach (see below). A way to consider the dielectric properties of the solvent
(in most cases water) is to re-scale the dielectric permittivity of free space ε0 by
a factor D. This damps the long-range electrostatic interactions according to ε =
Dε0. Using the macroscopic value for water (D = 78.0), ε then amounts to
78.0ε0. An alternative way is to introduce a distance dependence into the
electrostatic interactions by defining an effective dielectric constant ε = Dε0rij,
which modifies the Coulomb term in Eq. (16) according to
(22).
By using an effective, distance-dependent dielectric constant, the ability of
bulk water to shield electrostatic interactions can be mimicked without the
presence of explicit solvent molecules.[67] Nevertheless, these kinds of
simulations lack the influence of solvent molecules in terms of hydrogen bonds
between solvent molecules and the protein and the structural influence of the
solvent on the molecule of interest, e.g. a protein. Other approaches that do not
use explicit solvent molecules are continuum solvent models, like the
generalized Born (GB) method.[67]
21 1 04
N Ni j
couli j i ij
q qE
D rπ ε= = +
⎛ ⎞= ⎜ ⎟⎜ ⎟
⎝ ⎠∑ ∑
17
The use of explicit water molecules is the most common way to model the
bulk properties of the solvent correctly. The definition of a (in most cases
rectangular) box that contains the molecule(s) of interest (e.g. a protein) and
also a large number of water molecules causes the problem that water
molecules at the solvent/vacuum boundary would evaporate and therefore get
lost in the adjacent vacuum. This problem led to the definition of so-called
periodic boundaries. The box containing the system of interest and the solvent
molecules is surrounded by its periodic images and therefore no longer has a
border with the vacuum. A solvent molecule that leaves the box on one side
enters the box at the same time on the other side. One important condition of
periodic boundary conditions (PBC) is that an atom of the real molecule must
not interact with another real atom and its image at the same time (minimum
image convention). Therefore, spherical cutoffs for the non-bonding interactions
should be defined that must be smaller than half the smallest box dimension.[67]
The trajectories (series of structures) produced by MD simulations can
represent different statistical ensembles, depending on the variables that are
held constant during the simulation. In the so-called microcanonical ensemble
NVE, the number of particles N, the volume V and the energy is defined to be
constant and the temperature is calculated. In most cases, however, the
temperature, which corresponds to the kinetic energy of the particles, is held
constant. One example of constant temperature MD is the canonical ensemble
NVT. Temperature is controlled by adding or removing kinetic energy, which is
done by coupling the system weakly to an external heat bath. According to
Berendsen,[73] the rate of temperature change is proportional to the temperature
difference between bath (T0) and system:[67]
(23)
The velocities v are scaled by the factor λ to come closer to the desired
temperature T0:
(24)
[ ]0( ) 1 ( )
T
dT t T T tdt τ
= −
1 2
01 1( )T
TtT t
λτ
⎡ ⎤⎛ ⎞∆= + −⎢ ⎥⎜ ⎟⎝ ⎠⎣ ⎦
18
T(t) is the actual temperature at the time t, ∆t is the integration time step and
the relaxation time τT represents the strength of the coupling (smaller values
mean stronger coupling to the bath).[67]
It is common practice to sample an isothermal isobaric ensemble (NPT),
which represents normal laboratory conditions. Similarly to temperature control,
the system is coupled to an external bath with the desired target pressure P0.
By rescaling the dimensions of the periodic box and the atomic coordinates by
the factor µ at each integration step ∆t, the volume of the box and the forces of
the solvent molecules acting on the box walls are adjusted:
(25)
τP corresponds to a relaxation time that determines the coupling of the
modulated variable to the external bath. The pressure scaling can be isotropic,
which means that the factor is the same in all three spatial directions or – more
realistic – anisotropic, which means that the box dimensions change
independently during the simulation.[67]
The non-bonded interactions – van der Waals and Coulomb interactions –
require special attention because every atom interacts with every other atom in
the system. These interactions are expressed mathematically by the double
sum in Eq. (16). It is obvious that these sums require most of the computing
time in a simulation. As explained above, these interactions are in many cases
simply truncated at a certain cutoff distance (typically 8-10 Å). This approach is
only a minor problem for the van der Waals interactions, which are small at the
cutoff distance (because a [12-6]-potential is used to model these interactions).
The Coulomb interactions, however, are long-range interactions (two full
charges interact by more than 3 kcal mol-1 at a distance of 100 Å[74]), and using
a cutoff for these interactions can cause serious problems. A possible way to
avoid truncation of electrostatic interactions in periodic boundary conditions is to
apply the so-called particle-mesh Ewald (PME) method, which follows the
Ewald summation method of calculating the electrostatic energy for a number of
charges.[75] Originally devised by Ewald in 1921 to study the energetics of ionic
crystals,[76] it was used to simulate a crystal of bovine pancreatic trypsin inhibitor
[ ]1 3
01 ( )P
t P P tµτ
⎧ ⎫∆= − −⎨ ⎬⎩ ⎭
19
(BPTI) in an MD study by York and Darden in 1994.[77] They showed that the
use of PME dynamics produced results in much better agreement with X-ray
crystallography than MD simulations using a distance cutoff. In the case of
DNA, the use of PME electrostatics leads to more stable MD trajectories with
geometries closer to experimental data.[67, 78] An improved stability of the protein
structure is also observed in the case of the Tet-repressor if PME is used (see
section 3.2.2).
1.2.4 QM/MM Methods
Hybrid quantum mechanical/molecular mechanical (QM/MM) methods have
become a standard tool for describing large molecular systems, such as
enzymes (see [79-90] and the references therein). The fundamental concepts of
QM/MM techniques were given by Warshel and Levitt in 1976.[80] The basic idea
of all QM/MM approaches is to treat that part of the system that undergoes the
most important electronic changes during a reaction or upon binding a substrate
quantum mechanically and the rest of the system by molecular mechanics
(Figure 2).[80]
Figure 2. Partitioning of a molecular system into a QM- and an MM-part. The QM-part
is surrounded by the MM-part, and the latter can – depending on the implementation –
again be surrounded by a boundary region (BR).
20
The Hamiltonian describing the whole system consists of Hamiltonians for
the QM-and the MM-parts and of a Hamiltonian that describes the interaction
between the QM- and the MM-regions.[81]
(26)
Analogously, the total energy of the system consists of three terms:
(27)
In order to take the influence of the enzyme and the surrounding solvent on
the QM part into account, the electrostatic and the van der Waals interactions
between the QM- and the MM-part need to be calculated. To describe the
polarization of the QM wavefunction by the MM environment, the interactions
between the electrons/nuclei of the QM part with the MM point charges are
added to the QM Hamiltonian. If semiempirical methods are used, the integrals
describing the electrostatic interaction between the MM charges and the QM
electrons also need to be corrected.[81]
Various quantum mechanical methods, such as ab initio Hartree Fock,
density functional or semiempirical methods have been combined with force
fields like AMBER,[91] MM3,[92] CHARMM,[87, 88] GROMOS[93, 94] or the Tripos
Force Field.[83, 84, 89, 95-99] In the quantum mechanical part of the system, which is
usually polarized by the charges in the environment, electronic properties and
reaction mechanisms can be studied. The flexibility of the environment can be
simulated by an appropriate force field.
The system partitioning is easy if a complex between an enzyme and a non-
covalently bound substrate is to be computed. In this case, the ligand and the
active-site water molecules or metal ions can be chosen as the QM part and the
enzyme as the MM environment. If a substrate is bound covalently to the
enzyme or if the integration of some parts of the protein environment into the
QM region is desirable for other reasons, special techniques for the treatment of
the QM/MM border must be used (see also [80, 81, 83, 84, 89, 100] and the references
therein). One approach, which was proposed by Kollman and Singh in 1986[101]
and was taken up again by Bash and Karplus in 1990,[102] uses so-called link
/ˆ ˆ ˆ ˆ
QM MM QM MMH H H H= + +
/QM MM QM MME E E E= + +
21
atoms. In this concept, the free valences of the QM atoms bonded to MM atoms
are saturated by hydrogen atoms. These ‘dummies’ are considered explicitly in
the quantum mechanical calculation but are ‘invisible’ for the MM atoms.[81, 100]
Another approach, which allows bonds between the QM and the MM part and
does not use link atoms, is the local self-consistent field method (LSCF) by
Rivail and co-workers.[103] This method uses a combination of hybrid- and
atomic orbitals to represent the QM subsystem. The method takes up the idea
of Warshel and Levitt to use hybrid orbitals for the description of covalent bonds
across the QM/MM border. The LSCF-method is a development that arises from
the strictly localized molecular orbitals (SLMO)-approach by Rivail et al.[81, 100]
The different methods to treat the QM/MM border have been compared, but so
far none of them has proven to be better than the others.[4, 81]
In addition to biochemical questions, a wide range of applications of QM/MM
methods has been reported, e.g. solvation and reactivity of small molecules in
condensed phases or adsorption to zeolite surfaces.[104, 105]
Details regarding the QM/MM implementation in the semiempirical program
package VAMP are given in sections 2.1.3, 2.2, 3.1.7 and 3.2.3.
1.2.5 Molecular Docking
Molecular docking algorithms are intended to predict the structure of
intermolecular complexes formed between two or more molecules. In most
cases, docking techniques are used to determine the binding modes of protein
inhibitors bound to the active site of enzymes whose 3D-structure is known from
X-ray crystallography or homology modeling.[4, 100] Docking can be a powerful
tool for investigating unknown protein-inhibitor complexes or in the field of
computer aided drug design, where it is an important step in the in silico
screening of virtual libraries.[106-110] Many applications of docking techniques
have been reported (see e.g. [109-113] and the references therein). The abundant
supply of genomic data provides a rich basis for future computational studies,
e.g. protein-ligand docking studies using modeled protein structures as
targets.[114-116]
22
The simulation of the structural adaptation of an enzyme upon ligand binding
(induced fit)[117] requires a docking algorithm that treats both the ligand and the
protein environment as flexible molecules. Most commonly used docking
programs, however, do not consider the flexibility of the protein environment.
They simulate the binding of either a rigid[118] or a flexible ligand within a rigid
environment. Recently flexible ligand docking techniques were classified into
the following categories:[119, 120] fast shape matching (DOCK,[118, 121]
EUDOC),[122] incremental construction (FlexX,[123] Hammerhead),[124] Tabu
search (Pro_Leads,[125] SFDOCK),[126] genetic algorithms (Gold,[127] AutoDock
3.0,[128] Gambler),[129] evolutionary programming,[130] simulated annealing
(AutoDock 2.4),[131, 132] Monte Carlo simulations (MCDock,[133] QXP)[134] and
distance geometry (Dockit).[135] Several research groups have recently
developed algorithms for docking simulations accounting for a flexible protein
environment.[111, 127, 136]
As docking programs are intended to find the ‘correct’ binding mode of a
ligand within a very large conformational space, they must use efficient search
techniques in combination with fast scoring functions for evaluating the complex
geometries generated during the search.[4, 137] The combination of two or more
scoring functions (consensus scoring) has been suggested to perform better
than single scoring.[119, 129, 138]
More details regarding the docking algorithm used in this work are given in
sections 2.1.3 and 2.2.
23
2 QM/MM Docking: A New Protocol and its Evaluation for Known Test Systems
2.1 Introduction
2.1.1 Combining Docking Techniques with other Computational Methods
There are several reasons for combining docking techniques with other
computational methods: assessment of the quality of the scoring functions, re-
ranking the structures generated by docking, simulation of the structural
adaptations that occur in an enzyme upon ligand binding, a more detailed
description of the binding mode of the ligand and, in the case of QM methods,
information about reaction mechanisms and electronic properties. Several
studies using a combination of a rigid-environment docking algorithm with
subsequent molecular dynamics simulations and/or force field geometry
optimizations of the docked protein-inhibitor complexes have been
published.[139-143] In these studies, different degrees of protein flexibility were
used at different stages of the simulation protocol. In some cases, side chain
flexibility was simulated, while other studies used harmonic constraints or
allowed flexibility only for some residues. Some even used totally rigid
proteins.[144, 145] Semiempirical quantum mechanical/molecular mechanical
(QM/MM) studies of protein-bound ligands within rigid protein environments
have also been reported.[79, 100] Some excellent reviews concerning the
dynamics of molecular recognition have been published.[146-148]
In this work, a combined QM/MM docking approach was used to take the
flexibility of the protein environment into account.[95-99] AutoDock, which uses a
rigid protein environment, first explores the possible binding sites and the
conformational space of the ligand. The subsequent QM/MM geometry
optimization of the complex obtained by docking can refine the ligand-protein
interactions and gives a very good description of the electronic properties of the
region treated quantum mechanically. The use of a flexible environment allows
the simulation of the structural adaptation of an enzyme upon ligand binding.
24
2.1.2 The Protein-Inhibitor Systems Investigated
The QM/MM docking approach was validated using a set of structurally
known protein-inhibitor complexes of the enzymes HIV-1-protease (PDB entries
1hvr,[149] 1upj,[150] 7upj),[151] β-trypsin (3ptb),[152] thrombin (1bcu,[153] 1c5n,[154]
1bhx),[155] carboxypeptidase A (3cpa)[156] and dihydrofolate reductase (4dfr,[157]
1rg7).[158] The ligands investigated are shown in Figure 3. The crystallographic
data of the complexes were taken from the Protein Data Bank.[159-161] In order to
examine whether the QM/MM docking method also achieves good results if the
ligand was not crystallized with the enzyme and for investigating the induced fit,
compounds 1 and 6 were docked into the ligand-free HIV-1-protease and
thrombin structures, respectively (PDB entries 3phv[162] and 1c5l).[154]
Figure 3. Ligands 1-9 were used for the validation of the QM/MM docking approach.
For the investigation of the induced fit, ligands 1 and 6 were docked into ligand-free
crystallized structures of HIV-1-protease and thrombin, respectively. PDB codes of the
protein structures used are given in parentheses.
25
2.1.3 Overview of the Method
First, the binding position and the binding mode of a flexible ligand within a
rigid protein environment was explored. The structures obtained from this
docking with AutoDock were subsequently optimized using a combined QM/MM
approach. In contrast to the docking simulations, where protein flexibility cannot
be taken into account, both the ligand and the protein were treated as flexible
molecules in the QM/MM optimizations.
In the docking step of the protocol, AutoDock 3.0 was used.[128] The program
allows the docking of a flexible ligand within a rigid protein environment. Each of
the search algorithms implemented in AutoDock 3.0 was tested: the Lamarckian
genetic algorithm (LGA), which combines aspects of a global and a local search
procedure in a way similar to Lamarckian evolution, the ‘normal’ genetic
algorithm (GA) and simulated annealing (SA). The structures generated by the
docking algorithm were clustered and ranked according to their docked
energies. After a visual inspection of the docking results, the lowest energy
structures of the 3 ‘best’ clusters were selected for the QM/MM geometry
optimizations.
A ligand bound to an enzyme is polarized by the protein environment.
Although some pure molecular mechanics approaches consider polarization
effects,[4] a quantum mechanical optimization of the docked protein-inhibitor
structures gives a more accurate description of the electronic properties of the
ligand. As quantum mechanical calculations on whole enzymes are
computationally very demanding, a QM/MM approach for the optimization of the
protein-inhibitor structures obtained by docking was chosen. The QM/MM
optimizations were performed using the semiempirical program package VAMP
7.5.[163]
The original QM/MM implementation in VAMP used a rigid point charge
environment with associated van der Waals radii.[104] As the mechanical
response of the enzyme upon ligand binding was one of the main points of
interest, a flexible protein environment was used in this study. Our method for
hybrid QM/MM calculations has been described in detail elsewhere,[83, 84, 89, 104,
164] so only a rough overview is given here.
26
As the protein-inhibitor complexes generated by AutoDock contain non-
covalently bound ligands, the simplest strategy for partitioning the system into
QM and MM parts is to treat the ligand quantum mechanically and the protein
environment classically. In some cases, however, the QM/MM border crosses
covalent bonds, e.g. if active-site amino acids are included in the QM part. In
these cases, a modified link atom approach was used (see computational
details section).[84]
The QM and the MM parts are optimized alternately. The polarization of the
QM wavefunction by the point charges in the MM environment is taken into
account via an additional one-electron term in the Fock matrix.[104] The MM part
uses either the Cornell et al. 1995 AMBER[66] or the Tripos 5.2[165] force field.
Three major contributions are considered for the calculation of the QM/MM
interactions: van der Waals and Coulomb interactions and the polarization of
the QM part by the MM environment. The van der Waals parameters describing
the QM/MM interaction are taken from the Universal Force Field (UFF) by
Rappé et al.[166] No back-polarization of the MM part by the QM part is included.
In order to compare the results of the QM/MM optimizations with those
obtained from a pure MM optimization of the docked structures, also molecular
mechanics geometry optimizations of the docked structures of ligand-free HIV-
1-protease/ligand 1 (3phv) were also performed using the AMBER 5 program
suite[167] and Sybyl 6.7[168] (see computational details). This system was also
used to examine the influence of the force field implementation in the
semiempirical program package VAMP, particularly that of the low memory
Broyden-Fletcher-Goldfarb-Shanno (BFGS) optimizer[169, 170] used for the MM
part, by performing pure MM optimizations using VAMP.
2.2 Computational Details
The protein structures used in this study were retrieved from the Protein Data
Bank.[159-161] For the docking step, all water molecules in the X-ray structure of
the protein except those in the active site were deleted. Essential hydrogens,
i.e. hydrogens bound to electronegative elements like nitrogen and oxygen,
27
were added and Kollman united atom charges[64] were assigned using Sybyl
6.7.[168] The N- and C-termini of the protein peptide chains were used in their
zwitterionic forms. Active site water molecules were included in the environment
in the docking procedure. These water molecules were orientated in a
chemically reasonable manner in order to optimize the active site hydrogen-
bonding network. For the investigation of the role of active site water molecules,
additional docking runs were performed using environments without any water
molecules. The ligand was saturated with hydrogen atoms and the Coulson
partial atomic charges of the ligand atoms were taken from AM1[36, 171]
calculations. The AutoDock suite of tools (rem-lp, mol2topdbq, check-qs,
addsol, AutoTors, mkgpf3 and mkdpf3) was used to process the ligand and the
environment as described in the AutoDock manual.[111]
As AutoDock treats hydrogen bonds to pyramidal nitrogens with a [12-10]
Lennard-Jones type potential function and planar nitrogens in delocalized
systems with a [12-6] potential energy function, these two types of nitrogen had
to be distinguished. In cases where the geometry of a nitrogen atom is
presumably only slightly pyramidalized, e.g. in amino-substituted aromatic
systems,[172] both sp2 and sp3 hybridization of the nitrogen atom was tested. It
was also differentiated between aromatic and aliphatic carbon atoms for the
correct calculation of the desolvation term in the scoring function.[128]
The grid maps for docking were calculated using AutoGrid 3.0. If the binding
site of the ligand was known from X-ray crystallographic data, the grids were
centered on the ligand heavy atom closest to the ligand’s center of mass in the
X-ray structure and grid dimensions of 603 points were used. For structurally
unknown protein-inhibitor systems, maps of 1203 grid points were calculated,
centered on the center of mass of the protein. The standard grid point distance
of 0.3750 Å was used, but also grid point distances of 0.1875 and 0.5000 Å
were checked.
The standard AutoDock Lennard-Jones parameters are based on the Weiner
et al. 1984, 1986 AMBER force field,[64, 65] which does not contain some of the
elements (Zn, Br, I) present in the systems investigated in this study. In these
cases, the standard AutoDock Lennard Jones parameters were supplemented
with parameters for bromine and iodine taken from a combined set of AMBER
28
1986 and Merck Force Field (MMFF) parameters.[111, 173] The zinc non-bonded
parameters were taken from the literature (R* = 1.1 Å, ε = 0.0125 kcal mol-1).[174,
175] A zinc charge of +0.701, determined from a semiempirical QM/MM
optimization of the X-ray structure of the complex, was used in the docking runs
of the carboxypeptidase A complex ligand 8/3cpa.
Each of the docking algorithms provided by AutoDock, the Lamarckian
genetic algorithm, the genetic algorithm and simulated annealing, was tested
intensively. Each docking experiment consisted of 50 runs of docking, or, in
order to check the reproducibility of the experiment, of two independent
dockings of 20 runs each. A cluster analysis was performed in order to rank the
structures generated by the docking experiment. In the cluster analysis, a user-
defined RMSD tolerance is used to determine whether two docked ligand
conformations are similar enough to be placed in the same cluster. The clusters
are ranked in order of increasing energy, by the lowest energy in each
cluster.[128] Depending on the size and the flexibility of the ligand, different
RMSD values for the cluster analysis were chosen: in most cases a value of 0.5
Å leads to well-populated clusters, each representing a characteristic binding
position of the ligand (1/1hvr, 2/1upj, 4/3ptb, 5/1bcu, 6/1c5n, 6/1c5l). Larger
ligands or difficult docking problems require an RMSD cluster tolerance of 1.0 Å
(8/3cpa, 3/7upj) or even 2.0 Å (7/1bhx, 9/4dfr, 9/1rg7, 1/3phv). Every docking
experiment was run twice, using different maximum changes of the variables
describing the ligand-position. Besides the default maximum step sizes for the
translation of the ligand (2.0 Å), its orientation in space (50°) and its internal
torsions (50°) values of 0.2 Å and 5°, respectively, were also checked, as
recommended in the literature.[128]
For the docking runs with the Lamarckian genetic algorithm (LGA) and the
classical genetic algorithm (GA), a population of 50 individuals was used. The
maximum number of energy evaluations was set to 1.5 million, the maximum
number of generations was 27,000. The best individual of each simulation
survived and was passed to the next generation without any genetic changes
(elitism value 1). For the dockings with the Lamarckian genetic algorithm 6% of
the population were refined into the nearest local minimum at each generation.
For this local search the pseudo-Solis and Wets method implemented in
29
AutoDock 3.0 was used.[128] The LGA parameters were chosen as described in
the literature.[128]
In the simulated annealing experiments, the starting geometry and starting
position of the ligand were set by the random-number generator. For a
comparison, docking runs were also performed, in which the starting position of
the ligand was outside the binding pocket. The starting temperature for the
annealing was chosen to be RT = 616 cal mol-1 (T = 310 K). Values of 61,600
cal mol-1 (T = 31,000 K), 1,200 cal mol-1 (T = 604 K) and 300 cal mol-1 (T = 151
K) were also checked. The number of cooling cycles of each docking run was
50. A maximum of 25,000 accepted or rejected steps per cycle was allowed.
The van der Waals parameters for the calculation of the internal energy of
the ligand were chosen analogously to the ligand-protein interaction
parameters. For all other parameters, which are not listed here, the default
values given by the AutoDock tools mkgpf3 and mkdpf3 were used. As
AutoDock cannot consider the flexibility of ligand ring systems, different ligand-
ring conformations were generated where appropriate before the docking runs
and the different conformations were treated as different molecules. The grid-
map calculations and docking runs were performed on an SGI Power Challenge
(R10000).
The docking results were inspected visually and the best structures in the
three lowest energy clusters were then selected for the QM/MM part of the
protocol. In the case of the complex formed between ligand 2 and HIV-1-
protease (PDB code 1upj), all docked structures were selected for the QM/MM
optimizations to allow a detailed comparison of the results of the docking and
the QM/MM steps.
For the preparation of the QM and the MM parts used in the QM/MM
optimizations, again Sybyl 6.7 was used. Hydrogens were added to both the
ligand and the protein, and the protonation states of the ligand and the acidic
and basic protein residues were checked thoroughly. In order to avoid strong
monopole-monopole interactions between QM and MM parts, care was taken
that the overall charge of the MM part was zero. This was achieved by
modifying the protonation states of histidines, lysines, aspartates or glutamates
30
located far away from the active site at the bulk orientated surface of the
protein. Gasteiger-Marsili charges[176-178] were calculated for all protein atoms
using the Sybyl program. In order to investigate whether active-site water
molecules should better be included in the QM part or in the environment, both
possibilities were checked systematically.
In all cases where the protein backbone was cut for QM/MM partitioning, an
analogous approach as described in the literature was chosen.[84, 89, 96] An
overview of the cutting scheme is given in Figure 4.
Figure 4. Cutting the protein backbone for QM/MM input generation. 10, 20: force
constants [kcal mol-1 Å-1] of the harmonic constraints fixating the positions of carbonyl-
carbons and amide-nitrogens at the edge of the QM part. The hydrogen atoms added
in order to saturate the MM part have a reduced van-der-Waals radius of 0.1 Å. The
Cα-atoms of R’ and R’’ are deleted.[84, 89, 96]
31
The MM-QM optimization sequence was iterated 20 times. Each iteration
consisted of 250 MM cycles, followed by an optimization of the QM part to a
gradient norm of 0.4 kcal mol-1 Å-1. For the QM optimizations the AM1
Hamiltonian[36, 171] and the Eigenvector Following optimizer were used.[179, 180]
The MM part was optimized using the Tripos 5.2 force field[165] and a low
memory BFGS algorithm.[169, 170] A constant dielectric constant of 4.0 was used
for both the MM- and the QM/MM-Coulomb interactions. This has proven to be
appropriate in studies by Schürer et al.[83, 84, 89] Missing metal parameters were
taken from the Sybyl program.[168] Although these parameters have not been
validated, they are adequate for merely structural metal ions.[84, 89] Although it is
also possible to use these unvalidated parameters for active site metal ions,
these atoms can be included in the QM part for an exact description.
In the QM/MM optimizations of carboxypeptidase A (3cpa)/glycyltyrosine (8)
both calculations with the zinc ion included in the MM part and with the zinc as
part of the QM region were performed. For the MM description of the zinc ion
the QM/MM derived zinc charge of +0.701 was used (see above). For a
quantum mechanical analysis of the zinc-protein and zinc-ligand interactions,
the ligand, the zinc ion and its first coordination sphere (HIS69, GLU72,
HIS196) as well as the amino acids ARG145 and GLU270 were incorporated
into the QM part.
The AMBER geometry optimizations were performed with the AMBER 5
program package[167] and the Cornell et al. 1995 force field.[66] A maximum of
10,000 optimization cycles was allowed, after 100 initial steepest-descent
cycles the optimizer switched to conjugate gradients. Tripos 5.2 force field
optimizations with Sybyl 6.7[168] used 20 cycles of simplex minimization and
then switched to conjugate gradients, also allowing a maximum of 10,000
cycles. In both the AMBER and the Sybyl optimizations a distance-dependent
dielectric of 4 ε and a cutoff of 12 Å were used. The gradient convergence
criterion was 0.05 kcal mol-1 Å-1. For the pure MM optimizations with VAMP
8.0[181] the same parameters as for the MM part in the QM/MM procedure were
chosen (see above). Again a maximum of 10,000 optimization cycles was
allowed.
32
2.3 Results and Discussion
2.3.1 Validation of the QM/MM Docking Approach
Protein structures of crystallized protein-inhibitor complexes were used for
the docking runs designed to validate the QM/MM docking approach. A ‘good’
docked ligand structure has a low docked energy, shows a small RMS deviation
with respect to the X-ray structure, and is found in a well-populated cluster. In
terms of these criteria the best docking results were achieved using the
Lamarckian genetic algorithm implemented in AutoDock. The maximum step
size did not affect the quality of the results significantly. The best results were
obtained with the default grid point distance of 0.375 Å.
The QM/MM optimizations confirmed the binding positions and binding
modes obtained by AutoDock. The structures obtained from docking and from
the QM/MM optimizations were superimposed on the corresponding X-ray
structures using their α-carbons. As a measure for the agreement of the
docking results with the crystallographic data, the RMS deviations between the
docked ligand and the ligand in the X-ray structure were calculated for all heavy
atoms of the ligand (Table 1). Figure 5 shows the results of the QM/MM docking
protocol for the complex formed between carboxypeptidase A (3cpa) and ligand
8.
33
Table 1. RMS deviations in Å calculated for ligand heavy atoms after docking and
QM/MM optimization, respectively, with reference to the X-ray structure of the
crystallized complex.
enzyme
ligand
PDB-
code
RMSD ligand
AutoD./X-ray
RMSD ligand
VAMP/X-ray
crystallographic
resolution in Å
HIV-1-protease 1 1hvr 0.69 0.80 1.8
2 1upj 0.83 0.99 2.2
3 7upj 0.63 0.99 2.0
β-trypsin 4 3ptb 0.34 0.58 1.7
thrombin 5 1bcu 0.89 0.43 2.0
6 1c5n 0.87 0.74 1.5
7 1bhx 0.85 0.87 2.3
carboxypeptidase A 8 3cpa 1.10 0.51 2.0
dihydrofolate reductase 9 1rg7 1.14 1.43 2.0
9 4dfr 1.27 1.25 1.7
In the AutoDock runs, the active-site water molecules were treated as part of
the protein environment. In order to describe the positions of active-site water
molecules and the active-site hydrogen bonding network adequately in the
QM/MM calculations, the ligand, the active site water molecules and the protein
residues that participate in the active-site hydrogen bonding network must be
included in the quantum mechanical region. In test calculations, the water
molecules were also treated both as a part of the QM and the flexible MM
regions. Both procedures were able to reproduce the crystallographic positions
of the active-site water molecules approximately. Water molecules included in
the QM part show clear hydrogen bonds to suitable groups of the ligand. As
atom-centered point charges were used in the MM part, it was not possible to
simulate a directed interaction mediated by lone pairs between the MM
environment and QM water molecules. Therefore, no clear orientation of the
water molecules in the QM part with respect to suitable groups of the MM
environment is found. If the active site water molecules are included in the MM
34
part, hydrogen bonding can be observed neither to the QM nor the MM parts.
Consequently, the treatment of active-site water molecules in the QM part is the
best approach.
Figure 5. The QM region used in the QM/MM optimization of the complex formed
between carboxypeptidase A (3cpa) and ligand 8. CPK sticks: QM/MM docking result.
Green sticks: heavy atoms of the residues included in the QM part in their X-ray
conformation. Only the backbone atoms of residues 70 and 71 are included in the QM
part. Graphics were generated with ViewerLite 4.2.[182]
2.3.2 Simulating the Induced Fit
The systems HIV-1-protease/ligand 1 and thrombin/ligand 6 were used to
test how well the QM/MM docking protocol could reproduce the induced fit. The
calculations were run twice, once using the protein structures taken from the
crystallized protein-inhibitor complexes (1hvr, 1c5n) and again with ligand-free
protein structures (3phv, 1c5l) as receptors.
35
In contrast to the structure of the crystallized complex (1hvr), the ligand-free
crystallized structure of HIV-1-protease (3phv) exhibits a large cavity at the
active site. The so-called flaps, formed by the loop residues 42-58 and 42′-58′
in the other monomer of the homodimer,[183] are in an open conformation. These
are important for inhibitor binding and possibly also to exclude water from the
catalytic site.[183] Nevertheless, AutoDock successfully predicted the binding of
the ligand to the two ASP25/ASP25′ residues (Figure 6a, Figure 7a). Despite
the large movements required, the QM/MM optimization was able to close the
flaps of the ligand-free structure containing the docked ligand, so that a
geometry very similar to that of the crystallized complex was obtained (Figure
6b, Figure 7b). As a measure for the structural changes between two
conformations, the RMS deviations for the superimposition of the protein α-
carbons were calculated.
Thus, the Cα atoms of the non-flap residues (1-41, 59-99 and 1′-41′, 59-99′)
of the QM/MM optimized complex of 3phv/1 and the X-ray structure of the
ligand-free crystallized enzyme (3phv) were superimposed. The small RMS
deviation of 1.21 Å corresponds to the small changes in the non-flap part of the
protein. After this superimposition the RMSD for the Cα atoms in the flaps
(residues 42-58 and 42′-58′) was calculated, the value thus obtained is 2.74 Å.
The highly mobile loop residues 49-51 and 49′-51′ show an RMS deviation of
4.33 Å. This relatively high RMSD corresponds to the large movement of the
flaps in the QM/MM optimization. A superimposition of all Cα atoms of the
QM/MM optimized complex of 3phv/1 and the X-ray structure of the crystallized
enzyme-inhibitor complex (1hvr) gives an RMSD value of 1.69 Å.
36
Figure 6. Ligand 1 docked into ligand-free crystallized HIV-1-protease (PDB-code
3phv). a) Left: AutoDock result. b) Right: Result of the VAMP QM/MM optimization.
Graphics were generated with ViewerLite 4.2.[182]
Figure 7. Ligand 1 docked into ligand-free crystallized HIV-1-protease (PDB-code
3phv). a) Left: AutoDock result. b) Right: Result of the VAMP QM/MM optimization. For
the protein the solvent accessible surface was generated using a probe radius of 1.4 Å.
The electrostatic potential based on Gasteiger charges is mapped on the surface (a
color representation is given in the summary at the beginning of this thesis). Graphics
were generated with ViewerLite 4.2.[182]
An analogous procedure applied to thrombin also yields the correct binding
mode of ligand 6 within the active site of the enzyme. No significant structural
changes occur in the QM/MM optimization of the complex of the ligand-free
crystallized enzyme (1c5l) and ligand 6, which leads to the relatively small Cα
RMSD value of 1.14 Å. The structures obtained by QM/MM docking using the
37
ligand-free protein structure (1c5l) and the X-ray structure of the crystallized
complex (1c5n) are almost identical, corresponding to a Cα RMSD of 1.11 Å.
The calculations show that the structural adaptations of thrombin upon
binding the ligand are much smaller than in the case of HIV-1-protease. This is
in good agreement with the high substrate specifity of thrombin: as there are no
significant structural changes in thrombin upon ligand binding, the enzyme can
only bind a ligand that fits well into the active site.
In order to investigate whether the results of the QM/MM optimization of
ligand 1 docked into ligand-free crystallized HIV-1-protease (3phv) could be
reproduced using a pure MM optimization, geometry optimizations of the
protein-inhibitor complex obtained from docking were performed using the
AMBER 5 program package and Sybyl 6.7. In contrast to the VAMP QM/MM
optimizations, both the AMBER and the Tripos force field optimizations were not
able to close the flaps of HIV-1-protease. Again, the Cα atoms in the non-flap
regions (residues 1-41, 59-99 and 1′-41′, 59′-99′) of the MM optimized
complexes of 3phv/1 and the X-ray structure of the ligand-free crystallized
enzyme (3phv) were superimposed. The RMSD values thus obtained are 0.53
Å (AMBER) and 0.86 Å (Sybyl). Subsequently the RMSDs for the Cα atoms in
the flaps (residues 42-58 and 42′-58′) were calculated to be 0.72 Å for the
AMBER geometries and 1.18 Å for the Sybyl optimized geometries. The RMS
deviations calculated for the loop residues 49-51 and 49′-51′ are also small,
0.71 Å (AMBER) and 1.95 Å (Sybyl). These relatively small RMSDs indicate
negligible movements of the flaps in the MM optimizations. A superimposition of
all Cα atoms of the MM optimized complexes of 3phv/1 and the X-ray structure
of the crystallized enzyme-inhibitor complex (1hvr) gives RMSD values of 1.94
Å (AMBER) and 1.90 Å (Sybyl).
Pure MM optimization of the docked complex of 1 and 3phv with VAMP 8.0
closes the flaps almost as well as the QM/MM optimization. The Cα
superimposition of the non-flap residues (residues 1-41, 59-99 and 1′-41′, 59′-
99′) of the VAMP MM optimized complex of 3phv/1 with the X-ray structure of
the ligand-free crystallized enzyme (3phv) gives an RMSD of 1.17 Å. The Cα
RMS deviation calculated for the flap residues (42-58 and 42′-58′) is 2.75 Å.
38
The loop residues 49-51 and 49′-51′ have an RMS deviation of 3.95 Å, 0.38 Å
smaller than the value obtained by the QM/MM optimization, thus indicating that
the QM/MM optimization closes the flaps slightly better than the pure MM
optimization. A direct comparison between the VAMP QM/MM and pure MM
results can be made by superimposing all α-carbons of these structures, giving
an RMSD value of 0.97 Å. The superimposition of the VAMP MM optimized
complex of 3phv/1 and the X-ray structure of the crystallized enzyme-inhibitor
complex (1hvr) gives an RMS deviation of 1.71 Å. This value is very similar to
the corresponding value obtained using the QM/MM method. These results
show that, unlike the AMBER and Sybyl MM optimizations, the pure VAMP MM
optimization closes the flaps almost as well as the VAMP QM/MM optimization.
The distinct structural adaptations in the VAMP QM/MM case thus result
predominantly from using a different optimizer than in the AMBER and Sybyl
MM optimizations.
In order to assess the contribution of QM/MM polarization of the ligand, which
results in a stronger receptor-ligand interaction and also helps to close the flaps,
the Coulomb interaction energy between ligand and environment was
calculated for the QM/MM optimization with and without ligand polarization and
for the VAMP MM optimization of 3phv/1. QM polarization gives the strongest
QM/MM Coulomb interaction energy (-9.34 kcal mol-1) for the standard VAMP
QM/MM optimization with QM polarization. For comparison, the QM/MM
Coulomb energy was calculated without polarization of the QM part by using
gas phase single point charges for the QM part. The value thus obtained (-6.39
kcal mol-1) is significantly smaller than in the case with QM polarization,
emphasizing the importance of QM polarization for the correct description of
ligand-receptor interaction. For the pure VAMP MM optimization of 3phv/1, a
ligand-receptor Coulomb interaction energy of only –2.76 kcal mol-1 is found.
These results indicate that, in addition to the optimizer effects discussed above,
the enhanced Coulomb receptor-ligand interaction in the QM/MM optimizations
also helps to close the flaps by increasing the forces acting on the MM-treated
protein.
39
2.3.3 Polarization of the Ligand by the Environment
When a ligand is bound to the active site of an enzyme, it is polarized by the
protein environment. Unlike most of the commonly available force fields, a QM
description of the ligand within a classical protein environment can describe this
polarization effect, but not the reverse polarization of the protein by the ligand.
For the HIV-1-protease/ligand 1 system (3phv) the additional dipole moment
induced in the ligand results in a stronger Coulomb interaction (see above)
between the hydroxylic groups of the ligand and the ASP25/ASP25′
carboxylates and between the ligand’s carbonyl oxygen and the backbone
amide hydrogens of ILE50 and ILE50′. This enhanced Coulomb interaction
helps to close the flaps.
For all ligands discussed in this work, the additional dipole moment
contributions induced by the environment lead to a larger ligand dipole moment
in the presence of the environment (Table 2).
The effect of the polarization of the ligand by the protein environment can be
visualized by superimposing the ligand structure obtained from the QM/MM
optimization of the ligand/protein complex with the corresponding structure
calculated as a single point in the gas phase.[83] The electron densities of the
superimposed structures, represented as natural atomic orbital point charges
(NAO-PCs),[184] were subtracted from each other and the resulting difference
electrostatic potential was mapped on the common van der Waals surface of
the molecules (Figure 8).
40
Table 2. Ligand dipole moments in [D] with and without protein environment. For
charged molecules, the dipole moment was calculated using the center of mass as
origin. The total ligand dipole moments with environment were derived from the
QM/MM calculations of the protein-inhibitor complexes. The permanent ligand dipole
moments without environment were calculated from gas-phase single point calculations
on the protein-bound ligand conformation.
enzyme
ligand
PDB-code
total dipole moment with environment
[D]
permanent dipole moment
without environment [D]
HIV-1-protease 1 1hvr 4.12 2.40
1 3phv 4.66 2.93
2 1upj 7.70 6.18
3 7upj 11.78 10.29
β-trypsin 4 3ptb 6.68 6.07
thrombin 5 1bcu 2.75 1.87
6 1c5n 15.15 14.79
6 1c5l 14.99 14.58
7 1bhx 27.48 26.49
carboxypeptidase A 8 3cpa 11.24 10.18
dihydrofolate reductase 9 1rg7 67.34 64.53
9 4dfr 75.66 73.86
41
Figure 8. Superimposition the QM/MM optimized structure of ligand 1 within the active
site of HIV-1-protease (3phv) and the corresponding structure calculated in the gas
phase. The induced electrostatic difference potential was mapped on the common van
der Waals surface of the molecules. For clarity, only the extrema are shown by
excluding the range between –1 and +7 kcal mol-1. The resulting dipole moment vector
(1.64 D) is represented by the arrow. Graphics were generated with the molecule
visualizing and building software TRAMP 1.1b.[185]
Ligand 1 optimized with AM1[36, 171] in the gas phase shows a permanent
dipole moment (2.4 D) parallel to the C2 axis of the molecule. The direction of
the dipole moment vector is ideal for the interaction with the negatively charged
ASP25/ASP25′ carboxylates in the active site. When ligand 1 is bound to the
active site of HIV-1-protease (1hvr), an additional dipole moment (1.67 D) in
almost the same direction as the permanent one is induced: the positive
polarization of the ligand near its hydroxylic groups is enhanced by the
negatively charged carboxylic groups of ASP25/ASP25′. The negative
polarization of the ligand near its carbonyl oxygen, which is hydrogen-bonded to
the amide hydrogens of ILE50/ILE50′, is also enhanced. The additional dipole
42
moment induced by the protein environment is 1.67 D, giving an overall dipole
moment of ligand 1 within the active site of the enzyme of 4.12 D.
2.4 Conclusions and Outlook
The results described above suggest that the combination of a fast docking
technique, in this case AutoDock, with the subsequent QM/MM optimization of
the docked structures gives promising results for the systems investigated.
Although the flap-closing observed for HIV-1-protease turned out to be primarily
an effect of the optimizer used in the QM/MM program, polarization of the ligand
plays an important role in helping to promote such induced fits. The QM/MM
approach used here only allows polarization in one direction (i.e. the ligand is
polarized by the environment) but this should have exactly the same type of
effect as the polarization of the receptor by the ligand. Adjustment of the (in any
case fictitious) dielectric constant used to compensate for the fact that a point
charge model is used for the receptor can therefore simulate the effect of the
missing polarizability of the receptor.
For virtual screening, the protocol can be automated by using scripting
languages (e.g. Perl[186] or Python[187]). The computational costs, which are
moderate for the docking step (several minutes for 50 docking runs on an AMD
Athlon 900 MHz CPU, depending on the molecular system) and more
demanding for the QM optimizations (several hours on an AMD Athlon, system-
dependent), can be lowered by using fewer docking runs and applying less
restrictive convergence criteria in the QM optimization step. This is justifiable for
screening purposes.
In this study, the QM/MM binding energies were not related to the observed
binding constants, although a linear response technique based on components
of our QM/MM docking energy has been used successfully[188] and could also
be applied to the ligand/receptor systems reported above. Other promising HIV-
1 protease data sets, which could be used for a correlation of computed
QM/MM binding energies with Ki values are e.g. given in [143, 189, 190].
43
However, two problems require particular attention: Firstly, future studies
require a systematic protocol for predicting the presence or absence, the
positions and the hydrogen-bonding pattern of active-site water molecules. This
is important for many ligands that bind together with waters, such as in the HIV-
1-protease complex 7hvp.[191] Several approaches for the prediction of active
site water molecules have been published.[192-194]
Secondly, the static picture presented here cannot deal with induced fits that
are accompanied by a significant energy barrier for the protein movement. This
can be achieved by QM simulated annealing of the docked structures rather
than simple optimization, for example. Entropic binding effects will in any case
require extensive sampling of molecular dynamics trajectories or Monte Carlo
calculations or even free energy perturbation calculations. These can be
achieved relatively easily with a quantum mechanical ligand (plus waters) within
a classical environment.
44
3 Simulating FRET: A Combined Molecular Dynamics-QM/MM Approach
3.1 Introduction
3.1.1 The Tetracycline Class of Antibiotics
Tetracycline (Tc, 10, Figure 9) and its derivatives form a family of widely used
broad-spectrum antibiotics that show bacteriostatic activity against Gram-
positive and Gram-negative bacteria and other infective pathogens like
ricketsias, mycoplasms, viruses and protozoans. In addition to their medicinal
applications, tetracyclines are still used as nutritive agents in livestock farming
in some countries. In 1948, 7-chlorotetracycline (Aureomycin, 11) was isolated
from cultures of Streptomyces aureofaciens, in 1950 oxytetracycline
(Terramycin, 12) was described, and in 1953 the chemical structures of both
compounds were determined by Woodward, Conover and their colleagues.
Tetracycline was produced from Aureomycin by catalytic hydrogenation in
1953.[195, 196]
Tetracyclines are inhibitors of protein biosynthesis. The antibiotics bind to the
30S subunit of the bacterial 70S ribosome and inhibit the binding of aminoacyl-
tRNA to the acceptor A-site of the ribosome.[195, 196]
Under physiological conditions, tetracycline forms calcium complexes, which
are incorporated in bone and tooth tissue.[195, 196] Most interestingly, tetracycline
was isolated in human bones from several populations in ancient Nubia and
Egypt (350 AD – 550 AD). Tetracycline extracted from these bones retained its
antibiotic potential after 14 centuries. The proposed explanation for the
tetracycline enrichment in these ancient probes refers to the bread-baking and
beer-brewing processes used in these communities, which provide good
growing conditions for streptomycetes (bread-baking and brewing were closely
linked).[197]
45
Figure 9. Tetracycline (10), 7-chlorotetracycline (Aureomycin, 11), oxytetracycline
(Terramycin, 12) and doxycycline (13).[195, 196]
Tetracyclines are pharmacologically and chemically “promiscuous” agents. In
addition to their antibiotic properties, they also have agonist, antagonist and
physiological properties against a wide range of proteins, enzymes and
macromolecular targets in eukaryotic cells.[195] Tetracycline and its derivatives
have recently been proposed as a new class of antagonists in prion diseases as
they prevent the aggregation of prion protein peptides and their acquisition of
protease resistance in vitro and in vivo.[198-200]
The structural chemistry of tetracyclines in solution is complicated by several
factors. Most tetracyclines contain three functional groups that are subject to
protonation-deprotonation equilibria. Additionally, the molecules can exist in
several conformations, which again depend on the protonation state.[201-204]
Tetracyclines form complexes with divalent metal cations like Mg2+ or Ca2+. The
structure of these complexes depends on the protonation state and the
conformation of the drug. Several experimental and theoretical studies have
tried to characterize the conformations and tautomers of tetracyclines and their
complexes in solution.[201-205]
46
Othersen, Lanig and Clark used a systematic surface-scan technique using
semiempirical calculations to identify metal-binding sites of tetracycline.[205]
Furthermore, we have used density functional theory (DFT) to investigate the
conformations and tautomeric forms of neutral tetracycline and
anhydrotetracycline in aqueous solution.[202, 203] The results suggest that there
are two main conformations of tetracycline and anhydrotetracycline (‘extended’
and ‘twisted’, Figure 10). In the case of tetracycline in solution, the extended
conformation is 3-3.5 kcal mol-1 more stable than the twisted one for equivalent
tautomers. As many as six different tautomeric forms lie within 10 kcal mol-1 of
the most stable one (Figure 11). The energetic preference for the extended
conformation is due to a solvent effect.[202] Anhydrotetracycline in solution
adopts the twisted conformation in the unionized form and the extended one as
a zwitterion. In contrast to tetracycline, anhydrotetracycline is predicted to favor
zwitterionic structures in solution quite clearly.[203]
Figure 10. B3LYP/6-31G(d) optimized structures for the two ring conformations of
completely neutral tetracycline.[202]
extended conformation twisted conformation
47
Figure 11. Our best estimate of the relative stabilities of neutral tetracycline tautomers
in aqueous solution at 298 K.[202]
3.1.2 The Tet-Repressor/Operator (TetR/tetO) System
The Tet repressor/operator (TetR/tetO) system is a regulatory switch in the
most important resistance mechanism of Gram-negative bacteria against the
tetracycline class of antibiotics.
When tetracycline enters the bacterial cell, it forms a chelate complex with
divalent metal ions, mostly Mg2+. Under physiological conditions, the metal
cation is chelated by the deprotonated 1,3-keto-enol group O11/O12- of the
zwitteranionic tetracycline (Figure 12).[202, 205-208]
48
Figure 12. Structure of the physiologically active tetracycline-magnesium complex.[206-
208]
The expression of the protein predominantly responsible for the resistance,
the tetracycline antiporter (TetA), is under tight transcriptional control of TetR.
TetA is an intrinsic membrane transport protein that acts as antiporter across
the cell membrane and couples the efflux of the Tc/M2+ complex with the uptake
of H+.
In the absence of Tc, the homodimeric Tet repressor TetR binds specifically
to two operator sequences of the DNA (tetO) and thus prevents the expression
of the genes tetA and tetR, i.e. the genes coding itself and the TetA protein.
Two molecules of the [TcMg]+-complex bind to the TetR binding sites. This
induces a cascade of allosteric conformational changes in the TetR protein,
which result in an increased center-to-center separation of the DNA-recognition
helices (from 36.6 Å in the operator complex to 39.6 Å in after inducer-
binding).[206-208] These conformational changes lead to the dissociation of the
TetR/DNA complex (Figure 13). Thus, the genes tetA and tetR can be
transcribed, TetA and TetR are synthesized and TetA transports Tc out of the
bacterial cell.[206-208] Figure 13 shows the DNA-bound and the induced TetR
(complexed to two molecules of [TcMg]+). A schematic overview of the induction
mechanism is given in Figure 14.
49
Figure 13. Left: TetR bound to DNA (no Tc present, PDB[159-161] code 1qpi[209]). Right:
Induced TetR, [TcMg]+ docked into the binding sites (PDB[159-161] code 2trt[210]).
Graphics were generated with DeepView/Swiss-PdbViewer 3.7 SP5[211, 212] and POV-
Ray 3.5.[213]
Figure 14. The resistance mechanism of Gram-negative bacteria against tetracycline
(taken from [214]).
50
Thus, investigations of the TetR/tetO system are very important, not only in
order to elucidate the resistance mechanism but also because TetR/tetO is
used as a controllable switch of gene expression in a variety of organisms, even
in eukaryotic systems.[215-217] TetR regulatory systems may even hold promise
in gene therapy.[215]
3.1.3 Fluorescence Resonance Energy Transfer (FRET)
One of the most promising spectroscopic techniques for the analysis of
biopolymers in action and interaction uses fluorescence resonance energy
transfer (FRET).[218-225] FRET was first described by Theodor Förster in
1948.[218, 219] It is a quantum mechanical process that describes the nonradiative
energy transfer from a donor (D) to an acceptor (A) chromophore (Figure 15).
The energy transfer takes place without occurrence of a photon. The donor and
acceptor are coupled by a dipole-dipole interaction. The basic requirements for
FRET are sufficient overlap of the donor emission spectrum and the absorption
spectrum of the acceptor, proximity of donor and acceptor (D-A separation 10-
100 Å) and an almost parallel orientation of their transition dipoles. In 1967,
Lubert Stryer and Richard Haugland[225] confirmed an important prediction of
Förster theory. They showed that the efficiency of FRET is inversely
proportional to the sixth power of the donor-acceptor separation. A FRET signal
can either be detected from the donor-stimulated fluorescence of the acceptor
or from the quenching of the donor fluorescence. Because of its strong distance
dependency, FRET is frequently used for the examination of molecular
distances in biological macromolecules (“spectroscopic ruler”[225]), an important
technique for the observation of molecular interactions and conformational
changes (e.g. protein folding) with high spatial (1-10 nm) and time (< 1 ns)
resolutions. It can be combined with single-molecule detection techniques, even
in vivo.[222] In addition to these applications, FRET is also used in molecular
diagnostics, e.g. in “real-time polymerase chain reaction”.[222] “Molecular
beacons”[226, 227] and “scorpion oligonucleotides”,[228] are used for quantitative
nucleic acid detection.[222] The activity of proteins has also been monitored by
FRET techniques.[222] Apart from these spectroscopic applications, energy
51
transfer is also an important process in photosynthesis.[229, 230] More details
regarding FRET are given in the computational details section.
Figure 15. Fluorescence resonance energy transfer (FRET) from Trp43 (blue) to the
inducer tetracycline (green). Only one of the two monomers of TetR is shown (PDB[159-
161] code 2trt[210]). Graphics were generated with DeepView/Swiss-PdbViewer 3.7
SP5[211, 212] and POV-Ray 3.5.[213]
3.1.4 Tryptophan Fluorescence
Protein fluorescence is caused by the three aromatic amino acids tryptophan
(Trp, W), tyrosine (Tyr, Y) and phenylalanine (Phe, F). The dominating
fluorophore, tryptophan, displays complex spectroscopic properties. Its
emission spectrum is very sensitive to its local environment and fluorescence
quenchers. Complex time-resolved fluorescence decays are even found for
tryptophan itself, as well as for proteins that contain only a single tryptophan
residue. Such proteins show mostly multiexponential intensity decay curves.[221]
The emission maximum of tryptophan in water occurs near 350 nm and
depends strongly upon polarity and local environment. In most experiments,
protein fluorescence is excited at the absorption maximum near 280 nm or at
longer wavelengths. Therefore, in most cases phenylalanine is not excited.
Protein absorption at 280 nm is due to tyrosine and tryptophan residues. At 23°
C in neutral aqueous solution, the quantum yields of tyrosine and tryptophan
52
are about 0.14 and 0.13, respectively. At wavelengths longer than 295 nm, the
absorption is primarily due to tryptophan. Thus, many papers report the use of a
295 nm excitation, which avoids exciting tyrosine and phenylalanine.[221]
Another complicating factor is the existence of two nearly isoenergetic
excited singlet states called 1La and 1Lb. For most chromophores, the long-
wavelength absorption corresponds to a single electronic transition from the
electronic ground state (S0) to the first excited singlet state (S1). In the case of
tryptophan, indole and their derivatives, the long-wavelength absorption (240-
300 nm) consists of two overlapping transitions to the 1La and 1Lb excited states.
These states have similar energies, and, depending on the environment, either
of them can have the lower energy, although 1La is usually S1. According to
Kasha’s rule, the lowest-energy state is the emitting one. The transition dipoles
for the 1La and 1Lb transitions are oriented almost perpendicular to each other
(Figure 16). Tryptophan emission in solution and most proteins is unstructured
and takes place from the 1La state.[221]
The 1La transition is mainly an excitation from the highest occupied molecular
orbital (HOMO) to the lowest unoccupied orbital (LUMO). Additionally, there is a
smaller amount of a HOMO-1 to LUMO+1 transition. The 1Lb transition is an
almost equal superposition of the HOMO → LUMO+1 and HOMO-1 → LUMO
excitations. Compared to 1Lb, the 1La transition has a sizeable charge transfer
(CT) component, with the electron density primarily decreasing on atoms 1 and
3 while increasing on atoms 4, 7 and 9 (see Figure 16). The resulting increase
in permanent dipole moment is usually computed to be 3-5 D, with the pyrrole
ring becoming more positive.[231] This dipole change causes the tryptophan 1La
transition to be very solvent-sensitive.[231, 232]
The solvent-sensitivity of the indole chromophore can be illustrated by
examining the emission spectra of indole in cyclohexane, ethanol and their
mixtures. In a completely nonpolar solvent, 1Lb can be the lowest-energy excited
state, resulting in a structured emission. The presence of polar solvent
decreases the energy of 1La, so that its unstructured emission dominates and
the unstructured 1Lb emission cannot be observed.[221] A recent study describes
the inversion of the first two excited states of indole depending on solvent
polarity.[233] It analyzes indole solvatochromism in a wide range of different
53
solvents, solvent mixtures and the gas phase. Solvent polarity is measured
using the solvent polarity/polarizability scale (SPP) introduced by Catálan et
al.[234, 235] In solvents with a high solvent polarity (SPP > 0.8) 1La emission is
observed, while in solvents with SPP < 0.8 indole emits from its 1Lb state.[233]
For tryptophan in water, internal conversion between 1Lb and 1La should be
much faster than 400 fs, while solvent relaxation has a time constant of ~ 1.2
ps, as analyzed by femtosecond laser spectroscopy.[236]
In proteins, there is a similar effect of the environment. Tryptophan in a
completely apolar environment shows a blue-shifted structured emission, while
tryptophan near hydrogen-bonding groups or exposed to water emits at longer
wavelengths. In contrast to many other proteins, where 1La emission
dominates,[221, 232, 237] the single-tryptophan protein azurin Pfl from
Pseudomonas fluorescens, where the indole group is located in the
hydrophobic core of the protein, shows a structured emission spectrum
characteristic for the 1Lb state.[221]
Figure 16. Electronic transitions between S1 and S0 in tryptophan. (Transition dipoles
red and blue, respectively. Adapted from [221, 231]).
3.1.5 Time-Resolved Tryptophan Emission: The Rotamer Model and
Alternative Approaches
Time-resolved protein fluorescence spectroscopy provides considerable
information about protein structure and function. Intensity decay curves can be
54
sensitive to bound ligands or other binding partners that affect the mobility of
tryptophan, displace tryptophan residues to new locations or act as quenchers
or energy acceptors in an RET (resonance energy transfer) process.[221]
The fluorescence decay of tryptophan in neutral aqueous solution is known to
be double-exponential, with decay times of 3.1 and 0.51 ns.[221, 238] Early work
explains these two decay times by simultaneous emission from both the 1La and
the 1Lb excited states.[239] However, today there is consensus that the dual
exponential decay results from the presence of rotational isomers[221, 238, 240, 241]
and that tryptophan, unless in a completely nonpolar environment, only emits
from its 1La state.[221] The side chain of tryptophan in solution can exist in
various conformations by rotation about the Cα-Cβ (χ1 angle) and the Cβ-Cγ (χ2
angle) bonds. These conformations or rotamers can interchange on the
nanosecond timescale. Thus, a tryptophan solution can be regarded as a
mixture of these rotamers.[221]
In aqueous solution at pH 7, zwitterionic tryptophan can be self-quenched in
an intramolecular process involving the indole ring and the positively charged
ammonium group by charge-transfer from the indole ring to the adjacent
electrophile.[221] This quenching effect is quite general and occurs for a number
of tryptophan-containing peptides as well as tryptophan itself. The self-
quenching process can explain the biexponential decay of aqueous tryptophan
at pH 7.[221] The basic idea is that quenching is most efficient in one of the
rotamers (ammonium gauche to the indole moiety) and this species is
responsible for the shorter (0.51 ns) decay time.[221] A complete explanation,
however, is much more complex.[221] The assignment of lifetimes to three
rotamers obtained by rotation about the Cα-Cβ bond has been discussed
controversially in the literature.[221, 238, 240, 241] However, the situation becomes
even more complex when rotations about the Cβ-Cγ bond are also taken into
account.[221] The carboxylate group can also act as a quencher, although
weaker than the ammonium group.[238, 240] Nevertheless, a mechanism based on
ground-state rotamers can be regarded as the dominant origin of the
nonexponential decay of tryptophan in a first approximation (“rotamer
model”).[221, 238, 240, 241]
55
In a protein, the amino and carboxyl groups of tryptophan are neutral.
Therefore, the effect of the charged ammonium group is no longer relevant. It is
known that neutral tryptophan analogues have simpler decay kinetics. A model
compound frequently used, N-acetyl-L-tryptophanamide (NATA), has essentially
the same structure as a tryptophan residue in a protein. The fluorescence decay
of NATA is essentially single-exponential,[221, 240, 242] with a lifetime of ~ 3.0 ns in
water at pH 7.[221, 238, 240]
Fluorescence decays of tryptophan in proteins are typically complex and are
mostly fitted with two or more exponential functions. Multiexponential tryptophan
decays in proteins are not surprising for multi-tryptophan proteins. In order to
simplify the situation, single-tryptophan protein mutants are usually
examined.[221] However, even most single-tryptophan proteins have bi- or
triexponential decays. Known exceptions are e.g. apoazurin from Pseudomonas
fluorescens[243] and, under certain conditions, ribonuclease T1 (RNase T1) from
Aspergillus oryzae,[221] which are single-exponential.
The origin of multiexponentiality in single-tryptophan proteins can be
regarded as a result of protein structure. A possible reason for multiexponential
decays is the existence of multiple protein conformations. As amino acids near
the emitting tryptophan can quench its fluorescence, different conformations
result in different decay times.
However, multiexponentiality is even observed when a protein is assumed to
be in a single conformation.[221] In these cases, tryptophan multiexponential
decay in proteins is usually attributed to the existence of multiple ground-state
conformers of the tryptophan residue.[244] Similar to tryptophan in solution, the
“rotamer model” is also used to explain fluorescence decays of tryptophan
residues in proteins. By rotation about the Cα-Cβ and/or the Cβ-Cγ bond, the
tryptophan side chain can exist in different low-energy conformations. Each
conformation has a different microenvironment, leading to a distinct
nonradiative decay rate due to a different proximity and orientation to quencher
groups.[238, 240, 244] Possible tryptophan-quenching pathways that are sensitive to
environment are excited-state electron and proton transfers and solvent
quenching.[244] Charge transfer from the excited indole moiety to the carbonyl
group of the peptide bond is probably the main nonradiative quenching pathway
56
of tryptophan in proteins.[240, 244] Another possible quenching process is excited-
state proton transfer from a strong proton donor to an adjacent aromatic carbon
of the indole moiety. Histidine residues can quench tryptophan fluorescence by
charge-transfer (CT) or proton-transfer mechanisms.[244]
Many studies explain multiexponential tryptophan decays in proteins on the
basis of the rotamer model or using modifications of this model.[244-252]
Alternative models to the rotamer model invoke spectral migration due to
structural relaxation of the surrounding protein, perhaps also coupled with
changes in the arrangement of the indole chromophore relative to the backbone
of the surrounding protein domain or collisional fluorescence quenching. A
quantitative model describing such a situation has been developed by Toptygin
and Brand (“relaxation model”).[245, 253, 254]
FRET from an excited donor to two or more acceptor molecules can explain
multiexponential decays of proteins without invoking multiple protein
conformations. Bajzer and Prendergast showed that energy transfer from an
excited donor to N acceptor molecules in the surrounding protein matrix yields
2N exponential components in the decay function.[255]
Heme proteins like myoglobin are examples where very efficient RET can
occur, yielding very short (ps) lifetimes.[221] Another example of a protein with
FRET quenching of tryptophan is the Tet-repressor/tetracycline complex, where
the bound inducer tetracycline acts as energy acceptor (see following section).
In the present study, we show how a biexponential tryptophan fluorescence
decay is obtained from a molecular dynamics (MD) trajectory containing two
tryptophan rotamers by calculating the FRET rate constant for every MD
snapshot. FRET rate constants are calculated from the results of a QM/MM CI
calculation of the protein (see below).
An alternative approach to the traditional analysis of decay curves using a
small number of exponential components that correspond to particular protein
conformations assumes continuous lifetime distributions. These lifetimes
correlate to a large number of conformations that interconvert at a rate of the
same order of magnitude as the excited-state decay rate.[256] In these cases,
analysis can be performed using numerical methods.[257] Especially results for
57
proteins that contain a large number of tryptophans can be interpreted using
lifetime distributions.[221, 258]
Another complicating factor that leads to nonexponential decays arises from
donor-acceptor distance distributions in the case of FRET quenching. In the
ideal case, a relatively constant donor-acceptor distance yields a
monoexponential decay curve. If the donor-acceptor distance is a distribution
with a relatively broad FWHM (full width at half maximum), the donor decay
becomes more than single exponential. It is possible to recover the donor-
acceptor distance distribution from the nonexponential decay of the donor.[221]
3.1.6 Spectroscopic and Computational Investigations of the Tet
Repressor/Tetracycline System
Much of the experimental data available for the tetracycline repressor and its
complexes with tetracyclines is derived from time-resolved fluorescence
spectroscopy.[244-249, 259] Tryptophan 43 (Trp43, W43), an amino-acid residue
situated in the DNA-binding domain of the tetracycline repressor, is frequently
used as a probe for exploring the conformation of the protein in time-resolved
fluorescence measurements of TetR.[244, 246, 248, 249]
The fluorescence-decay curves obtained from these measurements are
usually fitted using 2 or 3 exponential functions, suggesting that species with 2
or 3 different fluorescence lifetimes are present.[244, 246, 248, 249] Usually, these
different tryptophan species are interpreted in the framework of the “rotamer
model”, which assigns the 2-3 lifetimes found for Trp to 2-3 discrete rotamers of
this residue.[244, 246, 248] Some of the published experimental studies employ
molecular modeling approaches based on molecular mechanics in order to
assign Trp rotamers to the lifetimes observed experimentally.[244, 247]
In our study, we have investigated the complex formed between class D TetR
and tetracycline.[210] The experimental analogue to our model is the complex
formed by the single-tryptophan TetR mutant W75F with phenylalanine
substituting tryptophan at position 75 in both subunits and therefore only Trp43
58
present, for which a biexponential Trp decay with lifetimes of 2.55 and 0.42 ns
is reported.[246]
In 1986, first evidence for fluorescence resonance energy transfer (FRET)[218-
221] from Trp43 to the inducer tetracycline was reported.[260] Since then, several
authors have described FRET from Trp43 to the inducers tetracycline and
anhydrotetracycline. The emission of tryptophan 43 at ~330-380 nm and the
absorption of tetracycline at ~360-400 nm show sufficient overlap for FRET.[244,
246, 259, 260] FRET from a tryptophan residue to the inducer is not restricted to
Trp43, it has also been reported for single Trp mutants of TetR with Trp in other
positions.[247, 259]
In 1993, Bajzer and Prendergast showed that energy transfer can explain
multiexponential decays of tryptophan in proteins.[255] As the existing studies of
the TetR/tetracycline complex have not investigated the effects of FRET from
different Trp rotamers to the inducer on Trp kinetics, we have used our
molecular dynamics-QM/MM approach in order to examine the effect of FRET
quenching in different Trp rotamers.[261-265]
3.1.7 Overview of the Method
The existing models for the interpretation of TetR spectroscopic data are
largely speculative, so that conclusions for the induced tetracycline repressor
(TetR) need to be validated by computer simulations in order to confirm the
interpretation of the experimental results.
Therefore, we have developed a combined molecular dynamics-QM/MM
approach, which allows us to simulate absorption and fluorescence spectra as
well as fluorophore quenching by FRET.[261-265]
A classical molecular dynamics simulation, for which we use the AMBER[167,
266, 267] program, gives “hot” geometries of a protein, which are the basis for
combined quantum mechanical/molecular mechanical (QM/MM) single point
configuration interaction (CI)-calculations using the semiempirical program
package VAMP.[50, 51]
59
The relevant chromophores are embedded in the protein environment and
the solvent using a quantum mechanics/molecular mechanics (QM/MM)-CI-
approach in which the environment is represented by a classical force field. The
QM/MM method uses a rigid point charge environment with associated van der
Waals radii.[83, 104, 164] The chromophores are treated quantum-mechanically
(QM part) and are polarized by the charges in the environment (MM part). The
polarization of the QM wavefunction by the point charges in the MM
environment is taken into account via an additional one-electron term in the
Fock matrix.[104, 164, 268]
The semiempirical CI calculations provide all the variables necessary to
calculate both the absorption and fluorescence spectra and the FRET energy-
transfer probabilities according to Förster theory,[218-221] as shown below.
In contrast to calculating UV/Vis spectra directly from relaxed minimum
geometries, which neglect the dynamics of the molecule modeled, our approach
allows the simulation of a protein system in solution at 298 K.
In our simulation, properties of a bulk system that consists of a large number
of molecules in reality are simulated by analyzing the temporal evolution of a
single-protein test system. This resembles the determination of bulk
fluorescence properties from single-molecule measurements, as described by
Hochstrasser, Lee and Tang.[269] Both approaches derive ensemble properties
from single-molecule data (“ensemble model”).
Figure 17. Overview of the combined molecular dynamics-QM/MM approach.
60
3.1.8 Molecular Dynamics Simulations for the Simulation of Protein Spectra
and FRET
A widely used approach for the simulation of spectroscopic data is the use of
a classical force-field molecular dynamics simulation, which generates the
chromophore geometries. In a first approximation, as described in a study of
fluorescence anisotropy decays,[270] the transition dipole of the chromophore
can be assumed to have a fixed geometry relative to the aromatic ring system,
e.g. an angle of 38° relative to the indole long axis for 1La in tryptophan.[270]
Several studies have combined classical methods with quantum mechanics
in order to simulate spectral properties of proteins. Combined quantum
mechanical and molecular mechanical (QM/MM) calculations and molecular
dynamics simulations of bacteriorhodopsin in the membrane matrix were used
to investigate the opsin shift.[271] Molecular dynamics simulations of cytochrome
c in explicit solvent were performed to investigate the structural distortions in the
protein, and semiempirical calculations were used to calculate the resultant
changes in the spectroscopic transitions.[272]
The fluorescence wavelengths of tryptophan in different proteins were
investigated using a combination of classical molecular dynamics with QM/MM
calculations.[232, 273] In this work, the quantum mechanics were used to assign
the charges to the Trp chromophore and to give the electronic transition energy
as a function of the reaction field produced by the MD force field. The Trp
dynamics and atomic coordinates were governed entirely by the MD with QM
modified charges on Trp. A trajectory for a solvated protein was started with the
Trp charge distribution and geometry of the ground state. Every 10 fs, the MD
trajectory was interrupted, an electronic structure calculation performed and the
geometry and charges of the Trp residue updated. Following 5 ps of ground-
state dynamics, electronic excitation was initiated by switching the charges and
geometry to those of the 1La state, and continuing the trajectory out to typically
30 ps, interrupting every 10 fs to perform a QM calculation as described
above.[232] An enhanced version of this protocol was used to examine
tryptophan fluorescence quantum yields in different proteins.[274]
In principle it is possible to use quantum mechanical molecular dynamics
simulations for the calculation of protein spectra. Such calculations are,
61
however, computationally very demanding and are therefore only applied to
small systems and for the simulation of short trajectories. An example is the
simulation of an infrared spectrum of p-benzoquinone in water by analyzing the
results of a QM/MM molecular dynamics simulation with Fourier transforms of
autocorrelation functions or instantaneous normal-mode analysis of
snapshots.[275]
Several groups have reported the use of computational methods to simulate
fluorescence quenching by FRET. A protein structure obtained by sequence
alignment was refined by force-field minimization and the structural data were
used to calculate Trp fluorescence quenching within the framework of Förster
theory and the rotamer model. The transition moment for 1La in tryptophan was
assumed to be fixed relative to the chromophore (in the indole plane, 39±7°
clockwise from the pseudosymmetry axis of indole).[276] Another study
investigated polymer folding using Brownian-dynamics simulations. Donor-
acceptor distances obtained from these simulations were used to calculate
FRET data using the Förster radius and the FRET rate constant (see below) as
parameters.[277] Monte-Carlo simulations and a statistical model were used to
model FRET in a system of actin filaments that includes interactions between
multiple donors and acceptors and photobleaching. As no transition dipoles
were used, a mean value of 2/3 for the orientation factor κ2 (see below) was
assumed.[278]
Classical molecular dynamics simulations have been used increasingly for
simulating FRET recently. The structural information obtained from MD
simulations of Poly(benzylphenyl ether) dendrimers was used[279] to calculate
rate constants for FRET within the Förster framework (see below for details).
Again, the orientation factor κ2 was assumed to be 2/3 and the donor-acceptor
distance was obtained from the MD simulation.[279] A study investigating the
effect of phosphorylation on the conformational states of a small peptide (a 10-
residue mitogen-activated protein kinase substrate found near the N-terminus of
tyrosine hydroxylase) also used MD simulations to determine donor-acceptor
distances.[280] Using a literature value for the Förster distance R0 of the donor-
acceptor pair, i.e. without considering the transition dipoles of the
62
chromophores, FRET efficiencies were calculated using a statistical mechanical
approach. The results agreed well with experimental values.[280]
Tryptophan-tryptophan FRET homotransfer as one fluorescence quenching
process among others in the HIV-1 integrase catalytic core has also been
investigated by molecular dynamics simulations.[281] Transition dipoles of a
model chromophore (3-methylindole) in a ground-state geometry calculated by
a single QM CI calculation were used to describe the S0 ↔ 1La transitions;
FRET data were calculated and the simulated fluorescence decays were
compared with experimental data.[281]
Conformational distributions of prion protein repeat systems were
investigated using a combination of FRET experiments and molecular dynamics
simulations.[282] The MD data were used to assist the interpretation of
experimental steady-state fluorescence FRET data and to characterize the
conformational distributions of the system. Typical transition moments for the
chromophores were obtained by a model QM calculation of the donor excited
state and the acceptor ground state. For each MD snapshot, the orientation
factor κ2 was calculated using these typical transition moments and the FRET
efficiency was calculated using the donor-acceptor distance from the MD
trajectory. An averaged energy transfer efficiency was calculated using the
Förster equations and these data were compared with the results from
experimental steady-state fluorescence measurements. The experimental
Förster distances were also compared with calculated Förster distances
obtained from the simulations.[282]
In all the above studies of FRET, the MD force field represented the
electronic ground state and the force field was not adapted to the excited state
(in contrast to the work of Callis and Vivian[232]). In contrast to our approach,[261-
265] the transition dipoles representing the absorption/emission process between
the electronic ground state and the first excited singlet state were not calculated
quantum mechanically for every MD snapshot in the above papers.
63
3.2 Computational Details
3.2.1 Validation: Glycyltryptophan Absorption and Emission Spectra
For a validation of our approach, we performed a short MD simulation of
zwitterionic glycyltryptophan (Gly-Trp, GW) in water. Zwitterionic
glycyltryptophan was generated and solvated in a water box of 734 water
molecules (exceeding the peptide dimensions by 10 Å in each direction) using
xLEaP from the AMBER 5 suite.[167] We used the Cornell et al. 1995 force field
for minimization and dynamics.[66] After 10 steps of steepest-descent energy
minimization 24,990 conjugate gradient minimization steps were performed
using Sander.[167] We used periodic boundary conditions and a cutoff of 12 Å for
electrostatic and van der Waals interactions. A constant dielectric of 1 was
used, the non-bonded pair list was updated every 25 steps. After minimization,
2 ns of constant pressure molecular dynamics were performed. The system was
heated to 300 K, after 3 ps the system was in equilibrium. The equilibration data
were not used for later processing. The step size was 2 fs and bonds involving
hydrogen were constrained using SHAKE. We used a time constant for heat-
bath coupling of 0.2 ps, temperature was kept constant using the Berendsen
coupling algorithm.[73] Snapshot geometries were dumped every ps, giving 1997
structures to be processed further.
The MD snapshot geometries were converted into VAMP[50] inputs by Perl[186]
scripts. We performed VAMP AM1[36] single point CI calculations with an active
space of 26 orbitals around the HOMO-LUMO border. Single and pair double
excitations were considered in the CI.[283] For comparison, we tested several
types of semiempirical calculations. QM/MM calculations were performed using
a rigid point-charge environment. The QM/MM implementation in VAMP is
described below. Glycyltryptophan was chosen as QM-part, while the water
molecules formed the MM-part. MM charges were taken from the force field
used in the AMBER MD simulation before.[66] Another set of QM-calculations
considered the influence of the surrounding solvent using the self-consistent
reaction field (SCRF)-technique as described by Rauhut, Clark and Steinke.[268]
Data extraction and calculation of absorption and fluorescence data were
performed using Perl[186] scripts. The emission wavelength for S1 → S0 was
64
calculated for an electronically and solvent relaxed excited state. For the
simulation of steady-state and time-resolved fluorescence spectra, we
calculated the Einstein coefficient of spontaneous emission Anm for each
snapshot from the transition moment mnMuuv
and the wavelength of the S1 → S0
transition (n = 0, m = 1).[284, 285]
(28)
If the frequency of the electronic transition ν is given in [s-1], the transition
dipole mnMuuv
in [Cm], the Planck constant h in [Js], the permittivity of free space
ε0 in [AsV-1m-1], and the velocity of light in vacuum c in [ms-1], Anm has the
dimension [s-1], i.e. it represents the rate constant for the intrinsic fluorescence
deactivation of the first excited singlet state S1.
A steady-state fluorescence spectrum of the whole ensemble/the whole MD
trajectory with a resolution of 0.1 nm was obtained as follows: The Einstein
coefficients Anm were binned in wavelength intervals of 0.1 nm. The spectrum
was constructed by plotting the sums of the Anm values in each bin vs. the
emission wavelengths λ.
In order to obtain a time-resolved fluorescence decay that represents the
intrinsic fluorescence emission of the system, the concentration of the first
excited singlet state is calculated by summation over all MD snapshots,
assuming for each snapshot a first-order decay of S1 with a rate constant Anm,i
(Eq. (29), t: relative time, Anm,i: Einstein coefficient of spontaneous emission, i:
snapshot). The concentration of S1 is proportional to the emission and is
therefore used as a measure of fluorescence intensity.[221]
(29)
In order to interpret the results in terms of a mono-, bi- or three-exponential
model with lifetimes τj, we fitted the calculated fluorescence decays with sums
of exponentials using the Levenberg-Marquardt algorithm as implemented in
Origin 7.0.[286] (αj: preexponential factors/amplitudes; j = 1,2,3, …)
,
1( ) nm iA t
Si
c t e−=∑
3 3 2
30
163
mnnmA Mh cπ νε
= ⋅uur
65
(30)
Alternatively, we also performed fits that allowed an additional constant C.
This is not physically correct for our simulated data, because the fluorescence
intensity should be 0 for t → ∞. In an experimental study, such a constant could
result from background noise.
(31)
The number of exponential components was increased progressively until the
fit was judged adequate from the reduced χ2 value and the squared correlation
coefficient r2. The fitting curve and the weighted residuals, which represent the
difference between the fitting function and the actual data points, each divided
by the y-value of the data point, were inspected visually.
In the fitting procedure, the parameters of the fitting function are varied until
the reduced χ2 value, which describes the goodness of the fit and should be
near unity, is at its minimum. Values of 2rχ much smaller than one indicate that
an insufficient number of data points is used in the analysis. (xi: independent
variable; yi: dependent variable in the data set to be fitted; f(xi, p1, p2,…): fitting
function that represents the theoretical model believed to explain the process
that produced the data set; p1, p2, …: parameters of the fitting function; n: total
number of data points; p: number of adjustable parameters, n-p: number of
degrees of freedom) For 2rχ calculation, each data point was weighted with a
quantity of yi-1.[221, 287]
(32)
Another measure to describe the quality of the fit is the squared correlation
coefficient r2, which should also be near unity. It is calculated as follows.[4]
(33)
1
( ) jt
jj
I t e τα−
=∑
( ) 22 21 1 11 2 1 2( , ,...) ; , ,...
ir i in p n p yi
p p y f x p pχ χ− −= = −⎡ ⎤⎣ ⎦∑
1
( ) jt
jj
I t e Cτα−
= +∑
2 2 2 2, , ,
22 2 2
( ) ( ) ( ) ( )1
( ) ( ) ( )
calc i i i calc i i calc ii i i i
i i ii i i
y y y y y y y yr
y y y y y y
− − − − −= = = −
− − −
∑ ∑ ∑ ∑∑ ∑ ∑
66
In this initial study, we used the nonlinear least-squares analysis described
above for the analysis of the calculated decay curves. However, there are other
models available for the analysis of complex decay data, e.g. the method of
moments, Laplace transformation, the maximum entropy method, Prony’s
method, sine transforms, phase plane and analytical methods.[221, 257]
Absorption spectra were generated analogously to the static fluorescence
spectra by plotting the sums of the oscillator strengths for a given wavelength
vs. the wavelength. All transitions S0 → Si in all snapshots were considered.
Wavelength data were calculated for the vertical transition between S0 and the
target state.
3.2.2 Simulating FRET in TetR: Molecular Dynamics Simulations
As a spatial basis for the model generation, we used an X-ray structure of the
induced class D tetracycline repressor (TetRD) dimer (PDB[159-161] code 2trt[210]).
The missing residues of the unresolved loop in the X-ray structure (residues
156-164) were inserted into the protein.[288] The model contains two identical
monomers generated by symmetry transformation. Each monomer consists of
the repressor protein (Ala2 to Ile207) co-crystallized with one tetracycline (Tc)
molecule and one Mg2+-ion. Crystal water molecules were deleted before
further processing, with the exception of the three active-site water molecules in
each monomer coordinating the active-site magnesium ion. All molecular
mechanics calculations were performed using the AMBER 1994 force field
parameters (parm94) of Cornell et al.,[66] as implemented in the AMBER 5.0
software package.[167] Charges for the ligand in the zwitteranionic form were
derived after geometry optimization with the AM1 Hamiltonian[36] implemented in
the semiempirical program package VAMP 7.0a[289] using the VAMP
electrostatic potential fit procedure VESPA.[290]
The system was solvated in a water box of 9,337 TIP3P water molecules
(exceeding the peptide dimensions by 8 Å in each direction) using xLEaP from
the AMBER 5 suite.[167] Two Na+-ions were added to achieve neutrality of the
system, these ions replaced two TIP3P water molecules. The whole system
67
consisted of 34,641 atoms, 6,616 of which were atoms of the dimeric TetR/Tc
complex (Figure 18).
The system was subjected to 10 steps of steepest descent optimization,
followed by 5,000 steps of conjugate-gradient minimization, leading to an RMS
gradient of 0.13 kcal mol-1 Å-1, which is sufficient as a starting point for the MD.
During minimization, weak harmonic restraints with force constants of 5.0, 1.0
and 1.0 kcal mol-1 Å-2 on the main chain atoms, carbon atoms of Tc and oxygen
atoms of the active-site water molecules were applied.
Constant pressure periodic boundary water box MD simulations were
performed using the particle-mesh Ewald (PME) method.[75, 291] Taking the
energy-minimized coordinates as starting geometry, the system was
equilibrated for 2.4 ns, raising the temperature gradually to 300 K by coupling to
a temperature bath with a time constant of 0.2 ps, according to Berendsen.[73]
During the equilibration, positional restraints on main chain atoms, carbon
atoms of Tc and oxygen atoms of the active-site water molecules were reduced
stepwise. (0-300 ps: force constants of 5.0, 1.0 and 1.0 kcal mol-1 Å-2 on main
chain atoms, carbon atoms of Tc and oxygen atoms of the active-site waters;
300-600 ps: 1.0, 0.2, 0,2 kcal mol-1 Å-2; 600-900 ps:0.5, 0.1, 0.1 kcal mol-1 Å-2;
900-1200 ps: 0.1, 0.02, 0.02 kcal mol-1 Å-2; 1200-2400 ps: 0.02, 0.0, 0.0 kcal
mol-1 Å-2) In the production phase (after 2.4 ns) no restraints were applied. The
equilibration part of the MD trajectories was not used for data collection.
The final production run was carried out for 7,094.3 ps at 300 K. We used an
integration step size of 2.0 fs. Atomic coordinates, temperatures, and energies
were recorded every ps and used for subsequent analyses. SHAKE was used
to constrain bonds involving hydrogen. A cutoff of 10 Å applied to van der
Waals interactions was used during all simulations, the 1-4 electrostatic and 1-4
van der Waals interactions were scaled by 1/1.2 and 1/2, respectively. The non-
bonded pair list was updated every 25 steps. The calculations were run using a
constant dielectric constant (ε) of 1 because explicit solvent molecules were
present and periodic boundary conditions were applied. The dimensions of the
PME charge grid (upon which the reciprocal sums are interpolated) were
adjusted to be approximately equal to the box dimensions after each restart,
68
ensuring a constant grid spacing of approximately 1.0 Å. The order of the B-
spline interpolation for PME was 4, which implies a cubic spline approximation.
The MD simulations were performed using an MPI-parallel executable of
AMBER on an IA32/EM64T Xeon cluster.[292]
The trajectory segments obtained from various MD restarts were processed
and analyzed with Ptraj 6.5.[293] In several restart files one of the monomers was
found in a neighboring box. As no restraints were applied in the production
phase, a protein rotation of about 40° occurred in the second half of the
production phase. Due to the use of periodic boundary conditions, these effects
are only computational artifacts and do not affect the MD simulation. As we
needed a “proper” protein structure for analyzing the MD results and for QM/MM
input generation, the protein monomers were moved back into the box using the
“image” command of Ptraj. In order to remove protein rotation, the whole
trajectory was fitted on the minimized X-ray structure (only Cα atoms were used
for fitting and RMSD calculation). For assessing protein stability in the MD
production phase, an average structure of all production phase snapshots was
calculated using Ptraj and compared with the X-ray structure. For QM/MM input
generation, we processed all snapshot geometries with Ptraj so that the water
box was centered on Asp53. This was necessary in order to obtain structures
with the QM parts (Trp43 and Tc) roughly in the center of the QM/MM system.
As it was not necessary to take FRET between the two monomers into account
(see below), the monomers could be computed separately in the QM/MM part of
the protocol. Therefore, centering on Asp53/Asp53’ was performed separately
for each monomer.
69
Figure 18. The system used in the MD simulations of the TetR/Tc complex. Graphics
were generated with DeepView/Swiss-PdbViewer 3.7 SP5[211, 212] and POV-Ray 3.5.[213]
3.2.3 Simulating FRET in TetR: Quantum Mechanical Simulations
During the MD production phase geometries were dumped every ps, giving
7,093 snapshots to be processed for QM/MM input generation. These MD
snapshot geometries were converted into VAMP[50] inputs by Perl[186] scripts.
Initial QM calculations did not use any MM environment. In these preliminary
calculations only the chromophores tetracycline and tryptophan were extracted
from the protein and single-point QM calculations for each snapshot were
performed in the gas phase. Because of the high dependence of the Trp
transition dipole on the local environment,[231] we improved our protocol and
used QM/MM single point CI calculations of the snapshot geometries extracted
from the MD trajectory. Two different sets of QM/MM-CI calculations with VAMP
8.1[50] were performed for each monomer of TetR, one of them with the donor
and the other with the acceptor in the QM part. The MM environment was
represented by a rigid point-charge environment, as described by Clark et al.[104,
70
164] MM charges were taken from the AMBER 1994 force field.[66] The
polarization of the QM wavefunction by the point charges in the MM
environment was taken into account via an additional one-electron term in the
Fock matrix.[104, 268] Three major contributions were considered for the
calculation of the QM/MM interactions: van der Waals and Coulomb interactions
and the polarization of the QM part by the MM environment. The van der Waals
parameters describing the QM/MM interaction were taken from the Universal
Force Field (UFF) by Rappé et al.[166] No back-polarization of the MM part by
the QM part was included. As explicit solvent molecules were present, we used
a constant dielectric constant of 1.0 for both the MM- and the QM/MM-Coulomb
interactions. For more details on our QM/MM implementation, see [83, 84, 95, 104,
164].
As model calculations had shown that the high charge of the active-site Mg2+
ion in direct neighborhood of tetracycline disturbs the QM part too much, we
defined tetracycline together with the magnesium ion as the QM part in the set
of calculations for the FRET acceptor tetracycline, giving a good representation
of the electronic properties of the chromophore. The QM part with Tc and Mg2+
consisted of 65 atoms, which were surrounded by 34,576 MM atoms.
In the case of the CI calculations of the FRET donor tryptophan, the protein
backbone was cut for QM input generation. In our initial calculations without
environment, the amide bonds between tryptophan and its neighbors were cut.
This cutting scheme caused the problem that the Trp S0 → S1 transition did not
contain the π → π* transition characteristic for the indole chromophore, but
consisted practically exclusively of an excitation from several backbone σ-
orbitals to an antibonding backbone carbonyl-π* orbital at the C-terminal end of
the QM part. In order to derive the ideal cutting scheme for the QM/MM
calculations, we tested several patterns for QM/MM cutting. We tested the 3-
methylindole fragment of Trp43, tryptophan itself, the tripeptide Gly-Trp43-Gly
(Gly instead of Tyr42 and Hip44) and Tyr42-Trp43-Hip44 (which is the real
situation in the protein) as QM parts. For all these patterns we conducted
QM/MM-CI calculations with coordinates of all snapshots of the MD trajectory,
and additionally, model TORQUE CI calculations of the chromophores with
VAMP 8.1.[50] In the TORQUE CI calculations, the angles χ1 and χ2 were
71
changed systematically in steps of 30° and the molecules were optimized with
exception of the χ1 and χ2 angles. Both the QM/MM CI calculations of the
snapshots of the complete MD trajectory and the TORQUE model calculations
showed that the chromophore is best described if a tripeptide Gly-Trp-Gly with
the coordinates of Tyr42, Trp43 and Hip44 is used as QM part in the QM/MM-CI
calculations of the MD snapshots. The cutting scheme finally used is shown in
Figure 19. For more details regarding the QM/MM system partitioning, see [84,
96]. In the CI calculations for the FRET donor chromophore, the 44 QM atoms of
the Gly-Trp-Gly fragment were embedded in 34,597 MM atoms.
In the AM1[36] QM/MM single-point CI calculations of both the FRET donor
and the acceptor (7,093 calculations for each monomer) geometries, we chose
an active space of 26 orbitals symmetrically distributed around the HOMO-
LUMO border. Single and pair double excitations were considered in the CI.[283]
For all donor and acceptor snapshots, the transition dipoles for all S0 → Si
transitions were calculated.
The 28,372 QM/MM-CI calculations were performed on an IA32/EM64T Xeon
cluster.[292]
72
Figure 19. Cutting the protein backbone for the QM/MM calculations. The van der
Waals radii of the encircled hydrogen atoms were reduced to 0.1 Å, in order to avoid
van der Waals clashes with the hydrogens added to saturate the QM part. See [84, 96] for
details on the treatment of the QM/MM border.
3.2.4 Calculation of FRET Data
A fluorescing molecule in an excited state emits radiation with a rate constant 1Dτ − . In the presence of an energy acceptor, and if donor emission and acceptor
73
absorption show a sufficient overlap, energy can be transferred from the donor
D to the acceptor A via dipole-dipole interaction.[221]
The efficiency, E, of energy transfer is the fraction of photons absorbed by
the donor that are transferred to the acceptor. This fraction is given by
(34),
which is the ratio of the transfer decay rate kT to the total decay rate of the
donor.[221] The transfer efficiency is typically measured using the relative
fluorescence intensity of the donor, in the absence (FD) and the presence (FDA)
of the acceptor. The FRET efficiency can also be calculated from the lifetimes
under these conditions (τD and τDA).[221]
(35)
(36)
The total rate constant for donor deactivation is given by the sum of the rate
constant for intrinsic donor fluorescence (in absence of an acceptor), the FRET
rate constant (kT) and the sum of the rate constants of other nonradiative
processes (knr), like internal conversion or intersystem crossing.
(37)
The Förster rate of excitation transfer depends on the lifetime of the donor in
absence of an acceptor 1Dτ − and strongly on the Förster distance R0 and the
donor-acceptor distance r.[221]
(38)
The Förster distance R0 is the donor-acceptor (D-A) distance at which FRET
occurs with 50 % efficiency and at which the FRET rate is equal to the decay
rate of the donor in the absence of the acceptor. R0 depends on the donor
quantum yield in the absence of an acceptor QD, the orientation factor κ2, the
601( )T
D
Rk rrτ
⎛ ⎞= ⎜ ⎟⎝ ⎠
1D T nrk k kτ −= + +
1T
D T
kEkτ −=
+
1 DA
D
E ττ
= −
1 DA
D
FEF
= −
74
normalized overlap integral describing the overlap of donor emission and
acceptor absorption 4
0
( ) ( )D AF dλ ε λ λ λ∞
∫ , Avogadro‘s number N and the refractive
index of the medium n, which is typically assumed to be 1.4.[221]
(39)
The orientation factor κ2 describes the relative orientation of the donor and
the acceptor transition dipoles.[221]
(40)
θT is the angle between the emission transition dipole of the donor and the
absorption transition dipole of the acceptor, θD and θA are the angles between
these dipoles and the line connecting donor and acceptor (Figure 20).
Depending on the relative orientation of donor and acceptor, κ2 can range from
0 (orthogonal dipoles) to 4 (collinear and parallel dipoles). For parallel dipoles,
κ2 = 1 (Figure 21). If donor and acceptor rotate freely, the mean value of κ2 is
2/3, a value that is assumed in many studies, as mentioned above.[221]
In general, the transition-dipole directions are the directions of light
polarization that give maximum absorbance. The transition dipole is not a
measure of the change in dipole accompanying the change in state. Its direction
is the direction that electron charge is displaced by mixing the ground and
excited states during the transition and involves the product of the ground- and
excited-state wavefunctions (transition density). The sign of the transition dipole
direction, i.e. the phase of the transition density, is not physically significant
when considering a single molecule because the absorption and emission rates
depend on its square.[231] Therefore, there is no difference between a transition
dipole mnMuuv
and mnM−uuv
calculated for a single chromophore in our model.
In order to check the model, we simulated the range of values of κ2 for two
freely movable dipoles with a Perl[186]-script. As expected, the range of values
was 0 ≤ κ2 ≤ 4 and the mean value of κ2 was 0.667898.
2 2(cos 3cos cos )T D Aκ = Θ − Θ Θ
26 40 5 4
0
9000(ln10) ( ) ( )128
DD A
QR F dNn
κ λ ε λ λ λπ
∞
= ∫
75
Figure 20. Calculation of the orientation factor κ2 from the donor (D) and the acceptor
(A) transition dipoles and the line connecting the donor and the acceptor chromophore
centers.
Figure 21. Typical dipole-dipole situations.
For calculating FRET rate constants from QM/MM-CI data, we only evaluated
FRET from the donor Trp to the Tc acceptor within the same monomer of TetR.
Because of the high distance-dependence of FRET (~r-6), FRET between the
two monomers of TetR should be orders of magnitude (factor 2.5 ⋅ 107) weaker
than FRET within one monomer (Trp-Tc distance within a monomer: 32.13 Å,
Trp-Tc distance between two monomers: 49.25 Å; both distances calculated for
Trp43 CE2 – Tc C1A). Therefore, we calculated FRET rate constants for the
donor-acceptor pairs in the two monomers independently from each other.
According to Kasha’s rule, Trp emission always occurs from the S1 state.
Therefore, the transition dipole for this transition in each snapshot was used for
further analysis. In the range of energetically accessible acceptor (Tc)
2 1κ =2 4κ = 2 0κ =
2 2(cos 3cos cos )T D Aκ = Θ − Θ Θ
76
transitions of a given snapshot (higher wavelength value than the corresponding
Trp emission wavelength) the transition with the highest oscillator strength was
used for the calculation of the orientation factor κ2. The orientation factor for
energy transfer in each monomer was calculated for each MD snapshot using
Eq. (40).
According to Eq. (39), the Förster radius only depends on the orientation
factor κ2 if the donor quantum yield, the refractive index of the medium and the
D-A spectral overlap are assumed to be constant. Such an assumption is
justifiable for our system because there is no reason why these values should
differ significantly between two snapshots.
(41)
Therefore, in our model a relative FRET rate constant can be calculated from
the orientation factor and the donor-acceptor distance r. The relative FRET rate
constant kT,rel. was calculated for each MD snapshot.
(42)
(43)
We constructed the decay of the Trp S1 state from these relative FRET rate
constants, assuming that no other deactivation processes than FRET occur in a
first approximation. The range of FRET rate constants represents a multi-
exponential decay in the ensemble of structures in the MD trajectory. Each of
the snapshots decays according to a first-order law with a relative rate constant
kT,i.. Therefore, we can express the concentration of Trp molecules in the S1
state 1( )Sc t by summation over all snapshots (Eq. (44), kT,i: relative FRET rate
constant for snapshot i).
(44)
The decay curves thus calculated were fitted with sums of exponential
functions in order to obtain relative lifetimes as described in the validation
26
1~ Tk rκ ⋅
6 20 ~ R κ
2, . 6
1 T relkr
κ= ⋅
,
1( ) T ik t
Si
c t e−=∑
77
section. The intrinsic fluorescence decays (without FRET quenching) of Trp43
were calculated analogously using the Einstein coefficients calculated for each
snapshot (Eq. (29)). Donor and acceptor emission/absorption spectra were
calculated as shown in the validation section.
For visual inspection of the transition dipoles calculated for the donor and the
acceptor, we extracted the relevant transition dipoles for each donor and
acceptor snapshot from the VAMP result files. Data processing was performed
using Perl.[186] For better visualization, the donor S1 → S0 vector was multiplied
by 25 and the starting point of the vector representative was chosen so that the
transition dipole had its midpoint in the center of the indole chromophore. The
xyz-coordinate files thus obtained were fitted on the first snapshot structure of
the MD production phase considering the heavy atoms of the indole
chromophore with Quatfit.[294] The fitted transition dipoles were superimposed
together with the first snapshot structure for a graphical representation of the
variation of transition dipoles relative to the indole chromophore. Additionally,
we produced a multiple xyz-file for an animation of snapshots together with
transition dipoles. We used Jmol 5[295] for visualization of xyz-files with transition
dipoles.
As there is no difference between a transition dipole mnMuuv
and mnM−uuv
in our
model, we multiplied all transition dipoles pointing to the “wrong” direction
(which was discriminated in terms of neighborhood to atom 7 in Figure 16) by -1
prior to cluster analysis. As data analysis with TSAR 3.3,[296] which only allows
cluster analyses to be performed for two or more clusters, indicated the
presence of only one cluster of tryptophan transition dipoles, we determined the
structure representing the cluster center by calculating the mean transition
dipole and selecting the structure with the smallest RMSD with respect to the
mean transition dipole. This structure was used for a detailed analysis of the
molecular orbitals of the donor chromophore with MS Modeling 3.01.[51]
For the FRET acceptor tetracycline, an analogous approach for visualization
of transition dipoles was chosen. In this case, the transition with the largest
oscillator strength in the range of energetically accessible transitions was
selected for each snapshot. Additionally, all energetically accessible transition
dipoles were inspected visually for representative snapshots.
78
3.3 Results and Discussion
3.3.1 Validation: Glycyltryptophan Absorption and Emission Spectra
The steady-state absorption and emission spectra simulated for zwitterionic
glycyltryptophan are shown in Figure 22. The semiempirical calculations with
the continuum solvent model SCRF gave better values for the fluorescence
wavelengths than the other set of calculations, in which the solvent was
modeled as the MM environment. We therefore used the SCRF results to
generate glycyltryptophan spectra. The band shapes of the experimental
spectra (Figure 23) are reproduced well. However, both the simulated
absorption and the emission spectrum are blue-shifted, by about 60 and 20 nm,
respectively. Additionally, the FWHM of the fluorescence peak is too small
compared to the experimental values. The blue shifts of the calculated spectra
relative to experiment can be attributed to the neglect of dispersion shifts in the
theory. We have now developed a calculational technique to treat these
dispersion interactions within a QM/MM framework.[297]
Figure 22. Calculated absorption spectrum (left) and fluorescence spectrum (right) of
zwitterionic glycyltryptophan. The absorption and emission data were smoothed by
averaging 50 (absorption) and 20 (emission) adjacent data points.
79
Figure 23. Experimental absorption (left) and fluorescence spectrum (right) of Trp in
aqueous solution (pH7). (Adapted from [221])
In order to simulate the time-resolved fluorescence decay of glycyltryptophan
in water, we analyzed the distribution of Einstein coefficients Anm (Figure 24)
using Eq. (29). The decay curve thus obtained was fitted with one, two and
three exponentials, according to Eq. (30). A reduced χ2 value of 6.51 and a
squared correlation coefficient r2 of 0.96 indicate that one exponential function
is not sufficient to model the fluorescence decay of Gly-Trp. A good fit is
obtained if biexponential behavior is assumed. Although the reduced χ2 value of
0.15 indicates that the number of data points is not ideal, the r2 value (1.00) and
the visual inspection of the fit (Figure 25) suggest that in our simulation the Gly-
Trp decay is biexponential with fluorescence lifetimes of 4.81·10-7 and 1.42·10-7
seconds. The small value of 2rχ , which would ideally be 1, results from the
limited number of data points in our model study. As model calculations had
shown that the number of data points required to obtain a 2rχ value near 1
would have been much higher than the number of structures we used
(especially in the FRET simulations), we did not increase the number of
snapshots in our simulations. Nevertheless, the quality of the fit is good. The
amplitudes of the two components are 226 and 1686. The ratio of the two
lifetimes (3.39 : 1) agrees astonishingly well with the experimental value of
4.34 : 1 (1.26 ns; 0.29 ns),[240] although the absolute values are off by two
orders of magnitude. Visual inspection of the fitting curve, and the fact that
fitting with three exponentials gives too small values for 2rχ ( 2
rχ = 0.00;
80
r2 = 1.00) indicates that fitting with two exponentials is sufficient for a description
of the fluorescence decay.
Figure 24. Distribution of Einstein coefficients of spontaneous fluorescence emission,
calculated for zwitterionic glycyltryptophan in water.
Figure 25. Simulated fluorescence decay for zwitterionic glycyltryptophan in water. The
decay curve was fitted using two exponentials. Upper part of diagram: grey: simulated
fluorescence decay, black: fitting curve)
81
3.3.2 Simulating FRET in TetR: Molecular Dynamics Simulations
Unlike gas-phase simulations and water-box simulations without using the
PME method (not shown), our PME MD simulations of the TetR/Tc complex
gave a remarkably stable protein structure over the whole MD trajectory. The
RMS deviation calculated for all Cα atoms after a Cα-superimposition of the
average structure of the MD production phase (7,094.3 ps) and the minimized
X-ray structure is 1.60 Å. The superimposition of the average structure of the
MD production phase and the minimized X-ray structure (Figure 26) shows that
differences are restricted to very few residues at the amino and carboxylic
termini of the monomers and in some loop regions (see below). The RMS
calculated for all Cα atoms of monomer 1 after superimposition of all Cα atoms
of both monomers was 1.46 Å (1.44 Å if the Cα atoms of monomer 1 were
superimposed), the corresponding value for monomer 2 was 1.74 Å (1.73 Å if
the Cα atoms of monomer 2 were superimposed). This suggests that monomer
2 undergoes a small, but significant structural displacement in the MD
simulation (see below).
82
Figure 26. Superimposition (Cα fit, RMSD 1.60 Å) of the average structure of the MD
production phase (7,094.3 ps, colors: red: large RMSD value, blue: small RMSD value)
and the minimized X-ray structure (CPK colors). Side chain and hydrogen atoms were
omitted for clarity (with exception of Trp43/Trp43’). Monomer 1 is in the upper part of
the image, monomer 2 in the lower part.
In order to obtain information about flexible regions in the MD, we analyzed
the RMS fluctuations of the two monomers over the MD production phase
(Figure 27, Figure 28). Fluctuations with respect to the minimized X-ray
structure were recorded using Ptraj 6.5[293] for the Cα atoms and for the heavy
atoms of each residue (in the latter case average, mass-weighted fluctuations
were obtained). Analysis of the fluctuations shows a high mobility of the N- and
C-terminal residues and of the loop regions. Especially the unresolved loop in
the X-ray structure (residues 156-164 and 156’-164’) shows high fluctuations. It
seems remarkable that also the helices α2/α3 and α2’/α3’ in the DNA binding
domain (helices α1/α2/α3 and α1’/α2’/α3’) show a relatively high mobility.
83
Figure 27. Fluctuations with respect to the minimized X-ray structure for the C� atoms
(green) and for the heavy atoms of each residue (red, mass-weighted). Data shown for
monomer 1. Helical regions are indicated with a blue line on the x-axis. Fluctuations
were recorded using Ptraj 6.5.[293] Please note that residue 1 (Ala1) in our numbering is
Ala2 in the original PDB entry 2trt.[210]
Figure 28. Fluctuations with respect to the minimized X-ray structure for the Cα atoms
(green) and for the heavy atoms of each residue (red, mass-weighted). Data shown for
monomer 2. Helical regions are indicated with a blue line on the x-axis. Fluctuations
were recorded using Ptraj 6.5.[293] Please note that residue 1 (Ala1) in our numbering is
Ala2 in the original PDB entry 2trt.[210]
84
As we were mainly interested in the parameters that influence FRET from
Trp43 to Tc, we analyzed the χ1/ χ2 side chain angle distribution of Trp43 and
Trp43’ and the donor-acceptor distance between Trp42/Trp43’ and the
tetracyclines in each monomer.
Our water box simulations reveal flips to different Trp43 χ1/χ2 mean values in
addition to fast fluctuations in the side-chain angles (Figure 29). Due to the fast
fluctuations, using the term “rotamer” to describe the major Trp43 side chain
angle flips is somehow problematic, as each “rotamer” resembles an ensemble
of structures with χ1/χ2 fluctuating around a certain mean value. Nevertheless,
we use “rotamer” to address the structural classes of Trp generated by major
χ1/χ2 flips in addition to fast fluctuations of the side-chain angles.
A cluster analysis of the χ1 and χ2 side chain angles with TSAR 3.3[296]
indicated 2 clusters for monomer 1 and 5 clusters for monomer 2. The
structures representing the cluster centers for both monomers are shown in
Figure 30. In contrast to the indole chromophores of Trp 43/Trp43’, which
moved in the surrounding solvent (variation of χ1/χ2 angles), tetracycline did not
change its conformation and binding mode over the whole trajectory.
The analysis of the donor-acceptor distance (given by the distance between
the geometric center of the heavy atoms in the indole chromophore (donor) and
the geometric center of the heavy atoms in the BCD chromophore of
tetracycline (acceptor)) in monomer 1 shows an almost Gauss-shaped distance
distribution centered at 33.22 Å, whereas in monomer 2 the distance between
Trp43’ and Tc exhibits a larger variation (Figure 31). One reason for the large
variation of the donor-acceptor distance in monomer 2 is the presence of more
Trp43 rotamers in monomer 2 than in monomer 1 (because the geometric
centers of the chromophores are analyzed for FRET and not the Cα positions in
the case of Trp). In order to eliminate this influence, we also analyzed the
distances between the Trp43/Trp43’ Cα atoms and the C1A atom of the
corresponding tetracycline (not shown). These data (and the RMSD data, as
explained above) show that monomer 1 stays almost perfectly in its X-ray
conformation while in monomer 2 the DNA binding domain (helices α1’/α2’/α3’)
undergoes a small, but significant structural displacement of the protein
85
backbone. The secondary and tertiary structure of the protein, however, is not
modified in this process.
Figure 29. χ1/ χ2 side-chain angle distribution of Trp43 (monomer 1, left) and Trp43’
(monomer 2, right). Black: χ1, red: χ2. Side-chain angle clusters (Figure 30) are marked
by the colored lines at the bottom of the diagrams.
Figure 30. Trp43 structures representing the cluster centers for monomer 1 (left) and
monomer 2 (right) obtained from a cluster analysis of the χ1/ χ2 side chain angles. Left
(mon. 1): red (χ1 = -61.29°, χ2 = -209.27°), orange (χ1 = -59.14°, χ2 = -63.55°). Right
(mon. 2): red (χ1 = -179.07°, χ2 = -155.72°), orange (χ1 = -176.78°, χ2 = -110.76°), blue
(χ1 = -250.24°, χ2 = -74.15°), green (χ1 = -172.38°, χ2 = -288.17°), yellow (χ1 = -62.59°,
χ2 = -52.07°).
86
Figure 31. Distribution of distances between the geometric center of the heavy atoms
in the indole chromophore (donor) and the geometric center of the heavy atoms in the
BCD chromophore of tetracycline (acceptor). Left: monomer 1, right: monomer 2.
Because of the strong donor-acceptor distance-dependence of FRET, we
used only monomer 1 for the simulation of decay curves, because the almost
ideal Gauss-distribution of donor-acceptor distances in this case allows us to
investigate the connection between Trp side chain conformation and
fluorescence quenching by FRET in the absence of other factors affecting
FRET. Additionally, a test system with two instead of five rotamers seems much
more appropriate to understand experiments that give two lifetimes for Trp43 in
the TetR/Tc complex.
3.3.3 Simulating FRET in TetR: Quantum Mechanical Simulations
The transition dipoles calculated for the tryptophan S1 → S0 transition are
shown in Figure 32 for all MD snapshots. The directions of the transition dipoles
indicate that tryptophan emits from its 1La state. Almost no transition dipole is
found in the area where 1Lb transition dipoles would be expected, which would
be perpendicular to those of 1La. All transition dipoles are in-plane with the
indole ring system. The ring distortion of the indole chromophores in the “hot
87
geometries” results in the considerable distribution of Trp transition dipoles
around the ideal Trp 1La transition.
Figure 32. Calculated Trp transition dipoles. Left: All dipoles of the complete MD
trajectory after fitting the indole chromophore of each snapshot on the first structure in
the MD production phase. Middle, left: Only every 10th dipole shown. Middle, right:
Every 10th dipole, side-view. Right: Structure representing the cluster center (snapshot
at 3674 ps). Graphics were generated with Jmol 5.[295]
The molecular orbitals involved in the Trp S1 ↔ S0 transition are shown in
Figure 33 for the structure representing the cluster center of all Trp transition
dipoles (snapshot at 3674 ps). All orbitals affected by the CI excitations are
shown and the CI excitations are indicated by red arrows. All microstates
contributing to the Trp S1 ↔ S0 transition result from single excitations (although
pair double excitations were included in the CI). In the literature 1La is reported
to consist of a HOMO → LUMO and, to a lesser extent, of a HOMO-1 →
LUMO+1 transition, and 1Lb from a HOMO → LUMO+1 and HOMO-1 → LUMO
transition.[231] The analysis of the molecular orbitals of the “hot structure” that
represents the cluster center of all Trp S1 ↔ S0 transition dipoles in our
simulation shows that four π-π* microstates in the indole chromophore are
involved. Single excitations from both HOMO-2 and HOMO to LUMO and
LUMO+1 contribute to the Trp S1 ↔ S0 transition. HOMO-1 is not an orbital of
the indole ring system and is not involved in the CI, therefore our HOMO-2
corresponds to HOMO-1 in the literature.[231]
88
Figure 33. Molecular orbitals (left) involved in the Trp S1 ↔ S0 transition, calculated for
the Trp structure representing the cluster center of all snapshots (right, snapshot at
3674 ps). Energies and CI-coefficients of the CI microstates are printed near the
arrows indicating the CI excitation. HOMO-1 is not an orbital of the indole ring system
and is not involved in the CI.
The emission spectrum calculated for Trp43 in the protein/water box
environment is shown in Figure 34. Compared with the experimental spectrum,
there is a considerable blue-shift and the FWHM is still too small, although
larger than that calculated for the model system Gly-Trp.
89
Figure 34. Calculated fluorescence spectrum of Trp43, obtained from QM/MM
calculations in the protein/water box environment.
The time-resolved decays of the intrinsic fluorescence (i.e. without FRET) of
Trp43 in the protein environment/water box were calculated for the whole MD
trajectory (containing two rotamers) and for the two rotamers themselves. The
decay curves were calculated from the Einstein coefficients Anm (Figure 35,
Figure 36) of the snapshots of the MD trajectory, according to Eq. (29). They
were fitted with one, two and three exponentials, according to Eq. (30). We first
discuss the fluorescence decay obtained for the whole MD trajectory (Figure
37). A reduced χ2 value of 62.94 and a squared correlation coefficient r2 of 0.90
shows that one exponential function is not sufficient to model the intrinsic
fluorescence of Trp43. A good fit is obtained if biexponential behavior is
assumed. Although the reduced χ2 value of 0.42 is not ideal (and indicates that
the number of data points is too small, as discussed above), the r2 value (1.00)
and the visual inspection of the fit (Figure 37) suggest that the intrinsic Trp43
fluorescence decay is biexponential with fluorescence lifetimes of 5.33·10-7 and
8.62·10-8 seconds. The amplitudes of these two components are 654 and 6190.
The ratio of the two lifetimes (6.19 : 1) is larger than the experimental ratio
found for Trp 43 in TetR without inducer (2.34 : 1 for τ1 = 5.53 ns, τ2 = 2.36
ns[246] or 3.77 : 1 for τ1 = 4.9 ns, τ2 = 1.3 ns[249]) or other single tryptophan
proteins with solvent-exposed tryptophan, like glucagon (3.27 : 1; τ1 = 3.6 ns, τ2
90
= 1.1 ns).[221] As in our case, where two Trp rotamers are present, these
experimental studies of Trp43 in non-induced TetR assume two “classes”, local
environments or conformational states of Trp43.[246, 249] After visual inspection of
the fitting curve, we considered the fit with two exponentials to be sufficient for a
description of the fluorescence decay, because fitting with three exponentials
results in an abnormally small value of 2rχ ( 2
rχ = 0.00; r2 = 1.00).
Figure 35. Distribution of Einstein coefficients of spontaneous fluorescence emission,
calculated for Trp43 in TetR. Data represents the whole trajectory (2 rotamers).
If we assume that the detected radiation does not converge to zero for t → ∞,
i.e. that there is some background noise present, we can fit the intrinsic Trp43
decay curve using an offset, as shown in Eq. (31). Although there is no physical
basis for such an assumption, we checked this fitting option because it is also
used to fit experimental data. If a background noise of 114.64 relative units is
assumed, the ratio of the two lifetimes (2.88 : 1; lifetimes 2.07·10-7 and 7.17·10-8
seconds) is very similar to the ratios observed experimentally (see above).[246,
249] We thus consider our results to be consistent with the experiment within the
resolution of the fitting procedure.
91
Figure 36. Distribution of Einstein coefficients of spontaneous fluorescence emission,
calculated for Trp43 in TetR. Data represents rotamer 1 (left) and rotamer 2 (right) of
Trp43.
Figure 37. Simulated fluorescence decay for Trp43. Data represents both rotamers of
Trp43. The decay curve was fitted using two exponentials. Upper part of diagram: grey:
simulated fluorescence decay, black: fitting curve)
92
Two experimental studies report three lifetimes for Trp43 in TetR without
inducer. These are explained by assuming three “species” or rotamers of
Trp.[244, 248] The ratios found for the lifetimes (~ 30 : 12 : 1) disagree with the
ratio obtained from a triexponential fit of our simulated data (15.98 : 2.20 : 1).
However, if the short-lived component in these experimental studies is assumed
to be an instrumental artifact and is neglected, a ratio of ~ 2.5 :1 is found, which
agrees with our simulated data.
In order to investigate whether decay curves individually calculated for the
two rotamers are mono- or biexponential, we fitted these curves with one, two
and three exponential functions. Most interestingly, one exponential function
cannot model the decay of the single rotamers satisfactorily, as shown by large 2rχ values and small correlation coefficients (rotamer 1, snapshots 1005-
5350: 2rχ = 39.56, r2 = 0.89; rotamer 2, snapshots 1-1004, 5351-7093: 2
rχ =
21.96, r2 = 0.91). Good fits are obtained if two exponentials are used, although
the reduced χ2 values indicate a small number of data points (rotamer 1: 2rχ =
0.24, r2 = 1.00; rotamer 2: 2rχ = 0.17, r2 = 1.00). Three exponentials are
definitely not justifiable (rotamer 1: 2rχ = 0.00, r2 = 1.00; rotamer 2: 2
rχ = 0.00, r2
= 1.00).
93
Figure 38. Biexponential fits to the calculated Trp43 fluorescence decays. Red: whole
MD trajectory (contains both rotamers of Trp43), green: rotamer 1, blue: rotamer 2 of
Trp43.
The fluorescence lifetimes and their ratios found for the rotamers (rotamer 1:
ratio 7.89 : 1, lifetimes 6.25·10-7 and 7.92·10-8 seconds; rotamer 2: ratio 5.23 : 1,
lifetimes 5.36·10-7 and 1.03·10-7 seconds) are similar to those obtained if the
whole MD trajectory (2 rotamers) is analyzed. A comparison of the
biexponential fitting curves (Figure 38) shows that the lifetimes of the decay
curves calculated for the whole MD trajectory and the two rotamers do not differ
significantly. This agrees well with the observation that the distributions of
Einstein coefficients of the whole trajectory (2 rotamers, Figure 35) and of the
rotamers themselves (Figure 36) have almost identical shapes. Therefore, the
existence of rotamers, which have been defined as major conformational flips in
addition to fast, continuous fluctuations of the Trp χ1/χ2 angles in this work, does
not have a significant influence on the fluorescence decay of Trp43 in the TetR
protein in the absence of a quencher. The distribution of Einstein coefficients
(Figure 35 and Figure 36), which can be transformed into a distribution of
94
lifetimes, results from a large number if individual tryptophan geometries with
different microenvironments, each of these geometries being one snapshot of
the MD trajectory. Alcala, Gratton and Prendergast have shown that an
ensemble consisting of a large number of conformations with an interconversion
rate of the same order of magnitude as the excited-state rate can be analyzed
using a continuous distribution of lifetime values.[256] Multiexponential fits as
performed above are only one possible approximation among others in such
cases.[221, 256] Therefore, the biexponential model applied above to the intrinsic
fluorescence of Trp43 in TetR is only an approximate view of a much more
complex situation, which arises from the fast, continuous fluctuations of the
χ1/χ2 angles of tryptophan in our simulation.
For each snapshot, we selected the Tc transition with the largest oscillator
strength from the energetically accessible transitions for calculating the
orientation factors (see computational details section). This is justifiable
because the energetically accessible transitions not selected for κ2 calculation
were randomly orientated and had much smaller transition dipoles and
consequently, also much lower oscillator strengths. The most intense transition
dipole lies in the plane of the BCD chromophore of tetracycline for each
snapshot (Figure 39). Analysis of the molecular orbitals revealed that π-π*
transitions in the aromatic system of ring D, the enolate system of ring B and C
and the conjugated system of ring A were involved in this transition.
95
Figure 39. Transition dipole representing the electronic transition with the maximum
oscillator strength of all energetically accessible tetracycline transitions in a typical
snapshot.
The absorption spectrum of tetracycline in the protein binding site calculated
from the QM/MM results (Figure 40) agrees very well with the experimental
spectrum of [TcMg]+ in water (Figure 41).
Figure 40. Calculated absorption spectrum of [TcMg]+ in the protein binding site
(QM/MM calculations).
96
Figure 41. Experimental absorption spectrum of [TcMg]+/[TcCa]+ in water (taken from [201]).
The calculated Trp emission and Tc absorption spectra overlap very well, so
that an important requirement for FRET is fulfilled. Although the absolute peak
position and FWHM of the simulated Trp emission are still not perfect (see
above), our model is not affected by this problem as it only uses the orientation
factor and donor-acceptor distance as parameters for FRET.
3.3.4 Simulating FRET in TetR: Tryptophan Quenching by FRET
In order to investigate how the presence of two Trp rotamers affects FRET
quenching of Trp43, we analyzed the FRET data calculated for the whole MD
trajectory (two rotamers) and compared these with the data calculated for the
two individual rotamers.
In most studies of FRET, it is assumed that donor and acceptor can undergo
unrestricted isotropic motion. In these cases the mean value for κ2 averaged
over the ensemble is 2/3. However, an isotropic distribution of the donor and
acceptor moieties is in most cases only a first approximation that represents the
real situation more or less correctly. “For those who have not enjoyed the
practice of FRET spectroscopy, the numerical value of κ2 is at best a nuisance
and at worst an insurmountable problem.”[224] Our approach allows the
calculation of κ2 for every snapshot of our simulation. For all snapshots of the
97
MD trajectory (containing two rotamers), the mean value of κ2 was 0.37,
indicating that donor and acceptor cannot rotate freely in the solvent. This is
plausible because donor and acceptor are embedded in the protein and can
therefore adopt only a limited number of orientations relative to each other. The
calculation of κ2 mean values for rotamer 1 (snapshots 1005-5350) gives a
lower (0.28) and for rotamer 2 a higher value (0.52) than that obtained for the
whole trajectory.
The histogram of the orientation factors of the whole trajectory is a
superimposition of the κ2 histograms of the two rotamers (Figure 42).
In the next step, we used Eq. (43) in order to calculate relative FRET rate
constants for each snapshot. Modulation with r-6 does not change the form of
the histograms significantly. Again, the histogram of the relative FRET rate
constants of all snapshots is a superimposition of the FRET rate constant
histograms of the two rotamers (Figure 43). Rotamer 1 has smaller and rotamer
2 has larger relative FRET rate constants.
98
Figure 42. Calculated orientation factors. Sum of histogram bins represents the whole
MD trajectory (snapshots 1-7093) with both rotamers of Trp43 (κ2 mean value: 0.37).
Green: rotamer 1 (snapshots 1005-5350, κ2 mean value: 0.28), blue: rotamer 2
(snapshots 1-1004 and 5351-7093, κ2 mean value: 0.52) of Trp43.
Figure 43. Calculated relative FRET rate constants. Sum of histogram bins represents
the whole MD trajectory (snapshots 1-7093) with both rotamers of Trp43. Green:
rotamer 1 (snapshots 1005-5350), blue: rotamer 2 (snapshots 1-1004 and 5351-7093)
of Trp43.
99
Decay curves for all MD snapshots (Figure 44) and the two rotamers (Figure
45) were obtained using Eq. (44). The decay curves were fitted with one, two
and three exponential functions, according to Eq. (30). Additionally, we also
used Eq. (31), which allows an additional offset in the exponential fitting.
A fit with one exponential function is definitely not sufficient for the decay
curve calculated from all MD snapshots ( 2rχ = 77.18, r2 = 0.75). A correlation
coefficient of 0.99 and the visual inspection of the fitting curve and the residuals
convinced us that a fit with two exponentials is adequate to model the decay
curve, although the reduced χ2 value of 1.81 is not ideal. The weighted
residuals are very small. The form of the residual curve, however, indicates that
the exponential model is not absolutely correct for this fitting problem of
simulated data. Fitting with three exponentials gives a very small reduced χ2
value (0.02) and does not improve the correlation coefficient significantly (1.00).
We therefore consider two exponentials to be sufficient to model the decay
curve computed from all MD snapshots.
The relative Trp fluorescence lifetimes obtained from the biexponential fit for
the whole MD trajectory are 7.56 and 0.45 with amplitudes of 1892 and 4477.
The ratio of these two lifetimes (16.93 : 1) is larger than the value for the
experimental analogue, the W75F TetR/Tc complex (ratio 6.07 : 1), which is
also reported to be biexponential.[246] A lifetime ratio much closer to the
experimental results (7.36 : 1; relative lifetimes 2.16, 0.29) is obtained if an
additional constant (741.82) is allowed in the exponential fitting (according to
Eq. (31)).
100
Figure 44. Calculated Trp fluorescence decay. Only S1 quenching by FRET is
considered. Upper part: grey: calculated intensities, red line: fit with 2 exponentials.
Data represents both rotamers of Trp43.
Figure 45. Biexponential fits to the calculated Trp fluorescence decays. Only S1
quenching by FRET is considered. Red: whole MD trajectory (contains both rotamers
of Trp43), green: rotamer 1, blue: rotamer 2 of Trp43.
101
As in the case of the whole trajectory, the decay curves calculated for the two
rotamers cannot be fitted satisfactorily with one exponential function (rotamer 1,
snapshots 1005-5350: 2rχ = 38.19, r2 = 0.79; rotamer 2, snapshots 1-1004,
5351-7093: 2rχ = 44.33, r2 = 0.64). Two exponentials give good fits, as judged
from the reduced χ2 values (rotamer 1: 2rχ = 1.06; rotamer 2: 2
rχ = 0.87) and
the correlation coefficients (rotamer 1: r2 = 0.99; rotamer 2: r2 = 0.99). The
biexponential fits of the decay curves calculated for the two rotamers and the
whole MD trajectory (containing two rotamers) are shown in Figure 45. Fitting
with three exponentials gives very small reduced χ2 values (rotamer 1: 0.02;
rotamer 2: 0.01) and does not improve the correlation coefficients significantly
(1.00 and 1.00). We therefore consider two exponentials to be sufficient to
model the decay curves computed for the two rotamers.
The ratio of the two relative lifetimes is 16.29 : 1 for rotamer 1 (snapshots
1005-5350; relative lifetimes 9.97, 0.61) and 15.22 : 1 for rotamer 2 (snapshots
1-1004, 5351-7093; relative lifetimes 5.36, 0.35). Both ratios are higher than in
the experiment (which, of course, does not discriminate between rotamers), and
become more similar to the experimental value if an additional constant (Eq.
(31)) is allowed in the fit (rotamer 1: 7.66 : 1, rotamer 2: 6.28 : 1).
The fitting analysis of the FRET decay curves both of the whole MD trajectory
(containing two rotamers) and of the rotamers themselves requires two
exponential functions. Therefore, two lifetimes are found, even if only one
tryptophan ‘rotamer’ (which of course represents numerous conformations
interchanging rapidly by variation of the tryptophan χ1 and χ2 side chain angles,
as shown in the MD results section) is analyzed. This observation disagrees
with the classical rotamer model, which assumes that a specific tryptophan
conformation is responsible for the existence of a corresponding lifetime. In
contrast to the classical model, the numerous relative orientations of the donor
and acceptor dipoles that are caused by continuous χ1 and χ2 side-chain angle
fluctuations are the reason for the biexponential decay curves. This is plausible
because FRET is extremely sensitive to the relative orientation of the
donor/acceptor transition dipoles. Another factor that complicates the situation
102
is that even a perfect, monomodal, Gauss-shaped donor/acceptor distance
distribution may lead to complex, nonexponential decay curves.[221]
The short-lived lifetime components of the individual rotamers, 0.61 for
rotamer 1 (0.36 if Eq. (31) is used) and 0.35 for rotamer 2 (0.26 if Eq. (31) is
used) do not differ significantly from the value found for the whole trajectory
(0.45 or 0.29 if Eq. (31) is applied), although their values follow the same trend
as those of the long-lived lifetime component (see below). This short-lived
component probably results from the continuous fluctuations of the Trp43 side
chain angles, and is not analyzed further.
As shown above (Figure 42-Figure 43), rotamer 1 has smaller κ2 values
(mean value: 0.28) and, consequently, smaller relative FRET rate constants
than rotamer 2 (κ2 mean value: 0.52). Therefore, the long-lived lifetime
component for rotamer 1 (9.97 or 2.74, if Eq. (31) is used) is larger than that of
rotamer 2 (5.36 or 1.65), while the long-lived lifetime calculated for the whole
trajectory is between these values (7.56 or 2.16). The long-lived lifetime
component contains the information about the ‘rotamers’ of Trp43, which result
from major conformational flips of the tryptophan χ1 and χ2 side chain angles,
as discussed in the MD results section. The influence of the Trp ‘rotamer’ on the
long-lived Trp lifetime component shows that the idea of the classical rotamer
model is not completely wrong, it just oversimplifies the real situation (at least in
the case of the TetR/Tc system), where a combination of fast, continuous
tryptophan χ1 and χ2 side-chain angle fluctuations and additional flips of the
side-chain conformation leads to a complex situation. Computer simulations like
the present work provide a valuable tool for investigating such complex data.
3.4 Conclusions
Our results show that a combination of a classical force-field molecular
dynamics simulation and subsequent quantum mechanical single point
calculations of the “hot” snapshots provided by the MD simulation can be used
to simulate protein emission and absorption spectra and is also capable of
simulating complex phenomena like FRET.
103
Our simulation protocol enables us to examine different tryptophan rotamers
in a protein. As we are able to investigate distinct rotamers and their properties
independently from each other, we can deduce the influence of individual
rotamers on spectroscopic properties instead of analyzing the whole ensemble
consisting of all rotamers, which is the case in most experiments.
The analysis of the intrinsic fluorescence (without FRET quenching) of the
Trp43 residue located in the DNA binding domain of the TetR protein shows a
biexponential fluorescence decay. Most interestingly, this biexponentiality is
also observed if the two rotamers of Trp43 are analyzed independently.
Therefore, the intrinsic, unquenched fluorescence of Trp43 is – independent of
the existence of rotamers – biexponential. However, this “biexponentiality” is
enhanced by FRET quenching (see below). The fast Trp side-chain angle
fluctuations, which are observed in addition to infrequent, major flips of the Trp
χ1 and χ2 angles, lead to a large number of different tryptophan conformations
with distinctive microenvironments, which would be best described by a
distribution of lifetimes, as described by Alcala et al.[256] Thus, the exponential
model is an approximation.
The analysis of FRET quenching in the TetR/tetracycline complex also gives
biexponential decays, even if the rotamers are analyzed independently. The
biexponentiality of the FRET decay curves results from a large ensemble of
snapshots with different Trp side-chain angles and, therefore, different relative
donor-acceptor orientations. Additionally, the Gauss-shaped donor-acceptor
distance distribution enhances the nonexponential character of the decay
curves.
The individual rotamers have different probabilities for Förster energy transfer
(FRET). Thus, the relative lifetimes of the two rotamers are influenced in a
different way by FRET. Therefore, the classical rotamer model, which explains
multiexponential tryptophan fluorescence decays as a result of the existence of
multiple conformations, is not completely wrong but it oversimplifies the real
situation, at least in the case of the tetracycline/TetR complex. Our study shows
that a combination of fast Trp side chain fluctuations with additional
conformational flips of the Trp side chain results in the complex decay spectra
found experimentally. In contrast to other experimental and computational
104
studies, which do not calculate the orientation factor κ2 explicitly, our work gives
an accurate description of FRET and shows that future interpretations of FRET
data require a detailed view based on accurate molecular data. This result is of
immense importance for spectroscopic studies on enzymes and must result in
the reinterpretation of many measurements. Instead of a simple qualitative
interpretation of fluorescence decay curves, very extensive MD-QM/MM-CI
simulations, which are extremely compute-intensive, are required. Our work
emphasizes that in this area only the combination of simulations validated by
comparison with experimental results (or vice versa, according to the point of
view) can provide information about protein conformations and dynamics.
3.5 Outlook
In this initial study, we have presented only relative FRET rate constants. In
order to obtain a complete kinetic model for Trp S1 deactivation, we need
absolute FRET rate constants, which require calculation of the spectroscopic
donor-acceptor overlap and the donor quantum yield calculated for each
snapshot. A complete model also contains the intrinsic Trp fluorescence and all
non-radiative pathways. Future studies also require the parameterization of a
force field for the tryptophan S1 state and the development of a simulation
protocol that switches to excited Trp during the MD simulation (as in the work of
Callis and Vivian[232]), which should be used instead of the ground-state MD
used in this work to give more accurate charges and geometries for the excited
Trp chromophores.
105
4 References
[1] J. A. Pople, in Encyclopedia of Computational Chemistry (Eds.: P. v. R.
Schleyer, N. L. Allinger, T. Clark, J. Gasteiger, P. A. Kollman, H. F.
Schaefer III, P. R. Schreiner), Wiley, Chichester, 1998.
[2] P. v. R. Schleyer, in Encyclopedia of Computational Chemistry (Eds.: P.
v. R. Schleyer, N. L. Allinger, T. Clark, J. Gasteiger, P. A. Kollman, H. F.
Schaefer III, P. R. Schreiner), Wiley, Chichester, 1998.
[3] C. J. Cramer, Essentials of Computational Chemistry, John Wiley &
Sons, West Sussex, 2002.
[4] A. R. Leach, Molecular Modelling: Principles and Applications, 2nd ed.,
Prentice-Hall, Harlow, 2001.
[5] T. Clark, Anwendung von semiempirischen MO-Methoden in der
Organischen Chemie, lecture, Universität Erlangen-Nürnberg, Erlangen,
Germany, 2004.
[6] Science Citiation Index Expanded web site:
http://isi15.isiknowledge.com/portal.cgi/portal.cgi?DestApp=WOS&Func=
Frame
[7] Nobel Foundation web site: http://nobelprize.org/chemistry/
[8] L. Terfloth, in Chemoinformatics - A Textbook (Eds.: J. Gasteiger, T.
Engel), WILEY-VCH, Weinheim, 2003, pp. 597-622.
[9] P. Tollman, P. Guy, J. Altshuler, A. Flanagan, M. Steiner, A Revolution in
R&D, The Boston Consulting Group, Boston, MA, 2001.
[10] Boston Consulting Group web site: http://www.bcg.com
[11] W. A. Warr, in Handbook of Chemoinformatics (Ed.: J. Gasteiger),
WILEY-VCH, Weinheim, 2003, pp. 1604-1639.
[12] T. Clark, in Chemoinformatics - A Textbook (Eds.: J. Gasteiger, T.
Engel), WILEY-VCH, Weinheim, 2003, pp. 376-397.
[13] F. Jensen, Introduction to Computational Chemistry, John Wiley & Sons,
Chichester, 1998.
[14] T. Herz, Ph.D. thesis, Universität Erlangen-Nürnberg, Erlangen,
Germany, 1999.
106
[15] G. Wedler, Lehrbuch der Physikalischen Chemie, 3rd ed., VCH,
Weinheim, New York, 1987.
[16] T. Clark, A Handbook of Computational Chemistry, Wiley & Sons, New
York, 1985.
[17] T. Clark, in Handbook of Chemoinformatics (Ed.: J. Gasteiger), WILEY-
VCH, Weinheim, 2003, pp. 947-975.
[18] E. Schrödinger, Ann. Physik 1926, 79, 489.
[19] M. Born, J. R. Oppenheimer, Ann. Physik 1927, 84, 457-484.
[20] V. Fock, Z. Physik 1930, 61, 126.
[21] J. C. Slater, Phys. Rev. 1930, 36, 57-64.
[22] C. C. J. Roothaan, Rev. Mod. Phys. 1951, 23, 69.
[23] G. G. Hall, Proc. Roy. Soc. (London) 1951, A205, 541.
[24] J. A. Pople, R. K. Nesbet, J. Chem. Phys. 1954, 22, 571-572.
[25] A. Szabo, N. S. Ostlund, Modern Quantum Chemistry, McGraw Hill
Publishing Company, New York, 1989.
[26] J. A. Pople, D. P. Santry, G. A. Segal, J. Chem. Phys. 1965, 43, S129-
S135.
[27] J. A. Pople, G. A. Segal, J. Chem. Phys. 1965, 43, S136-S149.
[28] M. J. S. Dewar, J. Am. Chem. Soc. 1952, 74, 3341.
[29] J. A. Pople, Trans. Farad. Soc. 1953, 49, 1375.
[30] R. Pariser, R. G. Parr, J. Chem. Phys. 1953, 21, 466, 767.
[31] J. A. Pople, D. L. Beveridge, P. A. Dobosh, J. Chem. Phys. 1967, 47,
2026-2033.
[32] R. C. Bingham, M. J. S. Dewar, D. H. Lo, J. Am. Chem. Soc. 1975, 97,
1285-1310.
[33] M. J. S. Dewar, W. Thiel, J. Am. Chem. Soc. 1977, 99, 4899-4917.
[34] J. J. P. Stewart, Prog. 455, Quantum Chem. Prog. Exch., 1983.
[35] Quantum Chemistry Program Exchange web site:
http://qcpe.chem.indiana.edu/
[36] M. J. S. Dewar, E. G. Zoebisch, E. F. Healy, J. J. P. Stewart, J. Am.
Chem. Soc. 1985, 107, 3902-3909.
[37] J. J. P. Stewart, J. Comput. Chem. 1989, 10, 209-264.
[38] M. J. S. Dewar, C. Lie, J. Yu, Tetrahedron 1993, 49, 5003-5038.
[39] A. J. Holder, R. D. Dennington II, C. Lie, Tetrahedron 1994, 50, 627.
107
[40] W. Thiel, A. A. Voityuk, J. Mol. Struct. (THEOCHEM) 1994, 313, 141-
154.
[41] W. Thiel, A. A. Voityuk, J. Phys. Chem. 1996, 100, 616-626.
[42] W. Thiel, Adv. Chem. Phys. 1996, 93, 703.
[43] W. Thiel, A. A. Voityuk, Theor. Chim. Acta 1992, 81, 391-404.
[44] W. Thiel, A. A. Voityuk, Theor. Chim. Acta 1996, 93, 315-315.
[45] W. Thiel, A. A. Voityuk, Int. J. Quantum Chem. 1994, 44, 807.
[46] W. Thiel, in Encyclopedia of Computational Chemistry, Vol. 3 (Eds.: P. v.
R. Schleyer, N. L. Allinger, T. Clark, J. Gasteiger, P. A. Kollman, H. F.
Schaefer III, P. R. Schreiner), Wiley, Chichester, 1998, pp. 1604-1606.
[47] P. Winget, A. H. C. Horn, C. Selçuki, B. Martin, T. Clark, J. Mol. Model.
2003, 9, 408-414.
[48] P. Winget, T. Clark, J. Mol. Model. 2005, in the press.
[49] A. A. Voityuk, N. Rösch, J. Phys. Chem. A 2000, 104, 4089-4094.
[50] T. Clark, A. Alex, B. Beck, F. Burkhardt, J. Chandrasekhar, P. Gedeck,
A. H. C. Horn, M. Hutter, B. Martin, G. Rauhut, W. Sauer, T. Schindler, T.
Steinke, VAMP 8.1, Computer-Chemie-Centrum, Universität Erlangen-
Nürnberg, Erlangen, Germany, 2002.
[51] MS Modeling 3.0.1, Accelrys Inc., San Diego, CA, 2003.
[52] MOPAC 2002, CAChe Group, Fujitsu America Inc., Beaverton, OR,
2002.
[53] AMPAC 8, Semichem Inc., Shawnee Mission, KS, 2004.
[54] D. N. Nanda, K. Jug, Theor. Chim. Acta 1980, 57, 95.
[55] J. Ridley, M. C. Zerner, Theor. Chim. Acta 1973, 32, 111-134.
[56] A. D. Bacon, M. C. Zerner, Theor. Chim. Acta 1979, 53, 21.
[57] W. D. Edwards, M. C. Zerner, Theor. Chim. Acta 1987, 72, 347.
[58] J. J. P. Stewart, P. Csaszar, P. Pulay, J. Comput. Chem. 1982, 3, 227-
228.
[59] W. Pan, T.-S. Lee, W. Yang, J. Comput. Chem. 1998, 19, 1101-1109.
[60] A. van der Vaart, V. Gogonea, S. L. Dixon, K. M. Merz Jr., J. Comput.
Chem. 2000, 21, 1494-1504.
[61] H. Lanig, in Chemoinformatics - A Textbook (Eds.: J. Gasteiger, T.
Engel), WILEY-VCH, Weinheim, 2003, pp. 338-358.
[62] N. L. Allinger, J. Am. Chem. Soc. 1977, 99, 8127-8134.
108
[63] N. L. Allinger, Y. H. Yuh, J.-J. Lii, J. Am. Chem. Soc. 1989, 111, 8551-
9556.
[64] S. J. Weiner, P. A. Kollman, D. A. Case, D. T. Nguyen, J. Comput.
Chem. 1986, 7, 230-252.
[65] S. J. Weiner, P. A. Kollman, D. A. Case, U. C. Singh, C. Ghio, G.
Alagona, S. Profeta Jr., P. Weiner, J. Am. Chem. Soc. 1984, 106, 765-
784.
[66] W. D. Cornell, P. Cieplak, C. I. Bayly, I. R. Gould, K. M. Merz Jr., D. M.
Ferguson, D. C. Spellmeyer, T. Fox, J. W. Caldwell, P. A. Kollman, J.
Am. Chem. Soc. 1995, 117, 5179-5197.
[67] H. Lanig, in Chemoinformatics - A Textbook (Eds.: J. Gasteiger, T.
Engel), WILEY-VCH, Weinheim, 2003, pp. 359-375.
[68] H. Lanig, Protein Modelling II: Computeranwendungen, lecture,
Universität Erlangen-Nürnberg, Erlangen, Germany, 2005.
[69] R. W. Hockney, Methods Comput. Phys. 1970, 9, 136-211.
[70] L. Verlet, Phys. Rev. 1967, 159, 98-103.
[71] W. F. van Gunsteren, H. J. C. Berendsen, Mol. Phys. 1977, 34, 1311-
1327.
[72] W. L. Jorgensen, J. Chandrasekhar, J. D. Madura, R. W. Impey, M. L.
Klein, J. Chem. Phys. 1983, 79, 926-935.
[73] H. J. C. Berendsen, J. P. M. Postma, W. F. van Gunsteren, A. Di Nola, J.
R. Haak, J. Chem. Phys. 1984, 81, 3684-3690.
[74] L. Pedersen, T. Darden, in Encyclopedia of Computational Chemistry,
Vol. 3 (Eds.: P. v. R. Schleyer, N. L. Allinger, T. Clark, J. Gasteiger, P. A.
Kollman, H. F. Schaefer III, P. R. Schreiner), Wiley, Chichester, 1998,
pp. 1650-1659.
[75] T. Darden, D. York, L. Pedersen, J. Chem. Phys. 1993, 98, 10089-
10092.
[76] P. Ewald, Ann. Physik 1921, 64, 253-287.
[77] D. M. York, A. Wlodawer, L. G. Pedersen, T. A. Darden, Proc. Natl.
Acad. Sci. U. S. A. 1994, 91, 8715-8718.
[78] T. E. Cheatham III, J. L. Miller, T. Fox, T. A. Darden, P. A. Kollman, J.
Am. Chem. Soc. 1995, 117, 4193-4194.
109
[79] M. Meyer, G. Wohlfahrt, J. Knäblein, D. Schomburg, J. Comput.-Aided
Mol. Des. 1998, 12, 425-440.
[80] A. Warshel, M. Levitt, J. Mol. Biol. 1976, 103, 227-249.
[81] G. Monard, K. M. Merz Jr., Acc. Chem. Res. 1999, 32, 904-911.
[82] J. Åqvist, A. Warshel, Chem. Rev. (Washington, DC, U. S.) 1993, 93,
2523-2544.
[83] G. Schürer, H. Lanig, T. Clark, J. Phys. Chem. B 2000, 104, 1349-1361.
[84] G. Schürer, A. H. C. Horn, P. Gedeck, T. Clark, J. Phys. Chem. B 2002,
106, 8815-8830.
[85] R. Castillo, J. Andrés, V. Moliner, J. Am. Chem. Soc. 1999, 121, 12140-
12147.
[86] G. Mlinsek, M. Novic, M. Hodoscek, T. Solmajer, J. Chem. Inf. Comput.
Sci. 2001, 41, 1286-1294.
[87] A. J. Mulholland, P. D. Lyne, M. Karplus, J. Am. Chem. Soc. 2000, 122,
534-535.
[88] L. Ridder, A. J. Mulholland, I. M. C. M. Rietjens, J. Vervoort, J. Am.
Chem. Soc. 2000, 122, 8728-8738.
[89] G. Schürer, Ph.D. thesis, Universität Erlangen-Nürnberg, Erlangen,
Germany, 2003.
[90] D. Bakowies, W. Thiel, J. Phys. Chem. 1996, 100, 10580-10594.
[91] P. L. Cummins, J. E. Gready, J. Comput. Chem. 1998, 19, 977-988.
[92] G. Ujaque, F. Maseras, A. Lledós, J. Am. Chem. Soc. 1999, 121, 1317-
1323.
[93] W. F. van Gunsteren, H. Liu, F. Müller-Plathe, THEOCHEM 1998, 432,
9-14.
[94] H. Hu, L. H., Y. Shi, Proteins: Struct., Funct., Genet. 1997, 27, 545-555.
[95] F. Beierlein, H. Lanig, G. Schürer, A. H. C. Horn, T. Clark, Mol. Phys.
2003, 101, 2469-2480.
[96] F. Beierlein, Diploma thesis, Universität Erlangen-Nürnberg, Erlangen,
Germany, 2001.
[97] F. Beierlein, T. Clark, Simulating the Induced Fit: A Combined QM/MM
Docking Approach, MGMS Young Modellers' Forum in Conjunction with
the RSC MMG, London, 2001.
110
[98] F. Beierlein, T. Clark, A QM/MM Docking Approach with a Flexible
Protein Environment, 16. Darmstädter Molecular Modelling Workshop,
Darmstadt, 2002.
[99] A. H. C. Horn, F. Beierlein, H. Lanig, G. Schürer, T. Clark, Quantum
mechanical/molecular mechanical (QM/MM) docking: A new protocol and
its evaluation for known test systems, 227th ACS National Meeting,
Anaheim, CA, March 28-April 1, 2004.
[100] B. Beck, T. Clark, in 3D QSAR in Drug Design, Vol. 2 (Eds.: H. Kubinyi,
G. Folkers, Y. C. Martin), Kluwer/Escom, Dordrecht, 1998, pp. 131-159.
[101] U. C. Singh, P. A. Kollman, J. Comput. Chem. 1986, 7, 718-730.
[102] M. J. Field, P. A. Bash, M. Karplus, J. Comput. Chem. 1990, 11, 700-
733.
[103] V. Théry, D. Rinaldi, J.-L. Rivail, B. Maigret, G. J. Ferenczy, J. Comput.
Chem. 1994, 15, 269-282.
[104] T. Clark, A. Alex, B. Beck, P. Gedeck, H. Lanig, J. Mol. Model. 1999, 5,
1-7.
[105] J. Gao, in Reviews in Computational Chemistry, Vol. 7, VCH Publishers,
New York, 1996, pp. 119-185.
[106] I. J. Enyedy, Y. Ling, K. Nacro, Y. Tomita, X. Wu, Y. Cao, R. Guo, B. Li,
X. Zhu, Y. Huang, Y.-Q. Long, P. P. Roller, D. Yang, S. Wang, J. Med.
Chem. 2001, 44, 4313-4324.
[107] S. Grüneberg, B. Wendt, G. Klebe, Angew. Chem. 2001, 113, 404-408.
[108] S. Grüneberg, B. Wendt, G. Klebe, Angew. Chem., Int. Ed. 2001, 40,
389-393.
[109] S. Grüneberg, M. T. Stubbs, G. Klebe, J. Med. Chem. 2002, 45, 3588-
3602.
[110] M. L. Lamb, K. W. Burdick, S. Toba, M. M. Young, A. G. Skillman, X.
Zou, J. R. Arnold, I. D. Kuntz, Proteins: Struct., Funct., Genet. 2001, 42,
296-318.
[111] AutoDock web page: http://www.scripps.edu/pub/olson-
web/doc/autodock
[112] FlexX web page: http://www.biosolveit.de/FlexX/
[113] M. Böhm, A. Mitsch, P. Wissner, I. Sattler, M. Schlitzer, J. Med. Chem.
2001, 44, 3117-3124.
111
[114] J. Drews, Science (Washington, DC, U. S.) 2000, 287, 1960-1964.
[115] A. H. Elcock, D. Sept, J. A. McCammon, J. Phys. Chem. B 2001, 105,
1504-1518.
[116] M. S. Lesney, Modern Drug Discovery 2001, 4(9), 34-41.
[117] D. E. Koshland, Proc. Natl. Acad. Sci. U. S. A. 1958, 44, 98-104.
[118] I. D. Kuntz, J. M. Blaney, S. J. Oatley, R. Langridge, T. E. Ferrin, J. Mol.
Biol. 1982, 161, 269-288.
[119] C. Bissantz, G. Folkers, D. Rognan, J. Med. Chem. 2000, 43, 4759-4767.
[120] R. C. Willis, Modern Drug Discovery 2001, 4(9), 26-32.
[121] T. J. A. Ewing, S. Makino, A. G. Skillman, I. D. Kuntz, J. Comput.-Aided
Mol. Des. 2001, 15, 411-428.
[122] E. Perola, K. Xu, T. M. Kollmeyer, S. H. Kaufmann, F. G. Prendergast, Y.
P. Pang, J. Med. Chem. 2000, 43, 401-408.
[123] M. Rarey, B. Kramer, T. Lengauer, G. Klebe, J. Mol. Biol. 1996, 261,
470-489.
[124] W. Welch, J. Ruppert, A. N. Jain, Chem. Biol. 1996, 3, 449-462.
[125] C. A. Baxter, C. W. Murray, D. E. Clark, D. R. Westhead, M. D. Eldridge,
Proteins: Struct., Funct., Genet. 1998, 33, 367-382.
[126] T. Hou, J. Wang, L. Chen, X. Xu, Protein Eng. 1999, 12, 639-647.
[127] G. Jones, P. Willett, R. C. Glen, A. R. Leach, R. Taylor, J. Mol. Biol.
1997, 267, 727-748.
[128] G. M. Morris, D. S. Goodsell, R. S. Halliday, R. Huey, W. E. Hart, R. K.
Belew, A. J. Olson, J. Comput. Chem. 1998, 19, 1639-1662.
[129] P. S. Charifson, J. J. Corkery, M. A. Murcko, W. P. Walters, J. Med.
Chem. 1999, 42, 5100-5109.
[130] D. K. Gelhaar, G. M. Verkhivker, P. A. Rejto, C. J. Sherman, D. B. Fogel,
L. J. Fogel, S. T. Freer, Chem. Biol. 1995, 2, 317-324.
[131] D. S. Goodsell, A. J. Olson, Proteins: Struct., Funct., Genet. 1990, 8,
195-202.
[132] G. M. Morris, D. S. Goodsell, R. Huey, A. J. Olson, J. Comput.-Aided
Mol. Des. 1996, 10, 293-304.
[133] M. Liu, S. Wang, J. Comput.-Aided Mol. Des. 1999, 13, 435-451.
[134] C. McMartin, R. S. Bohacek, J. Comput.-Aided Mol. Des. 1997, 11, 333-
344.
112
[135] Dockit, Metaphorics LLC, Mission Viejo, CA 92691.
[136] H. Claußen, C. Büning, M. Rarey, T. Lengauer, J. Mol. Biol. 2001, 308,
377-395.
[137] J. R. H. Tame, J. Comput.-Aided Mol. Des. 1999, 13, 99-108.
[138] H. Gohlke, G. Klebe, Curr. Opin. Struct. Biol. 2001, 11, 231-235.
[139] C. Pilger, C. Bartolucci, D. Lamba, A. Tropsha, G. Fels, J. Mol. Graphics
2001, 19, 288-296.
[140] S. W. Cho, S. Lee, W. Shin, J. Mol. Biol. 2001, 314, 863-878.
[141] C. A. Marhefka, B. M. Moore II, T. C. Bishop, L. Kirkovsky, A. Mukherjee,
J. T. Dalton, D. D. Miller, J. Med. Chem. 2001, 44, 1729-1740.
[142] L. Schaffer, G. M. Verkhivker, Proteins: Struct., Funct., Genet. 1998, 33,
295-310.
[143] E. Jenwitheesuk, R. Samudrala, BMC Struct. Biol. 2003, 3, 2.
[144] D. Hoffmann, B. Kramer, T. Washio, T. Steinmetzer, M. Rarey, T.
Lengauer, J. Med. Chem. 1999, 42, 4422-4433.
[145] J. Wang, P. A. Kollman, I. D. Kuntz, Proteins: Struct., Funct., Genet.
1999, 36, 1-19.
[146] J. Durup, Actual. Chim. Thér. 1993, 20, 309-326.
[147] R. C. Wade, S. Lüdemann, in Structure-Based Drug Design (Ed.: P. W.
Codding), Kluwer/Escom, Dordrecht, 1998, pp. 41-52.
[148] R. C. Wade, V. Sobolev, A. R. Ortiz, G. Peters, in Structure-Based Drug
Design (Ed.: P. W. Codding), Kluwer/Escom, Dordrecht, 1998, pp. 223-
232.
[149] P. Y. S. Lam, P. K. Jadhav, C. J. Eyermann, C. N. Hodge, Y. Ru, L. T.
Bacheler, J. L. Meek, M. J. Otto, M. M. Rayner, Y. N. Wong, C.-H.
Chang, P. C. Weber, D. A. Jackson, T. R. Sharpe, S. Erickson-Viitanen,
Science (Washington, DC, U. S.) 1994, 263, 380-384.
[150] S. Thaisrivongs, K. D. Watenpaugh, W. J. Howe, P. K. Tomich, L. A.
Dolak, K.-T. Chong, C.-S. C. Tomich, A. G. Tomasselli, S. R. Turner, J.
W. Strohbach, A. M. Mulichak, M. N. Janakiraman, J. B. Moon, J. C.
Lynn, M.-M. Horng, R. R. Hinshaw, K. A. Curry, D. J. Rothrock, J. Med.
Chem. 1995, 38, 3624-3637.
[151] H. I. Skulnick, P. D. Johnson, P. A. Aristoff, J. K. Morris, K. D. Lovasz, W.
J. Howe, K. D. Watenpaugh, M. N. Janakiraman, D. J. Anderson, R. J.
113
Reischer, T. M. Schwartz, L. S. Banitt, P. K. Tomich, J. C. Lynn, M.-M.
Horng, K.-T. Chong, R. R. Hinswaw, L. A. Dolak, E. P. Seest, F. J.
Schwende, B. D. Rush, G. M. Howard, L. N. Toth, K. R. Wilkinson, T. J.
Kakuk, C. W. Johnson, S. L. Cole, R. M. Zaya, G. L. Zipp, P. L. Possert,
R. J. Dalga, W.-Z. Zhong, M. G. Williams, K. R. Romines, J. Med. Chem.
1997, 40, 1149-1164.
[152] M. Marquart, J. Walter, J. Deisenhofer, W. Bode, R. Huber, Acta
Crystallogr., Sect. B 1983, 39, 480-490.
[153] E. Conti, C. Rivetti, A. Wonacott, P. Brick, FEBS Lett. 1998, 425, 229-
233.
[154] B. A. Katz, R. Mackman, C. Luong, K. Radika, A. Martelli, P. A.
Sprengeler, J. Wang, H. Chan, L. Wong, Chem. Biol. 2000, 7, 299-312.
[155] J. Wagner, J. Kallen, C. Ehrhardt, J.-P. Evenou, D. Wagner, J. Med.
Chem. 1998, 41, 3664-3674.
[156] D. W. Christianson, W. N. Lipscomb, Proc. Natl. Acad. Sci. U. S. A.
1986, 83, 7568-7572.
[157] J. T. Bolin, D. J. Filman, D. A. Matthews, R. C. Hamlin, J. Kraut, J. Biol.
Chem. 1982, 257, 13650-13662.
[158] M. R. Sawaya, J. Kraut, Biochemistry 1997, 36, 586-603.
[159] F. C. Bernstein, T. F. Koetzle, G. J. B. Williams, E. F. Meyer Jr., M. D.
Brice, J. R. Rodgers, O. Kennard, T. Shimanouchi, M. Tasumi, J. Mol.
Biol. 1977, 112, 535-542.
[160] H. M. Berman, J. Westbrook, Z. Feng, G. Gilliland, T. N. Bhat, H.
Weissig, I. N. Shindyalov, P. E. Bourne, Nucleic Acids Res. 2000, 28,
235-242.
[161] RCSB PDB web page: http://www.rcsb.org/pdb/
[162] R. Lapatto, T. Blundell, A. Hemmings, J. Overington, A. Wilderspin, S.
Wood, J. R. Merson, P. J. Whittle, D. E. Danley, K. F. Geoghegan, S. J.
Hawrylik, S. E. Lee, K. G. Scheld, P. M. Hobart, Nature (London, U. K.)
1989, 342, 299-302.
[163] T. Clark, A. Alex, B. Beck, J. Chandrasekhar, P. Gedeck, A. Horn, M.
Hutter, B. Martin, G. Rauhut, W. Sauer, T. Schindler, T. Steinke, VAMP
7.5 Build 18, Computer-Chemie-Centrum, Universität Erlangen-
Nürnberg, Erlangen, Germany, 2001.
114
[164] P. Gedeck, T. Clark, unpublished results.
[165] M. Clark, R. D. Cramer III, N. Van Opdenbosch, J. Comput. Chem. 1989,
10, 982-1012.
[166] A. K. Rappé, C. J. Casewit, K. S. Colwell, W. A. Goddard III, W. M. Stiff,
J. Am. Chem. Soc. 1992, 114, 10024-10034.
[167] D. A. Case, D. A. Pearlman, J. W. Caldwell, T. E. Cheatham III, W. S.
Ross, C. L. Simmerling, T. A. Darden, K. M. Merz, R. V. Stanton, A. L.
Cheng, J. J. Vincent, M. Crowley, D. M. Ferguson, R. J. Radmer, G. L.
Seibel, U. C. Singh, P. K. Weiner, P. A. Kollman, AMBER 5, University of
California, San Francisco, 1997.
[168] Sybyl 6.7, Tripos Associates, St. Louis, MO, 2000.
[169] D. C. Liu, J. Nocedal, Math. Programming 1985, 45, 503-528.
[170] J. Nocedal, Math. Comp. 1980, 24, 773-782.
[171] A. J. Holder, in Encyclopedia of Computational Chemistry, Vol. 1 (Eds.:
P. v. R. Schleyer, N. L. Allinger, T. Clark, J. Gasteiger, P. A. Kollman, H.
F. Schaefer III, P. R. Schreiner), Wiley, Chichester, 1998, pp. 8-11.
[172] K. C. Gross, P. G. Seybold, Int. J. Quantum Chem. 2000, 80, 1107-1115.
[173] T. A. Halgren, J. Am. Chem. Soc. 1992, 114, 7827-7843.
[174] K. M. Merz Jr., J. Am. Chem. Soc. 1991, 113, 406-411.
[175] S. C. Hoops, K. W. Anderson, K. M. Merz Jr., J. Am. Chem. Soc. 1991,
113, 8262-8270.
[176] J. Gasteiger, M. Marsili, Tetrahedron 1980, 36, 3219-3288.
[177] J. Gasteiger, H. Saller, Angew. Chem. 1985, 97, 699-701.
[178] J. Gasteiger, H. Saller, Angew. Chem., Int. Ed. 1985, 24, 687-689.
[179] A. Banerjee, N. Adams, J. Simons, J. Phys. Chem. 1985, 89, 52-57.
[180] J. Baker, J. Comput. Chem. 1982, 3, 214-218.
[181] T. Clark, A. Alex, B. Beck, F. Burkhardt, J. Chandrasekhar, P. Gedeck,
A. Horn, M. Hutter, B. Martin, G. Rauhut, W. Sauer, T. Schindler, T.
Steinke, VAMP 8.0 Build 27, Computer-Chemie-Centrum, Universität
Erlangen-Nürnberg, Erlangen, Germany, 2002.
[182] ViewerLite 4.2, Accelrys Inc., San Diego, CA, 2001.
[183] A. Gustchina, I. T. Weber, FEBS Lett. 1990, 269, 269-272.
[184] B. Beck, G. Rauhut, T. Clark, J. Comput. Chem. 1994, 15, 1064-1073.
115
[185] H. Lanig, R. König, T. Clark, TRAMP 1.1b, Computer-Chemie-Centrum,
Universität Erlangen-Nürnberg, Erlangen, Germany, 1995.
[186] Perl website: http://www.perl.com/
[187] Python website: http://www.python.org/
[188] A. Alex, P. Finn, THEOCHEM 1997, 398-399, 551-554.
[189] T. Brinck, P. Jin, J. Ma, J. S. Murray, P. Politzer, J. Mol. Model. 2003, 9,
77–83.
[190] S. Avram, C. Bologa, M.-L. Flonta, J. Mol. Model. 2005, 11, 105-115.
[191] A. L. Swain, M. M. Miller, J. Green, D. H. Rich, J. Schneider, S. B. H.
Kent, A. Wlodawer, Proc. Natl. Acad. Sci. U. S. A. 1990, 87, 8805-8809.
[192] M. L. Raymer, P. C. Sanschagrin, W. F. Punch, S. Venkataraman, E. D.
Goodman, L. A. Kuhn, J. Mol. Biol. 1997, 265, 445-464.
[193] A. Palomer, J. J. Pérez, S. Navea, O. Llorens, J. Pascual, L. García, D.
Mauleón, J. Med. Chem. 2000, 43, 2280-2284.
[194] M. Rarey, B. Kramer, T. Lengauer, Proteins: Struct., Funct., Genet. 1999,
34, 17-28.
[195] M. L. Nelson, in Tetracyclines in Biology, Chemistry and Medicine (Eds.:
M. Nelson, W. Hillen, R. A. Greenwald), Birkhäuser Verlag, Basel,
Boston, Berlin, 2001, pp. 3-63.
[196] J. Falbe, M. Regitz, CD Römpp Chemie Lexikon, 9th ed., Georg Thieme
Verlag, Stuttgart, New York, 1995.
[197] G. J. Armelagos, K. Kolbacher, K. Collins, J. Cook, M. Krafeld-
Daugherty, in Tetracyclines in Biology, Chemistry and Medicine (Eds.: M.
Nelson, W. Hillen, R. A. Greenwald), Birkhäuser Verlag, Basel, Boston,
Berlin, 2001, pp. 219-236.
[198] G. Forloni, N. Angeretti, R. Chiesa, E. Monzani, M. Salmona, O. Bugiani,
F. Tagliavini, Nature (London, U. K.) 1993, 362, 543-545.
[199] G. Forloni, M. Vari, L. Colombo, O. Bugiani, F. Tagliavini, M. Salmona,
Curr. Med. Chem.: Immunol., Endocr. Metab. Agents 2003, 3, 185-197.
[200] U. Cosentino, M. R. Varì, A. A. G. Saracino, D. Pitea, G. Moro, M.
Salmona, J. Mol. Model. 2005, 11, 17-25.
[201] S. Schneider, in Tetracyclines in Biology, Chemistry and Medicine (Eds.:
M. Nelson, W. Hillen, R. A. Greenwald), Birkhäuser Verlag, Basel,
Boston, Berlin, 2001, pp. 65-104.
116
[202] O. G. Othersen, F. Beierlein, H. Lanig, T. Clark, J. Phys. Chem. B 2003,
107, 13743-13749.
[203] K. Meindl, T. Clark, J. Phys. Chem. B 2005, 109, 4279-4284.
[204] H. Lanig, M. Gottschalk, S. Schneider, T. Clark, J. Mol. Model. 1999, 5,
46-62.
[205] O. G. Othersen, H. Lanig, T. Clark, J. Med. Chem. 2003, 46, 5571-5574.
[206] W. Hinrichs, C. Fenske, in Tetracyclines in Biology, Chemistry and
Medicine (Eds.: M. Nelson, W. Hillen, R. A. Greenwald), Birkhäuser
Verlag, Basel, Boston, Berlin, 2001, pp. 107-123.
[207] W. Saenger, P. Orth, C. Kisker, W. Hillen, W. Hinrichs, Angew. Chem.
2000, 112, 2122-2133.
[208] W. Saenger, P. Orth, C. Kisker, W. Hillen, W. Hinrichs, Angew. Chem.,
Int. Ed. 2000, 39, 2042-2052.
[209] P. Orth, D. Schnappinger, W. Hillen, W. Saenger, W. Hinrichs, Nat.
Struct. Biol. 2000, 7, 215-219.
[210] W. Hinrichs, C. Kisker, M. Düvel, A. Müller, K. Tovar, W. Hillen, W.
Saenger, Science (Washington, DC, U. S.) 1994, 264(5157), 418-420.
[211] N. Guex, M. C. Peitsch, Electrophoresis 1997, 18, 2714-2723.
[212] DeepView/Swiss PdbViewer website: http://www.expasy.org/spdbv/
[213] POV-Ray website: http://www.povray.org/
[214] H. Meiselbach, Diploma thesis, Universität Erlangen-Nürnberg, Erlangen,
Germany, 2003.
[215] M. Gossen, H. Bujard, in Tetracyclines in Biology, Chemistry and
Medicine (Eds.: M. Nelson, W. Hillen, R. A. Greenwald), Birkhäuser
Verlag, Basel, Boston, Berlin, 2001, pp. 139-157.
[216] J. E. Kudlow, in Tetracyclines in Biology, Chemistry and Medicine (Eds.:
M. Nelson, W. Hillen, R. A. Greenwald), Birkhäuser Verlag, Basel,
Boston, Berlin, 2001, pp. 159-175.
[217] E. Pook, in Tetracyclines in Biology, Chemistry and Medicine (Eds.: M.
Nelson, W. Hillen, R. A. Greenwald), Birkhäuser Verlag, Basel, Boston,
Berlin, 2001, pp. 125-137.
[218] T. Förster, Z. Naturforsch. 1949, 4a(5), 321-327.
[219] T. Förster, Ann. Phys. (Leipzig) 1948, 2, 55-75.
[220] L. Stryer, Annu. Rev. Biochem. 1978, 47, 819-846.
117
[221] J. R. Lakowicz, Principles of Fluorescence Spectroscopy, 2nd ed.,
Kluwer Academic/Plenum Publishers, New York, 1999.
[222] S. Brakmann, N. Nöbel, Nachr. Chem. 2003, 51, 319-323.
[223] R. E. Dale, J. Eisinger, Proc. Natl. Acad. Sci. U. S. A. 1976, 73, 271-273.
[224] C. G. Dos Remedios, P. D. J. Moens, J. Struct. Biol. 1995, 115, 175-185.
[225] L. Stryer, R. P. Haugland, Proc. Natl. Acad. Sci. U. S. A. 1967, 58, 719-
726.
[226] D. Summerer, A. Marx, Angew. Chem. 2002, 114, 3778-3780.
[227] D. Summerer, A. Marx, Angew. Chem., Int. Ed. 2002, 41, 3620-3622.
[228] N. Thelwell, S. Millington, A. Solinas, J. Booth, T. Brown, Nucleic Acids
Res. 2000, 28, 3752-3761.
[229] A. Damjanovic, T. Ritz, K. Schulten, Phys. Rev. E: Stat., Nonlinear, Soft
Matter Phys. 1999, 59, 3293-3311.
[230] M. Byrdin, P. Jordan, N. Krauss, P. Fromme, D. Stehlik, E. Schlodder,
Biophys. J. 2002, 83, 433-457.
[231] P. R. Callis, in Methods in Enzymology, Vol. 278 (Eds.: L. Brand, M. L.
Johnson), Academic Press, San Diego, New York, Boston, London,
Sydney, Tokyo, Toronto, 1997, pp. 113-150.
[232] J. T. Vivian, P. R. Callis, Biophys. J. 2001, 80, 2093-2109.
[233] J. Catalán, C. Díaz, Chem. Phys. Lett. 2003, 368, 717-723.
[234] J. Catalán, V. López, P. Pérez, R. Martín-Villamil, J. G. Rodríguez,
Liebigs Ann. 1995, 241-252.
[235] J. Catalán, V. López, P. Pérez, Liebigs Ann. 1995, 793-795.
[236] X. Shen, J. R. Knutson, J. Phys. Chem. B 2001, 105, 6260-6265.
[237] M. Chabbert, E. Piémont, F. G. Prendergast, H. Lami, Arch. Biochem.
Biophys. 1995, 322, 429-436.
[238] A. G. Szabo, D. M. Rayner, J. Am. Chem. Soc. 1980, 102, 554-563.
[239] D. M. Rayner, A. G. Szabo, Can. J. Chem. 1978, 56, 743-745.
[240] J. W. Petrich, M. C. Chang, D. B. McDonald, G. R. Fleming, J. Am.
Chem. Soc. 1983, 105, 3824-3832.
[241] D. Creed, Photochem. Photobiol. 1984, 39, 537-562.
[242] B. Zelent, J. Kusba, I. Gryczynski, M. L. Johnson, J. R. Lakowicz,
Biophys. Chem. 1998, 73, 53-75.
[243] C. M. Hutnik, A. G. Szabo, Biochemistry 1989, 28, 3935-3939.
118
[244] P. S. Antonini, W. Hillen, N. Ettner, W. Hinrichs, P. Fantucci, S. M.
Doglia, J.-A. Bousquet, M. Chabbert, Biophys. J. 1997, 72, 1800-1811.
[245] D.-T. Marian, C. Leypold, M. Kunz, S. Schneider, P. Schubert, O. Scholz,
W. Hillen, Photochem. Photobiol. Sci. 2002, 1, 841-851.
[246] P. Kaszycki, A. Guz, M. Drwiega, Z. Wasylewski, J. Protein Chem. 1996,
15, 607-619.
[247] M. Kunz, M. Kintrup, W. Hillen, S. Schneider, Photochem. Photobiol.
2000, 72, 35-48.
[248] C. Peviani, W. Hillen, N. Ettner, H. Lami, S. M. Doglia, E. Piémont, C.
Ellouze, M. Chabbert, Biochemistry 1995, 34, 13007-13015.
[249] M. Chabbert, W. Hillen, D. Hansen, M. Takahashi, J.-A. Bousquet,
Biochemistry 1992, 31, 1951-1960.
[250] A. Sillen, J. Hennecke, D. Roethlisberger, R. Glockshuber, Y.
Engelborghs, Proteins: Struct., Funct., Genet. 1999, 37, 253-263.
[251] A. Sillen, J. F. Díaz, Y. Engelborghs, Protein Sci. 2000, 9, 158-169.
[252] J. M. G. Martinho, A. M. Santos, A. Fedorov, R. P. Baptista, M. A. Taipa,
J. M. S. Cabral, Photochem. Photobiol. 2003, 78, 15-22.
[253] D. Toptygin, R. S. Savtchenko, N. D. Meadow, L. Brand, J. Phys. Chem.
B 2001, 105, 2043-2055.
[254] D. Toptygin, L. Brand, Chem. Phys. Lett. 2000, 322, 496-502.
[255] Z. Bajzer, F. G. Prendergast, Biophys. J. 1993, 65, 2313-2323.
[256] J. R. Alcala, E. Gratton, F. G. Prendergast, Biophys. J. 1987, 51, 925-
936.
[257] J.-C. Brochon, in Methods in Enzymology, Vol. 240 (Eds.: M. L. Johnson,
L. Brand), Academic Press, San Diego, New York, Boston, London,
Sydney, Tokyo, Toronto, 1994, pp. 262-311.
[258] E. Bismuto, G. Irace, S. D'Auria, M. Rossi, R. Nucci, Eur. J. Biochem.
1997, 244, 53-58.
[259] M. Kintrup, P. Schubert, M. Kunz, M. Chabbert, P. Alberti, E. Bombarda,
S. Schneider, W. Hillen, Eur. J. Biochem. 2000, 267, 821-829.
[260] M. Takahashi, L. Altschmied, W. Hillen, J. Mol. Biol. 1986, 187(3), 341-
348.
[261] F. Beierlein, T. Clark, in High Performance Computing in Science and
Engineering, Munich 2004 - Transactions of the Second Joint HLRB and
119
KONWIHR Status and Result Workshop, March 2-3, 2004, Technical
University of Munich and Leibnitz-Rechenzentrum Munich (Eds.: S.
Wagner, W. Hanke, A. Bode, F. Durst), Springer Verlag, Berlin,
Heidelberg, New York, 2004, pp. 245-260.
[262] F. Beierlein, H. Lanig, O. Othersen, S. Schneider, T. Clark, An MD/CI-
approach simulating FRET in proteins, 227th ACS National Meeting,
Anaheim, CA, United States, 2004.
[263] F. Beierlein, H. Lanig, O. Othersen, S. Schneider, T. Clark, An MD/CI
Approach for the Investigation of Fluorescence Resonance Energy
Transfer in Proteins, 17. Darmstädter Molecular Modelling Workshop,
Erlangen, Germany, 2003.
[264] H. Lanig, F. Beierlein, O. Othersen, S. Schneider, T. Clark, Combining
Molecular Dynamics Simulations with Semiempirical CI-Calculations to
Investigate Fluorescence Resonance Energy Transfer (FRET) within the
Tetracycline Repressor, 43rd Sanibel Symposium, St. Augustine, Florida,
February 22-March 1, 2003.
[265] F. Beierlein, O. Othersen, T. Clark, H. Lanig, A Molecular
Dynamics/Configuration Interaction (MD/CI)-Approach Simulating FRET
in Proteins, eCheminfo 2004 "Applications of Cheminformatics and
Modelling to Drug Discovery", http://echeminfo.colayer.net, 2004.
[266] D. A. Case, T. A. Darden, T. E. Cheatham III, C. L. Simmerling, J. Wang,
R. E. Duke, R. Luo, K. M. Merz, B. Wang, D. A. Pearlman, M. Crowley,
S. Brozell, V. Tsui, H. Gohlke, J. Mongan, V. Hornak, G. Cui, P. Beroza,
C. Schafmeister, J. W. Caldwell, W. S. Ross, P. A. Kollman, AMBER 8,
University of California, San Francisco, 2004.
[267] D. A. Case, D. A. Pearlman, J. W. Caldwell, T. E. Cheatham III, J. Wang,
W. S. Ross, C. L. Simmerling, T. A. Darden, K. M. Merz, R. V. Stanton,
A. L. Cheng, J. J. Vincent, M. Crowley, V. Tsui, H. Gohlke, R. J. Radmer,
Y. Duan, J. Pitera, I. Massova, G. L. Seibel, U. C. Singh, P. K. Weiner, P.
A. Kollman, AMBER 7, University of California, San Francisco, 2002.
[268] G. Rauhut, T. Clark, T. Steinke, J. Am. Chem. Soc. 1993, 115, 9174-
9181.
[269] M. Lee, J. Tang, R. M. Hochstrasser, Chem. Phys. Lett. 2001, 344, 501-
508.
120
[270] S. Ringhofer, J. Kallen, R. Dutzler, A. Billich, A. J. W. G. Visser, D.
Scholz, O. Steinhauser, H. Schreiber, M. Auer, A. J. Kungl, J. Mol. Biol.
1999, 286, 1147-1159.
[271] R. Rajamani, J. Gao, J. Comput. Chem. 2002, 23, 96-105.
[272] N. V. Prabhu, S. D. Dalosto, K. A. Sharp, W. W. Wright, J. M.
Vanderkooi, J. Phys. Chem. B 2002, 106, 5561-5571.
[273] P. R. Callis, B. K. Burgess, J. Phys. Chem. B 1997, 101, 9429-9432.
[274] P. R. Callis, J. T. Vivian, Chem. Phys. Lett. 2003, 369, 409-414.
[275] M. Nonella, G. Mathias, P. Tavan, J. Phys. Chem. A 2003, 107, 8638-
8647.
[276] C. Bosch Cabral, H. Imasato, J. C. Rosa, H. J. Laure, C. H. Tomich de
Paula da Silva, M. Tabak, R. C. Garratt, L. J. Greene, Biophys. Chem.
2002, 97, 139-157.
[277] G. Srinivas, A. Yethiraj, B. Bagchi, J. Phys. Chem. B 2001, 105, 2475-
2478.
[278] P. L. T. M. Frederix, E. L. de Beer, J. Phys. Chem. B 2002, 106, 6793-
6801.
[279] W. Ortiz, A. E. Roitberg, J. L. Krause, J. Phys. Chem. B 2004, 108, 8218-
8225.
[280] C. M. Stultz, A. D. Levin, E. R. Edelman, J. Biol. Chem. 2002, 277,
47653-47661.
[281] C. Laboulais, E. Deprez, H. Leh, J.-F. Mouscadet, J.-C. Brochon, M. Le
Bret, Biophys. J. 2001, 473-489.
[282] M. Gustiananda, J. R. Liggins, P. L. Cummins, J. E. Gready, Biophys. J.
2004, 86, 2467-2483.
[283] T. Clark, J. Chandrasekhar, Isr. J. Chem. 1993, 33, 435-448.
[284] Thomas Hebbeker web page: http://www-eep.physik.hu-
berlin.de/~hebbeker/lectures/i400_ind.htm
[285] J. Kurzawa, S. Schneider, personal communication.
[286] Origin 7G SR4, OriginLab Corporation, Northampton, MA.
[287] Origin 7.0 manual: http://www.originlab.com
[288] W. Hinrichs, W. Hillen, personal communication.
121
[289] T. Clark, A. Alex, B. Beck, J. Chandrasekhar, P. Gedeck, A. Horn, M.
Hutter, B. Martin, G. Rauhut, W. Sauer, T. Schindler, T. Steinke, VAMP
7.0, Oxford Molecular Group Plc, Oxford, UK, 1998.
[290] B. Beck, T. Clark, R. C. Glen, J. Comput. Chem. 1997, 18, 744-756.
[291] U. Essmann, L. Perera, M. L. Berkowitz, T. Darden, H. Lee, L. G.
Pedersen, J. Chem. Phys. 1995, 103, 8577-8593.
[292] RRZE web page: http://www.rrze.uni-erlangen.de/
[293] T. E. Cheatham III, Ptraj 6.5, University of Utah, Salt Lake City, Utah,
2003.
[294] Quatfit, http://ccl.net/cca/software/SOURCES/C/quaternion-mol-
fit/index.shtml.
[295] Jmol 5, http://jmol.sourceforge.net, 2001.
[296] TSAR 3.3, Oxford Molecular Ltd., Oxford, England, 2000.
[297] B. Martin, Ph.D. thesis, Universität Erlangen-Nürnberg, Erlangen,
Germany, 2004.
Publications
PUBLICATIONS
• F. Beierlein, T. Clark, Computer Simulations of Enzyme Reaction Mechanisms: Simulation of Protein Spectra, in High Performance Computing in Science and Engineering, Munich 2004 - Transactions of the Second Joint HLRB and KONWIHR Status and Result Workshop, March 2-3, 2004, Technical University of Munich and Leibnitz-Rechenzentrum Munich, Germany (Eds.: S. Wagner, W. Hanke, A. Bode, F. Durst), Springer-Verlag, Berlin, Heidelberg, New York, 2004, pp. 245-259.
• F. Beierlein, H. Lanig, G. Schürer, A. H. C. Horn, T. Clark, Quantum Mechanical/Molecular Mechanical (QM/MM) Docking: An Evaluation for Known Test Systems. Mol. Phys. 2003, 101, 2469-2480.
• O. G. Othersen, F. Beierlein, H. Lanig, T. Clark, Conformations and Tautomers of Tetracycline, J. Phys. Chem. B. 2003, 107, 13743-13749.
• F. R. Beierlein, O. G. Othersen, H. Lanig, S. Schneider, T. Clark, Simulating FRET: Is the Rotamer Model Correct?, in preparation.
TALKS AT INTERNATIONAL CONFERENCES
• F. Beierlein, H. Lanig, O. Othersen, S. Schneider, T. Clark, An MD/CI-Approach Simulating FRET in Proteins, 227th ACS National Meeting, Anaheim, CA, March 28-April 1, 2004.
• F. Beierlein, H. Lanig, O. Othersen, S. Schneider, T. Clark, An MD/CI Approach for the Investigation of Fluorescence Resonance Energy Transfer in Proteins, 17. Darmstädter Molecular Modelling Workshop, Erlangen, May 27-28, 2003.
• F. Beierlein, T. Clark, A QM/MM Docking Approach with a Flexible Protein Environment, 16. Darmstädter Molecular Modelling Workshop, Darmstadt, May 7-8, 2002.
• F. Beierlein, T. Clark, Simulating the Induced Fit: A Combined QM/MM Docking Approach, MGMS Young Modellers' Forum in Conjunction with the RSC MMG, London, November 30, 2001.
POSTERS
• F. Beierlein, O. G. Othersen, H. Lanig, S. Schneider, T. Clark, Simulating FRET: A Combined Molecular Dynamics-QM/MM Approach, 19. Darmstädter Molecular Modelling Workshop, Erlangen, May 3-4, 2005.
• F. Beierlein, H. Lanig, O. Othersen, S. Schneider, T. Clark, Simulating Fluorescence Resonance Energy Transfer: A Combined Molecular Dynamics-QM/MM-CI Approach, 45th Sanibel Symposium, St. Simons Island, GA, March 5-11, 2005.
• J. Bulitta, F. Schönfeld, F. Beierlein, H. Lanig, T. Clark, R. Troschütz, Chiral Structure Separation for Drug Synthesis: Docking of Methotrexate and Derivatives into the Dihydrofolate Reductase (DHFR) Binding Site, 17. Darmstädter Molecular Modelling Workshop, Erlangen, May 27-28, 2003.
• F. Beierlein, H. Lanig, G. Schürer, A. Horn, T. Clark, Simulating the Induced Fit: A Combined QM/MM Docking Approach, Model(l)ing 2001, The Annual International Meeting of the Molecular Graphics and Modelling Society, Erlangen, September 17-21, 2001.
Curriculum Vitae
FAMILY STATUS
• Name: Frank Rainer Beierlein
• Born: January 24, 1973, Erlangen, Germany
• Married, no children
EDUCATION
• Since Sept 2001 Ph.D. student, Ph.D. thesis “QM/MM Docking and Simulations of FRET”, supervisor: Prof. Dr. Tim Clark, Computer-Chemistry-Center, University of Erlangen-Nuremberg
• Aug 2001 Diplom-Chemiker Univ.
• Jan – Aug 2001 Diploma thesis “QM/MM Docking Techniken“, supervisor: Prof. Dr. Tim Clark, Computer-Chemistry-Center, University of Erlangen-Nuremberg
• March 1997 Vordiplom
• April 1995 Change of university courses: chemistry (diploma course) (1 semester credited with previous studies)
• March 1995 Zwischenprüfung
• Oct 1992 – March 1995 Teacher training (biology and chemistry), University of Erlangen-Nuremberg
• July 1992 Abitur, intensive courses in mathematics and biology
• Sept 1983 – July 1992 Emil-von-Behring-Gymnasium, Spardorf
EMPLOYMENT
• July 2002 – July 2003; since Nov 2003 scientific assistant, Computer-Chemistry-Center, University of Erlangen-Nuremberg, Prof. Dr. Clark
• Jan – June 2002 scientific assistant, Institute of Physical and Theoretical Chemistry, University of Erlangen-Nuremberg, Prof. Dr. Schneider
• Sept – Dec 2001 scientific assistant, Computer-Chemistry-Center, University of Erlangen-Nuremberg, Prof. Dr. Clark
• Jan – June 1998 Student co-worker, Institute of Physical and Theoretical Chemistry, University of Erlangen-Nuremberg, Prof. Dr. Löhmannsröben
• Sept – Oct 1993 Practical teacher training, Emil-von-Behring-Gymnasium, Spardorf
• July – Aug 1990 Student trainee, Siemens Medical Solutions, Erlangen
SCHOLARSHIP
• March 1998 Fellow of the Studienstiftung des deutschen Volkes (German National Academic Foundation), financial funding until Sept 2000