Tekst źródłowy strony 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 <ref name="dill1995">{{cite journal |author=Dill K.A. |title=Principles of protein folding - A perspective from simple exact models |journal=Protein science |volume=4 |issue=4 |year=1995 |id={{Entrez Pubmed|7613459}} |pages=561–602 |pmid=7613459}}</ref>. 
__TOC__


== Wstęp ==
[[Image:Hp2d2 1.png|thumb|upright|300px|(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.]]
[[Image:Hp2d2 inter.png|thumb|upright|300px|(2) Obrót wokół szóstego aminokwasu (zaznaczono na zielono) skutkuje utworzeniem kontaktu H-H między ósmym i piątym aminokwasem.]]
[[Image:Hp2d2 2.png|thumb|upright|300px|(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:
: <math>\mathbf e_x = (1, 0),\; \mathbf e_y = (0, 1) </math>
będą wektorami bazowymi w przypadku dwuwymiarowym, zaś:
: <math>\mathbf e_x = (1, 0, 0),\; \mathbf e_y = (0, 1, 0),\; \mathbf e_z = (0, 0, 1)</math>
wektorami bazowymi w przypadku trójwymiarowym. Siecią kwadratową nazywać będziemy zbiór:
: <math> LATTICE_{2D} = \{ x \mathbf e_x + y \mathbf e_y \mid x,y\in \mathbb Z  \} </math>
zaś zbiór:
: <math> LATTICE_{3D} = \{ x \mathbf e_x + y \mathbf e_y + z \mathbf e_z \mid x,y,z\in \mathbb Z  \} </math>
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:
: <math>\mathbf a = \mathbf b + \mathbf e \quad \or \quad \mathbf b = \mathbf a + \mathbf e </math>.

Niech <math>CHAIN_n=\{1,...,n\}</math> będzie zbiorem aminokwasów tworzących peptyd, gdzie <math> n </math> - długość peptydu. Wówczas strukturę przestrzenną wyrażać będziemy przez funkcję:
: <math>\mathbf s \colon CHAIN_n \to LATTICE </math>
spełniającą warunki:
: <math>\mathbf s ( 1 ) = ( 0,0,0 ) </math>
: <math>\forall_{i<n} \mathbf s (i+1) \sim \mathbf s(i) </math>
: <math>\forall_{i \not= j}\mathbf s(i) \not= \mathbf s(j) </math>
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 <math> Pat \colon CHAIN_n \to \{H,P\} </math>. 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 
:<math> 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 \} </math>
będzie liczbą kontaktów H-H w peptydzie o strukturze '''s'''. Energia układu wyraża się przez:
: <math> E ( \mathbf s )= -\varepsilon \cdot K_{HH}( \mathbf s )    </math>
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 <math> \mathbf x = (x_1, \ldots , x_n) </math>, gdzie ''n'' jest liczbą stopni swobody. Przyjmujemy, że własność ''A'' objawia się jako średnia po próbce pewnej przestrzeni mikrostanów, tzn.:
: <math> \langle A \rangle = Z^{-1} \int_{\Omega}{A( \mathbf x ) f( \mathcal{H}( \mathbf x ) ) d \mathbf x}  </math>
gdzie ''f'' jest funkcją rozkładu gęstości prawdopodobieństwa, <math>\Omega </math> jest przestrzenią dostępnych stanów układu (nazywana również w szerszym kontekście: przestrzenią fazową), zaś:
: <math> Z = \int_{\Omega} f( \mathcal{H} ( \mathbf x ) ) d \mathbf x </math>
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:
: <math> \langle A \rangle = \sum_{i=1}^{N} A_i \cdot p_i    </math>
gdzie <math> p_i </math> 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 <math> E_i </math>, dane jest rozkładem Boltzmanna:
: <math> p_i = \frac{e^{-E_i/kT}}{\sum_{j=1}^{N}e^{-E_j/kT}} </math>
gdzie ''k'' - stała Boltzmanna. Suma w mianowniku zapewnia normalizację rozkładu <math> p_i </math>:
: <math> \sum_{j=1}^{N} p_j = 1 </math>

== Metoda Monte Carlo - algorytm Metropolisa ==
W celu wyznaczenia <math> \langle A \rangle </math> 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: 
:<math> \pi (a) = e^{-E_a/kT}  </math> 
gdzie ''a'' jest mikrostanem o energii <math> E_a  </math>. 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: 
<br>1. Zainicjuj ciąg mikrostanów, tworząc pierwszy, dowolny mikrostan ''X''. 
<br>2. Oblicz energię <math> E_X </math>. 
<br>3. Dokonaj dozwolonej transformacji peptydu (transformacje opisano dalej, dla ustalenia uwagi - dokonujemy obrotu części peptydu wokół losowo wybranego aminokwasu). 
<br>4. Wyznacz energię <math> E_Y </math> uzyskanego w wyniku transformacji mikrostanu ''Y''. 
<br>5. Zaakceptuj nowy mikrostan (X:=Y) z prawdopodobieństwem <math>p(X,Y)=min \left\{ 1, \frac{\pi (X) }{\pi (Y)} \right\} </math> i wróć do 3. albo zakończ, jeśli wygenerowano ciąg o długości ''M''.

Dysponując ciągiem <math> (\mathbf x_n ) </math> mikrostanów uzyskanych w algorytmie Metropolisa, możemy wyznaczyć średnią wartość ''A'':
: <math> \langle A \rangle \approx \frac{1}{M} \sum_{i=1}^{M} A( \mathbf x_i )   </math>

== 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 <math> T_i </math> jest w mikrostanie <math>\mathbf x_i </math> o energii <math> E( \mathbf x_i ) </math>, zaś ''j''-ta replika w odpowiednio: temperaturze <math> T_j </math>, mikrostanie <math> \mathbf x_j </math> i energii <math> E( \mathbf x_j ) </math>. Z rozkładu jednostajnego losujemy parę kolejnych replik (i,j) , które z prawdopodobieństwem <math> p_s </math> zostaną zamienione miejscami:
:<math> p_s = min \{ 1, e^{-\Delta} \}, </math>
gdzie
:<math> \Delta = \left( \frac{1}{kT_j} - \frac{1}{kT_i} \right) ( E( \mathbf x_i ) - E( \mathbf x_j )  )   </math>
Po zamianie ''i''-ta replika symulowana jest w temperaturze <math>T_j</math> , a ''j''-ta w temperaturze <math>T_i</math>. 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 <math> T_i </math> może zostać oszacowana wzorem:
: <math> \langle A \rangle \approx \frac{1}{M} \sum_{i=1}^{M} A( \mathbf x_i )   </math>

== 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 <math> T_{max} =1</math>, a temperaturą końcową <math>T_{min} =0.1 </math>. Temperatura będzie w trakcie symulacji maleć o czynnik <math> \delta=0.05 </math>. Po przeprowadzeniu <math> x=10000 </math> transformacji (dokonując akceptacji/odrzuceń wygenerowanych konformacji) w danej temperaturze <math> T </math>, przeprowadzamy kolejnych <math> x </math> kroków w temperaturze <math> T - \delta </math> i tak dalej, aż do osiągnięcia <math> T_{min} </math>. 

W każdej temperaturze wygenerujemy <math> m(T) < x </math> 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 <math> C_v </math>, średni moment bezwładności <math> I </math> 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).'''
[[Image:Cv.png|thumb|upright|300px|Ciepło właściwe w funkcji temperatury, uzyskane w symulowanym wyżarzaniu.]]
[[Image:HistogramT0_3.png|thumb|upleft|300px|Histogram wystąpień stanów w zależności od liczby kontaktów.]]
[[Image:MeanI.jpg|thumb|upleft|300px|Średni moment bezwładności w funkcji temperatury.]]

Ciepło właściwe wyraża się wzorem:
: <math> C_V (T) = \frac{\langle E(T)^2 \rangle - \langle E(T) \rangle ^2 }{k T^2}   </math> 
Warto zauważyć, że w <math> C_V </math> 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:
: <math> I = \sum_{i=1}^{n} (\mathbf s_i - \mathbf s_0)^2    </math>
gdzie <math> s_0 </math> 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 <math> Cv(T) </math> oraz - analogiczny do wykresu <math> Cv </math> - wykres <math> \langle I(T) \rangle </math>.

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 ==
* [http://en.wikipedia.org/wiki/Hydrophobic-polar_protein_folding_model Model HP w Wikipedii]
* [http://en.wikipedia.org/wiki/Microstate_%28statistical_mechanics%29 Mikrostan w Wikipedii]
* [http://en.wikipedia.org/wiki/Stochastic_process Proces stochastyczny w Wikipedii]
* [http://en.wikipedia.org/wiki/Markov_Chain Łańcuch Markowa w Wikipedii]
* [http://en.wikipedia.org/wiki/Partition_function_%28statistical_mechanics%29 Suma statystyczna (funkcja podziału) w Wikipedii]
* [http://en.wikipedia.org/wiki/Metropolis_algorithm Algorytm Metropolisa]
* [http://www.ncbi.nlm.nih.gov/pmc/articles/PMC2143098/ Praca K. A. Dilla z 1995 roku]
* [http://www.pymol.org/ Strona WWW projektu PyMOL]