Od strony matematycznej prezentowane mo- dele należą do klasy ciągłych układów dynamicznych (patrz np

M A T E M A T Y K A S T O SO W A N A 1, 2 0 0 0

U^{r sz u l a} F^{o r y ś} (Warszawa)

M od ele m a te m a ty cz n e w epidem iologii i im m unologii

W stęp. Artykuł niniejszy jest poświęcony wybranym zagadnieniom mo- delowania infekcji na poziomie całej populacji oraz na poziomie pojedyn- czego osobnika w tej populacji. Od strony matematycznej prezentowane mo- dele należą do klasy ciągłych układów dynamicznych (patrz np. [Sk]). Jest to tylko jeden z możliwych sposobów podejścia do tych zagadnień, ale mimo swej klasyczności nie tracący na aktualności i pozwalający w jasny sposób opisać w języku symboli matematycznych procesy zachodzące w trakcie mo- delowanych zjawisk.

Na początku zajmiemy się podstawowymi modelami stosowanymi w epi- demiologii (por. [Co], [Ii], [KM], [HT], [My]). Będziemy chcieli przewidzieć, jaki może być efekt epidemii - czy liczba osób, które chorują, będzie ma-

lała, czy też ustali się na pewnym poziomie. Zaprezentujemy dwa proste modele, jeden bez uwzględnienia efektu nabywania odporności na daną in- fekcję i drugi z uwzględnieniem tego efektu (model Kermacka-McKendricka [KM]). Przedstawimy dokładną analizę asymptotyki obu modeli w celu zi- lustrowania używanych metod matematycznych.

Opisane modele są często uogólniane na dwa sposoby. W pierwszym z nich uwzględnia się wpływ historii. Układ równań zwyczajnych zamie- nia się wtedy na układ z opóźnionym argumentem (mogą to być modele z opóźnieniem dyskretnym, jak w przypadku omawianego w tym artykule modelu Marczuka [Mk], bądź też modele różniczkowo-całkowe, w których uwzględnia się wpływ całej historii na chwilę obecną, por. [HT]). Drugi spo- sób uogólnienia to modele z uwzględnieniem rozmieszczenia przestrzennego - powstaje wtedy układ równań cząstkowych (można np. uwzględnić dyfuzję, jak w [MM] czy ogólniej w [Bsl]). Podkreślić należy, że oba typy uogólnień prowadzą do układów dynamicznych nieskończenie wymiarowych, w odróż- nieniu od klasy modeli typu Kermacka-McKendricka, które są układami skończenie wymiarowymi. Można także dokonywać uogólnień na bazie teo- rii równań różniczkowo-stochastycznych, w których oprócz elementów deter-

[35]

ministycznych wprowadza się czynniki losowe. Przykłady tego typu modeli można znaleźć np. w [Gd], [BG], [Py].

Następnie przejdziemy do objaśnienia najprostszych modeli immunolo- gicznych (patrz np. [Mk]), w których uwzględnia się tylko dynamikę anty- genu (czyli czynnika wywołującego infekcję), a następnie dołącza się interak- cje z przeciwciałem (elementem systemu immunologicznego, który niszczy antygen).

Na koniec omówimy konstrukcję i własności modelu Marczuka ([Mkl, 2, 3], [Bh], [BF1, 2], [BS], [Fś4, 5, 7], [SY]), jednego z prostszych modeli mate- matycznych w immunologii, w którym uwzględnia się wszystkie podstawowe procesy tzw. humoralnej odpowiedzi odpornościowej (tzn. reakcji odporno- ściowej na antygen swobodny, nie związany z komórkami organizmu).

Model Marczuka był także na wiele sposobów uogólniany (porównaj np.

[Mkl, 2, 3], [Bh], [Bs2], [BF2, 3], [BSk], [Fśl, 2, 3, 6], [FŻ1, 2]), powsta- wały także modele o podobnej strukturze matematycznej, lecz uwzględnia- jące zjawisko ograniczonej możliwości reakcji systemu immunologicznego, co przejawia się użyciem w modelu tzw. funkcji przełączenia, znanych także w literaturze biomatematycznej jako czynnik typu Michaelisa-Mentena (np.

[PW], [KP], [KMn]). W celu zobrazowania tego typu modelowania przedsta- wimy przykładowe modele - zaproponowany w [KMn] układ równań opisu- jący reakcję systemu immunologicznego na rozwój komórek nowotworowych oraz podobny model [KP], w którym dodatkowo uwzględniona została rola interleukiny I L 2 w odpowiedzi odpornościowej skierowanej przeciw nowo- tworowi.

Warto podkreślić, że w modelu Marczuka po raz pierwszy uwzględniony został bardzo istotny element każdej odpowiedzi odpornościowej - opóźnie- nie reakcji w stosunku do momentu wtargnięcia antygenu do organizmu.

W związku z tym, mimo dość prostej postaci modelu, jego analiza matema- tyczna jest znacznie bardziej skomplikowana niż w przypadku modeli bez opóźnienia (porównaj np. [DG], [Gy], [He], [Kg], [KN]).

Modele epidemiologiczne. Najprostszy model epidemiologiczny za- kłada istnienie dwóch klas osobników: klasy osobników zdrowych podatnych na infekcję (I) i klasy osobników zainfekowanych (II). Każdy osobnik z klasy I jest podatny na infekcję, zatem może stać się osobnikiem z klasy II. Po przej- ściu infekcji osobnik wraca do klasy zdrowych lub umiera. Wobec tego model ten nie uwzględnia efektu nabycia odporności. Modele takie znane są pod nazwą SIS, od angielskiego susceptible - podatny i infected - zainfekowany.

Zakładamy dla uproszczenia, że liczebność populacji jest stała i wynosi N (śmiertelność jest równoważona w tym modelu przez rozrodczość). Niech S(t), I (t ) oznaczają gęstości osobników klasy I i II odpowiednio. Zatem N S (t) jest liczbą osobników zdrowych, a N I(t) liczbą osobników zainfeko-

wanych. Zakładamy także, że osobniki, które się rodzą, są zdrowe. Infekcja rozprzestrzenia się przez kontakt między osobnikiem zainfekowanym i zdro- wym. Liczba takich kontaktów jest proporcjonalna do N I(t)S(t). Liczba zachorowań jest zatem równa AN S(t)I(t), gdzie A oznacza średni? liczbę kontaktów powodujących infekcję. Po przebyciu choroby osobniki zdrowieją.

Współczynnik wyzdrowień jest w tym modelu stały i równy 7 . Niech n ozna- cza współczynnik urodzin i śmierci (z założenia współczynniki te są sobie równe). Wtedy 7 + fi jest frakcją osobników, którzy opuszczają klasę II (po- przez wyzdrowienie lub śmierć). Ostatecznie model będzie miał postać

( i )

gdzie:

$ (NS(t)) = nN - \NS(t)I(t) + 7JV/(t) - nNS(t),

£(NI(t)) = ANS(t)I(t) - yNI(t) -

• fiN = const jest liczbą nowo narodzonych osobników w dowolnej chwili,

• XN S(t)I(t) jest liczbą osobników opuszczających klasę I w wyniku zainfekowania i w związku z tym zasilających klasę II w chwili t,

• 7 -/V7 (t) jest liczbą osobników zdrowiejących, zatem opuszczających klasę II i zasilających klasę I w chwili t,

• fiN S(t) jest liczbą osobników zdrowych umierających w chwili t,

• fiN I(t) jest liczbą osobników chorych umierających w chwili t.

Po podzieleniu przez N dostaniemy

(2) f S = /i - ASI + l i - nS,

K } 1 1 = XSI - 7/ - fil.

Z założeń modelu wynika, że S = 1 — I. Zauważmy, że

a zatem S + I — 1 jest punktem stacjonarnym, czyli podstawowe założe- nie jest spełnione. Wobec tego układ (2) możemy zredukować do jednego równania różniczkowego postaci

(3) i = —XI2 + (A — (7 + aO)/-

Przy założeniu I / 0 przeprowadzimy zamianę zmiennych y — j . Wtedy y = —jz l i równanie (3) przyjmuje postać

(4) y = X - (A - (7 + y))y.

Równanie typu (4) jest znane pod nazwą równania Bernoulliego. Wpro- wadźmy współczynnik

(5 = —- — . 1 + M Wtedy (4) daje się zapisać jako

(5) y = (7 + /i)(ń- ( 5 - l)y).

Rozwiążemy równanie (5). Jeśli 6 ^ 1, to 5 i stąd

(6⁾

y(t) = y0e“(7+,‘)('5- 1)< + ^ ( 1 - e-h+MW-i)*)

ó — 1

e ( 7 + / i ) ( < 5 - l ) t

m = ih + — 1) ’

a dla <5 = 1 mamy y(t) = yo + At, czyli (7)

gdzie /o 7^ 0 jest gęstością klasy osobników zainfekowanych w początko- wej chwili to — 0. Możemy teraz pozbyć się założenia I ^ 0, przepisując równania (6) i (7) w postaci

Joe(7+X)(*-l)<

1 + h yri (eh+Z'K'*-1)* - 1) oraz

odpowiednio, ponieważ dla Jq = 0 mamy /(t) = 0 dla dowolnego t > 0.

Będziemy wszakże rozważać tylko sytuację Iq / 0, gdyż przypadek I(t) = 0 jest mało interesujący.

St w ie r d z e n ie 1. Jeśli 6 < 1, to lim^+oo 7(t) = 0; natomiast jeśli 6 > 1, to lim^+oo I(t) = 1 - j.

Widzimy zatem, że w pierwszym przypadku epidemia wygasa, natomiast w drugim - ustala się na pewnym poziomie. Takie zjawisko nazywamy pan- demią. Kiedy może dojść do pandemii? Nasze założenia mówią, że może tak być, jeśli A > 7 + fi, a zatem wtedy, gdy frakcja osobników napływają- cych do klasy zainfekowanych (czyli A) jest większa niż frakcja osobników opuszczających tę klasę (7 + /x).

Założenia, które poczyniliśmy w zaprezentowanym modelu, są mało re- alistyczne. W rzeczywistości założenie o stałej wielkości populacji lub stałej liczbie kontaktów A nie będzie spełnione. Jeśli np. rozpatrzymy populację osobników w wieku szkolnym w dużym mieście, to liczba kontaktów w czasie trwania roku szkolnego będzie większa niż w czasie wakacji. Możemy sobie zatem wyobrazić, że zamiast stałej liczby kontaktów będziemy rozpatrywać funkcję okresową A(t), opisującą kontakty w różnych sezonach. Okazuje się, że dla takiego modelu można udowodnić podobne twierdzenie (sformułowane przez Hethcote’a) dla średniej wartości

7 + M’

gdzie A = ^ \q A(s) ds jest średnią wartością funkcji A(t) na przedziale o dłu- gości okresu, tzn. T. Jeśli S < 1, to następuje zanik infekcji, natomiast dla 5 > 1 epidemia nie wygasa, rozwiązanie dąży do pewnej granicznej funkcji okresowej.

Przedstawimy teraz model zaproponowany przez Kermacka i McKen- dricka [KM], w którym osobniki po przejściu infekcji uodporniają się na nią (lub umierają). Takie modele oznaczane są symbolem SIR (susceptible - in- fec ted - resistant, czyli podatny - zainfekowany - odporny). Symbolika tego typu nie jest w literaturze jednoznaczna, np. Hethcote [HT] używa określe- nia R - removed, co oznacza osobnika, który wyzdrowiał, ale nie musi być odporny (choć może).

Równania modelu Kermacka-McKendricka mają następującą postać:

( i(N S ( t) ) = -A N S(t)I{t),

(8) j = AN S (t)I(t) - ■yNĄt),

gdzie N R (t) oznacza liczbę osobników wyleczonych i odpornych na daną in- fekcję, a pozostałe oznaczenia nie ulegają zmianie, przy czym jak poprzednio zakładamy, że liczebność populacji jest stała, więc

N S(t) + N I(t) + N R(t) = N ,

ale tym razem upraszczamy model i nie uwzględniamy w nim rozrodczo- ści i śmiertelności. Poszczególne składniki prawych stron modelu (8) mają następujące znaczenie:

• AN S (t)I(t) jest liczbą osobników zdrowych podatnych na infekcję, za- silających w wyniku zainfekowania klasę osobników chorych w chwili i,

• 7 N I(t) jest liczbą osobników zdrowiejących i nabywających odporności na infekcję w chwili t.

Przy naszych założeniach można zatem zredukować układ (8) do dwóch równań postaci

f 5 = -A S I,

\ i = X S I - 7 /,

z warunkami początkowymi 5(0) = So > 0, J(0) = Io > 0, 5 + I E [0,1].

Zauważmy, że rozwiązanie układu (9) istnieje i jest wyznaczone jedno- znacznie dla dowolnego warunku początkowego (prawa strona jest wielomia- nowa).

St w ie r d z e n ie 2. Jeśli warunek początkowy jest nieujemny, to rozwią- zania układu (9) są także nieujemne.

Dow ód. Dla S ^ 0 i I ^ 0 możemy przepisać układ (9) w postaci

= -A J, j/ = AS — 7 i wycałkować stronami, skąd dostaniemy

S(t) = S0e“M»/w<is oraz I(t) = I 0e ^ s ^-~i'>ds,

z czego łatwo zauważyć, że jeśli tylko 5o > 0 i io > 0, to rozwiązania pozostają dodatnie dla wszystkich t > 0. Jeśli natomiast Io — 0, to I(t) = 0 i S(t) = So dla wszystkich t > 0. Dla So — 0 dostajemy S(t) = 0 oraz I =

—7J, czyli I(t) = ioe-7 *, co jest oczywiście nieujemne dla nieujemnego Io. ■ St w ie r d z e n ie 3. Rozwiązania układu (9) dla nieujemnych warunków początkowych są określone dla wszystkich t > 0.

Dow ód. Dla nieujemnych warunków początkowych rozwiązania są nie- ujemne, zatem S(t) jest funkcją malejącą dla tych t > 0, dla których ist- nieje. Wobec tego jest ograniczona przez So- Stąd I < (ASo — 7 )i, czyli I(t) < /oe^5°-7 ^, a zatem dla dowolnego i > 0, jeśli rozwiązanie istnieje w przedziale [0, £), to jest w tym przedziale ograniczone (jeśli A5o < 7 , to I(t) jest malejąca, czyli ograniczona przez Io, natomiast gdy XSq > 7 , to I(t) jest rosnąca, zatem ograniczona przez /oe^50-7^). Z postaci układu (9) wynika, że także pochodne funkcji S(t) i I(t) są ograniczone w [0, i) dla dowolnego t > 0, a zatem istnieją granice obu funkcji przy t —> t i wtedy gra- nice te można potraktować jako nowe wartości początkowe dla początkowego momentu czasu t. Wobec tego rozwiązania przedłużają się do nieskończono- ści. ■

St w ie r d z e n ie 4. Jeśli 0 < S0 + Jo < 1, to 0 < S(t) + I(t) < 1, dla każdego rozwiązania układu (9) i wszystkich t > 0.

Dowód. Wynika z analizy obrazu fazowego przedstawionego na rys. 1. ■ Zdefiniujmy teraz S = K Będziemy w dalszym ciągu zakładać, że waru- nek początkowy jest nieujemny. Zachodzi wtedy następujące twierdzenie:

Tw ie r d z e n ie 1. Jeśli 5So < 1, to funkcja I(t) jest malejąca i limf_*+oo I{t) = 0. Jeśli 5Sq > 1, to I(t) najpierw rośnie do wartości 1 — Ro — | — ln(J5o), a potem maleje i limf_>+00J(t) = 0. Funkcja S(t) jest malejąca i dąży do wartości 5oo, która jest rozwiązaniem równania

(11) l - i ? 0 - S oo + i l n ^ = 0

w przedziale (0, |), przy czym Ro jest wartością początkową funkcji R(t).

Dow ód. Znajdziemy równanie orbit układu (9) i przeanalizujemy jego obraz fazowy.

Przy założeniu 5/ / 0 (jeśli któraś ze zmiennych jest równa 0, to po- stępujemy jak w dowodzie stwierdzenia 2) dzielimy równania układu (9) stronami i dostajemy

<H_ _ 1 1 dS ~ 5 S czyli

I(S ) = i l n S - S + C,O gdzie C jest pewną stałą.

Stąd funkcja

7 + S - - A n S = C ó

jest całką pierwszą układu (9), przy czym

C = Iq + Sq — - In Sq = l ~ Ro — - In Sq.

ó o

Izokliny zerowe dla zmiennej S są wyznaczone przez równanie S I = 0, a dla zmiennej I przez równanie I(XS — 7 ) = 0. Wobec tego, jeśli I = 0, to każdy taki punkt jest punktem krytycznym, a izoklina 5 = 0 jest jednocze- śnie rozwiązaniem układu. Interesujący nas obszar zamyka się w trójkącie o wierzchołkach (0,0), (0,1), (1,0). Widzimy, że S jest funkcją stale malejącą wewnątrz trójkąta, natomiast I maleje dla S < | i rośnie dla S > Dla 5 = 1 funkcja 7(5 ) przyjmuje wartość maksymalną, jeśli jest w tym punk- cie określona. Punkty postaci (5, 0) są punktami krytycznymi. Kiedy taki punkt może być stabilny? Macierz naszego układu ma postać

(12) 0 -A 5

0 A5 — 7

zatem wartościami własnymi są 0 i A(5 — |), a więc o stabilności można mówić tylko, gdy Ś <

Nasze orbity pozostają w obszarze zwartym, wewnątrz którego nie ma punktów krytycznych. Zatem muszą dążyć do któregoś z punktów krytycz- nych, znajdujących się na brzegu obszaru. Stąd I(t) —* 0 przy t —> -foo i jeśli 5So < 1, to I(t) stale maleje. Jeśli 5 So > 1, to I(t) najpierw rośnie, potem maleje i osiąga wartość maksymalną dla I(S ) = /(|). Granica funkcji S(t) musi spełniać równanie (11). ■

M odele im m unologiczne. Przejdziemy teraz do omówienia modelowa- nia działania systemu immunologicznego (odpornościowego). Przedstawiony w tym rozdziale sposób podejścia do modelowania działania układu odpor- nościowego pochodzi z wykładu doc. dr hab. Jacka Waniewskiego z Instytutu Biocybernetyki i Inżynierii Biomedycznej PAN.

Na początku sformułujemy proste modele, które w rezultacie doprowadzą nas do głównego modelu - modelu Marczuka.

Substancję, która wywołuje reakcję odpornościową, nazwiemy antyge- nem. W zależności od typu antygenu, jego dynamika jako gatunku (lub kinetyka jako substancji) może być różna. Oznaczmy przez V(t) stężenie antygenu w chwili t. Jeśli mamy do czynienia z antygenem aktywnym, ta- kim jak bakterie, wirusy, pasożyty, komórki nowotworowe, to podstawą do budowania modelu będzie równanie wzrostu (równanie Malthusa)

r > 0;

jeśli antygen jest pasywny, np. trucizny, jady, szczepionki, to równanie wyj- ściowe ma postać

dV_

dt = ~8V, s > 0.

Prócz tego możemy mieć do czynienia z antygenami produkowanymi w orga- nizmie, jak autoantygeny i antygeny transplantacyjne, i wtedy podstawowe równanie może mieć postać

dV

dt = p - s V , P,s > 0,

przy czym w tym przypadku istnieje dodatni punkt krytyczny V =

Drugim zasadniczym elementem systemu immunologicznego jest prze- ciwciało, które ma za zadanie związanie antygenu i usunięcie go z organi- zmu. Stężenie przeciwciał oznaczymy przez F(t). Podobnie jak w klasycz- nym modelu drapieżca-ofiara zakłada się, że szybkość usuwania antygenów z organizmu jest proporcjonalna do iloczynu V (t)F(t). Oznacza to, że im większa jest liczba zarówno antygenów, jak i przeciwciał, tym większe jest prawdopodobieństwo spotkania pomiędzy nimi i usunięcia antygenu przez

przeciwciało z organizmu. W najprostszym modelu, gdy uwzględnimy tylko reakcje zachodzące pomiędzy antygenami i przeciwciałami, szybkość produk- cji przeciwciał może być np. proporcjonalna do ilości antygenu. Po wejściu w reakcję z antygenem przeciwciało jest także niszczone (usuwany j 3St cały kompleks antygen-przeciwciało). W takim przypadku dla antygenu aktyw- nego dostaniemy układ

('13'j [ V = r l V - c l VF,

V ’ \ P = r2V - c2VF.

Oznacza to, że zarówno antygenów, jak i przeciwciał przybywa proporcjonal- nie do liczebności tych pierwszych, a ubywa ich proporcjonalnie do iloczynu liczebności obydwu. W przypadku modelu (13) decydujący o zachowaniu rozwiązań jest znak wyrażenia riC2 — r2c\. Jeśli r\c2 — r2c\ < 0, to sys- tem immunologiczny jest zawsze w stanie poradzić sobie z antygenem i jego stężenie maleje do 0. Natomiast stężenie przeciwciał ustala się na pewnym poziomie (zależnym od warunków początkowych, gdyż punkty krytyczne w tym modelu nie są izolowane). Jeśli natomiast riC2 — r2Ci > 0, to sys- tem immunologiczny, w zależności od warunków początkowych, albo nie radzi sobie z infekcją i stężenie antygenu rośnie, z matematycznego punktu widzenia może rosnąć do nieskończoności, natomiast w praktyce po przekro- czeniu pewnego poziomu progowego stężenia antygenu następuje całkowite porażenie systemu immunologicznego (podobnie jeśli stężenie antygenu jest mniejsze niż pewna wartość progowa, to uznaje się, że antygen został wy- eliminowany z organizmu - stężeń poniżej pewnego poziomu nie jesteśmy w stanie doświadczalnie zaobserwować), albo zachowanie jest takie jak w po- przednim przypadku (jeśli początkowe stężenie przeciwciał jest dostatecznie duże, a początkowe stężenie antygenu odpowiednio małe).

Podobnie jak w przypadku modelu Kermacka-McKendricka rozwiązania układu (13) istnieją i są jednoznaczne dla dowolnego warunku początko- wego. Kolejnym krokiem w analizie takiego modelu jest analiza nieujemności i przedłużalności rozwiązań. Pokażemy zatem:

St w ie r d z e n ie 5. Jeśli Vq,Fo > 0, gdzie Vb, Fo oznaczają wartości funk- cji w punkcie to = 0, to V (t),F (t) > 0 i rozwiązania układu (13) istnieją dla wszystkich t > 0.

Dowód. Jeśli Vo = 0, to V(t) — 0 dla wszystkich t, dla których ist- nieje, gdyż każdy punkt postaci (0, F ) dla dowolnego F jest punktem kry- tycznym. Wtedy także F (t) — Fo, zatem jest nieujemne. Wobec tego in- teresujący jest tylko przypadek Vq > 0. Z jednoznaczności rozwiązań do- stajemy V(t) > 0 dla wszystkich t, dla których rozwiązanie istnieje. Za- uważmy, że aby F(t) < 0 dla pewnego t > 0, musi istnieć pierwszy taki mo- ment to, w którym funkcja F(t) przekracza oś poziomą, a zatem F(to) = 0 i H to) < 0. Ale z drugiej strony F(io) = riV (to) > 0, co jest sprzeczne

z definicją punktu to. Zatem rozwiązania są nieujemne dla wszystkich £, dla których istnieją. Pokażemy teraz, że istnieją dla wszystkich t > 0. Za- uważmy, że jeśli funkcje są nieujemne, to V < r\V, czyli V(t) < Voerit, więc jest ograniczone w dowolnym skończonym przedziale czasu. Stąd wy- nika, że F < r2Voerit, co implikuje F(t) < FoVo^erit. Zatem rozwiązania układu (13) są ograniczone, ich pochodne także, czyli dowolne rozwiązanie przedłuża się do nieskończoności. ■

Dla modelu (13) można policzyć równania orbit i następnie wyznaczyć równanie, które musi spełniać granica rozwiązania. Po podzieleniu stronami równań układu (13) dostajemy

dV d F

c i a i - F

C2 <22 — F ’ gdzie aj = —, j = 1, 2,

Ci

dla F 7^ a2, a po scałkowaniu z warunkiem początkowym (Vq, Fq) (14)

(15)

V = V0 + - (F - F0 + (ai — a,) In --- i-_{c2 V} F0 - 0,2 / V = V ó + — { F - F o ) d l a o i = a 2.

C2

/ —& 7

Tw ie r d z e n ie 2. Jeśli Vo,Fo > O, to rozwiązania układu (13) mają na- stępującą asymptotykę:

• dla a\ = a2 każde rozwiązanie dąży do jednego z punktów krytycznych,

• dla oi < o,2 każde rozwiązanie dąży do jednego z punktów krytycznych postaci (O, F ),

• dla ai > 02, jeśli warunek początkowy spełnia nierówność Vo > — (f 0 - oi + (o2 - ai) ln —— — ' j ,

c2 V oi - o2 /

to każde rozwiązanie ma własność lim^+oo V (t) = +oo i lim^+oo F(t)

= 02, a w przeciwnym razie dąży do jednego z punktów krytycznych postaci (O,F ).

Dow ód. Niech o = ai = 02. Wtedy orbity układu (rys. 2) leżą na liniach prostych o równaniu (15). Każdy punkt postaci (O, F ), (K,o) jest punktem krytycznym. Jeśli F > o, to obie zmienne maleją, jeśli F < a, to obie rosną.

W obu przypadkach orbity są ograniczone (w pierwszym przypadku przez (Ko, Fo), w drugim przez (Vfo + ^ (o -F o ), a)), gdyż z jednoznaczności rozwią- zań F nie może przekroczyć wartości a. Wobec tego dla wartości (Vo, Fo) powyżej prostej V = ^ ( F — o) rozwiązania dążą do punktu krytycznego (O, Fo — §-Ko), natomiast poniżej tej prostej do punktu (Ko + §-(a — Fo), a).

Niech ai < 02- Wszystkie orbity (rys. 3) oprócz F = a2 leżą na krzywych o równaniu (14). Jeśli F0 > 02, to ^ > O, zatem K (F) jest funkcją rosnącą.

Natomiast obie zmienne maleją w tym obszarze względem czasu. Mamy

więc orbity monotoniczne w obszarze ograniczonym - zatem dążą do punktu krytycznego. Z jednoznaczności rozwiązań wynika, że lim^+oo V(t) = o, i wtedy granica lim^+oo F (t) musi spełniać równanie

(16) 0 = Vo + — ( lim F(t) - F0 + (a2 - a{) ln lim^ + °° ^ .

C2 V t—+~ł-oo — a 2 /

Jeśli ai < Fq < a2, to ^ < 0, czyli V (F) jest funkcją malejącą, natomiast w zależności od czasu V(t) maleje, a F(t) rośnie, przy czym F jest ogra- niczona przez a2. Zatem znów lim^+oo V(t) = 0 i limt_^+00 F(t) spełnia równanie (16).

Rys. 2. Obraz fazowy układu (13) w przypadku oi = a2 = a

Rys. 3. Obraz fazowy układu (13) w przypadku ai < a2

Jeśli Fo < ai, to V (F ), F(£), V(t) są rosnące i ograniczone w tym obszarze przez (V& + g (a i - F0 + (o2 - ai) In fF=SJ)» oi)- Gdy F osiąga wartość ai, to wchodzimy do omówionego wcześniej obszaru i V(t) zaczyna maleć.

Rys. 4. Obraz fazowy układu (13) w przypadku ai > <Z2

Dla a\ > <22 zmienia się dynamika (rys. 4) układu (13). Jeśli F = <22, to V = c\V{a\ — <22) > 0, czyli V (t) rośnie eksponencjalnie do nieskończoności.

W przypadku Fo < a2 obie zmienne rosną w czasie, V (t) jest nieograniczona, zatem lim*_>+0o V(t) = + 00.

Zauważmy, że każda z krzywych całkowych (14) przechodzi albo przez (V", 0), albo przez (0 ,F ), F < (Z2, co można przyjąć jako warunek począt- kowy. Mamy zatem

V = V + - (f + (o2- Ol) ln lub

V = — (f — F + (0.2 - <21) ln Z— .

C2 V F - a 2J

Graniczną krzywą jest

V = — (f + (a2 - ai) ln —— .

C2 V 02 /

Zauważmy także, że aby V —> +00 przy ograniczonym F, musimy mieć F —* a 2.

Jeśli Fo G (a2,ai], to V (t) w dalszym ciągu rośnie, ale F(t) maleje, czyli znów lim ^j-oo V'(t) = -foo i limf_»+00 F (t) = a2. Jeśli natomiast Fo > a\, to albo F(t) > a\ dla wszystkich t > 0 i wtedy obie zmienne maleją, za- tem mamy zbieżność do jednego z punktów krytycznych (0, F ), albo istnieje

to > O takie, że F(to) = oi, i wtedy V zaczyna rosnąć - przechodzimy do obszaru omówionego wcześniej. Całkowa krzywa graniczna przechodzi przez punkt (0, ai). ■

W przypadku antygenu pasywnego stężenie V (t) zawsze maleje. Podob- nie możemy rozpatrywać inne warianty układu antygen-przeciwciało, ale wiadomo, że równie istotnym elementem systemu immunologicznego są ko- mórki produkujące przeciwciała, jak również komórki stymulujące i pomoc- nicze, nazywane limfocytami.

Omówimy teraz (według [Mkl, 2, 3], więcej szczegółów można znaleźć np.

w [Mz]) podstawy działania systemu immunologicznego i rolę jego poszcze- gólnych składników. Większość komórek odpornościowych powstaje z jed- nego substratu S produkowanego w szpiku kostnym. W zależności od drogi, którą później substrat ten przechodzi, przekształca się w komórki spełnia- jące różne funkcje. Jedną z grup komórek powstających z substratu S są T-limfocyty, wśród których wyróżniamy trzy zasadnicze podgrupy: T e - limfocyty efektywne, które neutralizują i rozkładają antygeny, Tp - limfo- cyty pomocnicze współdziałające w rozpoznaniu antygenu i T s - supresory, które regulują proces immunologiczny (zmniejszenie ich poziomu może spo- wodować reakcje alergiczne, czyli potraktowanie przez organizm substancji nieszkodliwej jak antygen). Drugą grupą komórek powstających z substratu S są B-limfocyty. Po przekształceniu B-limfocyty są zdolne do tworzenia komórek plazmatycznych, które produkują przeciwciała (immunoglobuliny różnych typów, m.in. M, A, G, D, F, oznaczane IgM itd.). Antygen obecny w organizmie może mieć dwojaką postać antygenu swobodnego, krążą- cego we krwi, bądź antygenu związanego, gdy zaatakuje komórki organizmu.

W pierwszym przypadku zachodzące reakcje określa się mianem humoralnej odpowiedzi odpornościowej, w drugim - odpowiedzi komórkowej. W przy- padku większości antygenów zachodzą oba te procesy. Czasami zachodzi tylko jeden lub można je rozseparować. Schemat obu reakcji jest dość po- dobny, tyle że biorą w nich udział inne czynniki systemu immunologicznego.

W odpowiedzi humoralnej najistotniejszą rolę odgrywają przeciwciała, na- tomiast w komórkowej głównym czynnikiem są limfocyty.

Opiszemy dla przykładu odpowiedź humoralną. Po wtargnięciu do orga- nizmu antygen krąży we krwi, dopóki nie zwiąże się z tzw. przeciwciałem pre- zentującym antygen. Następnie wiązanie to jest absorbowane przez limfocyt Tp, który się w pewnym sensie “przylepia” do niego. Powoduje to utworzenie tzw. kompleksu VT, który wraz z receptorem doczepia się do makrofagu L.

Całość doczepia się do B-limfocytu. Dany B-limfocyt może w zależności od swojej struktury przyłączyć różną liczbę takich kompleksów. Na tym eta- pie rozpoznawana jest struktura antygenu. Powstają komórki pamięci, które zapamiętują cechy antygenu, oraz komórki plazmatyczne, które wytwarzają

immunoglobuliny z odpowiednią strukturą do wiązania danego rodzaju an- tygenu. Antygeny związane z odpowiednimi przeciwciałami są wychwyty- wane przez limfocyty T e i transportowane do zniszczenia. Na początku, gdy organizm jest w fazie rozpoznawania antygenu, pojawiają się immunoglobu- liny M, wolniejsze, słabo przystosowane do rodzaju antygenu (tzw. nieswo- iste), następnie produkowane są szybkie IgG (swoiste). Część procesu zwana rozpoznaniem antygenu (do momentu wytworzenia kompleksu B-limfocytu z makrofagiem L) zachodzi tylko raz. Przy następnych zachorowaniach re- akcja immunologiczna jest dużo szybsza - organizm zna już typ antygenu (posiada komórki pamięci) i produkuje od razu immunoglobuliny swoiste (przeciw temu konkretnemu antygenowi). Z przedstawionego opisu reakcji wynika jednak jasno, że zawsze istnieje pewne opóźnienie procesu tworzenia przeciwciał. Aby powstały nowe przeciwciała, musi zajść szereg reakcji od momentu wniknięcia antygenu do organizmu.

Wyróżniamy trzy podstawowe typy przebiegu infekcji:

• forma subkliniczna - występuje najczęściej, gdy antygen jest już znany (np. odra, ospa wietrzna); porcja bakterii jest od razu absorbowana i je- śli dawka początkowa była niewielka, to taka infekcja jest praktycznie niezauważalna dla organizmu;

• forma ostra - zwykle antygen jest nieznany, szybko rozmnaża się po dostaniu do organizmu; albo organizm jest dostatecznie silny, aby so- bie poradzić, albo nie i wtedy może dojść do wyniku letalnego, czyli śmierci;

• forma chroniczna - bakterie cały czas egzystują w organizmie.

Model M arczuka. Przejdziemy teraz do konstrukcji modelu Marczuka (patrz [Mkl, 2, 3], [Bh]), który opisuje przebieg infekcji z udziałem anty- genu aktywnego, tzn. uproszczoną wersję odpowiedzi odpornościowej, bez uwzględnienia podziału na odpowiedź humoralną i komórkową. W modelu występują cztery zmienne:

• V(t) - gęstość antygenu w chwili t,

• F (t) - gęstość przeciwciał w chwili t,

• C (t) - gęstość komórek plazmatycznych w chwili t,

• m(t) - charakterystyka porażenia zainfekowanego organu.

Poszczególne zmienne interpretujemy w następujący sposób. V (t) jest stę- żeniem antygenu aktywnego. F (t) oznacza stężenie przeciwciał, ale można też interpretować tę zmienną znacznie szerzej, uwzględniając nie tylko prze- ciwciała w sensie biologicznym, ale także T-limfocyty. C(t) jest stężeniem komórek produkujących przeciwciała. Charakterystykę m(t) liczymy w na- stępujący sposób. Niech M będzie charakterystyką zdrowego organu (zwykle masa lub powierzchnia), natomiast M '{t) - charakterystyką zdrowej części

tego organu w chwili t. Wtedy m{t) — 1 — , a zatem m oznacza udział części zainfekowanej w całości organu. Z definicji m(t) G [0,1], przy czym jeśli m(t) = 1, to organ jest całkowicie porażony i uznajemy to za jedno- znaczne z letalnym wynikiem infekcji (śmiercią osoby zainfekowanej).

Dynamika antygenu jest opisywana takim samym równaniem jak w mo- delu (13), czyli antygen rozmnaża się proporcjonalnie do swej liczebności i ginie w reakcji z przeciwciałami, proporcjonalnie do iloczynu V (t)F(t).

Mamy zatem

Komórki plazmatyczne produkowane są po wtargnięciu antygenu do or- ganizmu i utworzeniu tzw. kompleksów VT (czyli połączeń między antyge- nami i T-limfocytami). Zakładamy, że liczba takich kompleksów jest pro- porcjonalna do iloczynu V (t)F(t), gdyż mogą one powstać dopiero po utwo- rzeniu wiązania antygen-przeciwciało prezentujące. Utworzenie kompleksów VT jest sygnałem do produkcji komórek plazmatycznych, co odbywa się z pewnym opóźnieniem, potrzebnym na powstanie tych wiązań. Faza ta na- zywa się rozpoznaniem antygenu. Produkcja komórek plazmatycznych jest w sposób automatyczny wyrównywana do poziomu normalnego dla danego organizmu, zwanego poziomem fizjologicznym i oznaczanego C*. Prócz tego produkcja tych komórek zależy od kondycji organizmu, a więc także od charakterystyki zainfekowanego organu. Im większe porażenie, tym jest ona mniejsza. Omówione procesy można opisać za pomocą równania

^ { t ) = a£ (m )F (t - r)V (t - r) - fic{C{t) - C*).

Rys. 5. Wykres typowej funkcji £

Funkcja £ jest odpowiedzialna za zmniejszenie produkcji komórek pla- zmatycznych spowodowane postępującym porażeniem zainfekowanego or- ganu. W ogólnym przypadku zakłada się, że jest to funkcja nierosnąca,

£(0) = 1, £(1) = 0. W modelu wyjściowym rozpatrywana była konkretna postać funkcji £ - funkcja kawałkami liniowa, na początku stale równa 1, a po przekroczeniu pewnej progowej wartości m* malejąca liniowo do 0, jak na rys. 5.

Jeśli nie ma antygenu, to równanie na zmianę stężenia komórek plazma- tycznych przyjmuje postać

ć ( t ) = - c ' ) , czyli

(17) C(t) = C* + (C0 - C ) e ~ ^ ł ,

a zatem limt_.+00 C(t) = C* i stężenie wraca do poziomu fizjologicznego.

Stężenie przeciwciał zależy od stężenia komórek plazmatycznych, przez które są produkowane. Przeciwciała giną śmiercią naturalną lub w reakcji z antygenami. Otrzymujemy zatem równanie

~ { t ) = gC(t) - ( n f + rnV (t))F(t).

Ostatnie równanie opisuje zmiany charakterystyki porażonego organu.

Porażenie zwiększa się wraz ze wzrostem stężenia antygenu i maleje dzięki naturalnym zdolnościom regeneracyjnym. Stąd

= °V (t) -

przy czym równanie to obowiązuje przy założeniu m(t) < 1. Jeśli w którymś momencie m (t) = 1, to zakłada się, że nastąpiło całkowite porażenie organu i letalny efekt infekcji.

Ostatecznie rozważamy następujący układ równań:

(V (t) = (f3-'yF (t))V (t),

I Ó(t) = aĘ(m )V (t - r)F (t - r) - ^ (C (t) - C*), ] H t) = eC (t) - (uf + rpfV (t))F(t),

[ m (t) = <jV(t) — nmm(t) dla m < 1,

z dodatnimi współczynnikami. Interpretacja biologiczna współczynników jest następująca:

• (3 - współczynnik reprodukcji antygenu,

• 7 - współczynnik wyrażający prawdopodobieństwo związania antygenu z przeciwciałem i siłę ich reakcji,

• a - współczynnik stymulacji systemu immunologicznego (opisujący efektywność reakcji),

(18)

• r - średnie opóźnienie procesu produkcji komórek plazmatycznych w stosunku do powstania wiązań antygen-przeciwciało i kompleksów VT,

• fic - współczynnik śmiertelności komórek plazmatycznych, którego od- wrotność równa jest średniej długości życia takiej komórki,

• g - współczynnik produkcji przeciwciał przez jedną komórkę plazma- tyczną,

• t] - liczba przeciwciał potrzebna do związania jednego antygenu,

• n f - współczynnik śmiertelności przeciwciał, którego odwrotność rów- na jest średniej długości życia przeciwciała,

• o - współczynnik wyrażający siłę porażenia organu przez antygen,

• Hm - współczynnik odnawiania się organu, którego odwrotność równa się średniemu czasowi regeneracji organu.

Warunek początkowy dla autonomicznych układów z opóźnieniem ma postać (0, X °), gdzie X ° jest wektorem funkcyjnym ze współrzędnymi okre- ślonymi na przedziale [—r, 0]. W przypadku układu (18) X ° =

F°,m P), ale z postaci układu wynika, że możemy ograniczyć się do znajo- mości funkcji V° i F° na przedziale [—r, 0] (bo tylko te funkcje są w modelu uwzględnione z opóźnionym argumentem), natomiast w przypadku współ- rzędnych C i m wystarczy mam znajomość wartości w punkcie to — 0.

W ogólnej teorii równań z opóźnieniem zakłada się, że współrzędne wektora X ° są ciągłe. Natomiast w przypadku układu (18) standardowo zakłada się, że warunek początkowy odpowiada sytuacji zdrowego organizmu zara- żonego w momencie czasu to = 0 pewną dawką antygenu, czyli V°(ti) = 0 dla h < 0 oraz V'°(0) = Vq > 0. Widzimy więc, że ten warunek początkowy odpowiada zadaniu warunku początkowego dla modelu bez opóźnienia. Pod- kreślić trzeba, że ogólna teoria równań z opóźnieniem jest rozwijana dla cią- głego warunku początkowego, ale warunek początkowy z funkcją skokową także daje się analizować w taki sam sposób.

Analizę modeli ze stałym opóźnieniem dyskretnym możemy przeprowa- dzać tzw. metodą kroków. Idea tej metody jest następująca. Niech omawiany model ma postać

(19) X (t) = f ( X ( t ) ,X ( t - r ) ) ,

gdzie X (t) € W1, n > 1 oraz / : R2n —> Rn jest funkcją ciągłą spełnia- jącą lokalnie warunek Lipschitza. Jeśli X ° jest ciągłą funkcją początkową określoną na [—r, 0], to dla t € [0, r] równanie (19) przyjmuje postać

(20) X (t) = f( X ( t ) ,X ° ( t - T ) ) ,

jest zatem równaniem zwyczajnym. Jeśli rozwiązanie równania (20) istnieje na całym przedziale [0,r], to możemy dalej postąpić analogicznie, przecho- dząc na przedział [r, 2r] itd. Ogólnie - jeśli wiemy, że rozwiązanie istnieje

na przedziale [kr, (k + l)r] dla dowolnej liczby naturalnej k, to możemy wyznaczyć rozwiązanie na następnym przedziale [(k + 1 )r, (k + 2)r].

Łatwo zatem zauważyć, że rozwiązania układu (18) istnieją i są wyzna- czone jednoznacznie, gdyż prawa strona jest funkcją ciągłą lokalnie speł- niającą warunek Lipschitza. Pokażemy, że rozwiązania są nieujemne dla nieujemnych warunków początkowych. Zastosujemy metodę z [Fśl], gdyż pierwotne dowody z [Mkl] oraz [Bh] posiadają pewne braki.

St w ie r d z e n ie 6. Jeśli współrzędne wektora początkowego ^Xq są nie- ujemne, to rozwiązanie układu (18) ma współrzędne nieujemne.

Dow ód. Zauważmy, ze pierwsze równanie układu (18) możemy zapisać w postaci całkowej

V(t) = VbeS°(/3“7jF(s))<is,

gdzie Vq = P°(0), i stąd jeśli tylko Vo > 0, to także V(t) > 0 dla t > 0 takich, dla których rozwiązanie istnieje.

Z drugiego równania mamy

C(t) = C* + (Cq - C*)e~^ct + ae~^ct \ t(m (s))V (s - t t)F(s - r )e ^ s ds o

i jeśli tylko funkcja F jest nieujemna, to C jest także nieujemna. Warunek początkowy jest nieujemny. Załóżmy zatem, że istnieje pierwszy taki punkt to > 0, w którym funkcja F staje się ujemna, czyli F(to) = 0 oraz F(to) < 0.

Oznacza to, że F(to) = gC(to) < 0. Jeśli C(to) > 0, to nierówności są sprzeczne. Ponieważ to jest pierwszym punktem o podanych własnościach, więc

(7(t0) > C* + (C0 - C *)e-/ic< = C*( 1 - e~»ct) + C0e ~ ^ > 0,

o ile Co > 0 lub C * ,p c > 0. Ale nawet gdyby te współczynniki mogły się zerować, to na odcinku [0, to) funkcja C musiała rosnąć, a zatem nie może być C(to) = 0. Wynika stąd, że obie funkcje C i F są nieujemne.

Z ostatniego równania wynika, że

m (t) = moe_/im* + ae~^mt ^V(s)e^mSds > 0,t o

wobec nieujemności warunku początkowego i funkcji V. •

UWAGA 1. Wystarczy oczywiście założyć nieujemność funkcji V° i F ° oraz wartości C°(0) i m°(0), natomiast dla standardowego warunku począt- kowego - tylko nieujemność wartości funkcji początkowych w to = 0.

Pokażemy teraz, że rozwiązania układu (18) są ograniczone na każdym odcinku skończonym [0,t*], t* > 0, a zatem przedłużają się do nieskończo- ności. Podobnie jak w przypadku stwierdzenia 6 dowód pochodzi z [Fśl].

St w ie r d z e n ie 7. Jeśli współrzędne początkowego wektora funkcyjnego X 0 są nieujemne, to rozwiązanie układu (18) jest ograniczone dla dowolnego skończonego momentu czasu t* > 0.

Dowód. Z pierwszego równania układu (18) wynika, że V < (3V, a zatem V (t) < Vbe^, czyli V(t) jest ograniczone dla każdego skończonego przedziału czasu. Podobnie m < a\oe^ i stąd m(t) < mo + ^ e ^ . Z drugiego równania modelu (18) dostajemy

Ć < aV (t - r )F (t - r) - fic(C - C*), więc

t

C{t) < C* + (C0 - C *)e-Mc* + ae-/Xct \ V{s - t)F {s - r )e ^ s ds.

o

Załóżmy teraz, że F (t) jest maksymalną wartością funkcji F na przedziale [0,t\. Jeśli nie jest, to mamy ograniczenie wartości funkcji przez poprzednią wartość maksymalną. Przy naszym założeniu

t

C(t) < C * + (C0 - C*)e~Mc< + a e - ^ 1 \ V0e^sF {t)e ^ s ds o

i wtedy

C (t) < C* + (Co - C *)e -'*‘ + aV oF(t)— — r.e/3t

Mc I P

Podstawiamy powyższą nierówność do trzeciego równania układu (18) i do- stajemy kolejną nierówność różniczkową, która prowadzi do ograniczoności F(t) dla dowolnego skończonego przedziału, a po powrocie do równania na C (t), także ograniczoność C(t).

Podobnie badamy zachowanie układu (18) w przypadku (3 = 0 (patrz [BS], [FŚ6]). .

Znajdziemy teraz punkty krytyczne układu (18) w przypadku £(ra) = 1 i zbadamy ich stabilność (dla £(m) < 1 można także przeprowadzić ana- logiczną analizę [Bh], [Fśl], z tym że odpowiednie wzory będą znacznie bardziej skomplikowane). Punkt krytyczny X = (V",Ć, F,fri) spełnia na- stępujący układ równań:

(21)

0 = ((3 - 7F)V,

0 = aVF - Mc(C -_C*),

0 = qC - (fJLf + 7 ]') V ) F ,

l 0 = oV —

W zależności od parametrów mogą istnieć dwa punkty krytyczne. Oznaczmy je A i B . Widzimy, że zawsze istnieje punkt krytyczny A = (0 ,C *,F *,0 ), gdzie F* = A4/, który odpowiada stanowi zdrowia organizmu - nie ma anty-

genu, organizm jest zdrowy, przeciwciała i komórki plazmatyczne pozostają na swoich poziomach fizjologicznych F*, C*.

Dla drugiego punktu krytycznego B mamy F = Po rozwiązaniu układu (21) z tym warunkiem dostajemy

- _ ~ 7 * 1*) omz g _ a 0 - r n i 2HcC m ^ P(aQ - 777/xc) j( a g - r\7/ic) '

Oczywiście fh = — . Punkt krytyczny B odpowiada chorobie chronicznej organizmu. W organizmie występuje niezerowe stężenie antygenu i rożne od fizjologicznych poziomy przeciwciał i komórek plazmatycznych oraz pewne, stosunkowo niewielkie, porażenie zainfekowanego organu, tzn. fh < m*.

Zauważyć należy, że przy założeniu £(m) = 1 ostatnie równanie nie wpływa na przebieg rozwiązań trzech pierwszych równań układu (18), a za- tem układ ten możemy zredukować do trzech równań

( V = ( p - j F ) V ,

(23) l Ć = aV (t - r)F (t - r) - fj,c{C - C % [ F = gC - (fif + r]'yV)F}

z punktami krytycznymi A = (0,(7*, F *), B = (V ,Ć ,& ), gdzie V, C są zdefiniowane wzorem (22). Zachowanie zmiennej m zależy tylko od funkcji V.

Punkt krytyczny B istnieje, jeśli spełnione są nierówności 1) (3 > 7 F * i a g > rpf[i,c lub 2) (3 < 7F* i ag < i^ync.

Oznacza to, że choroba chroniczna może wystąpić tylko wtedy, gdy organizm jest silny (duże współczynniki a i £>), ale antygen wykazuje wysoki współ- czynnik reprodukcji w stosunku do poziomu fizjologicznego przeciwciał, lub odwrotnie - organizm jest słaby, ale antygen jest także mało aktywny.

Zbadamy teraz zachowanie układu (23) w pobliżu punktów krytycznych.

Niech

V = v + V , C = c + Ć, F = / + F . W nowych zmiennych układ (23) ma postać

( v = { 0 -' y f- ' y F )v - ' y V f_ ,

(24) < ć = av(t — r ) f( t — r) + aF v (t — r) + a V f(t — r) — /Licc, { f = gc - (/if + rryV)f - rjyFv - r^vf.

Załóżmy teraz, że nowe zmienne mają normę nie większą niż e, zatem można odrzucić wyrazy wyższego rzędu (twierdzenie o linearyzacji dla rów- nań z opóźnionym argumentem ma podobną postać do analogicznego twier- dzenia dla równań bez opóźnienia, patrz np. [He]), co prowadzi do liniowego

układu równań z opóźnieniem postaci ( v = (/?_- 7 F )v - 7 V j,

(25) < ć = aFv(t — r) + aVf(t — r) — pcc, l / = - (M/ + 77^)/ ~ t? 7 ^ -

Tw ie r d z e n ie 3. Jeśli (3 < 7F *, to rozwiązanie układu (25) dla punktu krytycznego (0,C * ,F * ) dąży przy t —> +00 do punktu (0,0,0).

Dow ód. Dla V = 0 i F = F * pierwsze równanie układu (25) przyjmuje postać v = {(3 — 7 F*)v, czyli v(t) = voe~^,F*~^t —> 0 przy t —> + 00. Stąd i z drugiego równania

t

c(t) = Coe- * * + e_f*ct 5 aF*woe'‘' se -(l'F' - /S)(s“T) ds, 0

a zatem

p—(yF*—P)t _ p -pct

c(t) = coe-^* + aF*«oe<1,JP' - « T--- - 5 ---=;---* 0.

Mc + (3 - 7 ^*

Podobnie pokazujemy, że / —> 0 przy t —» +oo. ■

Twierdzenie to oznacza, że punkt A jest lokalnie asymptotycznie sta- bilny przy założeniu (3 < 7 F *, co można również sprawdzić przez znalezie- nie wartości własnych. Będziemy szukać rozwiązań układu (25) w postaci {voext, coext, fo e xt). Dostaniemy wtedy

( Avo = (j3_- 7F )v0 - 7^/0,

(26) < Aco = aF voe~ Xr + aV Joe~ Xr - p cc0, { A/o = QCo - (pLf + T]jy)fo ~ rnFvo.

Wartości A spełniające równanie charakterystyczne det M — 0, gdzie

/ (3 — j F — A 0 ~ jV

(27) M = a F e ~ XT - p c - A aV e~ Xr

\ —rjyF g —fif — r)jV — A

nazywamy wartościami własnymi układu (23). Macierz M powstaje, ana- logicznie jak w przypadku równań bez opóźnienia, przez zróżniczkowanie układu (26) po uq, co i /o odpowiednio. Tak jak w przypadku bez opóźnienia, jeśli wartości własne mają części rzeczywiste ujemne, to punkt krytyczny, któremu odpowiada macierz M, jest lokalnie asymptotycznie stabilny, a jeśli istnieje wartość własna o dodatniej części rzeczywistej, to taki punkt jest niestabilny.

Dla punktu A macierz M ma postać

/ ( 3 - j F * - \ 0 0 \

« F * e - Ar -Mc - A 0 V -nyF* e - m/ - A j

(28) M =

a równanie charakterystyczne

(/3 - 7 ir * _ A )(Mc + A)(M/ + A) = 0.

Łatwo zatem widać, że jeśli parametry modelu są dodatnie, to warunkiem koniecznym asymptotycznej stabilności punktu A jest fi < 7F *, natomiast jeśli zachodzi nierówność przeciwna, to punkt A jest niestabilny.

Rozpatrzmy teraz przypadek standardowego warunku początkowego, czyli V(h) = 0 dla h < 0, U(0) = VQ > 0 oraz C(h) = C*, F {h ) = F*

i m (h) = 0 dla h G [—r, 0].

UWAGA 2. Dla standardowego warunku początkowego rozwiązanie ukła- du (18) ma własność C(t) > C* dla t > 0.

Dowód. W takim przypadku równanie (17) zamienia się na nierówność postaci

C(t) > C* + (C0 - C*)e“M , a ponieważ Co = C*, to C(t) > C*. •

Dla takiego warunku początkowego zachodzi następujące twierdzenie:

Tw ie r d z e n ie 4. Jeśli

(3 < 7 F* oraz * ^ H f{lF * - P ) v P i

to rozwiązanie układu (18) dąży do punktu krytycznego A przy t —» + 00, przy czym funkcja V jest malejąca i V(t) < Voe~at, gdzie a = — P > 0.

Dowód. Z nierówności fi < 7 F* wynika, że na początku funkcja V ma- leje. Załóżmy zatem, że istnieje pierwszy taki punkt, w którym V przestaje maleć, czyli takie to> że V(to) = 0 oraz V(t) > 0 w pewnym prawostronnym otoczeniu punktu to. Stąd F (to) = ^ i F(to) < 0.

Z trzeciego równania układu (18) i podstawienia gC* — U fF* wynika, że

F (t0) = L>C(t(i) - (/.i/ + rriV(t0))^

> qC * - (M/ + rrV {t0))~ > u ff" - (n, + = o.

Otrzymujemy zatem sprzeczność z definicją punktu to- Wynika z tego, że V jest malejąca. Łatwo zauważyć (por. [Bh]), że zachodzi oszacowanie ekspo- nencjalne.

Dalsza część dowodu jest uznana w [Mkl] i [Bh] za oczywisty wniosek i w związku z tym pominięta. Wydaje się jednak, że nie jest to aż tak oczywiste.

Przedstawimy zatem dowód z [Fśl]. Zajmiemy się najpierw zmienną m. Dla dowolnego e > 0 zachodzi nierówność V (t) < e^ 1 dla dostatecznie dużych t.

Wobec tego m < e/im — fimm , skąd m(t) < e(l — e /Xmt) < e dla dostatecznie dużych t. Zatem m (t) —» 0 przy t —> +oo.

Pokażemy teraz, że funkcja F jest ograniczona. Załóżmy, że F (t) = max {F (s) : s < t } . Wtedy

Ć < aVoe~a^~T^F(t) - fic(C - C *).

Ponieważ Cq = G*, więc

e ~at _ e - p ct

C(t) < C* + aV0eaTe-^ ct \ e ^ se~asF {s) ds < C* + aV0F{t)e'

o Mc — a

Stosując tę nierówność do trzeciego równania układu (18) i ponownie pod- stawiając gC* = U fF *, dostajemy

F < MfF* + ( ^ [ * _ e aT( e - at - e '^ *) -

\ Mc ^ /

i stąd F(t) < G(£), gdzie

1 — e~a* 1 — e~^ct

TP* I ar

Gm = F exp ---e Mc - a

4- m/F*(1 — e~^ft) exp I — -

\Mc - a

a Mc

/ — 1 — e~a*

- l i f t

a

Łatwo zauważyć, że lim^+oo G(£) < (1+m/)F* exp(afjf^ ajeQT), zatem funk- cja F(t) jest ograniczona. Jeśli opuścimy założenie, że F(t) jest maksymalną wartością funkcji F na przedziale [0, t], to oczywiście także dostaniemy ogra- niczenie funkcji F (albo przez G, albo przez poprzednią wartość maksy- malną). Wobec tego lim^+oo F(£)F(£) = 0. Stosując takie same oszacowa- nia jak w przypadku funkcji m, dostajemy tezę twierdzenia. ■

W przypadku modelu (23) można udowodnić znacznie silniejsze twier- dzenie (patrz [FŚ4, 7]). Okazuje się, że nie tylko dla Vb < F * mamy zbieżność do punktu krytycznego A. Przy spełnieniu pewnych dodatkowych warunków punkt A okazuje się być globalnie stabilny. Nie będziemy tego faktu dowo- dzić, sformułujemy tylko odpowiednie twierdzenie.

Tw ie r d z e n ie 5. Jeśli ag > r)j(iic + p)e^T i (3 < jF * , to każde rozwiąza- nie modelu (23) ze standardowym warunkiem początkowym dąży do punktu krytycznego A.

W przypadku punktu B macierz M ma postać

/ -A 0

(29) M = a^e-XT ^{- A}

V -riP 6 ^{-M / -} rnV - X

- y V aV e~ Xr