QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv),...

143
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

Transcript of QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv),...

Page 1: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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

Page 2: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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

Page 3: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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.

Page 4: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit
Page 5: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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)

Page 6: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

“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

Page 7: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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

Page 8: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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.

Page 9: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

Ü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

Page 10: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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.

Page 11: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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.

Page 12: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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.

Page 13: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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

Page 14: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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.

Page 15: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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)

Page 16: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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.

Page 17: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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

Page 18: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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

Page 19: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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]

Page 20: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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,

Page 21: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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.

Page 22: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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Ψ = Ψ

Page 23: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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 χ ε χ=∑

Page 24: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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ν νν

ψ φ=

=∑

Page 25: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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

Page 26: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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Ψ = Ψ + Ψ + Ψ +

Page 27: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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

Page 28: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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

Page 29: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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

Page 30: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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 −= + + +

Page 31: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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)

Page 32: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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

∂ ∂− = = =∂ ∂

Page 33: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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=∑

Page 34: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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π ε= = +

⎛ ⎞= ⎜ ⎟⎜ ⎟

⎝ ⎠∑ ∑

Page 35: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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

λτ

⎡ ⎤⎛ ⎞∆= + −⎢ ⎥⎜ ⎟⎝ ⎠⎣ ⎦

Page 36: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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µτ

⎧ ⎫∆= − −⎨ ⎬⎩ ⎭

Page 37: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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).

Page 38: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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= + +

Page 39: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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]

Page 40: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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.

Page 41: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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.

Page 42: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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.

Page 43: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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.

Page 44: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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,

Page 45: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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

Page 46: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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

Page 47: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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

Page 48: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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]

Page 49: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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.

Page 50: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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.

Page 51: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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

Page 52: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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.

Page 53: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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 Å.

Page 54: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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

Page 55: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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 Å.

Page 56: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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.

Page 57: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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).

Page 58: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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

Page 59: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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

Page 60: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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].

Page 61: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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.

Page 62: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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]

Page 63: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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]

Page 64: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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

Page 65: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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]

Page 66: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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.

Page 67: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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]).

Page 68: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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

Page 69: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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

Page 70: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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

Page 71: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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

Page 72: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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]

Page 73: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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

Page 74: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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

Page 75: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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

Page 76: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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]

Page 77: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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.

Page 78: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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,

Page 79: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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

Page 80: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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.

Page 81: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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

Page 82: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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

Page 83: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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

− − − − −= = = −

− − −

∑ ∑ ∑ ∑∑ ∑ ∑

Page 84: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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

Page 85: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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,

Page 86: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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.

Page 87: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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,

Page 88: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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

Page 89: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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]

Page 90: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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

Page 91: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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

= −

Page 92: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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

κ λ ε λ λ λπ

= ∫

Page 93: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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κ = Θ − Θ Θ

Page 94: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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−=∑

Page 95: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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.

Page 96: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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.

Page 97: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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;

Page 98: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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)

Page 99: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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).

Page 100: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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.

Page 101: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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]

Page 102: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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

Page 103: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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°).

Page 104: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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

Page 105: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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]

Page 106: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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.

Page 107: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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

Page 108: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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.

Page 109: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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)

Page 110: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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).

Page 111: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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

Page 112: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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.

Page 113: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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).

Page 114: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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

Page 115: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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.

Page 116: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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.

Page 117: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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)).

Page 118: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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.

Page 119: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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

Page 120: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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.

Page 121: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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

Page 122: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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.

Page 123: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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.

Page 124: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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.

Page 125: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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.

Page 126: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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.

Page 127: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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.

Page 128: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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.

Page 129: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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.

Page 130: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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.

Page 131: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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.

Page 132: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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.

Page 133: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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.

Page 134: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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.

Page 135: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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.

Page 136: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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

Page 137: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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.

Page 138: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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.

Page 139: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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.

Page 140: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit
Page 141: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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.

Page 142: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

• 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.

Page 143: QM/MM Docking and Simulations of FRET - OPUS 4 · Ein Inhibitor der HIV-1-Protease (PDB-Code 3phv), gedockt in das ligandfrei kristallisierte Enzym. Links: Ergebnis des Dockings mit

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