Zagadnienia

13.1 Prezentacja modelu Marczuka
13.2 Podstawowe własności układu ()

13. Model Marczuka humoralnej odpowiedzi odpornościowej

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.

13.1. Prezentacja modelu Marczuka

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

$V(t)$ — zagęszczenie populacji bakterii (antygenu) w organizmie w chwili $t$ ,
$C(t)$ — zagęszczenie komórek plazmatycznych w chwili $t$ ,
$F(t)$ — stężenie przeciwciał w chwili $t$ ,

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 $\beta>0$ 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 $\gamma>0$ , który odzwierciedla prawdopodobieństwo związania antygenu z przeciwciałem i siłę ich reakcji. Stąd

$\dot{V}(t)=\beta V(t)-\gamma V(t)F(t).$

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 $\varrho>0$ , czyli

$\dot{F}(t)=\varrho C(t)-\eta\gamma V(t)F(t)-\mu _{F}F(t),$

gdzie $\eta>0$ opisuje frakcję związanych przeciwciał ginących w reakcji odpornościowej, a $1/\mu _{F}$ odzwierciedla średni czas przetrwania przeciwciała.

Zauważyliśmy wcześniej, że produkcja komórek plazmatycznych uruchamiana jest przez powstanie kompleksów $VT$ i jest opóźniona w stosunku do mechanizmu uruchamiającego. Dla uproszczenia zakładamy, że opóźnienie to jest stałe, równe $\tau$ . W organizmie utrzymuje się też stały fizjologiczny poziom przeciwciał $C^{*}$ . Skoro tak, to jeśli nie ma antygenu, wewnętrzną dynamikę komórek plazmatycznych możemy opisać równaniem

$\dot{C}(t)=\mu _{C}(C^{*}-C(t)),$

gdzie $1/\mu _{C}$ oznacza średnią długość życia komórki plazmatycznej. Wtedy $C(t)=C^{*}+(C_{0}-C^{*})\text{e}^{{-\mu _{C}t}}\to C^{*}$ przy $t\to\infty$ , 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

$\dot{C}(t)=\alpha V(t-\tau)F(t-\tau)-\mu _{C}(C(t)-C^{*}),$

przy czym $\alpha>0$ 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 $M$ oznacza charakterystykę (masę lub powierzchnię) zdrowego organu, a $M(t)$ — zdrową część organu uszkodzonego przez antygen. Definiujemy

$m(t)=\frac{M-M(t)}{M}.$

Mamy zatem $m\in[0,1]$ , co opisuje stosunek części uszkodzonej do całości. Oczywiście $m=0$ dla organu zdrowego oraz $m=1$ dla organu całkowicie zniszczonego. Dynamikę $m(t)$ opisujemy za pomocą prostego równania liniowego

$\dot{m}(t)=\sigma V(t)-\mu _{m}m(t)\ \ \mbox{dla}\ \ m(t)\leq 1,$

gdzie $\sigma>0$ odzwierciedla zdolność antygenu do niszczenia organu, a $1/\mu _{m}$ oznacza średni czas regeneracji organu.

Wielkość $m(t)$ osłabia odpowiedź odpornościową poprzez zmniejszenie stymulacji produkcji komórek plazmatycznych, zatem zamiast współczynnika immunogenności $\alpha$ rozpatrujemy w modelu współczynnik $\alpha\xi(m)$ , gdzie $\xi(m)$ jest nierosnącą funkcją zmiennej $m$ , $\xi(0)=1$ i $\xi(1)=0$ .

Ostatecznie model Marczuka humoralnej odpowiedzi odpornościowej opisujemy układem 4. równań różniczkowych z opóźnionym argumentem (rrzoa)

$\begin{array}[]{l c l}\dot{V}(t)&=&(\beta-\gamma F(t))V(t),\\ \dot{C}(t)&=&\alpha\xi(m(t))V(t-\tau)F(t-\tau)-\mu _{C}(C(t)-C^{{*}}),\\ \dot{F}(t)&=&\varrho C(t)-(\mu _{F}+\eta\gamma V(t))F(t),\\ \dot{m}(t)&=&\sigma V(t)-\mu _{m}m(t)\ \ \mbox{dla}\ \ m(t)\leq 1.\end{array}$

(13.1)

W większości przypadków istnieje pewna progowa wielkość $m^{*}$ , do której $m$ nie wpływa znacząco na reakcję odpornościową i można założyć, że $\xi(m)=1$ dla $m<m^{*}$ .

$Przebieg typowej funkcji $\xi(m)$, czyli funkcja tożsamościowo równa 1 na odcinku $[0,m^{*}]$ oraz malejąca liniowo do 0 na $[m^{*},1]$.$

Rys. 13.1. Przebieg typowej funkcji $\xi(m)$ , czyli funkcja tożsamościowo równa 1 na odcinku $[0,m^{*}]$ oraz malejąca liniowo do 0 na $[m^{*},1]$ .

Typowo w analizie matematycznej układu (13.1) zakładamy, że $m<m^{*}$ i wtedy możemy zredukować układ do 3 rrzoa, choć dla konkretnej postaci funkcji $\xi$ można też badać np. stabilność całego układu (13.1).

13.2. Podstawowe własności układu (13.1)

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 $[-\tau,0]$ o wartościach w $(\mathbb{R}^{+})^{3}\times[0,1]$ . Ponieważ jednak tylko współrzędne $V(t)$ i $F(t)$ występują z opóźnionym argumentem, więc możemy zadać tylko ich wartości na całym przedziale $[-\tau,0]$ , ograniczając się do wartości $C_{0}$ i $m_{0}$ dla zmiennych $C$ i $m$ . 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.

Istnienie lokalnego rozwiązania układu (13.1).

Niech $V_{0}(s)$ , $F_{0}(s):[-\tau,0]\to\mathbb{R}^{+}$ . Wtedy dla $t\in[0,\tau]$ układ (13.1) sprowadza się do niejednorodnego układu rrz postaci

$\begin{array}[]{l c l}\dot{V}(t)&=&(\beta-\gamma F(t))V(t),\\ \dot{C}(t)&=&\alpha\xi(m(t))G_{0}(t)-\mu _{C}(C(t)-C^{{*}}),\\ \dot{F}(t)&=&\varrho C(t)-(\mu _{F}+\eta\gamma V(t))F(t),\\ \dot{m}(t)&=&\sigma V(t)-\mu _{m}m(t)\ \ \mbox{dla}\ \ m(t)\leq 1,\end{array}$

(13.2)

gdzie $G_{0}(t)=V_{0}(t-\tau)F_{0}(t-\tau)$ jest znaną funkcją, gdyż dla $t\in[0,\tau]$ argument $t-\tau\in[-\tau,0]$ . Funkcje opisujące współrzędne prawej strony układu (13.2) są wielomianami kwadratowym względem współrzędnych, czyli są klasy $\mathbf{C}^{1}$ (faktycznie $\mathbf{C}^{{\infty}}$ ), więc spełniają założenia twierdzeń o lokalnym istnieniu i jednoznaczności rozwiązań równań różniczkowych zwyczajnych. Istnieje zatem odcinek $[0,\bar{t})\subset[0,\tau]$ , na którym mamy jednoznaczne rozwiązanie układu.

Nieujemność lokalnego rozwiązania układu (13.1).

Pokażemy najpierw, że rozwiązanie to jest nieujemne. Pierwsze równanie układu (13.2) zapisujemy w równoważnej postaci całkowej

$\frac{\dot{V}}{V}=\beta-\gamma F\ \Longrightarrow\ \ln\left|\frac{V}{V_{0}}\right|=\int\limits _{0}^{t}(\beta-\gamma F(s))ds\ \Longrightarrow\ V(t)=V_{0}\text{e}^{{\beta t-\int\limits _{0}^{t}\gamma F(s)ds}},$

gdzie w powyższym wzorze $V_{0}=V_{0}(0)$ . Zatem $V(t)\geq 0$ dla $V_{0}\geq 0$ , co więcej, $V(t)>0$ dla $V_{0}>0$ , 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 $\xi(m)=1$ drugie równanie układu (13.2) jest równaniem liniowym niejednorodnym. Rozwiązujemy je stosując metodę uzmienniania stałej i otrzymujemy

$C(t)=C^{*}+(C_{0}-C^{*})\text{e}^{{-\mu _{C}t}}+\int\limits _{0}^{t}G_{0}(s)\text{e}^{{-\mu _{C}(t-s)}}ds.$

(13.3)

Ponieważ $G(s)\geq 0$ oraz $C^{*}+(C_{0}-C^{*})\text{e}^{{-\mu _{C}t}}=C^{*}\left(1-\text{e}^{{-\mu _{C}t}}\right)+C_{0}\text{e}^{{-\mu _{C}t}}>0$ dla $t>0$ , to $C(t)>0$ dla $t\in[0,\tau]$ . Zauważmy też, że w tym przypadku współrzędna $C$ jest zdefiniowana za pomocą wzoru (13.3) na całym przedziale $[0,\tau]$ .

W ogólnym przypadku, gdy $\xi(m)<1$ , oszacujemy wartość pochodnej $\dot{C}$ z dołu przez $-\mu _{C}C(t)$ , więc

$\dot{C}\geq-\mu _{C}C\ \Longrightarrow\ C(t)\geq C_{0}\text{e}^{{-\mu _{C}t}}.$

Zatem $C(t)\geq 0$ dla $C_{0}\geq 0$ . Zauważmy także, że jeśli $C_{0}=C^{*}$ , to z powyższych nierówności wynika, że $C(t)>C^{*}$ dla dowolnego $t>0$ .

W analogiczny sposób dowodzimy nieujemności zmiennej $m(t)$ .

Skoro $C(t)\geq 0$ dla $t\in[0,\bar{t})$ , to podobnie oszacujemy $\dot{F}$ , czyli

$F(t)\geq F_{0}\text{e}^{{-\int\limits _{0}^{t}(\eta\gamma V(s)+\mu _{F})ds}}.$

Zatem $F(t)\geq 0$ dla $F_{0}\geq 0$ .

Istnienie rozwiązania układu (13.1) na przedziale $[0,\tau]$ .

Pokażemy teraz, że z nieujemności rozwiązań wynika ich istnienie na całym przedziale $[0,\tau]$ . W tym celu oszacujemy nasze rozwiązania z góry. Dla pierwszej zmiennej mamy

$\dot{V}\leq\beta V\ \Longrightarrow\ V(t)\leq V_{0}\text{e}^{{\beta t}},$

zatem wzrost $V(t)$ jest co najwyżej wykładniczy.

Dla zmiennej $C(t)$ mamy

$\dot{C}\leq\alpha G_{0}\ \Longrightarrow\ C(t)\leq C_{0}+\alpha\int\limits _{0}^{t}G_{0}(s)ds,$

czyli $C(t)$ jest również ograniczona na każdym podprzedziale odcinka $[0,\tau]$ .

Z ograniczoności $C(t)$ i $V(t)$ wynika ograniczoność $F(t)$ i $m(t)$ . 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ę $\lim _{{t\to{\bar{t}}^{-}}}(V(t),C(t),F(t),m(t))$ (istnienie granicy wynika bezpośrednio z ograniczoności pochodnej) oraz punkt $\bar{t}$ i stosując twierdzenie o istnieniu rozwiązań. Wobec dowolności $\bar{t}$ wnioskujemy, że rozwiązanie istnieje na całym $[0,\tau]$ .

Globalne istnienie rozwiązań układu (13.1).

Stosując metodę kroków zakładamy, że istnieje nieujemne rozwiązanie układu (13.1) na odcinku $[(n-1)\tau,n\tau]$ i powtarzamy powyższe rozumowanie w celu wykazania istnienia nieujemnego rozwiązania na kolejnym odcinku $[n\tau,(n+1)\tau]$ . Na mocy zasady indukcji matematycznej jednoznaczne, nieujemne rozwiązanie układu (13.1) istnieje dla dowolnego $t\geq 0$ .

Punkty krytyczne układu (13.1).

Znajdziemy teraz punkty krytyczne badanego układu. Rozwiązujemy następujący układ równań algebraicznych

$\begin{array}[]{l c l}0&=&(\beta-\gamma\bar{F})\bar{V},\\ 0&=&\alpha\xi(\bar{m})\bar{V}\bar{F}-\mu _{C}(\bar{C}-C^{{*}}),\\ 0&=&\varrho\bar{C}-(\mu _{F}+\eta\gamma\bar{V})\bar{F},\\ 0&=&\sigma\bar{V}-\mu _{m}\bar{m}.\end{array}$

(13.4)

W ogólnym przypadku, aby rozwiązać układ (13.4), musimy zadać konkretną postać funkcji $\xi$ . Dla uproszczenia załóżmy więc, że $m<m^{*}$ , czyli $\xi(m)=1$ . 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 $V=0$ i $C=C^{*}$ . Wtedy $F=\tfrac{\varrho}{\mu _{F}}C^{*}=F^{*}$ i oczywiście $m=0$ . $F^{*}$ nazwiemy poziomem fizjologicznym przeciwciał.

Niech $X_{1}=(0,C^{*},F^{*})$ odzwierciedla stan zdrowia organizmu.

Jeśli $V\ne 0$ , to z pierwszego równania układu (13.4) mamy $\bar{F}=\frac{\beta}{\gamma}$ . Wstawiając tę zależność do drugiego i trzeciego równania obliczmy

$\bar{V}=\frac{\mu _{C}(\mu _{F}\beta-\gamma\varrho C^{*})}{\beta(\alpha\varrho-\eta\gamma\mu _{C})},\ \bar{C}=\frac{\alpha\mu _{F}\beta-\eta\gamma^{2}\mu _{C}C^{*}}{\gamma(\alpha\varrho-\eta\gamma\mu _{C})}.$

Punkt krytyczny $X_{2}=(\bar{V},\bar{F},\bar{C})$ odzwierciedla chorobę przewlekłą, oczywiście pod warunkiem, że wszystkie jego współrzędne są dodatnie. Współrzędne $\bar{V}$ i $\bar{C}$ możemy przepisać w nieco innej postaci

$\bar{V}=\frac{\mu _{C}\mu _{F}(\beta-\gamma F^{*})}{\beta(\alpha\varrho-\eta\gamma\mu _{C})},\ \bar{C}=\frac{\mu _{F}(\alpha\beta\varrho-\eta\gamma^{2}\mu _{C}F^{*})}{\gamma\varrho(\alpha\varrho-\eta\gamma\mu _{C})},$

z czego wynika, że jeśli

$\beta>\gamma F^{*}$ i $\alpha\varrho>\eta\gamma\mu _{C}$ albo
$\beta<\gamma F^{*}$ i $\alpha\varrho<\eta\gamma\mu _{C}$ ,

to $\bar{V}>0$ . Zauważmy dalej, że przy warunku A mamy:

$\alpha\beta\varrho-\eta\gamma^{2}\mu _{C}F^{*}>\gamma F^{*}(\alpha\varrho-\eta\gamma\mu _{C})>0,$

(13.5)

zatem $\bar{C}>0$ . Przy warunku B nierówność (13.5) zmienia znak, ale znak zmienia też mianownik $\bar{C}$ , czyli $\bar{C}>0$ . Widzimy więc, że warunki A, B są konieczne i dostateczne dla istnienia dodatniego punktu krytycznego. Zauważmy także, że w tym przypadku $\bar{m}=\frac{\sigma}{\mu _{m}}\bar{V}>0$ oraz $\bar{V}<\frac{\mu _{m}}{\sigma}m^{*}$ . Ostatnia nierówność oznacza, że w omawianym przypadku $m<m^{*}$ poziom stacjonarny antygenu $\bar{V}$ nie może być zbyt duży. W literaturze mówi się o lekkiej chorobie przewlekłej dla $m<m^{*}$ , natomiast gdy $m>m^{*}$ — o ciężkiej.

Lokalna stabilność punktów krytycznych układu (13.1).

Badanie lokalnej stabilności punktów krytycznych rrzoa nie różni się w początkowej fazie od analogicznej analizy dla rrz. Niech $(\tilde{V},\tilde{C},\tilde{F})$ oznacza punkt krytyczny. Zdefiniujmy nowe zmienne: $v=V-\tilde{V}$ , $c=C-\tilde{C}$ , $f=F-\tilde{F}$ . Układ w nowych zmiennych ma postać

$\begin{array}[]{l c l}\dot{v}(t)&=&(\beta-\gamma\tilde{F}-\gamma f(t))v(t)-\gamma\tilde{V}f(t),\\ \dot{c}(t)&=&\alpha(\tilde{V}f(t-\tau)+\tilde{F}v(t-\tau)+v(t-\tau)f(t-\tau))-\mu _{C}c(t),\\ \dot{f}(t)&=&\varrho c(t)-\mu _{F}f(t)-\eta\gamma(\tilde{F}v(t)+\tilde{V}f(t)+v(t)f(t)).\\ \end{array}$

(13.6)

Linearyzujemy układ (13.6) wokół punktu $(0,0,0)$ zakładając, że $|v(t)|<\varepsilon$ , $|c(t)|<\varepsilon$ , $|f(t)|<\varepsilon$ dla $t>0$ . Składniki wyższego rzędu, czyli $v(t)f(t)$ i $v(t-\tau)f(t-\tau)$ szacują się przez $\varepsilon^{2}$ , co możemy pominąć. Po zlinearyzowaniu układ (13.6) przyjmuje postać

$\begin{array}[]{l c l}\dot{v}(t)&=&(\beta-\gamma\tilde{F})v(t)-\gamma\tilde{V}f(t),\\ \dot{c}(t)&=&\alpha(\tilde{V}f(t-\tau)+\tilde{F}v(t-\tau))-\mu _{C}c(t),\\ \dot{f}(t)&=&\varrho c(t)-\mu _{F}f(t)-\eta\gamma(\tilde{F}v(t)+\tilde{V}f(t)).\\ \end{array}$

(13.7)

W kolejnym kroku znajdujemy rozwiązania układu (13.7) postaci $(v(t),c(t),f(t))=(v_{0},c_{0},f_{0})\text{e}^{{\lambda t}}$ , $\lambda\in\mathbb{C}$ , $\lambda\ne 0$ . Otrzymujemy więc następujące zależności

$\begin{array}[]{l c l}0&=&(-\lambda+\beta-\gamma\tilde{F})v_{0}-\gamma\tilde{V}f_{0},\\ 0&=&\alpha\tilde{F}v_{0}\text{e}^{{-\lambda\tau}}+\alpha\tilde{V}f_{0}\text{e}^{{-\lambda\tau}}-(\mu _{C}+\lambda)c_{0},\\ 0&=&-\eta\gamma\tilde{F}v_{0}+\varrho c_{0}-(\eta\gamma\tilde{V}+\mu _{F}+\lambda)f_{0}.\\ \end{array}$

(13.8)

Układ (13.8) ma nietrywialne rozwiązania, jeśli

$\det\left(\begin{array}[]{c c c}\beta-\gamma\tilde{F}-\lambda&0&-\gamma\tilde{V}\\ \alpha\tilde{F}\text{e}^{{-\lambda\tau}}&-\mu _{C}-\lambda&\alpha\tilde{V}\text{e}^{{-\lambda\tau}}\\ -\eta\gamma\tilde{F}&\varrho&-\eta\gamma\tilde{V}-\mu _{F}-\lambda\\ \end{array}\right)=0.$

(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 $(0,0,0)$ 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 $X_{1}$ . Wtedy macierz układu ma postać dolnietrójkątną, zatem łatwo wyznaczamy wartości własne, gdyż spełniają równanie

$(\beta-\gamma\tilde{F}-\lambda)(\mu _{C}+\lambda)(\eta\gamma\bar{V}+\mu _{F}+\lambda)=0,$

czyli wartości własne są równe: $\lambda _{1}=\beta-\gamma F^{*}$ , $\lambda _{2}=-\mu _{C}$ , $\lambda _{3}=-(\eta\gamma\bar{V}+\mu _{F})$ , zatem są rzeczywiste i ujemne dla $\beta<\gamma F^{*}$ .

Wniosek 13.1

Jeśli $\beta<\gamma F^{*}$ , to rozwiązanie $X_{1}$ jest lokalnie asymptotycznie stabilne. Natomiast jeśli $\beta>\gamma F^{*}$ , to $X_{1}$ jest niestabilne.

Widzimy, że stabilność rozwiązania $X_{1}$ , czyli wyzdrowienie jest możliwe tylko przy niewielkiej reproduktywność antygenu $\beta$ w stosunku do poziomu fizjologicznego przeciwciał $F^{*}$ .

Badając stabilność rozwiązania $X_{2}$ otrzymujemy następujący pseudowielomian charakterystyczny

$W(\lambda)=\lambda^{3}+a\lambda^{2}+b\lambda-d+(f-g\lambda)\text{e}^{{-\lambda\tau}},$

(13.10)

gdzie

$\begin{array}[]{c}a=\mu _{C}+\mu _{F}+\eta\gamma\bar{V},\ b=\mu _{C}(\mu _{F}+\eta\gamma\bar{V})-\eta\beta\gamma\bar{V}\\ d=\eta\gamma\mu _{C}\beta\bar{V},\ g=\alpha\varrho\bar{V},\ f=\beta g,\end{array}$

(13.11)

a wyznaczenie współczynników pozostawiamy jako ćwiczenie. Przypomnijmy, że dla $\tau>0$ 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 $\tau=0$ , to z ciągłej zależności od parametrów wynika, że albo rozwiązanie jest stabilne dla wszystkich $\tau>0$ , albo istnieje krytyczna wartość parametru $\tau _{c}$ , dla której może nastąpić destabilizacja. Dla $\tau=0$ możemy zbadać stabilność punktu $X_{2}$ za pomocą kryterium Routha – Hurwitza, skąd dostajemy warunek stabilności

$0<f-d<a(b-g).$

(13.12)

Warunek $f-d>0$ jest równoważny $\alpha\varrho>\eta\gamma\mu _{C}$ , zatem stabilność $X_{2}$ dla $\tau=0$ możliwa jest tylko w przypadku A. Natomiast w przypadku B mamy stabilność punktu $X_{1}$ .

Zauważmy, że zmiana stabilności może nastąpić tylko gdy pseudowielomian charakterystyczny (13.10) ma pierwiastki czysto urojone $\lambda=i\omega$ . Wtedy

$-i\omega^{3}-a\omega^{2}+ib\omega-d=-(f-ig\omega)\exp(-i\omega\tau _{c})\Longrightarrow||-i\omega^{3}-a\omega^{2}+ib\omega-d||=||f-ig\omega||,$

czyli

$\omega^{6}+(a^{2}-2b)\omega^{4}+(b^{2}+2ad-g^{2})\omega^{2}+d^{2}-f^{2}=0.$

Zdefiniujmy funkcję pomocniczą

$\Phi(z)=z^{3}+(a^{2}-2b)z^{2}+(b^{2}+2ad-g^{2})z+d^{2}-f^{2},\ z=\omega^{2}.$

Jej dodatnie miejsca zerowe determinują potencjalne zmiany stabilności punktu krytycznego $X_{2}$ . Z warunku (13.12) wynika, że wyraz wolny wielomianu $\Phi$ jest ujemny, zatem istnieje co najmniej jedno dodatnie miejsce zerowe. Warunek A lub B implikuje też $a^{2}-2b>0$ , czyli dodatnie miejsce zerowe $\Phi$ jest wyznaczone jednoznacznie. Niech $\omega _{0}=\sqrt{z_{0}}$ , gdzie $\Phi(z_{0})=0$ . Podstawiając $\lambda=i\omega _{0}$ do pseudowielomianu (13.10) dostajemy układ równań, z którego wyznaczamy sinus i kosinus wartości krytycznej $\omega _{0}\tau _{c}$ . 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 $\left(\tau _{{c_{n}}}\right)_{{n=1}}^{{\infty}}$ z pierwszą wartością $\tau _{{c_{1}}}$ , taką że $\omega _{0}\tau _{{c_{1}}}\in(0,2\pi)$ i $\tau _{{c_{n}}}=\tau _{{c_{1}}}+2(n-1)\pi/\omega _{0}$ .

Zmiana stabilności następuje przy takiej krytycznej wartości $\tau _{c}$ , dla której istnieje para czysto urojonych wartości własnych $\pm i\omega$ i przy przekraczaniu osi urojonej wartości własne przechodzą z lewej półpłaszczyzny zespolonej na prawą, czyli $\tfrac{d}{d\tau}\Re W(\lambda)\Big|_{{\tau=\tau _{c}}}>0$ . Przypomnijmy, że znak tej pochodnej jest taki sam jak znak pochodnej funkcji pomocniczej $\Phi$ w odpowiednim miejscu zerowym — tutaj jedynym miejscu zerowym $z_{0}$ . W tej sytuacji mamy oczywiście dodatnią pochodną $\Phi^{{\prime}}(z_{0})$ i wartości własne przechodzą z lewej półpłaszczyzny zespolonej na prawą.

Wniosek 13.2

Jeśli spełnione są nierówności (13.12) i istnieje dodatni punkt krytyczny $X_{2}$ , to punkt ten jest stabilny dla $\tau=0$ i istnieje krytyczna wartość opóźnienia $\tau _{c}$ , przy której $X_{2}$ traci stabilność i pozostaje niestabilny dla $\tau>\tau _{c}$ . W $\tau=\tau _{c}$ zachodzi bifurkacja Hopfa, czyli pojawiają się nietrywialne rozwiązania okresowe o okresie $\tfrac{2\pi}{\omega _{0}}$ .

Przypomnijmy, że z nierówności $f-d>0$ wynika, że stabilność punktu $X_{2}$ przy $\tau=0$ jest możliwa tylko w przypadku spełnienia warunku A.

$Hodograf Michajłowa dla parametrów: $a=f=1$, $d=g=\tau=0.5$, $b=2$.$

Rys. 13.2. Hodograf Michajłowa dla parametrów: $a=f=1$ , $d=g=\tau=0.5$ , $b=2$ .

Korzystając z kryterium Michajłowa, można udowodnione następujące twierdzenie dotyczące stabilności punktu $X_{2}$ , które daje nam oszacowanie wielkości opóźnienia krytycznego $\tau _{c}$

Twierdzenie 13.1

Niech $\mu _{C}\tau\leq 1$ oraz

$0<\frac{f-d}{a-g\tau}<b-g-f\tau.$

(13.13)

Jeśli dodatkowo spełnione są nierówności (13.12) i istnieje dodatni punkt krytyczny $X_{2}$ , to punkt ten jest lokalnie asymptotycznie stabilny.

Twierdzenie 13.1 zadaje warunki wystarczające stabilności punktu $X_{2}$ 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 $X_{2}$ przy założeniu $m<m^{*}$ .

$Rozwiązania układu~\eqref{UklMar1} ilustrujące dwa podstawowe przebiegi choroby: zakażenie z wyzdrowieniem (po prawej) oraz chorobę przewlekłą (po lewej).$

Rys. 13.3. Rozwiązania układu (13.1) ilustrujące dwa podstawowe przebiegi choroby: zakażenie z wyzdrowieniem (po prawej) oraz chorobę przewlekłą (po lewej).

Symulacja ilustrująca sytuację, w której poprzez zwiększenie dawki antygenu możemy wyprowadzić organizm z choroby przewlekłej.

Rys. 13.4. Symulacja ilustrująca sytuację, w której poprzez zwiększenie dawki antygenu możemy wyprowadzić organizm z choroby przewlekłej.

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.