MD Simulation

Latest revision as of 10:51, 15 January 2010

Als Grundlage für die Moleküldynamik Simulation mit dem Programm GROMACS dient ein PDB File des Peptids Vasopressin.

Vorbereitung

Loggen Sie sich mit dem Usernamen und Passwort Ihres HRZ Kontos ins Windows System ein.

Die MD Simulation wird auf einem Linux Server durchgeführt, der über den Icon BCC Linux Terminal-Service.nxs zugänglich ist. Der Username und das Passwort für den Linux Server sind die gleichen wie bei Windows. Oeffnen Sie ein Terminal Fenster.

Holen Sie die Input Daten für die MD Simulation in Ihr Home Directory und führen Sie die folgenden Setup Befehle aus:

cd
bash
tar zxf MD.tgz
cd MD

Alle Befehle werden in einem Terminalfenster über die Kommandozeile eingegeben. Neben dem Terminalfenster wird noch ein Text-Editor, das Programm xmgrace für die graphische Anzeige von Diagrammen und als Molekülbetrachter das Programm VMD benötigt.

Die Programme xmgrace und VMD sind nur auf dem 'blade41' Linux Server vorhanden. Um Sie zu verwenden, öffnen Sie ein weiteres Terminal Fenster und loggen Sie mit den folgenden Befehlen in den 'blade41' Server ein

ssh -X blade41.rz.uni-frankfurt.de
bash

Aus Effizienzgründen sollte die MD Rechnung mit GROMACS jedoch auf dem Server durchgeführt werden, in den Sie ursprünglich eingeloggt worden sind.

MD Simulation im Vakuum

Wir beginnen mit einer Simulation im Vakuum von 100 ps Länge und wechseln dazu ins Unterverzeichnis 'vakuum'.

cd vakuum

Als erster Schritt wird nun das PDB-File mit dem GROMACS Programm pdb2gmx in das von GROMACS verwendete Format (*.gro) umgewandelt und ein File mit einer Topologie (*.top) erstellt

pdb2gmx -f vaso.pdb -o vaso.gro -p vaso.top -ignh -ter

Das Programm fragt nach dem zu verwendenden Kraftfeld. Hier wird GROMOS96 (Auswahl mit Taste „0“) gewählt. Auch muss von Hand eingegeben werden, wie C- und N-Terminus aussehen. Im Fall des als Amid blockierten C-Terminus ist „none“ einzugeben. Für den N-Terminus „NH3+“. Als nächstes erstellen wir mit dem GROMACS Programm editconf die Box, in der die Simulation abläuft

editconf -bt cubic -f vaso.gro -o vaso_b4em.gro -c -d 1.0

Durch diesen Befehl wird eine kubische Box angelegt (-bt cubic), das Peptid wird in der Box zentriert (-c) und die Größe der Box wird so gewählt, dass von ihrer Begrenzung bis zum Molekül etwa 1.0 nm liegen (-d 1.0).

Nun müssen wir, bevor wir mit der Simulation beginnen, noch eine Energieminimierung durchführen. Hierfür wird eine Datei mit Parametern benötigt. Sie wird im Texteditor erstellt und sollte wie folgt aussehen und em.mdp benannt werden:

title = vaso_vac
cpp = /lib/cpp
defi ne = -DFLEX_SPC
constraints = none
integrator = steep
dt = 0.002
nsteps = 400
nstlist = 10
ns_type = grid
rlist = 0.9
coulombtype = PME
rcoulomb = 0.9
rvdw = 0.9
fourierspacing = 0.12
fourier_nx = 0
fourier_ny = 0
fourier_nz = 0
pme_order = 4
ewald_rtol = 1e-5
optimize_fft = yes
;
; Energy minimizing stuff
;
emtol = 1000.0
emstep = 0.01

Nun müssen wir die Files mit dem GROMACS Programm grompp prozessieren. grompp steht für „Gromacs pre-processor“.

grompp -f em.mdp -c vaso_b4em.gro -p vaso.top -o vaso_em.tpr

Nun führen wir die Energieminimierung mit dem GROMACS Programm mdrun durch:

mdrun -s vaso_em.tpr -o vaso_em.trr -c vaso_b4md.gro -v

Das Flag -v sorgt dafür, dass das Programm viele Ausgaben auf der Kommandozeile ausgibt („verbose“). Nun können wir die eigentliche Simulation (im Vakuum) starten. Dazu benötigen wir wieder ein Parameterfile, dass wir diesmal md.mdp nennen.

title = vaso_vac
cpp = /lib/cpp
constraints = all-bonds
integrator = md
dt = 0.002 ; ps !
nsteps = 50000 ; total 100 ps
nstcomm = 1
nstxout = 500 ; collect data every 1 ps
nstvout = 0
nstfout = 0
nstlist = 10
ns_type = grid
rlist = 0.9
coulombtype = PME
rcoulomb = 0.9
rvdw = 0.9
fourierspacing = 0.12
fourier_nx = 0
fourier_ny = 0
fourier_nz = 0
pme_order = 4
ewald_rtol = 1e-5
optimize_fft = yes
;Berendsen temperature coupling is on
Tcoupl = berendsen
tau_t = 0.1
tc-grps = protein
ref_t = 300
;Pressure coupling is on
Pcoupl = berendsen
tau_p = 0.5
compressibility = 4.5e-5
ref_p = 1.0
; Generate velocities is on at 300 K
gen_vel = yes
gen_temp = 300.0
gen_seed = 173529

Nun verwenden wir wieder grompp.

grompp -f md.mdp -c vaso_b4md.gro -r vaso_b4md.gro -p vaso.top -o vaso_md.tpr

Und starten die Simulation

mdrun -s vaso_md.tpr -o vaso_md.trr -c vaso_pmd.gro -g vaso_md.log -e vaso_md.edr -v

Den Verlauf der Simulation kann man mit dem Programm VMD visualisieren (Einführung durch Assistenten).

MD Simulation in Wasser

Nachdem wir nun eine Simulation im Vakuum durchgeführt haben, wollen wir die Simulation noch einmal durchführen und wechseln dazu ins Unterverzeichnis 'loemi'.

cd ../loemi

Diesmal füllen wir die Simulationsbox allerdings mit Wassermolekülen und sorgen durch das Hinzufügen von Gegenionen auch dafür, dass die Ladung des Vasopressins ausgeglichen wird. Letzters ist vor allem notwendig, damit die Particle Mesh Ewald Methode richtig funktioniert. Nachdem wir auch diese Simulation durchgeführt haben, können wir die Ergebnisse vergleichen. Für die Simulation mit Lösungsmittel werden nur die Schritte in Detail erläutert, die sich von der Simulation im Vakuum unterscheiden.

Wir beginnen wieder mit der Umwandlung des pdb-Files und der Box:

pdb2gmx -f vaso.pdb -o vaso.gro -p vaso.top -ignh -ter
editconf -bt cubic -f vaso.gro -o vaso.gro -c -d 1.0

Nun müssen wir mit dem GROMACS Programm genbox die Box mit Wassermolekülen füllen. Wir verwenden das SPC Wassermodell(Single Point Charge).

genbox -cp vaso.gro -cs spc216.gro -o vaso_b4em.gro -p vaso.top

Nun führen wir wieder eine Energieminimierung durch. Dafür kann die em.mdp Datei aus der Vakuumsimulation erneut verwendet werden.

grompp -f em.mdp -c vaso_b4em.gro -p vaso.top -o vaso_b4ion.tpr

An dieser Stelle nutzen wir das erstellte *.tpr File, um mit dem GROMACS Programm genion Gegenionen einzufügen. Vasopressin besitzt eine positive Ladung, daher fügen wir Chloridionen ein.

genion -s vaso_b4ion.tpr -o vaso_b4em.gro -nname CL -nn 2 -g vaso_ion.log

Das Programm fragt nach einer Gruppe von Molekülen, von denen einige durch Chlorid ersetzt werden. Hier wählt man die Gruppe „SOL“.

Nun müssen wir auch noch Änderungen an der Datei vaso.top vornehmen. Nach dem include Befehl für das Kraftfeld ist folgende Zeile einzufügen:

#include “ions.itp“

In der Sektion [molecules] am Ende des Dokuments ist eine Zeile für die Chloridionen(CL) einzugeben. Die Anzahl der Chloridionen ist von der Zahl der Lösungsmittelmoleküle abzuziehen.

[ molecules ]
; Compound        #mols
Protein             1
SOL              2173
CL-                 2

Wir führen grompp noch einmal aus und starten dann die Energieminimierung.

grompp -f em.mdp -c vaso_b4em.gro -p vaso.top -o vaso_em.tpr
mdrun -s vaso_em.tpr -o vaso_em.trr -c vaso_b4pr.gro -v

Bevor wir mit dem eigentlichen Simulationslauf fortfahren, ist es notwendig, einen Lauf mit sog. Position-restrained Molecular Dynamics durchzuführen. Dabei wird die Bewegung des Proteins stark eingeschränkt („restrained“), während das Lösungsmittel sich frei bewegen kann. Dabei bildet es die energetisch optimale Struktur zum Protein aus. Die Relaxationszeit von Wasser beträgt etwa 10 ps. Daher simulieren wir für 20 ps, um auf der sicheren Seite zu sein. Auch für diesen Schritt benötigen wir eine Parameterdatei. Wir benennen sie pr.mdp

title = vaso_sol
cpp = /lib/cpp
defi ne = -DPOSRES
constraints = all-bonds
integrator = md
dt = 0.002 ; ps !
nsteps = 10000 ; total 20 ps
nstcomm = 1
nstxout = 250 ; collect data every 0.5 ps
nstvout = 1000
nstfout = 0
nstlog = 10
nstenergy = 10
nstlist = 10
ns_type = grid
rlist = 0.9
coulombtype = PME
rcoulomb = 0.9
rvdw = 0.9
fourierspacing = 0.12
fourier_nx = 0
fourier_ny = 0
fourier_nz = 0
pme_order = 4
ewald_rtol = 1e-5
optimize_fft = yes
;Berendsen temperature coupling is on
Tcoupl = berendsen
tau_t = 0.1 0.1 0.1
tc-grps = protein sol CL-
ref_t = 300 300 300
;Pressure coupling is on
Pcoupl = berendsen
tau_p = 0.5
compressibility = 4.5e-5
ref_p = 1.0
; Generate velocities is on at 300 K
gen_vel = yes
gen_temp = 300.0
gen_seed = 173529

Wir die Dateien vor und starten

grompp -f pr.mdp -c vaso_b4pr.gro -r vaso_b4pr.gro -p vaso.top -o vaso_pr.tpr
mdrun -s vaso_pr.tpr -o vaso_pr.trr -c vaso_b4md.gro -v

Nun können wir die eigentliche Simulation starten. Dafür können wir die md.mdp Datei aus der Vakuumsimulation verwenden, müssen aber dafür sorgen, dass die Sektion „temperature coupling“ genau so aussieht, wie in der Datei pr.mdp, also mit jeweils einem Eintrag für protein, sol und CL.

grompp -f md.mdp -c vaso_b4md.gro -r vaso_b4md.gro -p vaso.top -o vaso_md.tpr
mdrun -s vaso_md.tpr -o vaso_md.trr -c vaso_pmd.gro -g vaso_md.log -e vaso_md.edr -v

Auswertung

Nun haben wir beide Simulationen abgeschlossen und können uns die Trajektorie mit dem Programm VMD anschauen. Doch GROMACS enthält noch eine Menge an Hilfsprogrammen für die Auswertung.

Wir beginnen damit, die Struktur, die am Ende der Simulationen erhalten wurde, mit der Ausgangsstruktur zu vergleichen. Das GROMACS Programm g_confrms liefert den RMSD sowie ein File, in dem die beiden Strukturen überlagert sind:

g_confrms -f1 vaso.pdb -f2 vaso_pmd.gro -o fit.pdb

Um den Verlauf von Energie, Druck, Temperatur etc. zu visualisieren, verwenden wir as GROMACS Programm g_energy. Das Programm fragt, welche Daten es liefern soll

g_energy -f vaso_md.edr -o vaso_ener.xvg

Dateien des Typs *.xvg können entweder in Excel importiert werden oder (und das ist die bessere Methode) man betrachtet sie mit dem Programm xmgrace

xmgrace -nxy vaso_ener.xvg

Ein Maß für die Kompaktheit der Struktur ist der Trägheitsradius (radius of gyration), der mit dem GROMACS Programm g_gyrate berechnet werden kann.

g_gyrate -f vaso_md.trr -s vaso_md.tpr -o vaso_gyrate.xvg

Um sich ein Bild davon zu machen, ob und wie schnell die Simulation ein Minimum der Energiehyperfläche erreicht (also konvergiert), ist der Verlauf des RMSD mit der Zeit ein gutes Hilfsmittel. Mit dem GROMACS Programm g_rms kann man den RMSD für die gesamte Struktur bestimmen oder auch nur für das Rückgrat.

g_rms -s vaso_md.tpr -f vaso_md.trr -o vaso_rmsd.xvg -dt 10

(-dt 10 bedeutet, dass nur jeder 10. Frame herausgeschrieben wird, dadurch wird der Plot weniger „verrauscht“.)

Als letztes wollen mit dem GROMACS Programm g_hbond noch untersuchen, wie viele Wasserstoffbrücken das Peptid intramolekular und mit dem Lösungsmittel ausbildet (die zu wählenden Gruppen sind protein und protein sowie protein mit SOL)

g_hbond -f vaso_md.trr -s vaso_md.tpr -num vaso_hbond.xvg

Weitere Informationen und Literatur

• Praktikumsanleitung
• Wilfred F. van Gunsteren, Herman J. C. Berendsen: Moleküldynamik-Computersimulationen; Methodik, Anwendungen und Perspektiven in der Chemie (1990)
• Wilfred F. van Gunsteren et al.: Biomolekulare Modellierung: Ziele, Probleme, Perspektiven (2006)
• Christopher M. Dobson, Andrej Sali, Martin Karplus: Proteinfaltung aus theoretischer und experimenteller Sicht (1998)