Błony biologiczne są złożonymi układami wieloskładnikowymi, dlatego też w celu pozna-nia ich własności wykorzystywane są błony modelowe. Dzięki badaniom przeprowadzonym na uproszczonych modelach możliwe jest poznanie biofizycznych własności błon biologicznych oraz wpływu wybranych składników na ich strukturę i funkcję. Ze względu na ograniczenia metod eks-perymentalnych, badania domen cholesterolowych za pomocą metod obliczeniowych są bardzo ważną alternatywą.
Przedstawiana praca doktorska miała na celu zbadanie właściwości dwuwarstw o wysokim stężeniu cholesterolu, jak to ma miejsce w komórkach włókienkowych soczewki oka. Główne cele pracy:
1. Określenie wpływu wysokiego stężenia cholesterolu na dwuwarstwę fosfolipidową nasyconą cholesterolem, na przykładzie fosfolipidu POPC.
2. Uzyskanie charakterystyki domeny czysto cholesterolowej CBD na poziomie molekular-nym z rozdzielczością atomową i porównanie jej z domeną PCD w fazie ciekłokrystalicznej uporządkowanej lo.
3. Zbadanie procesu dyfuzji tlenu cząsteczkowego w dwuwarstwach lipidowych różniących się zawartością POPC i cholesterolu.
4. Poznanie organizacji lipidów w układach różniących się zawartością POPC i cholesterolu.
2
Materiały i metody
2.1 Symulacje dynamiki molekularnej
Symulacja dynamiki molekularnej (ang. molecular dynamics simulation) jest numeryczną metodą badania zachowania dynamicznego modeli komputerowych. W trakcie symulacji dy-namiki molekularnej na podstawie znanego położenia i sił działających na pojedyncze atomy generowana jest trajektoria układu symulacyjnego.
Ze względu na napotykane ograniczenia mocy obliczeniowych i możliwości zapisu trajektorii symulacyjnych, na początku badań nad dwuwarstwami lipidowymi, w latach 80-tych, koniecz-nością było stosowanie uproszczeń, takich jak modele uogólnionych atomów (ang. united atom model), gdzie grupa atomów zastępowana jest pseudoatomami, lub zastąpienie cząsteczkowej wo-dy ciągłym ośrodkiem modelującym jego polarne właściwości. Również całkowity maksymalny czas przeprowadzanych symulacji był ograniczony do około kilku nanosekund.
Wraz z potężnym rozwojem technicznym, wzrosły możliwości zastosowania symulacji dyna-miki molekularnej do badania układów większych i przez znacząco dłuższy czas. Dodatkowo, modele stały się bardziej realistyczne, dzięki zastosowaniu atomowych modeli cząsteczek i roz-puszczalnika. Obecnie przeprowadzane symulacje dynamiki molekularnej modeli błon trwają po kilkaset nanosekund, a błony liczą kilkaset cząsteczek lipidów.
2.1.1 Postać funkcji potencjału
W symulacjach dynamiki molekularnej, wyznaczenie położenia atomów w kolejnych krokach czasowych odbywa się poprzez przybliżenie rozwiązania równania ruchu Newtona dla każdego atomu z osobna. Równanie Newtona dla atomu i ma formułę:
Fi= d2ri
dt2 mi, (2.1)
35
36 2. Materiały i metody gdzie Fi jest siłą działającą na i-ty atom, mi jego masą, a ri jest wektorem położenia. Rów-nanie ruchu Newtona nie ma rozwiązania analitycznego dla układów wieloatomowych, dlatego też stosuje się metody numeryczne. Symulacje dynamiki molekularnej przedstawione w poniższej pracy przeprowadzono w oparciu o algorytm o nazwie leapfrog integration (Van Gunsteren i Be-rendsen, 1988), który jest modyfikacją algorytmu Verlet’a (Verlet, 1967). Algorytm ten położenie i prędkość po czasie δt wyznacza w oparciu o równania:
r(t + δt) = r(t) + δtv(t +1
2δt) (2.2)
v(t + δt) = v(t − 1
2δt) + δtFi
mi (2.3)
Warto zauważyć, że zarówno położenie jak i prędkość wyznaczane są z tym samym krokiem czasowym δt, jednakże wyznaczanie prędkości odbywa się o pół kroku czasowego wcześniej. Ten przeskok momentu wyznaczania położeń i prędkości jest wytłumaczeniem nazwy tego algorytmu.
Zaletami algorytmu Verlet’a i jego modyfikacji jest jednokrotne obliczanie sił w danym kroku czasowym oraz ograniczenie w tej implementacji zapisywanych danych do położenia, prędkości i siły, co przekłada się na prostą formułę matematyczną i szybkość algorytmów.
Jak widać we wzorze 2.3 do wyznaczenia prędkości konieczne jest obliczenie siły działającej na atom w danym kroku czasowym. Tę siłę wyznacza się ze wzoru:
Ð→Fi= −∇iEp, (2.4)
gdzie ∇i jest gradientem w miejscu położenia i-tego atomu funkcji energii potencjalnej Ep
układu. Dzięki przybliżeniu Borna-Oppenheimera, które jest podstawowym założeniem mecha-niki molekularnej, układ molekularny możemy opisać na gruncie fizyki klasycznej.
Do opisu energii układów wykorzystywane są pola siłowe, na które, oprócz opisanej wcześniej funkcji potencjału składa się zbiór parametrów (opisanych na stronie 37). Wartości parametrów są przypisane do danego typu pola siłowego i wyznaczane są w wyniku obliczeń kwantowo-mechanicznych lub eksperymentalnie.
Postać funkcji potencjału wykorzystana w pracy składa się z pięciu członów, opisujących oddziaływania wiążące (3 pierwsze człony) i niewiążące (2 ostatnie człony):
V = ∑
kb, kθ - stałe harmoniczne wiązań i kątów walencyjnych
b, b0 - aktualna i równowagowa długość wiązania kowalencyjnego
2.1. Symulacje dynamiki molekularnej 37
θ, θ0 - aktualna i równowagowa wartość kąta walencyjnego ψ - aktualna wartość kąta torsyjnego
qi, qj, rij - ładunki na atomach i i j oraz odległość między nimi
, r∗ - głębokość studni potencjału i promień van der Waalsa
Oddziaływania wiążące są odpowiedzialne za zachowanie równowagowych wartości długości wiązań kowalencyjnych, kątów walencyjnych i kątów torsyjnych. Dwa człony opisujące oddziały-wania niewiążące to: oddziałyoddziały-wania elektrostatyczne, poprzez potencjał kulombowski, i oddzia-ływania van der Waalsa, poprzez potencjał Lennarda-Jonesa “6-12”. W przypadku oddziaływań niewiążących, w celu zredukowania problemu obliczeniowego stosuje się promień odcięcia Rc
(ang. R cuttoff ), co uzasadniane jest dążeniem tych oddziaływań do zera wraz z odległością.
Jest to niezwykle ważne uproszczenie, gdyż obliczenie oddziaływań niewiążących zajmuje naj-więcej czasu. Na potrzeby pracy zarówno dla oddziaływań kulombowskich i van der Waalsa przyjęto promień odcięcia równy 1 nm.
2.1.2 Parametry funkcji potencjału
Cząsteczki lipidów zostały sparametryzowane w oparciu o zmodyfikowane parametry OPLS-AA (ang. Optimized Potentials for LiquidS All Atoms) (Jorgensen i Swenson, 1985). Parame-tryzacja OPLS-AA jest zoptymalizowana tak, by jak najlepiej odtwarzać własności układów rzeczywistych. Dodatkowa modyfikacja parmetryzacji OPLS-AA polegała na zmianie ładunków dla atomów węgla i atomów wodoru w łańcuchach węglowodorowych, tak by odtwarzały litera-turowy moment dipolowy (Pasenkiewicz-Gierula i inni, 2012).
Dla cząsteczek wody wykorzystano parametryzację z modelu TIP3P (ang. Three-point Trans-ferable Intermolecular Potential) (Jorgensen i inni, 1983), który jest kompatybilny z parame-tryzacją OPLS-AA. W modelu tym występują trzy centra oddziałyania obdarzone ładunkiem punktowym, a atom tlenu posiada dodatkowo komplet parametrów Lennard-Jonesa.
Dla cząsteczek tlenu użyto parametrów z pola siłowego CHARMM-22 (MacKerel Jr. i inni, 1998). Ładunek na atomach tlenu wynosił zero, a równowagowa długość wiązania kowalencyjnego wynosiła 1,208 Å.
2.1.3 Parametry kontrolujące przebieg symulacji
Wszystkie symulacje, przedstawione w poniższej pracy, zostały wykonane za pomocą pro-gramu GROMACS v 4.0 (Hess i inni, 2008). Symulacje przeprowadzono w stałej temperaturze 310 K i przy stałym ciśnieniu równym jeden bar. Kontrola tych dwóch parametrów odbywa-ła się poprzez sprzężenie z odbywa-łaźnią temperaturową Nose-Hoover (Hoover, 1985) i ciśnieniową Parrinello-Rahman (Parrinello i Rahman, 1981) z czasami relaksacji równymi 1 ps. Ze wzglę-du na anizotropową kontrolę ciśnienia, układy symulacyjne miały możliwość zmiany wymiarów pudełek symulacyjnych.
38 2. Materiały i metody
2.1.4 Metoda Umbrella Sampling
Metoda Umbrella Sampling jest techniką umożliwiającą próbkowanie stanów układów, które ze względu na istniejące bariery energetyczne nie mogą zostać zbadane za pomocą kla-sycznej symulacji dynamiki molekularnej (Lemkul i Bevan, 2010). Metodą Umbrella Sampling wyznaczany jest potencjał średniej siły PMF (ang. Potential of Mean Force). Potencjał średniej siły PMF opisuje energię swobodną Gibbsa umieszczenia cząsteczki w danej lokalizacji. Orygi-nalnie PMF został wprowadzony w celu opisu przejść energetycznych w układach o wysokich barierach energetycznych. Obecnie jest on wykorzystywany do badania procesu dyfuzji małych cząsteczek w białkach globularnych (Cohen i inni, 2006), białkach transbłonowych (Bernèche i Roux, 2003) lub dwuwarstwach lipidowych (Bauer i inni, 2011).
W przedstawianej pracy wykorzystano metodę Umbrella Sampling do badania układów zawierających uwodnione dwuwarstwy lipidowe oraz pojedynczą próbkującą cząsteczkę tlenu.
Każdy badany układ symulacyjny błon lipidowych podzielony został na warstwy wzdłuż osi Z, prostopadłej do dwuwarstwy. Grubość warstw wynosiła 1 Å. Dla każdej z nich generowana była osobna konformacja początkowa systemu, w której próbkująca cząsteczka tlenu umieszczana była w centrum wybranej warstwy. Opis uzyskania tych konformacji dla badanych układów mieści się w paragrafie 2.2.4.
Podczas serii symulacji Umbrella Sampling cząsteczka próbkująca, w tym przypadku cząsteczka tlenu, związana jest za pomocą dodatkowego potencjału harmonicznego. Dla układów błonowych potencjał harmoniczny nakładany jest tylko na oś Z, przez co cząsteczki tlenu mają swobodę dyfuzji w płaszczyźnie XY. Standardem w wyznaczaniu potencjału średniej siły PMF w oparciu o zbiór symulacji Umbrella Sampling jest technika ważonych histogramów WHAM (ang. Weighted Histogram Analysis Method). W pracy wykorzystano implementację dostępną w pakiecie programów GROMACS (Hub i inni, 2010). Błędy oszacowania zostały wyznaczone za pomocą próby typu Bootstrap, polegającej na szacowaniu błędów w oparciu o wielokrotne losowanie ze zwracaniem z próby, na podstawie 100 niezależnych profili.