Modelowanie biomolekularne

Modelowanie w ICM
Mod. biomolekularne
Projektowanie leków
Numeryka
Wizualizacja
Inna działalność

Spis treści

1. Dynamika molekularna

Metoda dynamiki molekularnej umożliwia analizę ruchów molekularnych [1,2,3]. Pozwala symulować trajektorię ruchu jąder atomowych w czasie wychodząc od klasycznych równań Newtona dla układu wielu cząstek:

\boldsymbol{F} = m_i \cdot \boldsymbol{a}_i

-\frac{dV}{d\boldsymbol{r}_i} = m_i \cdot \frac{d^2\boldsymbol{r}_i}{dt^2}

W fizyce klasycznej można obliczyć ścisłe rozwiązania tych równań dla jednej lub dwóch cząstek. Dla układów wieloatomowych trzeba stosować metody numeryczne. Zakłada się wtedy, że ruch pojedynczego atomu w krótkim czasie można przedstawić w postaci szeregu Taylora:

\boldsymbol{r}(t+\delta t) = \boldsymbol{r}(t) + \frac{d \boldsymbol{r}}{dt}\delta t + 1/2\frac{d^2\boldsymbol{r}}{dt^2}(\delta t)^2+... = \boldsymbol{r}(t) + \boldsymbol{v}(t) \delta t + 1/2 \boldsymbol{a}(t)(\delta t)^2 + ...

gdzie \boldsymbol{r}(t), \boldsymbol{v}(t), \boldsymbol{a}(t) to odpowiednio położenie, prędkość i przyspieszenie atomu w chwili t, a \boldsymbol{r}(t+\delta t) to położenie atomu po czasie \delta t \quad.

Początkowe położenia atomów pochodzą przeważnie z krystalografii lub danych NMR (Protein Data Bank). Początkowe prędkości sa natomiast losowo przypisywane na podstawie rozkładu Maxwella-Boltzmana dla danej temperatury. Do obliczeń położeń i prędkości po czasie \delta  t \quad stosowany jest na ogół algorytm Verleta, który zakłada liniowość przyspieszenia i prędkości w jednym kroku czasowym. Krok czasowy \delta t \quad musi być więc tak dobrany, aby spełniał powyższe założenie. Dla symulacji pełnoatomowych stosuje się krok czasowy rzędu femtosekundy, jest to czas o rząd wielkości krótszy niż najszybsze ruchy oscylacyjne w białku występujące dla wiązania C-H.

Wartość energii potencjalnej układu wyznacza się traktując jądra atomowe jako cząstki klasyczne poruszające się w odpowiednio zdefiniowanym polu jąder atomowych i elektronów. Energia molekuły przedstawiona jest jako suma prostych funkcji zależnych od współrzędnych jąder atomowych cząsteczki. Funkcje wraz z ich parametrami wyznaczonymi doświadczalnie lub teoretycznie dla różnych typów atomów tworzą tzw. pole siłowe. Najprostsze pola siłowe stosują model mechaniczny, który traktuje jądra atomowe uczestniczące w wiązaniu chemicznym jako kulki połączone harmoniczną spreżyną. Dla atomów nie związanych chemicznie stosuje się wyrażenia charakteryzujące tzw. oddziaływania niewiążące: potencjał kulombowski opisujący oddziaływania elektrostatyczne i często potencjał Lennarda-Jonesa opisujący oddziaływania Van der Waalsa.

Standardowe symulacje dynamiki molekularnej odpowiadają zespołowi statystycznemu mikrokanonicznemu (NVE), w ktorym liczba cząstek, objętość i całkowita energia są stałe. Taki zespół nie odpowiada dokładnie typowym warunkom doświadczalnym. Inne możliwe zespoły to NVT, w którym utrzymywana jest za pomocą termostatu stała temperatura lub NPT, w którym przy pomocy termostatu i barostatu utrzymywane są stała temperatura i ciśnienie [4].

Analiza trajektorii i danych pochodzących z symulacji dynamiki molekularnej [3,5]:

  • Wizualizacja - dobrze jest obejrzeć jak propagują się współrzędne w czasie, aby mieć pewność, że układ poprawnie ewoluuje. Istnieje wiele pakietów do wizualizacji trajektorii ruchu (VMD, RasMol, InsightII, PyMOL) ale najprościej jest wybrać taki, który czyta format pliku generowany przez wykorzystywany pakiet do dynamiki.
  • Własności ogólne trajektorii
    • wykres energii całkowitej, kinetycznej i potencjalnej w czasie, po zrównoważeniu układu te wartości powinny fluktuować wokół stałej wartości
    • wykres położenia środka masy w czasie
    • RMSD (root mean square deviation) - średnie odchylenie standardowe wybranych, bądź wszystkich, atomów molekuły w czasie. Jako strukturę odniesienia bierze sie przeważnie konformację poczatkową układu lub konformację po jego zrównoważeniu.
    • RMSF (root mean square fluctuation) - średnie odchylenie ruchu atomów od wartości średniej. Wartość RMSF dla danego atomu może być porównywana z izotropowymi czynnikami temperaturowymi B umieszczonymi w pliku PDB zgodnie z zależnością RMSF_i = \sqrt {\frac{3B_i}{8\pi^2}} i pod warunkiem, że dana struktura krystalograficzna została otrzymana z dobrą rozdzielczościa.
    • promień bezwładności - wielkość dająca miarę zwartości badanego układu w trakcie symulacji
    • analiza odległości między wybranymi atomami lub środkami ciężkości grup atomów, analiza zmian wartości wybranych kątów płaskich i dwuściennych
    • wiązania wodorowe - analiza liczby i stabilności wiązań wodorowych w trakcie symulacji, analiza odległości donor-akceptor lub kątów donor-wodór-akceptor
  • Analiza struktury drugorzędowej - dla białek można określać liczbę \alpha \quad-helis oraz \beta \quad-kartek, a także wykreślić wykres Ramachandrana ukazujący wartości kątów torsyjnych \phi \quad i \psi \quad dla reszt aminokwasowych
  • Principal Component Analysis - do analizy trajektorii dynamiki molekularnej można stosować metodę obliczania głównych kierunków ruchu układu, polegającej na określeniu kolektywnych drgań fragmentów układu o niskich częstościach i dużych amplitudach tzw. "principal component analysis" (PCA) lub "essential dynamics" (ED) [6,7].
  • Korelacje ruchów - ruchy nawet odległych segmentów w układach molekularnych są często skorelowane, co ma istotne znaczenie biologiczne. Analiza korelacji ruchów różnych części układu może pozwolić zaproponować ich ewentualne znaczenie funkcjonalne.
  • Analiza częstościowa - z trajektorii ruchu można otrzymać widmo oscylacyjne układu dokonując transformaty Fouriera funkcji autokorelacji prędkości. Dla dużych układów widmo jest trudne do analizy, więc lepiej jest obliczać widmo drgań samych środków ciężkości wybranych segmentów lub wybranych głównych kierunków ruchu otrzymanych z PCA.


Pakiety do symulacji i analizy dynamiki molekularnej dostępne w ICM to: Amber, Charmm, NAMD, Gromos, Gromacs, Sybyl. Inne możliwości to na przykład: DL_POLY, TINKER, MMTK.

2. Dynamika brownowska

Do symulacji asocjacji różnych molekuł w roztworze, np. ligandu i białka, w obecności pola elektrostatycznego generowanego przez białko stosowana jest metoda dynamiki brownowskiej. Trajektoria ruchu cząsteczki jest wyznaczana w oparciu o algorytm Ermaka i McCammona [8]:

\boldsymbol r(t+\delta t) = \boldsymbol r(t) + \frac{D \boldsymbol F \delta t}{kT} + \boldsymbol R

gdzie:

  • D oznacza współczynnik dyfuzji czasteczki,
  • \boldsymbol F oznacza położenie i siłę elektrostatyczną działającą na ligand,
  • T - temperaturę,
  • k - stałą Boltzmanna,
  • \delta t \quad - krok czasowy dynamiki,
  • \boldsymbol R - przesunięcie stochastyczne o własnościach: <\boldsymbol R> = 0 i <\boldsymbol R \cdot \boldsymbol R> = 6D\delta t. Przesunięcie \boldsymbol R jest dodawane, aby symulować zderzenia z cząsteczkami wody.

Symulacje dynamiki brownowskiej przeprowadzane są w następujący sposób [9]. Najpierw obliczany jest potencjał elektrostatyczny wewnątrz i na zewnątrz białka poprzez rozwiązanie równania Poissona-Boltzmanna (patrz podrozdział "Właściwości elektrostatyczne"). Siatka potencjału jest używana jako plik wejściowy do dalszych obliczeń. Następnie symulowany jest ruch ligandu w polu elektrostatycznym białka w oparciu o powyższe równanie; w każdym kroku czasowym obliczana jest siła elektrostatyczna działająca na ligand, dodawane jest do niej losowe przesunięcie i w ten sposób znajdowane jest nowe położenie ligandu. Aby uniknąć zderzeń białka z ligandem wszystkie kroki, które prowadzą do zachodzenia na siebie tych cząsteczek są odrzucane i dla takich przypadków jest losowane nowe \boldsymbol R.

Dynamika Brownowska

Rysunek: Dynamika brownowska. S oznacza sferę, na której losowo wybierane jest położenie startowe ligandu. Ligand porusza się po trajektorii wyznaczonej przez siły oddziaływania z polem elektrostatycznym oraz przez siły stochastyczne. R i U oznaczają odpowiednio sfery "reakcji" i "ucieczki", na których kończą się trajektorie klasyfikowane jako prowadzące albo nie prowadzące do związania ligandu z białkiem. W modelu tym ligand jest przeważnie reprezentowany jako sfera o całkowitym ładunku lub jako sztywny zbiór atomów z ładunkami cząstkowymi.

Symulacje, w których ligand startuje ze sfery S są przeprowadzane wielokrotnie tak, aby mieć wiarygodną informację statystyczną na temat częstości utworzenia kompleksu ligand-białko. Do obliczeń musi być także odpowiednio zdefiniowane tzw. kryterium reakcji, czyli obszar, w którym zakładamy, że ligand jest już związany z białkiem. Ponieważ w dynamice brownowskiej nie interesują nas ruchy wewnętrzne molekuł, a jedynie ich wzajemne oddziaływanie stosuje się krok czasowy \delta t \quad o wiele dłuższy niż w dynamice molekularnej (> 1 ps).

Zbiór niezależnych trajektorii pochodzących z dynamiki brownowskiej jest analizowany pod względem częstości utworzenia kompleksu ligand-białko. Trajektorie, które prowadzą zgodnie z kryterium reakcji do asocjacji molekuł mogą być dokładniej analizowane pod względem szlaku wędrówki ligandu w strone białka. Metoda dynamiki brownowskiej pozwala obliczyć stałe asocjacji molekuł, określić mechanizm asocjacji, a także zbadać czy siły elektrostatyczne mają istotny wpływ w procesie asocjacji molekuł.

Do symulacji metodą dynamiki brownowskiej może być wykorzystane oprogramowanie UHBD, SDA, a także MacroDox.

3. Właściwości elektrostatyczne

Istotnym elementem badań może być informacja na temat elektrostatycznych właściwości układu. Wiele pakietów pozwala wyznaczać pole elektrostatyczne wewnątrz i na zewnątrz różnego rodzaju molekuł. Molekuła reprezentowana jest jako zbiór ładunków punktowych umieszczonych w ośrodku ciągłym o niskiej stałej dielektrycznej przeważnie 2-4, chociaż wartości 10-20 są również stosowane przez niektórych autorów. Otoczenie, woda oraz rozpuszczone w niej jony, modelowane jest w sposób ciągły poprzez wysoką stałą dielektryczną równą 78-80 i odpowiednią siłę jonową.

Potencjał elektrostatyczny wokół molekuły jest obliczany na podstawie modelu Poissona-Boltzmanna (PB), którego podstawą jest jedno z równań Maxwella, równanie Poissona:

\nabla \cdot \epsilon(\boldsymbol{r}) \nabla \Phi (\boldsymbol{r}) = - 4 \pi \rho(\boldsymbol{r})

gdzie:

  • \Phi \quad - potencjał pola elektrostatycznego,
  • \epsilon \quad - wartość stałej dielektrycznej wewnątrz cząsteczki i w zewnętrznym środowisku wodnym,
  • \rho \quad - całkowita gęstość ładunku w badanym układzie uwzględniająca rozkład nieruchomych ładunków w molekule i gęstość ładunku jonów w roztworze.

Gęstość jonów w roztworze przybliżana jest rozkładem Boltzmanna i po wstawieniu do równania Poissona otrzymujemy nieliniowe równanie Poissona-Boltzmanna.

 \nabla \cdot \epsilon (\boldsymbol{r}) \nabla \Phi (\boldsymbol{r}) = - 4 \pi \rho_{wewn} (r) + \lambda (r) \kappa_0^2 sinh (\Phi /kT),

gdzie:

  • \rho_{wewn} \quad - rozkład ładunków w cząsteczce,
  • \kappa^2_0 \quad - stała Debyea-Hückela, związaną z siłą jonową roztworu,
  • \lambda \quad - parametr równy zeru wewnątrz cząsteczki i równy stałej dielektrycznej wody w obszarze dostępnym dla jonów rozpuszczonych w wodzie.

Równanie to dla białek i kwasów nukleinowych, a także ich kompleksów jest rozwiązywane numerycznie. W modelu tym potencjał elektrostatyczny wewnątrz i wokół molekuły pochodzący z rozwiązania tego równania wyznaczany jest na trójwymiarowej przeważnie regularnej siatce.

Poisson-Boltzmann

Rysunek: Model Poissona-Boltzmanna. Cząsteczka (przedstawiona jako czerwony obszar) jest scharakteryzowana rozkładem ładunku, (w praktyce wykorzystywany jest rozkład punktowych ładunków atomowych) oraz stałą dielektryczną, \epsilon \quad, o małej wartości. Środowisko wodne (niebieski obszar) opisane jest poprzez parametr \kappa \quad, związany z siłą jonową, oraz poprzez stałą dielektryczną o dużej wartości. Potencjał elektrostatyczny, jest obliczany w punktach trójwymiarowej siatki.

Zastosowanie modelu Poissona-Boltzmanna obejmuje między innymi wizualizacje potencjału elektrostatycznego wokól molekuły i określanie obszarów dodatnio i ujemnie naładowanych, badanie procesów dyfuzji ligandu do białka, teoretyczne miareczkowanie jonizowalnych reszt aminokwasowych w białkach, obliczanie elektrostatycznego wkładu do energii solwatacji oraz elektrostatycznego wkładu do energii wiązania cząsteczek. Więcej informacji na temat modelu Poissona-Boltzmanna można znaleźć w artykułach [10,11,12].

Pakiety służące do obliczeń potencjału elektrostatycznego to DelPhi (np. jako część pakietu InsightII firmy Accelrys), UHBD, APBS, GRASP.

4. Komputerowe miareczkowanie białek

Dane o strukturze białka pochodzące z krystalografii nie zawierają informacji o położeniach jąder atomów wodoru. Przygotowując cząsteczkę do symulacji metodą dynamiki molekularnej należy zatem uzupełnić jej strukturę o brakujące protony. Położenie większości z nich jest stabilne i dobrze ustalone względem sąsiednich ciężkich atomów. Co za tym idzie, dołączenie ich do symulowanych cząsteczek jest stosunkowo proste, gdyż opiera się na wykorzystaniu standardowych topologii aminokwasów, które są zdefiniowane w ramach pola siłowego.

Niektóre aminokwasy posiadają jednak zdolność przyłączania i odłączania protonów zależną od zmian pH otoczenia. Aminokwasy o charakterze kwasowym, takie jak kwas asparaginowy (ASP), kwas glutaminowy (GLU), tyrozyna (TYR) oraz grupa karboksylowa na C-końcu łańcucha polipeptydowego, w miarę wzrostu pH mają tendencję do odłączania protonu i pozostawania w postaci jonu ujemnego. Aminokwasy zasadowe - lizyna (LIS), arginina (ARG), histydyna (HIS) oraz grupa aminowa N-końca polipeptydu - wraz ze spadkiem pH przyłączają proton i stają się jonami dodatnimi. Stan uprotonowania poszczególnych aminokwasów w białku, w określonym pH warunkuje zatem przestrzenny rozkład ładunków, co ma istotny wpływ na wyniki symulacji. Poprawne wyznaczenie go jest trudnym i wciąż nie do końca rozwiązanym zagadnieniem. Jedną z możliwości jest wykorzystanie metody bazującej na badaniu właściwości elektrostatycznych białka w oparciu o równanie Poissona-Boltzmanna.

4.1 Teoretyczne podstawy miareczkowania

Każdy rodzaj cząsteczek wymieniających protony z otoczeniem jest scharakteryzowany przez wartość pK oznaczającą takie pH środowiska, przy którym średnio połowa molekuł tego rodzaju pozostaje w postaci zjonizowanej. Otoczenie białkowe powoduje, że pK_{a} aminokwasów w cząsteczce białka jest inne niż pK_{0} aminokwasów swobodnych. Składa się na to przede wszystkim oddziaływanie grupy miareczkowalnej z ładunkami białka oraz oddziaływanie między poszczególnymi grupami miareczkowalnymi. Oba te wkłady brane są pod uwagę w teoretycznym wyznaczaniu pK_a [18, 19].

Oddziaływanie grupy miareczkowalnej z białkiem uwzględnia się wyliczając wartość pK_a "wewnętrznego" - pK_w. Odpowiada ono wartości pK_a w białku, w którym wszystkie grupy miareczkowalne są w stanie neutralnym. Podstawą do jego obliczenia jest następujący cykl termodynamiczny (ilustrujący przypadek aminokwasu o charakterze kwasowym):

Cykl td.png

Wartość pK_0 jest znana z badań doświadczalnych, a wartości \Delta G_1 i \Delta G_2 odpowiadają entalpiom swobodnym przeniesienia odpowiednio nienaładowanego i naładowanego aminokwasu z wody do otoczenia białkowego. Na potrzeby obliczeń zakłada się, że wartości tych entalpii wynikają jedynie z oddziaływań elektrostatycznych. Pozwala to obliczyć wartość \Delta \Delta G = \Delta G_2 - \Delta G_1 na podstawie różnic odpowiednich energii elektrostatycznych:

\Delta \Delta G = (W_{prot}^- - W_{prot}^0) - (W_{wat}^- - W_{wat}^0) (1)

gdzie W_{prot/wat}^{-/0} oznaczają energię elektrostatyczną aminokwasu naładowanego lub obojętnego w białku bądź w wodzie. Do ich obliczenia można wykorzystać np. metodę Poissona-Boltzmana.

Znajomość pK_0 oraz \Delta \Delta G pozwala wyznaczyć pK_w na podstawie wzoru:

pK_w=pK_{0}-\gamma\Delta\Delta G/2.303RT (2)

gdzie \gamma wynosi +1 dla grup zasadowych oraz -1 dla kwasowych.

Wzajemne oddziaływanie grup miareczkowalnych oblicza się zakładając, że jest ono jedynie natury elektrostatycznej. Stąd niezerowy wkład powstaje tylko w przypadku, gdy obie grupy są naładowane. Jego wartość, \phi_{ij} można wyliczyć z następująco:

\phi_{ij}=W_{ic,jc}-W_{ic,jn}-W_{in,jc}+W_{in,jn} (3)

gdzie W_{ic/in,jc,jn} oznacza energię elektrostatyczną białka w przypadku, gdy grupy i oraz j są naładowane (ic, jc), bądź neutralne (in, jn). Wartości tych energii można również obliczyć korzystając z metody Poissona-Boltzmanna.

Stan protonacji białka. Rozważmy białko zawierające N grup mogących wymieniać protony z otoczeniem. Stan grupy i można opisać podając jej ładunek x_{i}. Dla grup kwasowych x\in \{-1,0\}, a dla zasadowych x\in \{0,1\}. Stan całego białka scharakteryzowany jest zatem przez wektor protonacji \vec x o wymiarze N. Różnica entalpii swobodnej białka w danym pH i w stanie \vec x w stosunku do białka, w którym wszystkie grupy pozostają neutralne wynosi:

\Delta G(pH,\vec x)=2.303RT \sum_{i=1}^{N} x_{i}(pH-pK_{w,i})+ \sum_{i=1}^{N}\sum_{j=i+1}^{N} x_{i}x_{j} \phi_{ij} (4)

Pierwszy człon opisuje zmianę entalpii swobodnej wynikającą z oddziaływania elektrostatycznego zjonizowanych grup z resztą białka oraz przejścia protonów z lub do środowiska. Drugi człon opisuje energię oddziaływania elektrostatycznego pomiędzy zjonizowanymi grupami.

Znalezienie najbardziej prawdopodobnego stanu protonacji białka w danym pH sprowadza się do wyznaczenia wektora protonacji odpowiadającego minimalnej entalpii swobodnej. Obliczenie entalpii dla każdego możliwego stanu protonacji jest jednak praktycznie nierealne. Dla białka, w którym jest np. 30 grup miareczkowalnych istnieje 2^{30} różnych wektorów protonacji. W praktyce zatem stanu o najniższej entalpii swobodnej poszukuje się np. metodą Monte Carlo.

4.2 Etapy komputerowego miareczkowania

  • Dodanie wszystkich możliwych protonów do struktury krystalograficznej białka. Jego wynikiem jest struktura, w której wszystkie reszty zasadowe (ARG, LYS, HIS, CYS, N-koniec) są w postaci naładowanej dodatnio, a wszystkie reszty kwasowe (ASP, GLU, TYR, C-koniec) są obojętne. W ten sposób wyznaczona zostaje lokalizacja protonów miareczkowalnych. Na tym etapie należy podjąć decyzję, do którego z dwóch atomów tlenu grup COO obecnych w białku dostawić proton. Omawiana metoda nie pozwala bowiem na jednoczesne miareczkowanie tych grup w obu miejscach.
    Krok ten można wykonać za pomocą dowolnego programu do modelowania białek. Należy zwrócić jednak uwagę, aby format, w jakim zapisana zostanie struktura wynikowa (w szczególności nazwy typów poszczególnych atomów wodoru) był rozpoznawany przez program wykonujący obliczenia Poissona-Boltzmanna.
  • Obliczenie potencjałów oddziaływania (\Delta \Delta G - wzór (1) oraz \phi_{ij} - wzór (3)) metodą Poissona-Boltzmanna. Do ich wyznaczenia niezbędne jest przypisanie ładunków cząstkowych oraz promieni poszczególnym atomom. W tym celu musi być znana parametryzacja dla obu stanów aminokwasów miareczkowalnych (uprotonowaego i zdeprotonowanego). Z reguły pliki zawierające odpowiednie parametry są częścią pakietów służących do miareczkowania.
    Przed wykonaniem tego kroku istotne jest podjęcie decyzji, który z dwóch protonów grup histydynowych uznać za miareczkowalny.
  • Wyznaczenie pK_a aminokwasów oraz wektorów protonacji. Znajomość pH otoczenia, pK_0 i \Delta\Delta G dla poszczególnych aminokwasów miareczkowalnych oraz potencjałów oddziaływania \phi_{ij} pozwala wyznaczyć wartość entalpii swobodnej dla dowolnego wektora protonacji (wzór 4). Dzięki temu, wykorzystując np. metodę Monte Carlo (w przypadku, gdy wstępne obliczenia były wykonane programem UHBD, można zastosować do tego celu program \verb+hybrid+), można uzyskać następujące wyniki:
    • wektory protonacji odpowiadające minimalnym wartościom entalpii swobodnej w zadanym pH
    • procentowy stan obsadzeń poszczególnych grup miareczkowalnych w danym pH
    • wartości pK_a dla poszczególnych aminokwasów (uzyskuje się je prowadząc obliczenia dla całego przedziału wartości pH w poszukiwaniu takiej, dla której procentowy stan obsadzeń bliski będzie 50%)

Często celem obliczeń jest uzyskanie stanu protonacji białka do dalszych symulacji np. metodą dynamiki molekularnej. W takiej sytuacji warto przyjąć stan protonacji białka opisany którymś ze znalezionych wektorów protonacji, a nie ustalać stan poszczególnych aminokwasów w oparciu o wyznaczone wartości pK_a.

4.3 Dostępne programy

Programy:

Serwery wykonujące automatyczne miareczkowanie:

Linki do stron związanych z miareczkowaniem:

5. Modelowanie nanoskalowych obiektów

Dla obiektów o rozmiarach powyżej kilkudziesięciu nanometrów symulacje pełnoatomowe są zbyt kosztowne obliczeniowo. Modelowanie tak dużych systemów wymaga często dostosowania do tych celów oprogramowania, a także wielokrotnego zredukowania stopni swobody badanego układu. Przykłady nanoskalowych obiektów biomolekularnych to rybosom, wirusy, polimeraza RNA czy włókna chromatyny.

Obliczenia dynamiki molekularnej dla takich układów, które zawierają kilkaset tysięcy atomów są zbyt czasochłonne i mogą być prowadzone tylko w krótkiej skali czasowej. Aby zmniejszyć ilość stopni swobody otoczenia w symulacjach dynamiki molekularnej można zastosować ciągły model wody na przykład stosując uogólniony model Borna [13] zaimplementowany w pakietach programowym Amber, Charmm. Można także zastosować techniki uproszczonej gruboziarnistej dynamiki molekularnej redukując znacznie liczbę stopni swobody molekuły i uzwzględniając środowisko wodne jedynie w postaci efektywnych potencjałów lub stosując dynamikę langevinowską. Model polega na zastąpieniu całej reszty aminokwasowej jednym bądź kilkoma pseudo-atomami. Przegląd takich metod i stosowanych dla nich pól siłowych dostępny jest na przykład w artykule [14].

Do obliczeń potencjału elektrostatycznego nanomolekuł [15] i do rozwiązywania równania Poissona-Boltzmanna dla dużych układów bardzo dobry jest program Adaptive Poisson-Boltzmann Solver [16] APBS, gdyż pozwala na obliczenia równoległe. Przy odpowiedniej liczbie procesorów umożliwia obliczanie siatki potencjału dla nanoskalowych układów w dobrej rozdzielczości (0.2-0.3\dot {A}). Siatki mogą być zapisane w formacie OpenDX, UHBD lub AVS. Mozna je oglądać np. używając darmowych programów do wizualizacji VMD lub OpenDX, PyMOL lub dostępnego w ICM UW pakietu AVS.

Potencjal podjednostki 30S

Rysunek: Potencjał elektrostatyczny na powierzchni dostępności wody podjednostki 30S rybosomu bakteryjnego zawierającej około 90 tysięcy atomów [17]. Skala potencjału w jednostkach kT/e.

Więcej informacji dostępnych jest na stronie www.trylska.com, a w sprawie konsultacji merytorycznych dotyczących modelowania maszyn molekularnych (doboru modelu, oprogramowania) prosimy pisać na adres: joanna@icm.edu.pl

6. Monte Carlo

Alternatywą dla dynamiki molekularnej może być prowadzenie symulacji metodami stochastycznymi, np. metodą Monte Carlo. Metoda ta doskonale się nadaje do przeszukiwania przestrzeni konfiguracyjnej układu.

6.1 Teoretyczne podstawy metody Monte Carlo

Celem metody Monte Carlo jest praktyczne obliczenie wielowymiarowych całek statystycznych, gdyż w przypadku wielowymiarowym rozkład losowy jest szybciej zbieżny numerycznie niż rozkład regularny, czyli np całkowanie metodą trapezów. Wykorzystanie metody Monte Carlo oznacza, że rozwiązujemy pewien problem przedstawiając go jako parametr jakiejś hipotetycznej populacji, której reprezentatywną próbkę generujemy za pomocą liczb losowych. Z otrzymanej próbki wartość parametru szacujemy statystycznie. Problemem samym w sobie jest jak najefektywniej stworzyć reprezentatywną próbkę, tzn. jaka najmniejsza próbka populacji może być uważana za reprezenatywną. Zazwyczaj rozwiązanie problemu przedstawia się w postaci całki oznaczonej z pewnej funkcji liczb losowych o rozkładzie płaskim. W przypadku jednowymiarowym można zapisać:

I\ =\ \frac{1}{a-b}\int _{a}^{b}\mathit{dx}f(x)

ale rzecz jasna powyższą całkę w prosty sposób można uogólnić do wielu wymiarów. W przypadku losowania liczb z rozkładu płaskiego, wartość całki I jest po prostu wartością oczekiwaną funkcji f, czyli średnią z próby, I = <f(x)>_{proba}, przy czym im gęściej próbkujemy, tym dokładniej ją wyznaczamy. Jednak w przypadku układów fizycznych, np: molekularnych, zadanych w pewnej przestrzeni fazowej, nawet bardzo gęste próbkowanie może prowadzić do artefaktów w oszacowaniu wartości całki. Problem polega na tym, by próbkować z tych rejonów przestrzeni fazowej, które dają największy wkład do całki, czyli losować liczby z innego niż jednostajny rozkładu. Stosunkowo popularną metodą zwiększającą dokładność w wyznaczaniu całek, czyli redukującą wariancję, i spełniającą powyższy warunek, jest metoda średniej ważonej. Polega ona na próbkowaniu przestrzeni wg arbitralnie wybranego rozkładu, \rho (x), gdzie \rho (x) jest dowolną unormowaną (\int _{a}^{b}\mathit{dx}\rho (x)=1) funkcją x. Wówczas wartość całki otrzymamy jako średnią wartość funkcji f(x) ważonej tym rozkładem, dokładniej:

I\ =\ \int _{a}^{b}\mathit{dx}\frac{f(x)}{\rho (x)}\rho (x)

czyli  I=<f(x)/\rho (x)>_{\rho}=N^{-1}\sum_{i=1}^{N} f(x)/\rho (x) , gdzie N jest liczbą prób, czyli losowań. Jeżeli \rho (x) będzie innym niż jednostajny rozkładem, to taką metodę całkowania nazwiemy próbkowaniem ważonym. Do budowy kanonicznego zespołu statystycznego, w którym próbkujemy przestrzeń zgodnie z mikroskopowym rozkładem Boltzmana, wykorzystał ją Metropolis.

6.2. Próbkowanie ważone, symulacje Monte Carlo w kanonicznym zespole statystycznym

Wartość średnią pewnej obserwabli Y w kanonicznym zespole statystycznym (czyli zespole o stałej liczbie cząstek N, stałej objętości v oraz stałej temperaturze T) wyznaczona jest przez całkę:

\langle Y\rangle =\int d\Gamma Y\rho _{\mathit{NvT}}(\Gamma ),

gdzie \rho_{\mathit{NvT}} jest rozkładem prawdopodobieństw w tym zespole. Numeryczne całkowanie metodą średniej ważonej polega na wyznaczeniu wartości Y\cdot \rho _{\mathit{NvT}} losując z dowolnego rozkładu \varrho . Przyjmując \rho =\rho_{\mathit{NvT}} mamy szansę dość dobrze oszacować Y, gdyż zawyczaj w takim wypadku średnia wartość Y po zespole jest równa średniej jego wartości w symulacji. Należy zatem wygenerować taki ciąg stanów, których rozkład będzie odpowiadał wymaganemu rozkładowi prawdopodobieństw.

W tym celu w symulacji generuje się skończony, regularny łańcuch Markowa, w którym prawdopodobieństwo przejścia do następnego stanu zależy tylko od stanu bezpośrednio go poprzedzającego. Dla takiego łańcucha rozład graniczny (na mocy twierdzenia o ergodyczności) jest jedynym rozkładem stacjonarnym. Dodatkowo, jeśli \vec{\rho } oznacza dowolny początkowy rozkład prawdopodobieństwa, \vec{\alpha } - graniczny, a \hat{P} (1-krokową) macierz przejść między stanami, to zachodzi: \lim _{n\to \infty }\vec{\rho
}\hat{P}^{n}=\vec{\alpha }.

W przypadku układów molekularnych (np: cieczy) elementy P_{ij} macierzy przejścia nie są znane, natomiast znany jest rozkład graniczny łańcucha. Chcąc wyznaczyć elementy macierzy przejścia można posłużyć się dodatkowym warunkiem (poza warunkiem stacjonarności) tzw. równowagi szczegółowej, który mówi, że jeśli \rho _{i} jest prawdopodobieństwem znalezienia układu w stanie i, a P_{ij} prawdopodobieństwem przejścia ze stanu i do stanu j, to \rho _{i}P_{\mathit{ij}}=\rho _{j}P_{\mathit{ji}}.

Rozwiązanie Metropolis'a na postać stochastycznej (  \sum_{j} P_{ij}=1 )macierzy przejść, spełniającej warunek równowagi szczegółowej, jest następujące:


P_{ij} =  \left\{ \begin{matrix}
\alpha_{ij} & jesli & \rho_i \ge \rho_j,i \ne
j \\ 
\alpha_{ij}\frac{\rho_j}{\rho_i} & jesli & \rho_i < \rho_j,i \ne j \\
1-\sum_{i\ne j}P_{ij} & jesli  & i=j
\end{matrix} \right.

gdzie \alpha jest symetryczną, stochastyczną macierzą graniczną. Uściślając powyższe dla zespołu kanonicznego, gdzie \rho
\sim \operatorname{e}^{-\beta V} mamy: \vec{\alpha }=\rho
_{\mathit{NvT}}, zaś warunek równowagi szczegółowej można przepisać w postaci:

P_{\mathit{ij}}\operatorname{e}^{-\beta
V_{i}}\ =\ P_{\mathit{ji}}\operatorname{e}^{-\beta V_{j}}

Zastosowanie algorytmu Metropolisa podczas symulacji sprowadza się do następującego mechanizmu(dla przejścia ze stanu i do j):

  • stan układu reprezentujemy jego energią wewnętrzną; generujemy nową konfigurację układu;
  • obliczamy różnicę energii między konfiguracjami, \Delta V_{\mathit{ij}},
  • jeśli energia nowej konfiguracji jest mniejsza bądź równa energii poprzedniej konfiguracji

(\Delta V_{\mathit{ij}}\le 0), układ przechodzi do nowej konfiguracji,

  • jeśli energia nowej konfiguracji jest większa (\Delta
V_{\mathit{ij}}>0), układ przechodzi do nowej konfiguracji z pewnym prawdopodobieństwem, wynikającym z rozkładu Boltzmana:

\frac{\rho _{j}}{\rho
_{i}}\ =\ \frac{Z_{\mathit{NvT}}^{-1}\operatorname{e}^{-\beta
V_{j}}}{Z_{\mathit{NvT}}^{-1}\operatorname{e}^{-\beta
V_{i}}}\ =\ \operatorname{e}^{-\beta
(V_{j}-V_{i})}\ =\ \operatorname{e}^{-\beta \Delta V_{\mathit{ij}}}

(czyli jeśli wylosowana z rozkładu płaskiego [0,1] liczba jest mniejsza niż \rho _{j}/\rho _{i}). Jak widać, próbkowanie przestrzeni nie wymaga wyznaczania konfiguracyjnej sumy statystycznej.

Nowe położenie cząstki i, otrzymuje się postępując zgodnie z poniższym wzorem:


{\vec{r}}_{i,\;\mathit{nowe}}\;=\;{\vec{r}}_{i,\;\mathit{stare}}\;+\;\Delta
r_{\mathit{MAX}}\cdot \;(2\;u-1)

gdzie u jest liczbą wylosowaną z rozkładu płaskiego [0,1], \Delta r_{\mathit{MAX}} jest największą dopuszczalną zmianą którejkolwiek ze współrzędnych wektora położenia .

6.3. Symulacje Monte Carlo w izotermiczno-izobarycznym zespole statystycznym

Korzystając z metody Monte Carlo dość łatwo można przejść do innych zespołów statystycznych. W przypadku zespołu izotermiczno-izobarycznego, czyli zespołu, który składa się z N cząstek, w którym utrzymywana jest stała temperatura T oraz stałe ciśnienie p, rozkład prawdopodobieństw ma postać:


\rho\ =\ Z_{\mathit{NpT}}^{-1}\operatorname{e}^{-\beta
(V(X)+\mathit{pv})},


Z_{\mathit{NpT}}=\int
\mathit{dv}\operatorname{e}^{-\beta \mathit{pv}}\ \int
\mathit{dX}\operatorname{e}^{-\beta V(X)},

gdzie p i v to ciśnienie i objętość układu, Z_{NpT} - konfiguracyjna suma statystyczna, X symbolizuje wektory położeń wszystkich cząstek układu. Nowa konfiguracja w tym zespole statystycznym powstaje poprzez ruch cząstek układu i zmianę jego objętości.

W szczególnym, a zarazem dość popularnym przypadku, kiedy układ jest reprezentowany przez pudło sześcienne wypelnione cząstkami, charakterystyczną wielkością jest długość krawędzi pudła, L. Wygodnie jest wówczas przejść do względnych położeń cząstek: {\vec{r}_{i}}^{*}=\vec{r}_{i}/L, gdzie L jest długością pudła, L=v^{1/3} . Rozkład graniczny jest wtedy proporcjonalny do wyrażenia:

\rho \ =\ Z_{\mathit{NpT}}^{-1}\operatorname{e}^{-\beta V(X)-\beta
\mathit{pv}+N\ln (v)}

Ostatni składnik w wykładniku wynika właśnie z przeskalowania współrzędnych, sumę statystyczną można bowiem zapisać:

Z_{\mathit{NpT}}\ =\ \int \mathit{dv}v^{N}\operatorname{e}^{-\beta
\mathit{pv}}\ \int \mathit{dX}^{*}\operatorname{e}^{-\beta
V(X^{*})}

gdzie “*” wskazuje na współrzędne przeskalowane. Zatem algorytm Metropolis'a dla przejścia ze stanu i do stanu j w tym zespole będzie wyglądał następująco:

  • stan układu reprezentujemy energią postaci: E_{i}=V_{i}(X^{*})+\mathit{pv}_{i}-\beta ^{-1}N\ln (v_{i}) ; generujemy nową konfigurację układu;
  • obliczamy różnicę wartości energii między konfiguracjami, \Delta E_{\mathit{ij}}=\Delta V_{\mathit{ij}}+p\Delta
v_{\mathit{ij}}-\beta ^{-1}N\ln (v_{j}/v_{i})
  • jeśli energia nowej konfiguracji jest mniejsza bądź równa energii poprzedniej konfiguracji

(\Delta E_{\mathit{ij}}\le 0), układ przechodzi do nowej konfiguracji,

  • jeśli energia nowej konfiguracji jest większa (\Delta
E_{\mathit{ij}}>0), układ przechodzi do nowej konfiguracji z prawdopodobieństwem wynikającym z rozkładu Boltzmana, \rho
_{j}/\rho _{i}=\operatorname{e}^{-\beta \Delta E_{\mathit{ij}}}.

Nową objętość układu uzyskuje się analogicznie jak położenie, czyli L_{n}=L_{s}+\Delta L_{\mathit{MAX}}\cdot
\;(2u-1).

6.4. Symulacje Monte Carlo w Wielkim Zespole Kanonicznym

W wielkim zespole statystycznym układ znajduje się w stałej temperaturze T, stałej objętości v oraz ma stały potencjał chemiczny, \mu. Rozkład prawdopodobieństw ma postać:


\rho \ =\ Q_{\mu
\mathit{vT}}^{-1}\operatorname{e}^{-\beta H(X,P)+\beta \mu N},


Q_{\mu
\mathit{vT}}\ =\ \sum _{N=0}^{\infty
}{\frac{1}{N!}}\frac{1}{h^{3N}}\;\operatorname{e}^{\beta \mu N}\;\int
\mathit{dX}\mathit{dP}\operatorname{e}^{-\beta H(X,P)}\ =\ \sum
_{N=0}^{\infty }{\frac{1}{N!}}\frac{1}{\Lambda
^{3N}}\;\operatorname{e}^{\beta \mu N}\;\int
\mathit{dX}\operatorname{e}^{-\beta V(X)}

dzie H(X,P) jest hamiltonianem układu, X i P to wektory położeń i pędów wszystkich cząstek układu (energia kinetyczna będzie zależała od N, dlatego pędy też należy uwzględnić), h jest stałą Planck'a, zaś \Lambda tzw. termiczną długością fali de Broglie'a, \Lambda =(\beta h^{2}/2\pi m)^{1/2}. Podobnie jak w poprzednim przypadku wygodnie jest wprowadzić współrzędne przeskalowane. Wówczas suma statystyczna ma postać Q_{\mu \mathit{vT}}=\sum _{N=0}^{\infty
}(N!)^{-1}(\operatorname{e}^{\beta \mu }\Lambda ^{-3})^{N}v^{N}\;\int
\mathit{dX}^{*}\operatorname{e}^{-\beta V(X^{*})}, zaś rozkład graniczny jest równy:


\rho \ =\ Q_{\mu \mathit{vT}}^{-1}\operatorname{e}^{-\beta
V(X^{*})+\beta \mu N-\ln (N!)-3N\ln (\Lambda )+N\ln (v)}

Jeden z przykładów algorytmu realizacji wielkiego zespołu kanonicznego [norman] jest następujący:

  • generujemy nową konfigurację przestrzenną układu (tj. przesuwamy minimum jedną z cząstek)
  • usuwamy jedną z cząstek
  • tworzymy jedną cząstkę w przypakowym położeniu w układzie

Zmiany położeń cząstek realizowane są jak w zespole kanonicznym. Jeśli cząstka zostaje usunięta z układu wówczas stosunek prawdopodobieństw wynosi:

\frac{\rho _{j}}{\rho _{i}}\ =\ \operatorname{e}^{-\beta \Delta
V_{\mathit{ij}}-\beta \mu +\ln (N)+\ln (\Lambda ^{3})-\ln
(v)}\ =\ \operatorname{e}^{-\beta \Delta V_{\mathit{ij}}-\beta \mu +\ln
(N\Lambda ^{3}/v)}\ =\ \operatorname{e}^{-\beta \Delta
D_{\mathit{ij}}}

gdzie D_{ij} można zdefiniować jako energię usunięcia (destrukcji) jakiejś jednej cząstki. Usunięcie cząstki następuje dalej zgodnie z algorytmem Metropolis'a, jeśli \Delta D_{\mathit{ij}}\le 0 jest akceptowane, jeśli \Delta D_{\mathit{ij}}>0 jest akceptowane z pewnym prawdopodobieństwem. Analogicznie, jeśli cząstka zostaje wprowadzona do układu, stosunek prawdopodobieństw wynosi:


\frac{\rho _{j}}{\rho _{i}}\ =\ \operatorname{e}^{-\beta \Delta
V_{\mathit{ij}}+\beta \mu -\ln (N+1)-\ln (\Lambda ^{3})+\ln
(v)}\ =\ \operatorname{e}^{-\beta \Delta V_{\mathit{ij}}+\beta \mu +\ln
(v/((N+1)\Lambda ^{3}))}\ =\ \operatorname{e}^{-\beta \Delta
C_{\mathit{ij}}}

gdzie tym razem C_{ij} można zdefiniować jako energię wprowadzenia (kreacji) jakiejś jednej cząstki. Kreacja cząstki następuje dalej zgodnie z algorytmem Metropolis'a: jeśli \Delta
C_{\mathit{ij}}\le 0 jest akceptowane, jeśli \Delta
C_{\mathit{ij}}>0 jest akceptowane z pewnym prawdopodobieństwem. Powyższe równania można przepisać wprowadzając pojęcie aktywności układu, a=\operatorname{e}^{\beta
\mu }/\Lambda ^{3}. Wtedy:


\Delta D_{\mathit{ij}}\ =\ \Delta
V_{\mathit{ij}}\;-\;\beta ^{-1}\ln (N/(a\;v))\;,


\Delta
C_{\mathit{ij}}\ =\ \Delta V_{\mathit{ij}}\;-\;\beta ^{-1}\ln
(a\;v/(N+1))\;.


Symulacje Monte Carlo można realizować w programie Charmm lub jego wersji zawartej w oprgramowaniu Accelrys - CHARMm, dostępnej w ICM.

6.5. Generatory liczb losowych

Podstawą realizacji metody Monte Carlo jest możliwość generowania długich ciągów prób losowych. Przy tym chcemy, by generator liczb losowych był dobry, tzn. nie był obciążony błędami aparaturowymi, przechodził jak najwięcej testów losowości, miał jak najdłuższy okres powtarzalności (najczęsciej są stosowane generatory liczb pseudolosowych, z samej ich definicji charakteryzujące się powtarzalnością). Jednym z najprostszych i dobrych, choć wcale nie najgorszym ciągiem liczb pseudolosowych z przedziału [0,1] jest ten określony następująco:


X_{n}= (reszta\ z\ dzielenia\ 16807\cdot A\cdot X_{n-1}\ przez\ A)\ /\ A

gdzie A = 2147483647,  X_{0} = C/A, a C jest dowolną liczbą z przedziału [1, A-1] [20, 21]. Ogólnie, generatory kongruencyjne liniowe [22] można zapisać w postaci:

X_{n}\ =\ reszta\ z\ dzielenia \ a\cdot X_{n-1}\ przez \ m

lub

X_{n}\ =\ reszta\ z\ dzielenia \ a\cdot X_{n-1}\ przez \ m ,

gdzie m=2^{k}, a k jest liczbą bitów z jaką reprezentowana jest liczba całkowita w komputerze. Takie generatory można znaleźć w bibliotece C (np popularny drand48) i Fortran (ran). Inne generatory:

6.6. Generatory rozkladu normalnego o zadanej średniej i wariancji

Pierwszy z przykladów generatora rozkladu normalnego opiera się na Centralnym Twierdzeniu Granicznym. W myśl CTG zmienna losowa, będąca nieskończona sumą niezależnych zmiennych losowych pochodzących z identycznych rozkladów, będzie miala rozklad normalny. W praktyce, nieskończoną sumę zastępujemy skończoną sumą N zmiennych, przy czym wystarczy już N rzędu kilkadziesiąt, czy kilkanaście. Czyli zmienną Y o standaryzowanym (średniej równej 0 i wariancji równej 1) rozkladzie normalnym tworzymy następująco:


Y=\frac{ \frac{1}{N} \sum_{i=1}^{N} X_{i} - \mu }{ \frac{\sigma}{\sqrt{N}} }

gdzie m i s są odpowiednio średnią i wariancją każdej ze zmiennych X_{i}. Zmienną o rozkladzie normalnym (ozn. Z) o dowolnej średniej \mu _{norm} i wariancji \sigma _{norm} latwo uzyskamy poprzez liniową tranformację zmiennej ze standaryzowanego rozkladu normalnego (czyli Y):


Z=\sigma _{norm} Y+ \mu _{norm}
.

Drugi z przykladów generowania zmiennych z rozkladu normalnego polega na wylosowaniu dwóch liczb U_{1} i U_{2} z rozkladu jednorodnego [0,1]. Poniższa definicja opisuje dwie zmienne niezależne Y_{1} i Y_{2} mające standaryzowany rozklad normalny:


Y_{1}=\sin (2\pi U_{1})\sqrt{-2 \ln U_{2} }


Y_{2}=\cos (2\pi U_{1})\sqrt{-2 \ln U_{2} }

7. Metody wyznaczania energii swobodnej

Opis makroskopowy układu dynamicznego wymaga podania kilku z parametrów stanu układu, takich jak ciśnienie, objętość, temperatura, energia wewnętrzna, energia swobodna, przy czym istotna jest nie ich wartość bieżąca lecz raczej średnia w czasie. Wyznaczenie zaś ich wartości średniej ściśle wiąże się ze znajomością sumy statystycznej dla badanego układu. Niestety, analitycznie wyznaczyć ją można tylko dla bardzo prostych, wyidealizowanych systemów jak gaz doskonały lub układy o bardzo małej gęstości. Dla bardziej realnych układów, z dość złożonym opisem oddziaływań międzycząsteczkowych nie jest to możliwe. Należy więc sięgnąć do innych metod - doświadczalnych lub obliczeniowych. Wśród tych ostatnich wyróżnia się dwie podstawowe metody prowadzenia ewolucji układu, deterministyczną (Dynamikę Molekularną) i stochastyczną (metody Monte Carlo, dynamika brownowska).
W każdym z przypadków (deterministycznym i stochastycznym), zazwyczaj zyskuje się wiedzę o średniej wewnętrznej energii układu, jego własnościach strukturalnych (funkcje rozkładu) i innych własnościach z nich wynikających. Jednakże do opisu pewnych zjawisk, takich, jak mechanizm reakcji chemicznych, wpływ rozpuszczalnika na stabilność konformacji, procesy solwatacyjne potrzebna jest wiedza o konfiguracyjnej swobodnej energii zawartej w układzie.
Wiele parametrów termodynamicznych można wyznaczyć jako średnie po zespole statystycznym, a odwołując się do [hipotezy ergodycznej], jako średnią z symulacji prowadzonej w odpowiednim zespole. Postępując analogicznie z energią swobodną, można ją także zapisać w postaci pewnej średniej po zespole:

 F\ =\ -\beta ^{-1}\ln Z_{\mathit{NvT}}\ =\ -\beta^{-1}\ln (v^{-1}\int \mathit{dX}\operatorname{e}^{-\beta V(X)})\ =\ \beta ^{-1}\ln \frac{1}{v^{-1}\int \mathit{dX}\operatorname{e}^{-\beta V(X)}}\ =\ \beta^{-1}\ln
   \frac{v^{-1}\int \mathit{dX}\operatorname{e}^{-\beta V(X)}\operatorname{e}^{\beta V(X)}}{v^{-1}\int \mathit{dX}\operatorname{e}^{-\beta V(X)}}\ =\ \beta ^{-1}\ln <\operatorname{e}^{\beta V(X)}>_{\mathit{NvT}}

(do pierwszego wyrażenia w drugiej linii wstawiona została jedynka w postaci: v^{-1}\int \mathit{dX}\operatorname{e}^{-\beta
V}\operatorname{e}^{\beta V} , gdzie Z_{\mathit{NvT}} to kanoniczna konfiguracyjna suma statystyczna, v to objętość przestrzeni konfiguracyjnej N cząstek.

Zatem, jak widać, energię swobodną Helmholtza można wyznaczyć jako średnią z symylacji wartość \operatorname{e}^{\beta V(X)}. Ponieważ jednak exp będzie szybko rosnącą funkcją potencjału, najistotniejszy wkład do średniej będą miały konfiguracje z dużą dodatnią jego wartością. Tymczasem, metoda Metropolis Monte Carlo próbkuje najgęściej konfiguracje dające jak najniższą wartość potencjału, co wynika z faktu, że prawdopodobieństwo znalezienia danej konfiguracji jest proporcjonalne do \operatorname{e}^{-\beta V}. W efekcie, powyższa średnia po zespole będzie numerycznie niezbieżna.
W praktyce jednak, często zamiast absolutnej wartości energii swobodnej, której de facto w praktyce wyznaczyć się nie da, interesuje nas i jest możliwa do wyznaczenia jej wartość względna, czasem określana jako praca potrzebna na przeniesienie układu z jednego stanu do drugiego: \Delta F=F_{1}-F_{0}, czyli \Delta F=-\beta ^{-1}\ln Z_{1}/Z_{0}, gdzie Z_i, i=0,1 jest konfiguracyjną sumą statystyczną dla układu odpowiednio w stanie 0 lub 1. Także i w tym wypadku powstały metody pozwalające na ominięcie wyznaczania bezpośrednio dowolnej z Z_i.
Zazwyczaj, na początku parametryzuje się hamiltonian - H=H(\lambda ),\;\lambda \in [0,1] - tak, że H(\lambda_{0})=H_{0} i H(\lambda _{1})=H_{1}, czyli tak, by przy przejściu od \lambda _{0} do \lambda _{1} układ gładko przechodził ze stanu 0 do stanu 1. \lambda określa się często jako parametr sprzężenia (coupling parameter).

Następnie, wyznaczenie zmiany energii swobodnej odbywa się poprzez [23]:

  • całkowanie termodynamiczne,
  • metodę zaburzeń
  • umbrella sampling
    • potencjał średniej siły
    • energię swobodną Landau'a

7.1. Całkowanie termodynamiczne

Istotą tego podejścia jest zróżniczkowanie energii swobodnej po parametrze sprzężenia:


\frac{\mathit{dF}}{d\lambda }\ =\ -\frac{d}{d\lambda }\beta ^{-1}\ln
Z\ =\ -\frac{\beta ^{-1}}{Z}\frac{dZ}{d\lambda }\ =\ -\frac{\beta
^{-1}}{Z}\int \mathit{dX}\frac{d}{d\lambda }\ \operatorname{e}^{-\beta
H(\lambda ,X)}\ =\ \frac{1}{Z}\int
\mathit{dX}\frac{\mathit{dH}}{d\lambda }\ \operatorname{e}^{-\beta
H(\lambda ,X)}

czyli:


\frac{\mathit{dF}}{d\lambda }\ =\ \langle \frac{\mathit{dH}}{d\lambda}\rangle

Zatem zmiana energii swobodnej przy przejściu od stanu 0 do 1 daje się wyrazić:


\Delta F\ =\ F_{1}-F_{0}\ =\ \int _{\lambda _{0}}^{\lambda
_{1}}d\lambda \frac{\mathit{dF}}{d\lambda }\ =\ \int _{\lambda
_{0}}^{\lambda _{1}}d\lambda \langle \frac{\mathit{dH}}{d\lambda
}\rangle

Najprostsza parametryzacja hamiltonianu to H=(1-\lambda )H_{0}+\lambda H_{1} wówczas pochodna pod całką wynosi H_1 - H_0. Serię symulacji prowadzi się z kolejnymi wartościami parametru sprzężenia, zbiera z nich kolejne średnie wartości pochodnej spod całki i w ostatnim etapie sumuje uzyskane wyniki ze wzorem na \Delta F.
Parametryzacji można dokonać oczywiście na wiele sposobów, w zależności od procesu, który się opisuje.

7.2. Metoda zaburzeń

Tę metodę stosuje się zazwyczaj dla przypadków, w których stan początkowy i stan końcowy nie różnią się znacznie, czyli zmianę energii można traktować jako zaburzenie. Zmianę energii swobodnej przy przejściu ze stanu 0 do stanu 1 możemy zapisać:


\Delta F\ =\ F_{1}-F_{0}\ =\ -\beta ^{-1}\ln \frac{Z_{1}}{Z_{0}}\ =\ -\beta ^{-1}\ln \frac{\int \mathit{dX}\operatorname{e}^{-\beta H_{1}}}{\int \mathit{dX}\operatorname{e}^{-\beta H_{0}}}

Wstawiając 1 postaci \operatorname{e}^{-\beta H_{0}}\operatorname{e}^{\beta H_{0}} do wyrażenia podcałkowego w liczniku, możemy całość przepisać:


\Delta F\ =\ -\beta ^{-1}\ln \frac{\int\mathit{dX}\operatorname{e}^{-\beta (H_{1}-H_{0})}\operatorname{e}^{-\beta H_{0}}}{\int\mathit{dX}\operatorname{e}^{-\beta H_{0}}}\ =\ -\beta ^{-1}\ln \langle\operatorname{e}^{-\beta \Delta H}\rangle _{0}.

7.3. Umbrella sampling [24]

7.3.1. Potencjał średniej siły

Jest to energia swobodna wyrażona jako funkcja wybranych stopni swobody układu (np odległości między wybranymi dwoma cząstkami, kąta torsyjnego). Nazwiemy ten wybrany stopień swobody \lambda.
Prawdopodobieństwo, że \lambda przyjmie wartość z przedziału <\lambda,\lambda+\Delta\lambda> przy dowolnej konfiguracji przestrzennej układu wynosi:


\rho (\lambda )\ =\ Z^{-1}\int\mathit{dX}_{(N-2)}\operatorname{e}^{-\beta V(X_{(N-2)};\lambda)}\ =\ Z^{-1}\operatorname{e}^{-\beta w(\lambda )},

dzie całkujemy po stopniach swobody N-2 cząstek, zaś dwie ustawiamy w “odległości” \lambda od siebie, w(\lambda) jest szukanym potencjałem średniej siły, Z jest konfiguracyjną sumą statystyczną pełnego układu:


Z\ =\ \int \mathit{dX}_{(N-2)}\mathit{dX}_{(1)}\mathit{dX}_{(2)}\operatorname{e}^{-\beta V(X_{(N-2)},X_{(1)},X_{(2)})}\ =\ \int\mathit{dX}\operatorname{e}^{-\beta V(X)}.

Jeśli przyjmiemy, że \lambda jest odległością między dwiema wybranymi cząsteczkami, to mnożąc \rho(\lambda) przez czynnik \rho /(N(N-1)) (\rho jest gęstością makroskopową układu) otrzymujemy funkcję rozkładu radialnego g(\lambda). Rozwiązując pierwsze równanie ze względu na w, otrzymujemy:


w(\lambda )\ =\ -\beta ^{-1}\ln g(\lambda )+const.

Określenie w(\lambda) potencjałem średniej siły bierze się z faktu, że jeśli oddziaływania w układzie są dwuciałowe, addytywne, to gradient potencjału w(\lambda) będzie opisywał uśrednioną siłę z jaką układ działa na jedną z wybranych cząstek.


7.3.2. Energia swobodna [25, 26]

W poprzednim wypadku liczba stopni swobody uległa zmniejszeniu. Można sobie łatwo wyobrazić sytuację, w której chcemy zwiększyć liczbę stopni swobody poprzez wprowadzenie dodatkowego parametru porządku \lambda. Może nim być np któryś z parametrów potencjału oddziaływania międzycząsteczkowego, mający w standardowych symulacjach ustaloną wartość. Suma statystyczna dla rozszerzonego o parametr \lambda układu ma postać:


Z_{\lambda }\;=\;\int \mathit{dX}d\lambda \;\operatorname{e}^{-\beta
V(X,\lambda )}.

Rozkład prawdopodobieństwa zmiennej \lambda wyraża się następująco:


\rho (\lambda )\ =\ Z_{\lambda }^{-1}\int\mathit{dX}\operatorname{e}^{-\beta V(X;\lambda )}.

Całka w powyższym wzorze jest sumą stanów pełnego układu (N cząstek) z ustaloną wartością parametru \lambda:Z\ =\int\mathit{{dX}}\operatorname{e}^{-\beta\mathit{{\sigma}}(\mathit{{X}}{;}\lambda )}. Na podstawie powyższych wzorów energię swobodną Helmholtza można zapisać:


F\ =\ -\beta ^{-1}\ln Z\ =\ -\beta ^{-1}\ln \rho (\lambda )-\beta^{-1}\ln Z_{\lambda }\ =\ -\beta ^{-1}\ln \rho (\lambda)\;+\ const\ =\ F(\lambda ).

Chcąc wyznaczyć F(\lambda) napotykamy jednak pewien problem. Dla danego \lambda, rozszerzona przestrzeń konfiguracyjna jest najgęściej próbkowana w rejonach z najbardziej prawdopodobnym \lambda, przez co zakres \rho(\lambda) jest dość wąski.
“Kontroli” rejonu próbkowania dokonuje się poprzez dodanie funkcji ważącej zależnej tylko od \lambda do energii potencjalnej układu V(X;\lambda):


V(X,\lambda )\;=\;V(X,\lambda)\;+\;V_{\lambda _{0}}^{w}(\lambda ),

gdzie V_{\lambda _{0}}^{w}(\lambda ) jest arbitralnie wybraną funkcją \lambda, centrowaną w \lambda_0. Prowadząc teraz serię symulacji wokół kolejnych \lambda_0, możemy “przepróbkować” w miarę dokładnie dowolny żądany zakres \lambda. W kolejnym kroku łączymy unormowane histogramy. Możemy po prostu skorzystać z zależności \Delta F_{\mathit{ij}}=F(\lambda _{j})-F(\lambda _{i})=-\beta^{-1}\ln (\rho _{j}(\lambda )/\rho _{i}(\lambda )), przy czym \rho_j i \rho_i muszą się częsciowo pokrywać lub dostosować do naszych potrzeb najbardziej optymalną do tych celów metodę WHAM [27].

Należy jednak pamiętać, że w ten sposób otrzymujemy rozład w zespole ważonym, a nie rozpatrywanym. Przejście z zespołu ważonego do oryginalnego wykonuje się następująco:

\rho (\lambda )\ =\ \frac{\operatorname{e}^{-\beta V}}{\int \mathit{dX}d\lambda \operatorname{e}^{-\beta V}}\ =\ \frac{Z_{\lambda }\operatorname{e}^{-\beta V}\operatorname{e}^{\beta V_{w}}}{Z_{\lambda}\int \mathit{dX}d\lambda \operatorname{e}^{-\beta V}\operatorname{e}^{\beta V_{w}}}\ =\ \rho_{w}(\lambda )\frac{\operatorname{e}^{\beta V_{w}}}{\langle\operatorname{e}^{\beta V_{w}}\rangle _{w}}.

Wstawiając teraz powyższy wynik do równania na F(\lambda) otrzymujemy:

F(\lambda )\ =\ -\beta^{-1}\ln \rho (\lambda\ )+const\ =\ -\beta ^{-1}\ln \rho _{w}(\lambda)\;-\;V_{w}(\lambda )\;+\;\beta ^{-1}\ln \langle\operatorname{e}^{\beta V_{w}}\rangle _{w}\ +\ const.

Podobnie, wszystkie średnie zebrane w zespole ważonym powinny zostać przepisane na średnie liczone z oryginalnego rozkładu. Postępując analogicznie jak w ostatnim przypadku, średnią wartość pewnej obserwabli Y w rozkładzie boltzmannowskim wyrazić można wzorem:

\langle Y\rangle \ =\ \frac{\langle Y\operatorname{e}^{\beta V_{w}}\rangle _{w}}{\langle \operatorname{e}^{\beta V_{w}}\rangle_{w}}.


Powyższe trzy metody można zrealizować przy użyciu programu gibbs pakietu Amber oraz w Charmmie, zarówno akademickim jak i w wersji Accelrys'a, CHARMm dostępnej w ICM.

Literatura

[1] J. A. McCammon, S. C. Harvey, Dynamics of proteins and nucleic acids, Cambridge UP, ISBN 0521307503, (1987)

[2] M. Karplus, Molecular dynamics of biological macromolecules: A brief history and perspective, Biopolymers, 68:350-358, (2003)

[3] T. Schlick, Molecular Modeling and Simulation, Springer, ISBN 038795404X, (2002)

[4] G. R. Kneller, T. Mülders, Nosé-Andersen dynamics of partially rigid molecules: Coupling of all degrees of freedom to heat and pressure baths, Phys. Rev. E 54:6825-6837, (1996)

[5] S. Pantano, M. Tyagi, M. Giacca, P. Carloni, Molecular dynamics simulations on HIV-1 Tat, Eur. Biophys. J., 33:344-351, (2004)

[6] A. E. Garcia, Large-amplitude nonlinear motions in proteins, Phys. Rev. Lett., 68:2696-2699, (1992)

[7] A. Amadei, A. B. M. Linssen, H. J. C. Berendsen, Essential dynamics of proteins, Proteins: Struct. Funct. Genet., 17:412-425, (1993)

[8] D. L. Ermak, J. A. McCammon, Brownian dynamics with hydrodynamic interactions, J. Chem. Phys., 69:1352-1360, (1978)

[9] R. R. Gabdoulline, R. C. Wade, Brownian Dynamics Simulation of Protein-Protein Diffusional Encounter, METHODS: A Companion to Methods in Enzymology, 14:329-341, (1998)

[10] B. Honig, A. Nicholls, Classical Electrostatics in Biology and Chemistry, Science, 268:1144-1149, (1995)

[11] F. Fogolari, A. Brigo, H. Molinari, The Poisson-Boltzmann equation for biomolecular electrostatics: a tool for structural biology, 15:377-392, (2002)

[12] N. A. Baker, J. A. McCammon, Electrostatic Interactions, In Structural Bioinformatics, Weissig M., Bourne P. E. Eds, Wiley, New York, (2003)

[13] D. Bashford D. A. Case, Generalized born models of macromolecular solvation effects, Ann. Rev. Phys. Chem., 51: 129-152, (2000)

[14] V. Tozzini, Coarse grained models for proteins, Curr. Opin. Struct. Biol., 15:144-150, (2005)

[15] T. Simonson, Macromolecular electrostatics: continuum models and their growing pains, Curr. Opin. Struct. Biol., 11:243-252, (2001)

[16] N. A. Baker, D. Sept, S. Joseph, M. J. Holst, J. A. McCammon, Electrostatics of nanosystems: Application to microtubules and the ribosome, 98:10037-10041, (2001)

[17] J. Trylska, R. Konecny, F. Tama, C. L. Brooks III, J. A. McCammon, Ribosome motions modulate electrostatic properties, Biopolymers, 74:423-431, (2004)

[18] J. M. Antosiewicz, J. A. McCammon, M. K. Gilson, Prediction of pH-dependent properties of proteins. J Mol Biol 1994;238:415-436.

[19] J. M. Antosiewicz, J. M. Briggs, A. H. Elcock, M. K. Glison, J. A. McCammon, Computing Ionization States of Proteins with a Detailed Charge Model. J Comp Chem. 1996;17(14):1633-1644.

[20] J. Koronacki, J. Mielniczuk, Statystyka, Wydawnictwa Naukowo-Techniczne, Warszawa 2001

[21] S. Park, K. Miller, Comm. ACM, 31, 1988, 1191-1201

[22] Zaproponowanl D. H. Lehmer w 1948, opublikowane w Proc. 2nd Symposium on Large-Scale Digital Computing Machinery, Cambridge: Harvard University Press, 1951, 141-146

[23] D. L. Beveridge, F. M. DiCapua, Free energy via molecular simulation: A primer, Computer Simulation of Biomolecular Systems (Vol 1): Theoretical and Experimental Applications, (W.F. van Gunsteren, P.K. Weiner eds, ESCOM, Leiden, 1989) pp 1-26

[24] M. P. Allen, D. J. Tildesley, Computer Simulation of Liquids, (Clarendon Press, Oxford, 1987)

[25] L. D. Landau, E. M. Lifshitz, Statistical Physics

[26] R.M. Lynden-Bell, Landau free energy, Landau entropy, phase transitions and limits of metastability in an analytical model with a variable number of degrees of freedom, Mol. Phys., 86, 1353-1374, (1995)

[27] S. Kumar, D. Bouzida, R. H. Swendsen, P. A. Kollman, J. M. Rosenberg, "The Weighted Histogram Analysis Method for Free Energy Calculations on Biomolecules. I. The Method, J.Comp.Chem., 13, 1011-1021, 1992