Twierdzenie 5.5.1 Niech układ o parametrze :
5.6 Matematyczny model wzrostu nowotworu
W ostatnim czasie teoria bifurkacji rr przyciągnęła dużą uwagę dając nowe możliwości w celu zbudowania różnorodnych modeli matematycznych. W niniejszym paragrafie rozpatrzymy jeden z takich model, który szczególnie dobrze ilustruje ich charakter. Dotyczy on reakcji immunologicznej osobnika na obcą tkankę
( w naszym przypadku będzie to nowotwór ). Model ten został zbudowany przez Rescigno i de Lisi’ego (1977), tutaj podamy jedynie jego krótki opis.
5.6.1 Budowa modelu.
Komórki nowotworu zawierają szczególny materiał ( antygeny ), który to wywołuje reakcje immunologiczną w wielu organizmach żywych. Reakcja ta polega na tym, że produkowane są specyficzne komórki – limfocyty, które atakują i zabijają komórki rakowe.
Model operuje następującymi zmiennymi ( we wszystkich przypadkach mamy na uwadze rozmiar populacji komórek ) : a) L – swobodne limfocyty na powierzchni nowotworu.
b) C – komórki rakowe ( wewnętrzne i na powierzchni nowotworu ) c) Cs - komórki rakowe na powierzchni nowotworu
d) C- - komórki rakowe na powierzchni nowotworu i niezwiązane z limfocytami.
e) Cf - komórki rakowe wewnątrz i na powierzchni nowotworu , nie związane z limfocytami.
Z tych definicji bezpośrednio wynika, że :
C = Cf – C- + Cs (5.88)
Zakłada się ponadto, że sferyczna forma nowotworu nie zmienia się, tak że :
Cs = K1C2/3 (5.89)
Gdzie : K1 – stała
Zakłada się również, że oddziaływanie komórek rakowych z limfocytami następuje tylko na powierzchni nowotworu.
Nie wszystkie komórki rakowe poddawane są działaniu limfocytów, dlatego tylko w części przypadków, kiedy spotyka się komórka rakowa z swobodnym limfocytem, następuje ich wiązanie. Przyjmiemy, że między ilością swobodnych i związanych limfocytów istnieje zależność o postaci :
Cs – C- = K2C-L (5.90)
Gdzie : K2 – stała
Z (5.88) i (5.89) wynika, że :
Cf = C – K1K2 LC2/3 / ( 1 + K2L ) (5.91)
C- = K1C2/3 / ( 1 + K2L ) (5.92)
Z równań tych wynika, że zmienne L, C można przyjąć jako podstawowe zmienne danego modelu.
Zakłada się, że wielkość L• / L dla populacji limfocytów ma dwie składowe : a) o stałym poziomie śmiertelności λ1
b) o poziomie śmiertelności α’1C- [ 1 – (L/LM ) ]
Wyrażenie b) pokazuje, że, kiedy L jest małe, stymulacja swobodnych limfocytów wzrasta liniowo ze wzrostem C- i że istnieje maksymalny rozmiar populacji LM, przy którym poziom stymulacji jest równy zeru.
Zatem, L spełnia równanie :
L• = - λ1L + α’1C- [ 1 – (L/LM ) ] (5.93)
Prędkość wzrostu populacji komórek nowotworowych zadana jest poprzez równanie :
C• = λ2Cf – α’2C- L (5.94) Pierwszy człon w prawej części równania (5.94) opisuje wzrost nowotworu, nie podlegającego działaniu limfocytów, drugi człon uwzględnia oddziaływanie swobodnych limfocytów z komórkami nowotworowymi na powierzchni nowotworu.
Podstawiając wartości C- i Cf z równań (5.91) i (5.92) możemy przepisać (5.93) i (5.94) w postaci :
x• = - λ1x + { α1xy2/3 [ 1 – (x/c ) ] / ( 1 + x ) } (5.95)
y• = λ2y – α2xy2/3 / ( 1 + x ) (5.95)
gdzie :
x = K2L , c = K2LM , y = K2C
λ1 , λ2 , α1, α2 – są stałymi dodatnimi parametrami.
Ponieważ x i y – są rozmiarami populacji, powinny być one nieujemne, a x nie może przewyższać c, ponieważ L jest ograniczone z góry przez wielkość LM.
Teraz przejdziemy do badania jakościowych następstw równań dynamicznych (5.95) i postaramy się zbudować dla tego układu portrety fazowe, pojawiające się dla różnych obszarów zmienności w/w parametrów.
5.6.2 Analiza dynamiki.
Łatwo zauważyć, że układ (5.95) posiada siodło w początku współrzędnych przy wszystkich wartościach parametrów, wynika to z twierdzenia o linearyzacji. Dodatnie półosie x i y są trajektoriami danego układu, są one separatysami siodła.
Aby zbadać nietrywialne punkty stałe, układ (5.95) zapiszemy w następującej postaci :
x• = xf(x, y) , y• = -y2/3g(x, y ) (5.96)
gdzie :
f(x, y) = - λ1x + { α1y2/3 [ 1 – (x/c ) ] / ( 1 + x ) } (5.97)
g(x, y) = λ2y1/3 – [ α2x / ( 1 + x ) ] (5.98)
Równania dla punktów stałych mają postać : f(x, y) = 0 , g(x, y) = 0
Z równań tych wynika, że :
y2/3 = ( λ1 /α1 ){ ( 1 + x ) / [ 1 – (x/c )] } = { (λ2 /α2 ) [ x / ( 1 + x ) ] }2 tak, że x-współrzędne nietrywialnych punktów stałych spełniają równanie :
ψ(x) = x2[ 1 – (x/c) ]/ ( 1 + x )3 = (λ1 /α1 ) ( λ2 /α2 )2 (5.99) Funkcja ψ(x) posiada jedno globalne maksimum przy x* = 2c/( c + 3 ), przy czym ψ(x*) = 4c2 /27( c + 1 )2 .
Podstawmy teraz :
µ1 = [ (λ1 /α1 ) ( λ2 /α2 )2 ] – [ 4c2 /27( c + 1 )2 ] (5.100) Równanie (5.99) :
a) nie posiada pierwiastków rzeczywistych przy µ1 > 0 b) posiada równo jeden pierwiastek x = x*, przy µ1 = 0
c) posiada dwa pierwiastki x*1, x*2 , 0 < x*1 < 2c/ ( c + 3 ) < x*2 < c przy µ1 < 0.
Geometrycznie przypadki a), b) i c) odpowiadają charakterowi przecięcia krzywych f(x, y) = 0 i g(x, y) = 0 pokazanych na rysunku 5.25.
Rys. 5.25 Równanie f(x, y) = 0 określa krzywą, wypukłą w dół i posiadającą asymptotykę x = c ; krzywa g(x, y) = 0 posiada asymptotykę y = ( α2 /λ2 )3
Możliwe położenia wzajemne tych krzywych pokazano na rysunkach : a) µ1 > 0 , b) µ1 = 0 , c) µ1< 0.
Rozpatrzmy teraz portret fazowy dla każdego z tych przypadków : a) µ1 > 0.
Nie istnieją nietrywialne punkty stałe. Jeśli przypomnimy, że początek współrzędnych to siodło i uwzględnimy znaki x• i y• , to można narysować portret fazowy przedstawiony na rysunku 5.26.
Zauważmy, że y → ∞ przy t → ∞ dla wszystkich stanów początkowych populacji. Odpowiada to niekontrolowanemu wzrostowi nowotworu.
Rys. 5.26 Portret fazowy układu (5.95) przy µ1 > 0. Wszystkie trajektorie asymptotycznie dążą do x = x przy t → ∞.
b) µ1 = 0.
Istnieje jeden nietrywialny punkt stały (x* , y* ), gdzie x* = 2c/ (c + 3 ).
Macierz współczynników zlinearyzowanego układu ma postać :
W = [ xfx xfy ] | (5.101)
[ y2/3gx y2/3gy ] | (x*, y* ) ponieważ f(x*, y* ) = g(x*, y* ) = 0.
Wykorzystaliśmy tutaj następujące oznaczenia : fx = ∂f/∂x , fy = ∂f/∂y , gx = ∂g/∂x , gy = ∂g/∂y Mamy dalej :
det W = x*y*2/3( fxgy – fygx ) | (x*, y* ) (5.102)
Jednak z rysunku 5.25 b widać, że przy µ1 = 0 nachylenia krzywych f(x, y) = 0 i g(x, y) = 0 w punkcie (x*, y* ) są jednakowe. Stąd wynika, ze :
fx /fy | (x*, y* ) = gx /gy | (x*, y* ) (5.103)
tak, że det W = 0. W tym przypadku (x*, y* ) – jest nieprostym punktem stałym i nie możemy określić jego charakteru za pomocą linearyzacji. Później powrócimy do dokładnego zbadania charakteru tego punktu stałego.
Globalny charakter portretu fazowego przy µ = 0 jest oczywiście taki, że dla większości stanów początkowych ( obszar stabilności dla punktu (x*, y* ) jest skończony ) pojawia się niekontrolowany wzrost nowotworu.
c) µ1 < 0.
Istnieją dwa nietrywialne punkty stałe : P1(x*1 , y*1 ) i P2(x*2 , y*2 ), przy czym 0 < x*1< x* < x*2 < c, gdzie x* = 2c/ ( c + 3 )
Jeśli Wi( i = 1, 2 ) – są macierzami współczynników zlinearyzowanego układu w punktach ( x*i , y*i ), to wzory (5.97) , (5.98) i (5.102) dają nam :
det Wi = {α1λ2xy2/3 [ (2c/x) – c – 3 ] / [ 3c ( 1 + x )2 ] }| (x*, y* ) (5.104) Stad wynika, że det W1 jest dodatni, a det W2 ujemny. Od razu możemy wnioskować, że punkt P2 jest siodłem, a stabilność lub niestabilność P1 określona jest poprzez znak Tr W1.
Jeśli obliczyć macierz (5.101) w punkcie ( x*1 , y*1 ), to otrzymamy :
Tr W1 = ( xfx + gyy2/3 ) | (x*1 , y*1 ) = 1/3 λ2 { 1 – 3(λ1/ λ2 ) ( 1 + c/ c ) [ x / ( 1 + x ) ( 1 – x/c )] } | x*1 = η(x*1 ) (5.105)
W celu wykluczenia y wykorzystaliśmy tutaj równość (5.97).
Funkcja :
x / ( 1 + x ) ( 1 – x/c ) ściśle wzrasta na odcinku (0, c) , a zatem η(x) jest funkcją x, ściśle malejącą na odcinku (0, x* ) ( przypomnijmy, że 0 < x* < c )
Dla Tr W1 = η( x*1) mogą pojawiać się różne możliwości, tak jak to widać na rysunku 5.27.
Rys. 5.27 Możliwe wykresy funkcji η(x). Zauważmy, że Tr W1 = η(x*1) jest dodatnie dla krzywych a) – c), równe zeru dla krzywej d) i ujemne dla krzywej e).
Zauważmy, że z η(x* ) ≥ 0 wynika η(x*1) > 0 ( krzywe a i b ), podczas gdy w przypadku η(x* ) < 0 wartość η(x*1) może być : dodatnia ( krzywa c ), zerowa ( krzywa d ) lub ujemna ( krzywa e ).
Otrzymany wynik możemy w następujący sposób wykorzystać z wykorzystaniem parametrów badanych modeli.
Jeśli x = x* = 2c/ ( c + 3 ), to :
η(x* ) = 1/3 λ2 { 1 – ( λ1/ λ2 ) [ 2( c + 3 ) / ( c + 1 )] } (5.106) Wielkość (5.106) zeruje się ( krzywa b na rys. 5.27 ) przy λ1/ λ2 = ( c + 1 ) / 2( c + 3 ).
Podstawimy teraz :
µ2 = ( λ1/ λ2 ) – [ ( c + 1 ) / 2( c + 3 )]
Wtedy otrzymamy, że :
a) jeśli µ2 ≤ 0, to η(x* ) ≥ 0 i Tr W1 = η( x*1 ) > 0, a zatem punkt P1 jest niestabilny.
b) jeśli µ2 > 0, to znak Tr W1 nie jest określony.
Jeśli µ2 > 0, to możemy twierdzić, że krzywa η(x ) posiada charakter krzywych c, d i e przy kolejnym wzroście parametru µ2, poczynając od zera. Odpowiednio do tego Tr W1 = η( x*1 ) jest na początku dodatni, następnie równy zeru i na koniec ujemny. Zatem, przy wzroście µ2 od 0 do ∞ ( przy stałej wartości µ1 ) punkt P1 jest niestabilny przy µ2 = 0 i przy wystarczająco dużych wartościach µ2 staje się stabilny.
Otrzymane wyniki podsumowuje rysunek 5.28. Pokazano na nim płaszczyznę µ1, µ2 podzieloną na następujące obszary : obszar A, gdzie istnieje tylko jeden punkt stały w początku współrzędnych, obszar B, gdzie punkt P1 jest niestabilny i obszar C, gdzie punkt P2 jest stabilny. Nie określiliśmy formy granicy między obszarami B i C ( rys. 5.28 ), dlatego też pokazuje go jedynie schematycznie linia punktowa.
Nie ma nietrywialnych punktów stałych
Rys. 5.28 Wynik badania stabilności względem przybliżenia liniowego dla układu (5.95). Charakter punktu osobliwego na prostej µ1 = 0 i granica między obszarami B i C nie mogą być określona z pomocą linearyzacji.
Lokalne zachowanie ukazane na rysunku 5.28, jest jakościowo równoważne zachowaniu prostszego dwuparametrowego układu :
x• = - ( µ-1 + y2 ) , y• = - ( x + µ-2 y + y2 ) (5.107)
Układ ten zbadał Takens (1974 ). Podstawowe wyniki globalnej analizy tego układu pokazano na rysunkach 5.29 i 5.30.
Zauważmy, że lokalne portrety fazowe w punktach stałych na rysunku 5.30b i c pokrywają się.
Zatem, z punktu widzenia zachowania lokalnego obszary B’1 i B’2 są równoważne i możemy ustanowić analogię między A i A’ ; B i B’ = B’1 ∪ B’2 ; C i C’.
Rys. 5.29 Wynik globalnej analizy układu (5.107). Płaszczyzna µ-1, µ-2 rozbita została na cztery obszary : A’, B’1 , B’2 , C’. Portrety fazowe w każdym z tych obszarów pokazano na rysunku 5.30.
Analogia między układami (5.95) i (5.107) rozciąga się nie tylko na zachowanie lokalne. W istocie wszystkie globalne portrety fazowe dla (5.95), otrzymane przez Swan’a (1977 ), posiadają jakościowe analogi w wynikach badań Takensa dla układu (5.107). Oczywiście, to nie dowodzi, że układy (5.95) i (5.107) są jakościowo równoważne, ponieważ badania Swan’a nie wyczerpują wszystkich możliwych portretów fazowych układu (5.95). Jednakże taka jakościowa
równoważność jest całkowicie możliwa.
Założenie o tym, że układy (5.95) i (5.107) są jakościowo równoważne, pociąga za sobą ważne następstwo dla układu (5.95). Takens podał pełną analizę układu (5.107) i może ona posłużyć jako wzorzec dla analizy układu (5.95).
Przykładowo, moglibyśmy założyć, ze granica między obszarami C i B na rysunku 5.28 określa bifurkacje Hopfa, tak jak ma to miejsce dla granicy między obszarami C’ i B’ na rysunku 5.29. Przypomnijmy, że linearyzacja układu (5.95) posiada czysto urojone wartości własne na granicy między C i B ( Tr W1 = 0 ) i zachodzi szybka zmiana stabilności punktu P1. Jest to „symptom” bifurkacji Hopfa. Cykl graniczny, charakterystyczny dla tej bifurkacji, również pojawia się w badaniach Swan’a. W jego pracy nie wspomniano jednak, jaka bifurkacja zachodzi na granicy między obszarami B’1 i B’2 ( rys. 5.29 ), jednakże z hipotezy o jakościowej równoważności wynika, że ma ona miejsce.
Rys. 5.30 Portrety fazowe dla układu (5.107), kiedy punkt (µ-1, µ-2 ) znajduje się w obszarach : a) C’ ; b) B’2 ; c) B’1 ; d) A’
Przy zmniejszaniu się µ2 ( przy ustalonej wartości µ1 < 0 ) możemy oczekiwać pojawienia się cyklu granicznego na granicy między obszarami B i C powodowanego przez bifurkacje Hopfa. Zatem, cykl graniczny powinien wzrastać do tej pory, aż nie dojdzie on do siodła P2, wtedy powinna nastąpić następna bifurkacja i ograniczone drgania pojawiające się w przypadku obecności cyklu granicznego, powinny zmienić niekontrolowany wzrost nowotworu ( zobacz rys. 5.30 c ).
Na koniec, jakościowa równoważność układów (5.95) i (5.107) pozwoliłaby nam zbadać charakter punktu stałego w początku współrzędnych układu (5.95) dla µ1= 0. Klasyfikacja Takensa zawiera wszystkie takie „zwyrodniałe”
przypadki dla układu (5.107). Przykładowo, przy µ-1= µ-2 = 0 nieprosty punkt stały układu (5.107) w początku współrzędnych posiada portret fazowy przedstawiony na rysunku 5.31. Portret ten jest równoważny portretowi otrzymanemu przez Swan’a dla nietrywialnego punktu stałego (x*, y* ) układu (5.95).
Rys. 5.31 Nieprosty punkt stały układu (5.107) Rys. 5.32 Portret fazowy dla układu (5.108) ze źródłem limfocytów.
przy µ-1= µ-2 = 0 Obszar stabilności dla punktu ( x0 , 0 ) zakreskowano.
Wszystkie cztery portrety fazowe pokazane na rysunku 5.30, realizują się i dla układu (5.95) przy pewnych wartościach parametrów ( µ1, µ2 ). Można się przekonać o tym, że najczęściej następuje niekontrolowany wzrost nowotworu.
Dla pozostałych obszarów ( µ1, µ2 ), gdzie portrety fazowe są analogiczne do rys. 5.30 a lub b, niekontrolowany wzrost może również nastąpić, jeśli stan początkowy populacji należy do określonych i ograniczonych obszarów płaszczyzny x, y.
Kluczową ideą dla danych modeli jest poszukiwanie jakichkolwiek realnych możliwości w celu rozszerzenia obszaru stabilności. Przykładowo, Rescigno i De Lisi (1977) rozpatrywali model ze źródłem limfocytów, które charakteryzuje się stałą prędkością ich dopływu. W takim przypadku równania dynamiczne (5.95) zamieniamy na równania o postaci : x• = - λ1( x – x0 ) + α1xy2/3 [ 1 – (x/c ) ] / ( 1 + x )] (5.108)
y• = λ2y – α2xy2/3 / ( 1 + x ) (5.108)
gdzie człon : λ1x0, x0 > 0 odpowiada źródłu limfocytów.
Układ (5.108) posiada stabilny punkt stały o współrzędnych ( x0, 0 ), który odpowiada całkowitej reemisji nowotworu.
Jeden z możliwych portretów fazowych przedstawiono na rysunku 5.32. Rozmiar obszaru stabilności punktu stałego ( x0, 0 ) zwiększa się przy zwiększaniu wielkości λ1λ22 / α1α22
Dalsze szczegóły omówionego modelu zainteresowany czytelnik znajdzie w artykule Swan’a (1977).
Ćwiczenia
Do punktu 5.4
1. Pokazać, ze funkcja V(x1, x2 ) = x12 + x22 jest silną funkcją Lapunowa w początku współrzędnych dla każdego z następujących układów :
2. Znaleźć obszar stabilności w początku współrzędnych dla wszystkich układów z ćwiczenia 1.
3. Pokazać, że funkcja V(x1, x2 ) = x12 + x22 jest silną funkcją Lapunowa w początku współrzędnych dla każdego z następujących układów :
Dla każdego z nich początek współrzędnych jest asymptotycznie stabilny.
4. Dowieść, że jeśli V jest silną funkcją Lapunowa dla równania x• = − X (x ) w pewnym otoczeniu początku współrzędnych, to x• = X (x ) posiada w początku współrzędnych niestabilny punkt stały. Wykorzystując ten wynik, pokazać, że układy :
a) x1• = x13 , x2• = x22 b) x1• = sin(x1) , x2• = sin(x2 )
c) x1• = - x13 + 2x12 sin(x1), x2• = x2 sin2(x2 ) w początku współrzędnych są niestabilne.
5. Dowieść, ze równania :
a) x•• + x• – 1/3 x•3 + x = 0 , b) x•• + x• sin(x•3 ) + x = 0 c) x•• + x• + x3 = 0 , c) x•• + x•3 + x3 = 0 posiadają asymptotycznie stabilne zerowe rozwiązania.
6. Dowieść, że funkcja V(x1, x2 ) = ax12 + 2bx1x2 + cx22 jest określona dodatnio wtedy i tylko wtedy, kiedy a > 0 i ac > b2. Z pomocą tego kryterium ( lub inaczej ) dowieść, że układ :
x1• = x2 , x2• = - x1 – x2 – ( x1 + 2x2 )( x22 – 1 ) jest asymptotycznie stabilny w początku współrzędnych.
Rozpatrzyć obszar | x2 | < 1. Znaleźć obszar stabilności.
7. Znaleźć obszar stabilności dla następujących układów, zbudować odpowiednie funkcje Lapunowa : a) x1• = x2 – x1( 1 – x12 – x22 )( x12 + x22 + 1 )
x2• = - x1 – x2( 1 – x12 – x22 )( x12 + x22 + 1 ) b) x1• = x2
x2• = - x2 + x23 + x15
8. Wykorzystać funkcje V(x1, x2 ) = (x1/a )2 + ( x2/b)2 i pokazać, że układ : x1• = x1( x1 – a ) , x2• = x2( x2 – b ) ; a, b > 0
posiada asymptotycznie stabilny początek współrzędnych.
Pokazać, że wszystkie trajektorie dążą do początku współrzędnych przy t → ∞ w obszarze : ( x1/a )2 + ( x2/b)2 < 1
9. Niech będzie zadany układ : x1• = x2
x2• = x2 – x13
pokazać, że można dobrać funkcję dodatnio określoną o postaci : V(x1, x2 ) = ax14 + bx12 + cx1x2 + dx22
taką, że V• (x1, x2 ) również będzie dodatnio określona.
Wyprowadzić z tego faktu, wniosek, że początek współrzędnych jest niestabilny.
10. Pokazać, że dla układu : x1• = x22 – x12 , x2• = 2x1x2
początek współrzędnych jest niestabilny. Wykorzystać funkcję : V(x1, x2 ) = 3x1x22 – x13
11. Pokazać, ze punkt stały w początku współrzędnych dla układu : x1• = x14 , x2• = 2x12x22 – x22
jest niestabilny, zastosować funkcję V(x1, x2 ) = αx1 + βx2 i dobrać w odpowiedni sposób stałe α, β.
Sprawdzić niestabilność początku współrzędnych, badając zachowanie separatys.