Theoretische Grundlagen molekulardynamischer Simulationen

Die molekulardynamische Simulation (MD-Simulation) ist eine der wichtigsten Methoden zur theoretischen Untersuchung biologischer Moleküle. Es ist äußerst schwierig, wenn nicht unmöglich, die Dynamik biologischer Moleküle direkt auf molekularer Ebene zu verfolgen.

Mittels MD-Simulationen wird versucht ein Modell des zu untersuchenden Moleküls zu erstellen, welches einerseits dessen dynamische Eigenschaften nachahmt und andererseits makroskopisch betrachtet werden kann. Diese computergestützte Methode berechnet das zeitabhängige Verhalten eines molekularen Systems und erlaubt es, Erkenntnisse über die komplexen und dynamischen Prozesse zu gewinnen, die in biologischen Systemen auftreten.

MD-Simulationen ermöglichen somit detaillierte Aussagen über die Struktur, Dynamik und Thermodynamik biologischer Moleküle und deren Wechselwirkungen beispielsweise mit

Substraten. Bei MD-Simulationen wird ein biologisches Molekül vereinfacht gesehen, als ein Netzwerk von untereinander durch Federn verbundene Kugeln betrachtet. Die Federn sind mit einer Kraftkonstante versehen und die Kugeln tragen Punktladungen in ihrem Mittelpunkt.

Der Radius einer Kugel entspricht dabei dem Van-der-Waals Radius des entsprechenden Atoms. Damit lassen sich gebundene sowie nicht gebundene Wechselwirkungen zwischen den Kugeln beschreiben, welche im Allgemeinen unter dem Begriff „Kraftfeld“ zusammengefasst sind. Die Positionen und Geschwindigkeiten der einzelnen Atome über einen bestimmten Simulationszeitraum bezeichnet man als Trajektorie. Sie wird durch Integration der Newton-Bewegungsgleichung erhalten.

8.1.5.1 Kraftfelder

Es stehen eine Vielzahl von Kraftfeldern für MD-Simulationen verschiedener Systeme zur Verfügung. Alle Kraftfelder berücksichtigen in Form von Kraftfeldgleichungen die Energie der bindenden sowie nicht bindenden Wechselwirkungen zwischen den Atomen. Unter den bindenden Wechselwirkungen werden dabei die Bindungslänge zwischen zwei Atomen, der Bindungswinkel zwischen drei Atomen und der Torsionswinkel zwischen zwei Ebenen zusammengefasst. Nicht bindende Wechselwirkungen setzen sich zusammen aus Van-der-Waals-Wechselwirkungen, die durch das Lennard-Jones-Potential beschrieben werden und elektrostatischen Kräften, denen ein Coulomb-Potential zu Grunde liegt. Die potentielle Energie des Systems ergibt sich aus der Summe der Einzelgleichungen des jeweiligen Kraftfeldes. In dieser Arbeit wurde das im AMBER Programmpaket enthaltene Kraftfeld parm99 verwendet,^462,463 welches für Proteine und Nukleinsäuren gut geeignet ist.

8.1.5.2 Integration der Newton-Bewegungsgleichung

Das zweite Newton’sche Gesetz der Bewegung (2) beschreibt den Zusammenhang zwischen der auf ein Atom mit der Masse mi wirkenden Kraft Fi und der daraus für das Atom resultierenden Beschleunigung ai.

i i

i ma

F = (2)

Die Kraft Fi entspricht darüber hinaus dem negativen Gradienten der potentiellen Energie in Bezug auf die Änderung der Atomposition (3).

i pot

i dr

F =−dE (3)

Gleichung (4) drückt die Beziehung zwischen der potentiellen Energie und der Positionsänderung der Atome als Funktion der Zeit aus:

2

Mit Hilfe der Atommassen und der aus den Kraftfeldgleichungen bekannten potentiellen Energie ist es nun möglich, die Position und Geschwindigkeit für jedes Atom über eine Reihe diskreter Zeitintervalle zu bestimmen. Die Zeitintervalle müssen dabei so eng gewählt werden, dass auch die Bewegungen mit der höchsten Frequenz abgebildet werden. Diese stellen in biologischen Systemen die Streckschwingungen der an Kohlenstoff oder Sauerstoff gebunden Wasserstoffe dar (fmax-1 ~ 10 fs).

8.1.5.3 Temperatur- und Druckkopplung

Biologische Vorgänge finden meistens bei konstanter Temperatur und konstantem Druck statt, weshalb MD-Simulationen von biologischen Makromolekülen ebenfalls unter diesen Bedingungen durchgeführt werden. Die Temperatur des Systems steht dabei über die kinetische Energie mit den Geschwindigkeiten der Atome in Beziehung. Weicht die Temperatur an einem bestimmten Zeitpunkt von der vorgegebenen Referenztemperatur ab, so wird in den nachfolgenden Zeitschritten die Geschwindigkeit der Atome durch einen in jedem Zeitschritt neu berechneten Faktor reskaliert und damit die Temperatur über eine gewisse Relaxationszeit wieder an den Referenzwert angeglichen. Die Regulation der Temperatur entspricht der Kopplung an ein virtuelles Wärmebad (Berendsen-Methode).⁴⁶⁴ Durch die Wahl der Relaxationszeit kann eine enge oder lose Kopplung an das Wärmebad erfolgen. Die Regulation des Drucks erfolgt analog zur Regulation der Temperatur. An Stelle der Geschwindigkeiten der Atome werden die interatomaren Atomabstände reskaliert, sodass eine gezielte Veränderung des Simulationsvolumens erreicht wird.

8.1.5.4 Periodische Randbedingungen

Die MD-Simulationen von P450 BM-3 wurden in der vorliegenden Arbeit in explizitem Wasser durchgeführt. Dazu wurde das Protein in einer Wasserbox mit der Form eines truncated octahedron (14-seitiger Polyeder mit acht sechseckigen und sechs viereckigen Seitenflächen) simuliert, welche periodische Randbedingungen aufwies. Alle Atome innerhalb der Wasserbox werden in alle drei Raumrichtungen unendlich oft repliziert. Sobald ein Atom die Grenze der Box überschreitet, tritt es auf der gegenüberliegenden Seite wieder in die Box ein. Dadurch kann das zu simulierende System klein gehalten werden und es wird Rechenzeit eingespart. Um elektrostatische Interaktionen des Proteins mit sich selbst in der benachbarten Wasserbox zu vermeiden, muss die Gesamtladung des Systems neutral sein und die Schichtdicke des Lösungsmittels um das Protein ausreichend groß gewählt werden.

8.1.5.5 Langreichweitige elektrostatische Wechselwirkungen (Particle-Mesh-Ewald-Algorithmus)

Bei der Berechnung des elektrostatischen Potentials bzw. der Coulomb-Energie werden alle Ladungen innerhalb der Zelle und zusätzlich die Ladungen in den periodischen Nachbarzellen berücksichtigt. Diese Berechnungen sind sehr aufwendig, deshalb verwendet man bei MD-Simulationen häufig Particle-Mesh-Ewald-Algorithmen (PME-Methode). Der Geschwindig-keitszuwachs der PME-Methode beruht auf einer Umformulierung der Ladungsverteilung und der Aufteilung der Coulomb-Summe in zwei schnell konvergierende Summen. Zu den realen Punktladungen wird jeweils eine gaußförmige, abschirmende Ladungsverteilung addiert und das elektrostatische Potential berechnet. Um die Abweichungen durch die abschirmenden Gaußladungen zu korrigieren, führt man eine weitere gaußförmige Ladungsverteilung mit entgegengesetztem Vorzeichen ein. Das elektrostatische Potential aus der Ladungsverteilung kann als Fourierreihe effizient mittels Fast-Fourier-Transformation berechnet werden.

8.1.6 Parametrisierung der Hämgruppe und des Substrates