Model HP

Model HP

Model HP (hydrophobic-polar protein folding model) to model polimeru wykorzystywany w badaniach nad ogólnymi zasadami rządzącymi procesem zwijania białek. Badania tego typu w przypadku modeli pełnoatomowych wiążą się ze znacznymi kosztami obliczeniowymi, podczas gdy w modelu HP, ze względu na uproszczoną charakterystykę układu, możliwe jest przeprowadzenie krótkiej symulacji (trwającej od kilku minut do kilku godzin), w trakcie której układ jest w stanie osiągnąć wszystkie możliwe mikrostany .

Wstęp

Plik:Hp2d2 1.png

(1) Dwuwymiarowy model HP o sekwencji: HPHPHHHHHPHP (wizualizacja w PyMOLu). Na czerwono zaznaczono aminokwasy hydrofobowe (H), na czarno aminokwasy polarne (P). Ponieważ w powyższym mikrostanie nie występują kontakty H-H, energia wynosi 0.

Plik:Hp2d2 inter.png

(2) Obrót wokół szóstego aminokwasu (zaznaczono na zielono) skutkuje utworzeniem kontaktu H-H między ósmym i piątym aminokwasem.

Plik:Hp2d2 2.png

(3) Transformacja została zaakceptowana, liczba kontaktów H-H wynosi 1, zatem nowa energia układu wynosi -ɛ.

Idea modelu HP opiera się na obserwacji, iż kluczową rolę w procesie zwijania białek pełni efekt hydrofobowy (w tym kontekście spotkać się można z terminem: "oddziaływania hydrofobowe"). W podstawowym modelu HP polimer zbudowny jest z monomerów H (hydrofobowych) oraz P (polarnych), przy czym wkład do energii pochodzi jedynie od H. Można więc myśleć o modelu HP jak o modelu białka, w którym alfabet aminokwasów ograniczony został do zbioru {H,P}. Aminokwasy znajdują się w węzłach sieci kwadratowej (square lattice) w przypadku modelu dwuwymiarowego (2D), bądź w węzłach sieci sześciennej (cubic lattice) w przypadku modelu trójwymiarowego (3D). Dwa aminokwasy nie mogą znajdować się w tym samym węźle. Natomiast jeśli dwa aminokwasy połączone są wiązaniem (przez analogię do wiązania peptydowego między aminokwasami w białkach), to muszą się one znajdować w sąsiednich węzłach. Mikrostan układu można określić poprzez: sekwencję peptydu oraz współrzędne poszczególnych aminokwasów.

Ewolucję układu w modelu HP zadaje zestaw dozwolonych transformacji struktury oraz rozkład prawdopodobieństwa przejść między mikrostanami. Przykładem dozwolonej transformacji może być obrót części białka o pewien kąt wokół wybranego aminokwasu (przykład przedstawiono po prawej). W przypadku modelu 2D istnieją trzy możliwe nietrywialne obroty. Jeżeli po dokonaniu obrotu żadne dwa aminokwasy nie zajmują tego samego punktu w przestrzeni, obrót uznajemy za dozwolony.

Po dokonaniu dozwolonej transformacji prawdopodobieństwo akceptacji nowego mikrostanu zależy jedynie od zmiany wartości energii. Innymi słowy: to, czy zaakceptujemy mikrostan uzyskany w wyniku transformacji zależy jedynie od mikrostanu przed transformacją; wcześniejsza historia układu nie ma tu znaczenia. Zatem ewolucja peptydu (ciąg mikrostanów wygenerowany w toku symulacji) jest realizacją procesu stochastycznego, w którym prawdopodobieństwo zdarzenia (akceptacja nowego mikrostanu) zależy jedynie od wyniku poprzedniego. Proces stochastyczny tego typu w przypadku dyskretnej przestrzeni stanów nazywany jest łańcuchem Markowa.

Sieć

Niech:

$\mathbf {e} _{x}=(1,0),\;\mathbf {e} _{y}=(0,1)$

będą wektorami bazowymi w przypadku dwuwymiarowym, zaś:

$\mathbf {e} _{x}=(1,0,0),\;\mathbf {e} _{y}=(0,1,0),\;\mathbf {e} _{z}=(0,0,1)$

wektorami bazowymi w przypadku trójwymiarowym. Siecią kwadratową nazywać będziemy zbiór:

$LATTICE_{2D}=\{x\mathbf {e} _{x}+y\mathbf {e} _{y}\mid x,y\in \mathbb {Z} \}$

zaś zbiór:

$LATTICE_{3D}=\{x\mathbf {e} _{x}+y\mathbf {e} _{y}+z\mathbf {e} _{z}\mid x,y,z\in \mathbb {Z} \}$

nazwiemy siecią sześcienną. Element sieci (węzeł) opisujemy przez podanie dwóch, bądź trzech liczb całkowitych (współrzędnych węzła), przykładowo dla sieci sześciennej: (0,1,-10).

Powiemy, że węzły a i b sąsiadują ze sobą na siatce (ozn. a ~ b), jeżeli istnieje wektor bazowy e taki, że:

$\mathbf {a} =\mathbf {b} +\mathbf {e} \quad \lor \quad \mathbf {b} =\mathbf {a} +\mathbf {e}$.

Niech $CHAIN_{n}=\{1,...,n\}$ będzie zbiorem aminokwasów tworzących peptyd, gdzie $n$ - długość peptydu. Wówczas strukturę przestrzenną wyrażać będziemy przez funkcję:

$\mathbf {s} \colon CHAIN_{n}\to LATTICE$

spełniającą warunki:

$\mathbf {s} (1)=(0,0,0)$

$\forall _{i<n}\mathbf {s} (i+1)\sim \mathbf {s} (i)$

$\forall _{i\not =j}\mathbf {s} (i)\not =\mathbf {s} (j)$

Podanie struktury (w postaci funkcji s) nie wystarcza do określenia miktrostanu układu, potrzebna jest jeszcze sekwencja.

Sekwencja

Sekwencja łańcucha określona jest przez wzorzec hydrofobowy $Pat\colon CHAIN_{n}\to \{H,P\}$. Rozważany model dzieli aminokwasy ze względu na właściwości oddziaływań dalekozasięgowych na dwie kategorie: hydrofobowe (H) oraz polarne (P). Dalekozasięgowość oddziaływań odnosi się do wzajemnych położeń aminokwasów w sekwencji, a nie w przestrzeni. Przykładowo: o obecności oddziaływań dalekozasięgowych możemy mówić w przypadku pary aminokwasów o numerach 1 i 4, bądź: 2 i 9, ale nie w przypadku par: 1 i 3, czy też 4 i 5. Szczegóły w poniższej sekcji Oddziaływania.

Oddziaływania

W najprostszym modelu HP rozważa się jedynie oddziaływania dalekozasięgowe pomiędzy aminokwasami hydrofobowymi. Energia danego mikrostanu zależy od liczby kontaktów występujących między aminokwasami H, niesąsiadującymi w peptydzie.

Niech

$K_{HH}(\mathbf {s} )=\#\{\{i,j\}\colon \mid i-j\mid >1,\quad \mathbf {s} (i)\sim \mathbf {s} (j),\quad Pat(i)=Pat(j)=H\}$

będzie liczbą kontaktów H-H w peptydzie o strukturze s. Energia układu wyraża się przez:

$E(\mathbf {s} )=-\varepsilon \cdot K_{HH}(\mathbf {s} )$

gdzie ɛ >0.

Interakcje pomiędzy aminokwasami hydrofobowymi odzwierciedlają ich tendencję do kierowania się do wewnątrz białka i tym samym unikania kontaktu z wodą. Należy podkreślić, że model HP uwydatnia jeden aspekt procesu zwijania białek (efekt hydrofobowy), ignoruje natomiast oddziaływania lokalne występujące w rzeczywistym białku - "sztywność" łańcucha (objawiająca się niedozwolonymi wartościami kątów φ-ψ na wykresie Ramachandrana) oraz wiązania wodorowe (istotne w α-helisach i β-kartkach). Proste modele, jak model HP, skłaniają do zadawania pytań: Które z własności białek udaje się odtworzyć pomimo poczynionych przybliżeń?

Średnia po zespole

Niech A będzie pewną własnością fizyczną badanego układu. Mikrostan układu oznaczymy przez $\mathbf {x} =(x_{1},\ldots ,x_{n})$, gdzie n jest liczbą stopni swobody. Przyjmujemy, że własność A objawia się jako średnia po próbce pewnej przestrzeni mikrostanów, tzn.:

$\langle A\rangle =Z^{-1}\int _{\Omega }{A(\mathbf {x} )f({\mathcal {H}}(\mathbf {x} ))d\mathbf {x} }$

gdzie f jest funkcją rozkładu gęstości prawdopodobieństwa, $\Omega$ jest przestrzenią dostępnych stanów układu (nazywana również w szerszym kontekście: przestrzenią fazową), zaś:

$Z=\int _{\Omega }f({\mathcal {H}}(\mathbf {x} ))d\mathbf {x}$

to sumą statystyczna nazywana również funkcją podziału. Rozkład f określa odpowiedni zespół statytyczny (mikrokanoniczny, kanoniczny,...).

W przypadku modelu HP liczba mikrostanów układu jest skończona (ozn. N), zaś średnią wartość A wyraża się w postaci sumy po dostępnych mikrostanach układu:

$\langle A\rangle =\sum _{i=1}^{N}A_{i}\cdot p_{i}$

gdzie $p_{i}$ jest prawdopodobieństwem uzyskania przez układ i-tego mikrostanu. Prawdopodobieństwo, że układ o określonej, stałej temperaturze T (używa się też określenia: w kontakcie z termostatem o temperaturze T) osiągnie i-ty mikrostan o energii $E_{i}$, dane jest rozkładem Boltzmanna:

$p_{i}={\frac {e^{-E_{i}/kT}}{\sum _{j=1}^{N}e^{-E_{j}/kT}}}$

gdzie k - stała Boltzmanna. Suma w mianowniku zapewnia normalizację rozkładu $p_{i}$:

$\sum _{j=1}^{N}p_{j}=1$

Metoda Monte Carlo - algorytm Metropolisa

W celu wyznaczenia $\langle A\rangle$ dla układu o temperaturze T wystarczy dysponować metodą do generowania mikrostanów zgodnie z rozkładem Boltzmanna. Metodą tego typu jest algorytm Metropolisa.

Wprowadźmy następujące oznaczenie:

$\pi (a)=e^{-E_{a}/kT}$

gdzie a jest mikrostanem o energii $E_{a}$. Istotą algorytmu Metropolisa jest stworzenie ciągu mikrostanów, będący realizacją łańcucha Markowa z prawdopodobieństwem przejść, zależącym od różnicy energii kolejnych mikrostanów. W przypadku modelu HP algorytm Metropolisa przebiega następująco:
1. Zainicjuj ciąg mikrostanów, tworząc pierwszy, dowolny mikrostan X.
2. Oblicz energię $E_{X}$.
3. Dokonaj dozwolonej transformacji peptydu (transformacje opisano dalej, dla ustalenia uwagi - dokonujemy obrotu części peptydu wokół losowo wybranego aminokwasu).
4. Wyznacz energię $E_{Y}$ uzyskanego w wyniku transformacji mikrostanu Y.
5. Zaakceptuj nowy mikrostan (X:=Y) z prawdopodobieństwem $p(X,Y)=min\left\{1,{\frac {\pi (X)}{\pi (Y)}}\right\}$ i wróć do 3. albo zakończ, jeśli wygenerowano ciąg o długości M.

Dysponując ciągiem $(\mathbf {x} _{n})$ mikrostanów uzyskanych w algorytmie Metropolisa, możemy wyznaczyć średnią wartość A:

$\langle A\rangle \approx {\frac {1}{M}}\sum _{i=1}^{M}A(\mathbf {x} _{i})$

Transformacje

...

Symulowane wyżarzanie (Simulated Annealing)

Najprostszym sposobem znajdowania konformacji o minimalnej energii jest systematyczne obniżanie temperatury podczas symulacji. Wadą tego rozwiązania jest to, że układ może łatwo zatrzymać się w lokalnym minimum energii, z którego wyjście przy obniżonej temperaturze okaże się niemożliwe (precyzyjniej: niezwykle mało prawdopodobne). Ponadto, zbieżność algorytmu przy niskich temperaturach jest dosyć wolna. Układ może stracić dużo czasu (kroków symulacji) w niecce reprezentującej lokalne minimum, bądź oscylując między stanami o tej samej energii.

Zamiana replik (Replica Exchange Monte Carlo)

W tym podejściu równolegle symuluje się wiele kopii układu, każdy w innej, stałej temperaturze. Załóżmy, że w pewnym momencie symulacji algorytmu Metropolisa i-ta replika o temperaturze $T_{i}$ jest w mikrostanie $\mathbf {x} _{i}$ o energii $E(\mathbf {x} _{i})$, zaś j-ta replika w odpowiednio: temperaturze $T_{j}$, mikrostanie $\mathbf {x} _{j}$ i energii $E(\mathbf {x} _{j})$. Z rozkładu jednostajnego losujemy parę kolejnych replik (i,j) , które z prawdopodobieństwem $p_{s}$ zostaną zamienione miejscami:

$p_{s}=min\{1,e^{-\Delta }\},$

gdzie

$\Delta =\left({\frac {1}{kT_{j}}}-{\frac {1}{kT_{i}}}\right)(E(\mathbf {x} _{i})-E(\mathbf {x} _{j}))$

Po zamianie i-ta replika symulowana jest w temperaturze $T_{j}$ , a j-ta w temperaturze $T_{i}$. Ponieważ prawdopodobieństwo zamiany maleje wykładniczo wraz ze wzrostem różnicy temperatur, rozważamy wyłącznie repliki sąsiednie.

Zamiana temperatur zmienia krajobraz energetyczny. W bardzo wysokich temperaturach bariery energetyczne znikają i można domniemywać, że prawdopodobieństwo odwiedzenia mikrostanu jest zadane rozkładem jednostajnym. Repliki, które utknęły w lokalnych minimach mogą zostać z nich wyzwolone przez przeniesienie do wyższej temperatury.

Wymiany nie powinny być zbyt częste. Po zmianie temperatury układ przez pewien czas się stabilizuje i przemieszcza w najbardziej prawdopodobny region krajobrazu energetycznego.

Po zakończeniu symulacji średnia A w temperaturze $T_{i}$ może zostać oszacowana wzorem:

$\langle A\rangle \approx {\frac {1}{M}}\sum _{i=1}^{M}A(\mathbf {x} _{i})$

Symulacje planowane na ćwiczeniach

Celem ćwiczeń jest zaimplementowanie modelu HP (domyślnie w języku Java) i sprawdzenie wyników przedstawionych w publikacji K. A. Dilla z 1995 roku. W pierwszym etapie przeprowadzimy symulowane wyżarzania modeli dwuwymiarowych trzech polimerów o sekwencjach:

• PHPPHPPHHPPHHPPHPPHP
• HPPHPPHPHPPHPHPHHH
• HPPPHHPPHPHHPHHH

(na ostatnich ćwiczeniach zawęziłem symulacje obowiązujące na 25. marca do jednej sekwencji, pierwszej na powyższej liście). W przypadku każdego peptydu temperaturą początkową będzie $T_{max}=1$, a temperaturą końcową $T_{min}=0.1$. Temperatura będzie w trakcie symulacji maleć o czynnik $\delta =0.05$. Po przeprowadzeniu $x=10000$ transformacji (dokonując akceptacji/odrzuceń wygenerowanych konformacji) w danej temperaturze $T$, przeprowadzamy kolejnych $x$ kroków w temperaturze $T-\delta$ i tak dalej, aż do osiągnięcia $T_{min}$.

W każdej temperaturze wygenerujemy $m(T)<x$ konformacji, które akceptować będziemy na drodze algorytmu Metropolisa i dla każdej takiej próby możemy wyznaczyć ciepło właściwe układu $C_{v}$, średni moment bezwładności $I$ oraz histogram wystąpień stanów w zależności od liczby kontaktów.

(z ostatnich ćwiczeń: na 25. marca proszę wykonać wykres Cv(T) oraz histogramy dla poszczególnych temperatur. Załączone rysunki przedstawiają czego spodziewać się można po wynikach).

Plik:Cv.png

Ciepło właściwe w funkcji temperatury, uzyskane w symulowanym wyżarzaniu.

Plik:HistogramT0 3.png

Histogram wystąpień stanów w zależności od liczby kontaktów.

Plik:MeanI.jpg

Średni moment bezwładności w funkcji temperatury.

Ciepło właściwe wyraża się wzorem:

$C_{V}(T)={\frac {\langle E(T)^{2}\rangle -\langle E(T)\rangle ^{2}}{kT^{2}}}$

Warto zauważyć, że w $C_{V}$ jest proporcjonalne do wariancji energii i można ją interpretować jako "miarę rozrzutu energii" stanów, jakie osiąga układ w danej temperaturze.

Moment bezwładności definiujemy następująco:

$I=\sum _{i=1}^{n}(\mathbf {s} _{i}-\mathbf {s} _{0})^{2}$

gdzie $s_{0}$ jest środkiem ciężkości, a n liczbą aminokwasów w strukturze. Natywna struktura większości białek jest globularna. Można więc przyjąć, że moment bezwładności jest dobrym przybliżeniem "stopnia zwinięcia białka".

Model HP "na piątkę"

Na ocenę bardzo dobrą z tej części ćwiczeń należy przeprowadzić dodatkowo (poza symulacją na zaliczenie, na 25. marca) dwie symulacje, dla białek o sekwencjach:

• HPPPHHPPHPHHHHHH
• HPPPHHPPHPHHPHHH

które różnią się sekwencyjnie tylko na jednej pozycji (praca Dilla, str. 20). Specyfikacje symulacji są te same, co poprzednio. Należy wykonać analogiczne histogramy wystąpień kontaktów, wykresy $Cv(T)$ oraz - analogiczny do wykresu $Cv$ - wykres $\langle I(T)\rangle$.

Czy różnica sekwencyjna na jednej pozycji powoduje zmianę minimalnej wartości energii dla tych dwóch białek? Czy obydwa białka "zwijają się" do minimum energii? A jak nie - to do minimum jakiego potencjału?

Linki zewnętrzne

• Model HP w Wikipedii
• Mikrostan w Wikipedii
• Proces stochastyczny w Wikipedii
• Łańcuch Markowa w Wikipedii
• Suma statystyczna (funkcja podziału) w Wikipedii
• Algorytm Metropolisa
• Praca K. A. Dilla z 1995 roku
• Strona WWW projektu PyMOL

Szablon:Cite journal