Kolejny wykład poświęcimy na zapoznanie się z nieco bardziej skomplikowanym modelem odpowiedzi odpornościowej. W 1975 r. Guri I. Marczuk (właściwie G. I. Marchuk według obowiązującej transkrypcji rosyjskiego nazwiska, Marczuk to spolszczona wersja tego nazwiska) zaproponował dość prosty model, który jest w stanie odzwierciedlić najważniejsze przebiegi odpowiedzi odpornościowej, w tym również chorobę przewlekłą (chroniczną). Model opisany jest za pomocą czterech równań z opóźnionym argumentem i z założenia odzwierciedla humoralną odpowiedź odpornościową, choć podobną strukturę modelu można także zaproponować w przypadku odpowiedzi komórkowej. Model ten będziemy w dalszym ciągu nazywać modelem Marczuka. W 1980 r. ukazała się książka Marczuka [14], w której oprócz tego najprostszego modelu przedstawił także kilka wersji bardziej rozbudowanych, w tym model opisujący równoległy przebieg odpowiedzi humoralnej i komórkowej. W późniejszych latach ukazały się także tłumaczenia tej książki na język angielski [11] i polski [12]. Wyniki doświadczalne dla pewnych chorób okazały się bardzo podobne do teoretycznych wyników uzyskanych na podstawie modelu Marczuka.
Z powstaniem modelu wiąże się ciekawa historia. Marczuk przez wiele lat był przewodniczącym Akademii Nauk ZSSR. Niestety zachorował na żółtaczkę zakaźną, a pobyt w szpitalu urozmaicał sobie poznawaniem mechanizmów odpornościowych. W rezultacie doprowadziło to do powstania modelu reakcji odpornościowej i zaowocowało wieloma pracami naukowymi w tej dziedzinie. Po powrocie do zdrowia Marczuk w dalszym ciągu prowadził działalność naukową, w ramach której odwiedził Wydział Matematyki, Informatyki i Mechaniki UW, gdzie zaprezentował model i wyniki symulacji wskazujące możliwości leczenia choroby przewlekłej za pomocą zaostrzenia infekcji. Zarówno sam model, jak i możliwości jego zastosowania, zafascynowały jednego z pionierów polskiej biomatematyki, prof. Wiesława Szlenka. Tematyka ta znalazła się później wśród zagadnień badawczych w powstałym po pewnym czasie Instytucie Matematyki Stosowanej i Mechaniki.
Zajmiemy się teraz omówieniem konstrukcji prostego modelu odpowiedzi odpornościowej zaproponowanego przez G. I. Marczuka, a następnie przedstawimy podstawowe własności układu.
Przypomnijmy, że reakcja odpornościowa następuje w wyniku wniknięcia do organizmu obcego czynnika, który nazywamy antygenem. W przypadku omawianego modelu antygenem są głównie bakterie, ponieważ np. wirusy wywołują znacznie bardziej skomplikowane reakcje, które wymagają bardziej skomplikowanego opisu matematycznego, w szczególności uwzględnienia także komórkowej odpowiedzi odpornościowej. Pojawienie się antygenu wywołuje kaskadę reakcji, które w efekcie doprowadzają do rozpoznania antygenu i uruchomienia produkcji przeciwciał swoistych przez komórki plazmatyczne. Mechanizmem uruchamiającym tę kaskadę jest powstanie tzw. kompleksu VT, czyli przyłączenie do antygenu przeciwciała prezentującego. Wobec tego zakładamy, że produkcja komórek plazmatycznych zależy od liczby wiązań antygen – przeciwciało, przy czym obserwujemy opóźnienie między powstaniem takiego wiązania a wyprodukowaniem przed układ odpornościowy komórek plazmatycznych produkujących przeciwciała swoiste.
W modelu występują następujące podstawowe zmienne
— zagęszczenie populacji bakterii (antygenu) w organizmie
w chwili
,
— zagęszczenie komórek plazmatycznych w chwili
,
— stężenie przeciwciał w chwili
,
przy czym (jak dotychczas) używamy określenia zagęszczenie w stosunku do osobników pewnej populacji (tu populacji komórkowych), natomiast stężenie dotyczy głównie substancji, a wiemy, że przeciwciała są białkami.
Dynamikę antygenu opisujemy w modelu Marczuka w taki sam sposób jak w najprostszym modelu odpowiedzi odpornościowej prezentowanym poprzednio. Zakładamy, że bakterie rozmnażają się za stałym współczynnikiem reprodukcji i giną w wyniku reakcji odpornościowych. Zatem śmiertelność spowodowana jest wiązaniem antygenów z przeciwciałami, proporcjonalnie do liczby tych
wiązań ze współczynnikiem
, który odzwierciedla prawdopodobieństwo związania antygenu z przeciwciałem i siłę ich reakcji.
Stąd
![]() |
Przeciwciała giną nie tylko w wyniku wiązania z antygenem, ale też w procesie naturalnej degradacji. Są one produkowane przez komórki plazmatyczne w stałym tempie , czyli
![]() |
gdzie opisuje frakcję związanych przeciwciał ginących w reakcji odpornościowej, a
odzwierciedla średni czas przetrwania
przeciwciała.
Zauważyliśmy wcześniej, że produkcja komórek plazmatycznych uruchamiana jest przez powstanie kompleksów
i jest opóźniona w stosunku do mechanizmu uruchamiającego.
Dla uproszczenia zakładamy, że opóźnienie to jest stałe, równe
.
W organizmie utrzymuje się też stały fizjologiczny poziom przeciwciał
. Skoro tak, to jeśli nie ma antygenu, wewnętrzną dynamikę komórek plazmatycznych możemy opisać równaniem
![]() |
gdzie oznacza średnią długość życia komórki plazmatycznej.
Wtedy
przy
, co oznacza, że w zdrowym organizmie po odchyleniu od poziomu fizjologicznego zawsze następuje samoczynny powrót do tego poziomu.
Jeśli antygen występuje w organizmie, to w powyższym równaniu uwzględniamy jeszcze składnik związany ze stymulacją odpowiedzi odpornościowej
![]() |
przy czym jest współczynnikiem immunogenności antygenu odzwierciedlającym siłę odpowiedzi odpornościowej na dany antygen.
Prócz tych trzech podstawowych zmiennych w modelu występuje zmienna dodatkowa, która opisuje stopień uszkodzenia organu zaatakowanego przez antygen, ponieważ Marczuk budował swój model z myślą o żółtaczce i stopniu uszkodzenia wątroby w wyniku tej choroby.
Niech oznacza charakterystykę (masę lub powierzchnię) zdrowego organu, a
— zdrową część organu uszkodzonego przez antygen. Definiujemy
![]() |
Mamy zatem , co opisuje stosunek części uszkodzonej do całości. Oczywiście
dla organu zdrowego oraz
dla organu całkowicie zniszczonego. Dynamikę
opisujemy za pomocą prostego równania liniowego
![]() |
gdzie odzwierciedla zdolność antygenu do niszczenia organu, a
oznacza średni czas regeneracji organu.
Wielkość osłabia odpowiedź odpornościową poprzez zmniejszenie stymulacji produkcji komórek plazmatycznych, zatem zamiast współczynnika immunogenności
rozpatrujemy w modelu współczynnik
, gdzie
jest nierosnącą funkcją zmiennej
,
i
.
Ostatecznie model Marczuka humoralnej odpowiedzi odpornościowej opisujemy układem 4. równań różniczkowych z opóźnionym argumentem (rrzoa)
![]() |
(13.1) |
W większości przypadków istnieje pewna progowa wielkość , do której
nie wpływa znacząco na reakcję odpornościową i można założyć, że
dla
.
Przejdziemy teraz do zbadania podstawowych własności układu (13.1). Jest on układem równań z opóźnionym argumentem, zatem jako warunek początkowy należy zadać funkcję określoną na odcinku o wartościach w
. Ponieważ jednak tylko współrzędne
i
występują z opóźnionym argumentem, więc możemy zadać tylko ich wartości na całym przedziale
, ograniczając się do wartości
i
dla zmiennych
i
. Należy przy tym pamiętać, że z
formalnego punktu widzenia potrzebne są wszystkie cztery współrzędne funkcji początkowej.
Istnienie, jednoznaczność i nieujemność rozwiązań dla nieujemnego warunku początkowego możemy wykazać stosując metodę kroków. Przypomnijmy, że opiera się ona na zasadzie indukcji matematycznej i pozwala zastosować znane techniki dotyczące równań różniczkowych zwyczajnych do analizy rrzoa.
Niech ,
. Wtedy dla
układ (13.1) sprowadza się do niejednorodnego układu rrz postaci
![]() |
(13.2) |
gdzie jest znaną funkcją, gdyż dla
argument
. Funkcje
opisujące współrzędne prawej strony układu (13.2)
są wielomianami kwadratowym względem współrzędnych, czyli są klasy
(faktycznie
), więc spełniają założenia twierdzeń o lokalnym istnieniu i jednoznaczności rozwiązań równań różniczkowych zwyczajnych. Istnieje zatem odcinek
, na którym mamy jednoznaczne rozwiązanie układu.
Pokażemy najpierw, że rozwiązanie to jest nieujemne. Pierwsze równanie układu (13.2) zapisujemy w równoważnej postaci całkowej
![]() |
gdzie w powyższym wzorze .
Zatem
dla
, co więcej,
dla
, czyli z matematycznego punktu widzenia wyzdrowienie może wystąpić tylko asymptotycznie. W praktyce zagęszczenie antygenu niższe od czułości urządzeń pomiarowych jest niewykrywalne i takie zagęszczenia uznajemy za równoważne wyzdrowieniu.
Dla drugie równanie układu (13.2) jest równaniem liniowym niejednorodnym. Rozwiązujemy je stosując metodę uzmienniania
stałej i otrzymujemy
![]() |
(13.3) |
Ponieważ
oraz
dla
, to
dla
.
Zauważmy też, że
w tym przypadku współrzędna
jest zdefiniowana za pomocą wzoru (13.3) na całym przedziale
.
W ogólnym przypadku, gdy ,
oszacujemy wartość pochodnej
z dołu przez
, więc
![]() |
Zatem dla
.
Zauważmy także, że jeśli
, to z powyższych nierówności
wynika, że
dla dowolnego
.
W analogiczny sposób dowodzimy nieujemności zmiennej .
Skoro dla
, to
podobnie oszacujemy
, czyli
![]() |
Zatem dla
.
Pokażemy teraz, że z nieujemności rozwiązań wynika ich istnienie na całym przedziale . W tym celu oszacujemy nasze rozwiązania z góry.
Dla pierwszej zmiennej mamy
![]() |
zatem wzrost jest co najwyżej wykładniczy.
Dla zmiennej mamy
![]() |
czyli jest również ograniczona na każdym podprzedziale odcinka
.
Z ograniczoności i
wynika ograniczoność
i
. Skoro funkcje te są ograniczone, więc ograniczone są także ich pochodne, czyli rozwiązanie można przedłużyć biorąc jako nowy warunek początkowy granicę
(istnienie granicy wynika bezpośrednio z ograniczoności pochodnej) oraz punkt
i stosując twierdzenie o istnieniu rozwiązań. Wobec dowolności
wnioskujemy, że rozwiązanie istnieje na całym
.
Stosując metodę kroków zakładamy, że istnieje nieujemne rozwiązanie
układu (13.1) na odcinku i powtarzamy powyższe rozumowanie w celu wykazania istnienia nieujemnego rozwiązania na kolejnym
odcinku
. Na mocy zasady indukcji matematycznej
jednoznaczne, nieujemne rozwiązanie układu (13.1)
istnieje dla dowolnego
.
Znajdziemy teraz punkty krytyczne badanego układu. Rozwiązujemy następujący układ równań algebraicznych
![]() |
(13.4) |
W ogólnym przypadku, aby rozwiązać układ (13.4), musimy zadać konkretną postać funkcji .
Dla uproszczenia załóżmy więc, że
, czyli
. W tym przypadku badanie możemy ograniczyć do układu trzech równań, a następnie sprawdzić zachowanie ostatniej zmiennej.
Łatwo widać, że układ (13.4) ma rozwiązanie opisujące stan zdrowia, czyli i
. Wtedy
i oczywiście
.
nazwiemy poziomem fizjologicznym przeciwciał.
Niech odzwierciedla stan zdrowia organizmu.
Jeśli , to z pierwszego równania układu (13.4) mamy
. Wstawiając tę zależność do drugiego i trzeciego równania obliczmy
![]() |
Punkt krytyczny odzwierciedla chorobę przewlekłą, oczywiście pod warunkiem, że wszystkie jego współrzędne są dodatnie.
Współrzędne
i
możemy przepisać w nieco innej postaci
![]() |
z czego wynika, że jeśli
i
albo
i
,
to .
Zauważmy dalej, że przy warunku A mamy:
![]() |
(13.5) |
zatem . Przy warunku B nierówność (13.5) zmienia znak, ale znak zmienia też mianownik
, czyli
. Widzimy więc, że warunki A, B są konieczne i dostateczne dla istnienia dodatniego punktu krytycznego. Zauważmy także, że w tym przypadku
oraz
.
Ostatnia nierówność oznacza, że w omawianym przypadku
poziom stacjonarny antygenu
nie może być zbyt duży. W literaturze mówi się o lekkiej chorobie przewlekłej dla
, natomiast gdy
— o ciężkiej.
Badanie lokalnej stabilności punktów krytycznych rrzoa nie różni się w początkowej fazie od analogicznej analizy dla rrz. Niech oznacza punkt krytyczny. Zdefiniujmy nowe zmienne:
,
,
. Układ w nowych zmiennych ma postać
![]() |
(13.6) |
Linearyzujemy układ (13.6) wokół punktu zakładając, że
,
,
dla
.
Składniki wyższego rzędu, czyli
i
szacują się przez
, co
możemy pominąć. Po zlinearyzowaniu układ (13.6) przyjmuje postać
![]() |
(13.7) |
W kolejnym kroku znajdujemy rozwiązania układu (13.7) postaci ,
,
. Otrzymujemy więc następujące zależności
![]() |
(13.8) |
Układ (13.8) ma nietrywialne rozwiązania, jeśli
![]() |
(13.9) |
Twierdzenie o linearyzacji gwarantuje, że jeśli wartości własne, czyli rozwiązania równania (13.8) mają części rzeczywiste ujemne, to rozwiązanie układu (13.6) jest lokalnie asymptotycznie stabilne, natomiast jeśli istnieje wartość własna o części rzeczywistej dodatniej, to rozwiązanie jest niestabilne.
Zbadamy wartości własne dla poszczególnych punktów krytycznych. Zaczniemy od . Wtedy macierz układu ma postać dolnietrójkątną, zatem łatwo wyznaczamy wartości własne, gdyż spełniają równanie
![]() |
czyli wartości własne są równe: ,
,
, zatem są rzeczywiste i ujemne dla
.
Jeśli , to rozwiązanie
jest lokalnie asymptotycznie stabilne. Natomiast jeśli
, to
jest niestabilne.
Widzimy, że stabilność rozwiązania , czyli wyzdrowienie jest możliwe tylko przy niewielkiej reproduktywność antygenu
w stosunku do poziomu fizjologicznego przeciwciał
.
Badając stabilność rozwiązania otrzymujemy następujący pseudowielomian charakterystyczny
![]() |
(13.10) |
gdzie
![]() |
(13.11) |
a wyznaczenie współczynników pozostawiamy jako ćwiczenie. Przypomnijmy, że dla pseudowielomian ma nieskończenie wiele pierwiastków, co utrudnia jego analizę (podczas gdy dla wielomianów możemy skorzystać np. z kryt. Routha – Hurwitza).
Wiemy, że jeśli rozwiązanie
jest stabilne dla , to z ciągłej zależności od parametrów
wynika, że albo rozwiązanie jest stabilne dla wszystkich
,
albo istnieje krytyczna wartość parametru
, dla której może nastąpić destabilizacja.
Dla
możemy zbadać stabilność punktu
za pomocą kryterium Routha – Hurwitza, skąd dostajemy warunek stabilności
![]() |
(13.12) |
Warunek jest równoważny
, zatem stabilność
dla
możliwa jest tylko w przypadku A. Natomiast w przypadku B mamy stabilność punktu
.
Zauważmy, że zmiana stabilności może nastąpić tylko gdy pseudowielomian charakterystyczny (13.10) ma pierwiastki czysto urojone . Wtedy
![]() |
czyli
![]() |
Zdefiniujmy funkcję pomocniczą
![]() |
Jej dodatnie miejsca zerowe determinują potencjalne zmiany stabilności punktu krytycznego . Z warunku (13.12) wynika, że wyraz wolny wielomianu
jest ujemny, zatem istnieje co najmniej jedno dodatnie miejsce zerowe.
Warunek A lub B implikuje też
, czyli dodatnie miejsce zerowe
jest wyznaczone jednoznacznie. Niech
, gdzie
. Podstawiając
do pseudowielomianu (13.10) dostajemy układ równań, z którego wyznaczamy sinus i kosinus wartości krytycznej
. Mamy zatem krytyczną wartość opóźnienia zdefiniowaną w sposób uwikłany za pomocą tych funkcji trygonometrycznych. Skoro jednak znamy tylko wartości sinusa i kosinusa, to oczywiście mamy nie jedną wartość krytyczną, a cały ciąg takich wartości
z pierwszą wartością
, taką że
i
.
Zmiana stabilności następuje przy takiej krytycznej wartości , dla której istnieje para czysto urojonych wartości własnych
i przy przekraczaniu osi urojonej wartości własne przechodzą z lewej półpłaszczyzny zespolonej na prawą, czyli
. Przypomnijmy, że znak tej pochodnej jest taki sam jak znak pochodnej funkcji pomocniczej
w odpowiednim miejscu zerowym — tutaj jedynym miejscu zerowym
. W tej sytuacji mamy oczywiście dodatnią pochodną
i wartości własne przechodzą z lewej półpłaszczyzny zespolonej na prawą.
Jeśli spełnione są nierówności (13.12) i istnieje dodatni punkt krytyczny , to punkt ten jest stabilny dla
i istnieje krytyczna wartość opóźnienia
, przy której
traci stabilność i pozostaje niestabilny dla
. W
zachodzi bifurkacja Hopfa, czyli pojawiają się nietrywialne rozwiązania okresowe o okresie
.
Przypomnijmy, że z nierówności wynika, że stabilność punktu
przy
jest możliwa tylko w przypadku spełnienia warunku A.
Korzystając z kryterium Michajłowa, można udowodnione następujące twierdzenie dotyczące stabilności punktu , które daje nam oszacowanie wielkości opóźnienia krytycznego
Niech oraz
![]() |
(13.13) |
Jeśli dodatkowo spełnione są nierówności (13.12) i istnieje dodatni punkt krytyczny , to punkt ten jest lokalnie asymptotycznie stabilny.
Twierdzenie 13.1 zadaje warunki wystarczające stabilności punktu przy ustalonych parametrach. Warunki te implikują, że hodograf Michajłowa ma kształt przedstawiony na rys. 13.2.
Widzimy więc, że układ (13.1) jest w stanie odzwierciedlić dwa podstawowe przebiegi choroby: zakażenie z wyzdrowieniem oraz chorobę przewlekłą (chroniczną), por. rys. 13.3. Ciekawą własnością symulacji tego modelu jest to, że zwiększając dawkę antygenu możemy wyprowadzić organizm z choroby chronicznej — mówimy o wyzdrowieniu poprzez zaostrzenie choroby, por. rys. 13.4. Trzeba jednak podkreślić, że taki efekt obserwuje się tylko w symulacjach, natomiast wszystkie wyniki analityczne wskazują na globalną stabilność punktu przy założeniu
.
Treść automatycznie generowana z plików źródłowych LaTeXa za pomocą oprogramowania wykorzystującego LaTeXML.
strona główna | webmaster | o portalu | pomoc
© Wydział Matematyki, Informatyki i
Mechaniki UW, 2009-2010. Niniejsze materiały są udostępnione bezpłatnie na licencji Creative Commons Uznanie autorstwa-Użycie niekomercyjne-Bez utworów zależnych 3.0 Polska.
Projekt współfinansowany przez Unię Europejską w ramach
Europejskiego Funduszu Społecznego.
Projekt współfinansowany przez Ministerstwo Nauki i Szkolnictwa Wyższego i przez Uniwersytet Warszawski.