żytnictwa. Aby możliwa była jednak stabilna koegzystencja pasożytów i replikaz w analizowanym systemie, efektywniejsza replikacja spowodowana przez lepsze powi-nowactwo pasożyta do replikazy musi być zrekompensowana aby nie doprowadzić do wymarcia replikaz i w efekcie całego systemu. Rekompensatę dla wzrostu wartości 𝑎𝑃 zapewnia wzrost wartości 𝑙𝑃 obserwowany na wykresie 7.2. Im wyższa wartość 𝑙𝑃, tym więcej czasu pasożyt spędza w stanie zwiniętym, w którym nie może być powiela-ny. Wzrost wartości 𝑎𝑃 jest więc rekompensowany przez wzrost wartości 𝑙𝑃, dzięki temu pasożyt nie doprowadza do wymarcia replikaz, a cały system może stabilnie koeg-zystować.
Tabela 7.1 Wartości parametrów używanych w symulacjach wykonanych na podstawie modelu opartego o równania różniczkowe.
Parametr Opis Wartość
𝑎𝑅 powinowactwo replikazy do replikazy 0,6
𝑎𝑃 powinowactwo pasożyta do replikazy (wartość użyta w równaniach 7.1 – 7.5)
0,55 𝑙𝑃 prawdopodobieństwo zwinięcia cząsteczki pasożyta
(war-tość użyta w równaniach 7.1 – 7.5) 0,2 𝑎𝑃′ powinowactwo pasożyta do replikazy (wartość użyta
w równaniach 7.6 – 7.12)
0,7 𝑙𝑃′ prawdopodobieństwo zwinięcia cząsteczki pasożyta
(war-tość użyta w równaniach 7.6 – 7.12) 0,7 𝑎𝑃𝑚𝑢𝑡 powinowactwo zmutowanego pasożyta do replikazy 0,75
𝑙𝑃𝑚𝑢𝑡 prawdopodobieństwo zwinięcia cząsteczki zmutowanego
pasożyta 0,6
K stała replikacji 1
𝑅0 początkowy udział cząsteczek replikazy w całym modelo-wanym systemie
0,1 𝑃0 początkowy udział cząsteczek pasożyta w całym
modelo-wanym systemie
0,1 𝑃𝑚𝑢𝑡0 początkowy udział cząsteczek zmutowanego pasożyta
w całym modelowanym systemie
0,1
𝐶𝑅𝑃0 początkowy udział cząsteczek kompleksów replikazy z pasożytem w całym modelowanym systemie
0 𝐶𝑅𝑅0 początkowy udział cząsteczek kompleksów replikazy
z replikazą w całym modelowanym systemie
0
𝐶𝑅𝑃𝑚𝑢𝑡
0 początkowy udział cząsteczek kompleksów replikazy ze
zmutowanym pasożytem w całym modelowanym systemie 0
𝑑 stała rozpadu 0,02
Rys. 7.1 Wykres przedstawiający wyniki numerycznej symulacji układu równań różnicz-kowych (równania od 7.1 – 7.5) przeprowadzonej z wykorzystaniem funkcji NDSolve
pa-kietu Mathematica.
Rys. 7.2 Wykres przedstawiający wyniki numerycznej analizy stanów równowagi układu równań różniczkowych (równania od 7.1 – 7.5) przeprowadzonej w pakiecie Mathematica
dla różnych wartości parametru 𝒂𝑷.
Warto jednak zauważyć, że z punktu widzenia ewolucji biologicznej najlepiej dostosowane cząsteczki:
powinny posiadać jak najwyższą wartość powinowactwa do replikaz, zatem maksymalizowana jest wartość parametru 𝑎
przebywają jak najwięcej w stanie rozwiniętym, co oznacza minimalizację 𝑙, dzięki czemu są częściej dostępne dla replikaz.
Gdyby zatem uwzględnić mutację, selekcja powinna faworyzować cząsteczki o takich właśnie parametrach i w efekcie powinniśmy zaobserwować przetrwanie najlepiej do-stosowanej populacji cząsteczek. Chociaż w równaniach różniczkowych modelujących zachowanie opisanego w niniejszej pracy systemu RP nie zawarto parametrów odpo-wiedzialnych za mutacje cząsteczek, to aby zbadać sytuację, w której w systemie poja-wia się lepiej dostosowana populacja, układ równań 7.1 – 7.5 zmodyfikowano tak, aby dodać do niego równania opisujące populację cząsteczek pasożytniczych o faworyzo-wanych ewolucyjnie parametrach, więc zwiększonej wartości 𝑎𝑃 i zmniejszonej warto-ści 𝑙𝑃. Układ równań uwzględniający dwie populacje pasożytów przedstawiony jest poniżej.
𝑹̇ = −𝟐𝒂𝑹𝑹𝟐+ [𝟐(𝟏 − 𝒂𝑹) + 𝟑K𝝎 + 𝟐𝒅]𝑪𝑹𝑹− 𝒂𝑷′(𝟏 − 𝒍𝑷′)𝑹
− 𝒂𝑷𝒎𝒖𝒕(𝟏 − 𝒍𝑷𝒎𝒖𝒕)𝑹𝑷𝒎𝒖𝒕+ [(𝟏 − 𝒂𝑷′) + K𝝎 + 𝒅]𝑪𝑹𝑷
+ [(𝟏 − 𝒂𝑷𝒎𝒖𝒕) + K𝝎 + 𝒅]𝑪𝑹𝑷𝒎𝒖𝒕− 𝒅𝑹 (7.6) 𝑷̇ = −𝒂𝑷′(𝟏 − 𝒍𝑷′)𝑹𝑷 + [(𝟏 − 𝒂𝑷′) + 2K𝝎 + 𝒅]𝑪𝑹𝑷− 𝒅𝑷 (7.7) 𝑷̇𝒎𝒖𝒕= −𝒂𝑷𝒎𝒖𝒕(𝟏 − 𝒍𝑷𝒎𝒖𝒕)𝑹𝑷𝒎𝒖𝒕+ [(𝟏 − 𝒂𝑷𝒎𝒖𝒕) + 2K𝝎 + 𝒅]𝑪𝑹𝑷𝒎𝒖𝒕
− 𝒅𝑷𝒎𝒖𝒕
(7.8)
𝑪̇𝑹𝑹 = 𝒂𝑹𝑹𝟐− [(𝟏 − 𝒂𝑹) + K𝝎]𝑪𝑹𝑹− 𝟐𝐝𝑪𝑹𝑹 (7.9) 𝑪̇𝑹𝑷= 𝒂𝑷′(𝟏 − 𝒍𝑷′)𝑹𝑷 − [(𝟏 − 𝒂𝑷′) + K𝝎]𝑪𝑹𝑷− 𝟐𝐝𝑪𝑹𝑷 (7.10) 𝑪̇𝑹𝑷𝒎𝒖𝒕 = 𝒂𝑷𝒎𝒖𝒕(𝟏 − 𝒍𝑷𝒎𝒖𝒕)𝑹𝑷𝒎𝒖𝒕− [(𝟏 − 𝒂𝑷𝒎𝒖𝒕) + K𝝎]𝑪𝑹𝑷𝒎𝒖𝒕− 𝟐𝐝𝑪𝑹𝑷𝒎𝒖𝒕 (7.11) gdzie:
𝝎 = 𝟏 − 𝑹 − 𝑷 − 𝑪𝑹𝑷− 𝑪𝑹𝑹− 𝑷𝒎𝒖𝒕− 𝑪𝑹𝑷𝒎𝒖𝒕 (7.12)
Układ równań 7.6 – 7.12 ponownie rozwiązany został numerycznie za pomocą programu Mathematica. Wykresy przedstawiające zmianę ilości cząsteczek poszczegól-nego typu w czasie pokazano na rysunku 7.3. Wyniki pokazują, że pasożyt z lepszym powinowactwem do replikazy, który więcej czasu jest dostępny jako matryca przewyż-sza liczebnie pasożyta o mniej korzystnych wartościach tych parametrów. Gdyby zatem dozwolona była mutacja rozważanych w tym podrozdziale parametrów, dobrze wymie-szany system opisywany przez ODE dążyłby do wymarcia na skutek ewolucji zbyt szkodliwych pasożytów.
Rys. 7.3 Wykres przedstawiający wyniki numerycznej symulacji układu równań różnicz-kowych (równania od 7.6 –7.12) przeprowadzonej z wykorzystaniem funkcji NDSolve
pa-kietu Mathematica.
W rozdziale tym zaprezentowany został model RP wykorzystujący równania różniczkowe. Model ten uwzględnia nie tylko tworzenie się cząsteczek, ale także bierze pod uwagę punkt pośredni w reakcji replikacji, jakim jest tworzenie kompleksów po-między cząsteczkami, co jest z kolei uzależnione od wartości parametrów opisujących modelowane cząsteczki. Ponadto model uwzględnia rozpad cząsteczek. W ramach prac nad doktoratem sformułowane zostały układy równań opisujące ilościowe zmiany czą-steczek w czasie. Określono także wartości analizowanych parametrów zapewniające stabilność opisywanego w tym rozdziale modelu oraz kierunek zmian wartości
parame-7.3 Podsumowanie
tru obrazującego dostępność cząsteczki dla replikazy 𝑙𝑃 wraz ze zmianą powinowactwa cząsteczki do replikazy 𝑎𝑃.
Z poprzedniego rozdziału opisującego model RP w oparciu o ETG wynikało, że zróżnicowanie tempa replikacji pasożyta i replikazy może teoretycznie umożliwić ich stabilną koegzystencję. Prezentowane w tym rozdziale wyniki pokazują, że rzeczywi-ście jest to możliwe (rysunki 7.1 i 7.2), jednak w przypadku gdy uwzględnilibyśmy mu-tacje cząsteczek RNA, tak jak to zachodzi ona w rzeczywistości, wtedy obserwujemy przetrwanie cząsteczki najlepiej dostosowanej (rysunek 7.3). Zatem uwzględniając mu-tacje obserwować będziemy ewolucję coraz bardziej szkodliwych pasożytów, o wyższym tempie replikacji, które sprawią, że replikazy będą coraz bardziej wykorzy-stywane do powielania pasożytów zamiast siebie samych. W efekcie replikazy wyginą, a w konsekwencji wymrze cały system.
Jednak w realnym świecie pierwotne życie oparte o cząsteczki RNA nie tylko przetrwało, ale także rozwinęło się w życie, które obserwujemy dzisiaj, a które charak-teryzuje się niezwykłym zróżnicowaniem i złożonością nawet w przypadku najprost-szych jego form. Równania różniczkowe chociaż niewątpliwie są bardzo użytecznym narzędziem pozwalającym na modelowanie obserwowanych zjawisk, to jednak zakłada-ją środowisko dobrze wymieszane oraz ciągłość analizowanych wartości. W praktyce może to doprowadzić do otrzymania rozwiązań, które są nierealistyczne. Dlatego w celu zbadania, czy ograniczenia modeli różniczkowych wpływają znacząco na wyniki, a tym samym możliwość przetrwania systemu złożonego z pasożytów i replikaz, kolejny roz-dział został poświęcony modelowaniu systemu RP za pomocą automatów komórko-wych. Metoda ta uwzględnia zarówno przestrzenny aspekt rzeczywistości, jak i pozwala rozważać dyskretne wartości liczebności populacji analizowanych cząsteczek. Otrzy-mane w ramach analiz ODE punkty bistabilne posłużyły jako punkt wyjścia do ustale-nia parametrów używanych w modelu CA przedstawiony w następnym rozdziale.
Model wykorzystujący automat komórkowy
W niniejszym rozdziale przedstawiony zostanie model RP wykorzystujący me-todę automatów komórkowych (CA), która pozwala na obejście ograniczeń równań różniczkowych. Model CA uwzględnia wszystkie aspekty systemu RP modelowane przez równania różniczkowe przedstawione w poprzednim rozdziale, a więc tworzenie się kompleksów replikacyjnych, ich dysocjację, powstawanie nowych cząsteczek oraz śmierć cząsteczek. Cząsteczkom, podobnie jak w przypadku modelu ODE, przypi-sane są atrybuty powinowactwa cząsteczki do replikazy, 𝑎, oraz prawdopodobieństwo przebywania przez cząsteczkę w stanie zwiniętym 𝑙.
W przeciwieństwo do modelu ODE, model CA w sposób bezpośredni pozwolił na uwzględnienie mutacji cząsteczek podczas replikacji. Mutowane są parametry opisu-jące cząsteczki, a więc 𝑎 oraz 𝑙. Dodatkowo model oparty o CA dzięki temu, że pozwa-la na modelowanie przestrzennych interakcji, pozwapozwa-la na modelowanie przemieszczania się cząsteczek, a więc dyfuzji. Poza tym, że metoda CA uwzględnia przestrzenny aspekt rzeczywistości, pozwala także rozważać dyskretne wartości liczebności populacji anali-zowanych cząsteczek. Model oparty o automat komórkowy jest dwuwymiarowy, co jest zgodne z pracami dowodzącymi, że umiejscowienie łańcuchów RNA na powierzchni minerałów umożliwia wydajniejsze tworzenie łańcuchów RNA (Ferris et al. 1996; Luther et al. 1998; Orgel 1998; von Kiedrowski i Szathmáry 2001).
Analiza wyników uzyskanych w tym rozdziale polega na obserwacji zmian śred-nich wartości parametrów powinowactwa cząsteczki do replikazy oraz prawdopodobieństwo przebywania przez cząsteczkę w stanie zwiniętym w miarę upływu czasu symulacji. Ponadto istotna jest obserwacja zachowania analizowanego systemu całościowo. Jest to możliwe dzięki wizualizacji, która pozwala
na obserwowanie form wyłaniających się pod postacią wzorów pojawiających się w czasie symulacji na planszy symulacyjnej.
Model wykorzystujący automat komórkowy analizowany w tym rozdziale oparty jest o prace Nobuto Takeuchiego oraz Paulien Hogeweg nad początkami życia na Ziemi (Hogeweg i Takeuchi 2003). Analiza stabilności systemu RP sformułowanego w oparciu o równania różniczkowe przedstawiona w rozdziale 7 pozwoliła na ograniczenie wartości analizowanych parametrów wpływających na replikację czą-steczek w przedstawianym w niniejszym rozdziale modelu opartym o automat komór-kowy. Prezentowane wyniki zostały w całości uzyskane za pomocą wykonanej przez autorkę pracy doktorskiej implementacji modelu RP wykorzystującego automat komór-kowy zaimplementowanej w języku NetLogo. Wyniki przedstawione w tym rozdziale zostały zaprezentowane na konferencji Polskiego Towarzystwa Bioinformatycznego w 2015.
Model systemu RP oparty o automat komórkowy składa się z siatki pól ze skończoną liczbą stanów dla każdego z tych pól. Jeśli zgodnie z definicją przedsta-wioną w rozdziale 2.3 przez 𝑆 oznaczymy skończony zbiór stanów, jaki może przyjąć komórka, wtedy w przypadku modelu RP, 𝑆 = {𝑅, 𝑃, 𝑅𝑅′, 𝑅𝑃, 𝑃𝑅, 𝜔}, gdzie 𝑅 to stan nieskompleksowanej replikazy, 𝑃 stan nieskompleksowanego pasożyta, a 𝜔 oznacza wolne miejsce, czyli modeluje zasoby potrzebne do replikacji. Stan 𝑅𝑅′ oznacza, że dane pole w stanie 𝑅 jest powiązane z innym polem typu 𝑅′ z jego sąsiedztwa, co mode-luje tworzenie się kompleksów dwóch replikaz. Analogicznie stan 𝑅𝑃 oznacza, że dane pole w stanie 𝑅 jest powiązane z innym polem typu 𝑃 (lub na odwrót) z jego sąsiedz-twa, co modeluje tworzenie się kompleksów replikaz z pasożytami. Wyodrębnienie sta-nów skompleksowanych jest istotne ze względu na istnienie konkretnych reguł właści-wych tylko dla stanu skompleksowanego. Rozważane w tym modelu jest sąsiedztwo Moore’a składające się z 8 komórek otaczających rozpatrywane pole. Dodatkowo w przypadku dyfuzji kompleksów rozpatrywane jest sąsiedztwo rozszerzone obejmują-ce także pole obejmują-centralne zajmowane przez fragment skompleksowanego agenta. Mode-lowana przestrzeń jest toroidalna, aby uniknąć efektu kumulacji agentów przy krawędzi siatki.