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

(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.

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

(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:

${𝐞}_x=(1,0),{𝐞}_y=(0,1)$

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

${𝐞}_x=(1,0,0),{𝐞}_y=(0,1,0),{𝐞}_z=(0,0,1)$

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

$LATTIC{E}_{2D}=\{x{𝐞}_x+y{𝐞}_y\mid x,y\in \mathbb{Z}\}$

zaś zbiór:

$LATTIC{E}_{3D}=\{x{𝐞}_x+y{𝐞}_y+z{𝐞}_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:

$𝐚=𝐛+𝐞\vee 𝐛=𝐚+𝐞$.

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

$𝐬:CHAI{N}_n\to LATTICE$

spełniającą warunki:

$𝐬(1)=(0,0,0)$

${\mathrm{\forall }}_{i<n}𝐬(i+1)\sim 𝐬(i)$

${\mathrm{\forall }}_{i≠j}𝐬(i)≠𝐬(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:CHAI{N}_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}(𝐬)=\mathrm{\#}\{\{i,j\}:\mid i-j\mid >1,𝐬(i)\sim 𝐬(j),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(𝐬)=-\epsilon \cdot K_{HH}(𝐬)$

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 $𝐱=(x_1,\dots ,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.:

$⟨A⟩=Z^{-1}{\int }_{\Omega }A(𝐱)f(\mathscr{H}(𝐱))d𝐱$

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(\mathscr{H}(𝐱))d𝐱$

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:

$⟨A⟩=\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 $⟨A⟩$ 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\{1,\frac{\pi (X)}{\pi (Y)}\}$ i wróć do 3. albo zakończ, jeśli wygenerowano ciąg o długości M.

Dysponując ciągiem $({𝐱}_n)$ mikrostanów uzyskanych w algorytmie Metropolisa, możemy wyznaczyć średnią wartość A:

$⟨A⟩\approx \frac{1}{M}\sum _{i=1}^{M}A({𝐱}_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 ${𝐱}_i$ o energii $E({𝐱}_i)$, zaś j-ta replika w odpowiednio: temperaturze $T_j$, mikrostanie ${𝐱}_j$ i energii $E({𝐱}_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 =(\frac{1}{k{T}_j}-\frac{1}{k{T}_i})(E({𝐱}_i)-E({𝐱}_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:

$⟨A⟩\approx \frac{1}{M}\sum _{i=1}^{M}A({𝐱}_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).

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

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

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

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

$C_V(T)=\frac{⟨E(T{)}^2⟩-⟨E(T){⟩}^2}{k{T}^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({𝐬}_i-{𝐬}_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 $⟨I(T)⟩$.

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