MD Simulation

3. Versuchsdurchführung und Auswertung

Als Grundlage für die Simulation dient ein PDB-File des Peptids Vasopressin. Entweder wird dies durch den Assistenten vorgegeben oder es wird das aus der Strukturrechnung gewonnene File verwendet (wenn die Qualität gut ist). In diesem Fall müssen noch einige Korrekturen vorgenommen werden (z.B. die Pseudoatome entfernt werden).

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

Wir beginnen mit einer Simulation im Vakuum von 100 ps Länge. Als erster Schritt wird nun das PDB-File 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 Cund 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 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 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 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 Parameterfi le, 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).

Nachdem wir nun eine Simulation im Vakuum durchgeführt haben, wollen wir die Simulation noch einmal durchführen. 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 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 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.

Nun führen wir 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

Nun können wir allerdings nicht sofort mit dem eigentlichen Simulationslauf fortfahren. Vorher 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 wird für 20 ps um auf der sicheren Seite zu sein. Auch für diesen Schritt benötigen wir eine Parameterdatei. Wir bennen 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

Nun bereiten 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

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 Programm g_confrms liefert den RMSD sowie ein pdb File, in dem die beiden Strukturen überlagert sind:

g_confrms -f1 vaso.pdb -f2 vaso_pmd.gro -o fi t.pdb

Um den Verlauf von Energie, Druck, Temperatur etc. zu visualisieren, verwenden wir 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 „Radius of Gyration“

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. Man kann 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 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