From Tomasz Żok

BioinformatykaStrukturalna2: WprowadzenieDoDynamikiMolekularnej

Wstęp

Pętle typu GNRA

• Pętla typu GNRA (gdzie N to dowolny nukleotyd, czyli A, C, G lub U, a R to puryna czyli A lub G) to najczęstszy motyw strukturalny RNA
• Cechą charakterystyczną wszystkich 8 możliwych układów jest niekanoniczne wiązanie G-A stabilizujące układ
• Dwa pozostałe nukleotydy N i R tworzą wiązania wodorowe z odległymi fragmentami jednoniciowymi lub białkami
• Utrzymują się w ten sposób przestrzenne interakcje wpływające na całościowe zwinięcie się struktury
• Biorą one udział w stabilizacji, szczególnie dla dużych struktur RNA (1000+ nukleotydów)
• Wizualizacja poniżej przedstawia strukturę drugorzędową dla 1ZIH (tutaj pętla to GCAA):

Symulacje dynamiki molekularnej

• Symulacja dynamiki molekularnej to iteracyjne powtarzanie rozwiązywania równań ruchu Newtona dla układu atomów (powiązanych odpowiednimi ograniczeniami wynikającymi ze struktury i właściwości własnych atomów)
• Punktem początkowym jest struktura (=współrzędne wszystkich jej atomów), najczęściej umieszczona w wodzie z jonami celem zobojętnienia
• Układ posiada też temperaturę (czyli de facto wektory prędkości dla każdego z atomów) i znajduje się pod określonym ciśnieniem
• Każdy krok symulacji odpowiada niewielkiej zmianie czasowej (często 1-2 fs [10^-15 s])
• Symuluje się różnej długości trajektorie w zależności od problemu badawczego. Przykładowo, białka zwijają się w czasie rzędu od kilkudziesięciu mikrosekund do nawet kilku sekund. Przeciętny czas trwania symulacji jest rzędu milisekund [10^-3 s], czyli kroków symulacji jest ~10¹²
• Dynamika molekularna pozwala też na:
  • sprawdzenie stabilności struktury po wprowadzeniu modyfikacji lub mutacji,
  • analizę interakcji czwartorzędowych (np. przeprowadzenie tzw. dokowania ligandu do miejsca aktywnego białka),
  • pełnoatomowe potwierdzenie eksperymentalnych wyników (np. temperatura topnienia RNA/DNA to temperatura, w której 50% składu próbki ulega denaturacji; przy pomocy symulacji można ten proces prześledzić dokładniej),
  • przewidzenie struktury trzeciorzędowej na podstawie jej sekwencji,
  • i wiele innych...
• Poniższy film przedstawia symulację pętli GNRA, w której autorom udało się z formy w pełni rozwiniętej (unfolded) otrzymać końcową postać znaną z wielu innych badań strukturalnych: (referencja)

(:html5video filename=GNRA:)

Przeprowadzenie własnej symulacji

Wymagane oprogramowanie

• NAMD
• VMD

Pliki potrzebne do przeprowadzenia symulacji

• Początkowa struktura (współrzędne atomów) w formacie PDB (w naszym przypadku: 1ZIH: pętla GCAA, 10 modeli NMR)
• Pliki parametrów i topologii: (źródło)
  • topologia dla kwasów nukleinowych,
  • parametry dla kwasów nukleinowych,
  • parametry dla wody i jonów
• Informacja o strukturze w formacie PSF (w dalszej sekcji)
• Plik konfiguracyjny opisujący przebieg symulacji (w dalszej sekcji)

Wizualizacja w VMD

1. Wczytaj strukturę do programu VMD (File → New molecule... → Browse...)
2. Klik oraz rolka od myszy pozwalają na obrót oraz zbliżenie/oddalenie
3. Domyślnie VMD wyświetla struktury w reprezentacji "Lines", zmień to następująco:
   1. Otwórz menu Graphics → Representations...
   2. Dla wybranej reprezentacji zmień Drawing Method na NewRibbons
   3. Dodaj reprezentację (przycisk Create Rep), ustaw Drawing Method na HBonds, Coloring Method na ColorID z wartością 1 red oraz Line Thickness równym 5
      
4. Żeby przedstawić wszystkie 10 modeli na jednej wizualizacji:
   1. Wyłącz rysowanie wiązań wodorowych (dwukrotne kliknięcie na liście reprezentacji w HBonds)
   2. Kliknięciem wybierz reprezentację NewRibbons i w zakładce Trajectory, w polu Draw Multiple Frames wpisz: 1:10 (w ogólności od:do gdzie od i do to numery modeli/ramek w trajektorii)
      
5. Żeby sprawdzić różnice pomiędzy modelami, wykorzystaj narzędzie Timeline:
   1. Wybierz menu Extensions → Analysis → Timeline
   2. W nowym okienku, wybierz menu Calculate → Calc. H-bonds... (z domyślnymi wartościami):
      • na osi X pokazane są numery modeli
      • na osi Y identyfikatory atomów
      • czarny kolor oznacza brak wiązania wodorowego, biały natomiast jego obecność
      • ta wizualizacja pozwala sprawdzić stabilność wiązań w warunkach eksperymentu NMR, ponieważ patrząc wierszami widać jak często wiązanie wystąpiło w różnych modelach
        
   3. Wybierz teraz menu Calculate → Calc. RMSD:
      • na osi X ponownie pokazane są numery modeli
      • na osi Y identyfikatory reszt
      • pola wyświetlane są w skali szarości (patrz: legenda w lewym dolnym rogu)
      • im jaśniejsze pole, tym bardziej dana reszta różni się od referencyjnej w pierwszym modelu (dlatego cała pierwsza kolumna jest czarna, RMSD wszędzie równe jest zero)
      • na rys. w punkcie 4 powyżej możesz zauważyć dwie zasady azotowe w innym ułożeniu niż pozostałe 8; to jest właśnie reszta nr 6 (cytydyna) i obserwację te potwierdza wizualizacja z RMSD (dwa wyraźnie najjaśniejsze pola)
        

Przygotowanie pliku PSF (Protein Structure File)

• Plik PDB zawiera współrzędne atomów, natomiast plik PSF ma zapisane informacje o strukturze: masy atomów, wiązania, itp.
• Plik PSF buduje się w oparciu o strukturę PDB oraz topologię, można wykorzystać do tego narzędzie dostępne w VMD w menu Extensions → Modeling → Automatic PSF Builder:
  1. W menu Options zaznacz Add solvation box oraz Add neutralizing ions
  2. Upewnij się, że Molecule ustawione jest na 1ZIH.pdb
  3. Usuń domyślny plik z topologią i wskaż na ten ściągnięty wcześniej (plik: top_all36_na.rtf)
     Uwaga! Zwróć uwagę na to, żeby był to plik z TOPologią, a nie PARametrami, czyli top_all36_na.rtf, a nie par_all36_na.prm
  4. Kliknij w Load input files
  5. W kroku drugim zaznacz Nucleic Acid i kliknij w Guess and split chains using selections (później wskaż lokalizację pliku PDB gdy o to zapyta VMD)
  6. W kroku trzecim wybierz Create chains; wtyczka AutoPSF właśnie utworzyla nowy plik PDB (z przemianowanymi resztami, z atomami wody oraz z jonami Na⁺ i Cl^-), a także plik PSF opisujący jego strukturę wewnętrzną
  7. Struktura umieszczona w wodzie z jonami jest automatycznie wczytana przez VMD, spróbuj zwizualizować strukturę ustawiając trzy aktywne reprezentacje w Graphical Representations:
     • Style = NewCartoon, Color = ColorID 0, Selection = nucleic
     • Style = VDW, Color = Name, Selection = ion
     • Style = Licorice, Color = Name, Selection = water, Bond Radius = 0.1

Symulacja z wykorzystaniem periodycznych warunków brzegowych

1. Ściągnij przedstawiony niżej szablon pliku konfiguracyjnego dla NAMD: (korzystając z linku poniżej GetCode). Nazwij go np. namd.conf
   
   1. #############################################################
   2. ## JOB DESCRIPTION                                         ##
   3. #############################################################
   4.  
   5. # Minimization and Equilibration of
   6. # 1ZIH in a Water Box
   7.  
   8.  
   9. #############################################################
   10. ## ADJUSTABLE PARAMETERS                                   ##
   11. #############################################################
   12.  
   13. structure          1ZIH_autopsf.psf
   14. coordinates        1ZIH_autopsf.pdb
   15.  
   16. set temperature    310
   17. set outputname     1ZIH_namd
   18.  
   19. firsttimestep      0
   20.  
   21.  
   22. #############################################################
   23. ## SIMULATION PARAMETERS                                   ##
   24. #############################################################
   25.  
   26. # Input
   27. paraTypeCharmm      on
   28. parameters          par_all36_na.prm
   29. parameters          par_water_ions.str
   30. temperature         $temperature
   31.  
   32.  
   33. # Force-Field Parameters
   34. exclude             scaled1-4
   35. 1-4scaling          1.0
   36. cutoff              12.0
   37. switching           on
   38. switchdist          10.0
   39. pairlistdist        14.0
   40.  
   41.  
   42. # Integrator Parameters
   43. timestep            2.0  ;# 2fs/step
   44. rigidBonds          all  ;# needed for 2fs steps
   45. nonbondedFreq       1
   46. fullElectFrequency  2  
   47. stepspercycle       10
   48.  
   49.  
   50. # Constant Temperature Control
   51. langevin            on    ;# do langevin dynamics
   52. langevinDamping     1     ;# damping coefficient (gamma) of 1/ps
   53. langevinTemp        $temperature
   54. langevinHydrogen    off    ;# don't couple langevin bath to hydrogens
   55.  
   56.  
   57. # Periodic Boundary Conditions
   58. cellBasisVector1    48.0    0.   0.0
   59. cellBasisVector2     0.0  44.0   0.0
   60. cellBasisVector3     0.0    0   44.0
   61. cellOrigin           0.25  1.48  0.3
   62.  
   63. wrapAll             on
   64.  
   65.  
   66. # PME (for full-system periodic electrostatics)
   67. PME                 yes
   68. PMEGridSpacing      1.0
   69.  
   70. #manual grid definition
   71. #PMEGridSizeX        45
   72. #PMEGridSizeY        45
   73. #PMEGridSizeZ        48
   74.  
   75.  
   76. # Constant Pressure Control (variable volume)
   77. useGroupPressure      yes ;# needed for rigidBonds
   78. useFlexibleCell       no
   79. useConstantArea       no
   80.  
   81. langevinPiston        on
   82. langevinPistonTarget  1.01325 ;#  in bar -> 1 atm
   83. langevinPistonPeriod  100.0
   84. langevinPistonDecay   50.0
   85. langevinPistonTemp    $temperature
   86.  
   87.  
   88. # Output
   89. outputName          $outputname
   90.  
   91. restartfreq         500     ;# 500steps = every 1ps
   92. dcdfreq             250
   93. xstFreq             250
   94. outputEnergies      100
   95. outputPressure      100
   96.  
   97.  
   98. #############################################################
   99. ## EXTRA PARAMETERS                                        ##
   100. #############################################################
   101.  
   102.  
   103. #############################################################
   104. ## EXECUTION SCRIPT                                        ##
   105. #############################################################
   106.  
   107. # Minimization
   108. minimize            100
   109. reinitvels          $temperature
   110.  
   111. run 2500 ;# 5ps

   [$[Get Code]]
2. Zwróć uwagę na zaznaczone linie:
   1. Przez 100 kroków odbywać się będzie minimalizacja energii układu, jest to obowiązkowy początek każdej symulacji dynamiki molekularnej
   2. Przez następne 2500 kroków, czyli przez 5 ps odbywać się będzie swobodna symulacja
   3. Umieść wskazane pliki w jednym katalogu wraz z namd.conf:
      • 1ZIH_autopsf.pdb oraz 1ZIH_autopsf.psf (znajdują się one w katalogu VMD, czyli prawdopobnie: C:\Program Files (x86)\University of Illinois\VMD)
      • par_all36_na.prm oraz par_water_ions.str (do ściągnięcia były wcześniej, linki na początku strony)
   4. Popraw wymiary i położenie cząsteczki wraz z wodą i jonami:
      1. W VMD otwórz menu Extensions → Tk Console
      2. Zaznacz wszystkie atomy i wyznacz min-max:
         set sel [atomselect top all]
         set m [measure minmax $sel]
         foreach {j1 j2} $m {}
         foreach {x2 y2 z2} $j2 {}
         foreach {x1 y1 z1} $j1 {}
      3. Wyznacz poniższą wartość i wpisz ją jako pierwsze pole wektora cellBasisVector1 (pozostałe pola zostaw równe zero):
         expr $x2 - $x1
      4. Wyznacz poniższą wartość i wpisz ją jako drugie pole wektora cellBasisVector2 (pozostałe pola zostaw równe zero):
         expr $y2 - $y1
      5. Wyznacz poniższą wartość i wpisz ją jako trzecie pole wektora cellBasisVector3 (pozostałe pola zostaw równe zero):
         expr $z2 - $z1
      6. Wyznacz środek i wpisz wynik jako cellOrigin:
         measure center $sel
3. Przygotuj plik uruchomieniowy do przeprowadzenia symulacji:
   "C:\Ścieżka Do NAMD\namd2.exe" namd.conf > simulation.log

   [$[Get Code]]

Analiza przebiegu symulacji

• Namd generuje zestaw plików wyjściowych, nas interesować będzie simulation.log (tekstowy log z przebiegu symulacji) oraz 1ZIH_namd.dcd (binarnie zapisana trajektoria = pozycje atomów zapisywane co pewną liczbę iteracji)
• Aby otworzyć trajektorię w VMD:
  1. Wybierz z menu File → New Molecule...
  2. Upewnij się, że pole Load files for: ma wartość New Molecule
  3. Wybierz plik 1ZIH_autopsf.psf (Uwaga! Wybierz plik PSF, a nie PDB. Dla tego drugiego wczytanie trajektorii też zadziała, ale wówczas VMD nie skorzysta z pełni informacji o strukturze jaka jest tylko w PSF)
  4. Po naciśnięciu Load okienko wczytywania powinno zostać aktywne, ale pole Load files for zawierać teraz będzie nazwę 1ZIH_autopsf.psf
  5. Wybierz plik 1ZIH_namd.dcd i kliknij w Load
• Aby przeanalizować plik logu w VMD:
  1. Wybierz z menu Extensions → Analysis → NAMD Plot
  2. W nowym okienku wybierz File → Select NAMD Log File i wskaż na plik simulation.log
  3. Zaznacz przykładowo BOND, ANGLE oraz DIHED odpowiadające trzem komponentom funkcji energii w modelu
     
  4. Wybierz z menu File → Plot Selected Data, na wykresie wyraźnie widać moment przejścia symulacji z trybu minimalizacji energii (optymalizacji) do swobodnej symulacji cząsteczki
     
• Prześledź interakcję między resztami G-A (skrajne nukleotydy w pętli GNRA):
  1. Wybierz menu Graphics → Representations...
  2. Ustaw następującą wizualizację: (pogrubienie interesujących nas reszt oraz wyświetlenie atomów wody wokół nich)
     • Style = Lines, Color = name, Selection = water and same resid as within 3.5 of (nucleic and resid 5 or resid 8)
     • Style = Licorice, Color = name, Selection = nucleic and resid 5 or resid 8
       
  3. W głównym oknie VMD przesuwaj po kolejnych zapisanych ramkach i obserwuj wzajemne ułożenie się reszt oraz cząsteczek wody
     

Analiza gotowej, dłuższej trajektorii

• Pełniejszą i trójetapową symulację (minimalizacja, podgrzewanie, swobodny ruch) możemy przeanalizować dzięki uprzejmości dr Sarzyńskiej:
• Aby wczytać pełną trajektorię:
  1. Pobierz i rozpakuj archiwum z wynikami oraz pliki PDB i PSF
  2. Wybierz w VMD File → New Molecule
  3. Upewnij się, że Load files for: ma wartość New Molecule
  4. Wskaż na plik 1ZIH_js.psf (ponownie, zwróć uwagę na to by użyć pliku PSF, a nie PDB) i kliknij w Load
  5. We wciąż otwartym oknie Molecule File Browser, z polem Load files for: ustawionym na 1ZIH_autopsf.psf wskaż na plik min.dcd z pobranego archiwum i kliknij w Load
  6. Powtórz poprzedni krok dla heat.dcd i następnie run1.dcd (w taki właśnie sposób "doczytuje" się kolejne ramki z zapisanej trajektorii, VMD dokłada je na końcu tworząc spójny, jednolity obraz symulacji)
• Uruchom całą symulację (przycisk "▶" w oknie VMD Main), możesz sterować suwakiem speed by spowolnić wizualizację kolejnych kroków
• Wykorzystaj narzędzie Timeline opisane wcześniej by zaobserwować otwieranie/zamykanie ostatniej pary zasad G1-U12
  
  
• Sprawdź przebieg zmian strukturalnych w czasie trwania symulacji:
  1. Wybierz z menu głównego okna VMD Extension → Analysis → RMSD Trajectory Tool
  2. W polu edycyjnym w lewym górnym rogu (jest to zaznaczenie fragmentu, którego ma dotyczyć analiza) wpisz nucleic
  3. Kliknij najpierw w przycisk ALIGN, następnie w RMSD
  4. Wybierz w menu File → Plot Data by zobaczyć wykres zmieniającego się RMSD w zależności od numeru ramki (Uwaga! Różne etapy symulacji miały różny odstęp czasowy pomiędzy zapisem współrzędnych do pliku trajektorii. Oznacza to, że między sąsiednimi ramkami nie zawsze jest ta sama różnica w czasie)