Zagadnienia

7.1 Model drapieżnik – ofiara z ograniczoną pojemnością środowiska dla ofiar
7.2 Model z kryjówkami dla ofiar

7. Modele oddziaływań między dwiema populacjami II

7.1. Model drapieżnik – ofiara z ograniczoną pojemnością środowiska dla ofiar

Jak już wspomnieliśmy, w modelu Lotki – Volterry zostały pominięte różne istotne czynniki, gdyż w zamyśle Volterry taki model powinien być możliwie najprostszy. W szczególności wewnętrzna dynamika populacji ofiar została zbudowana w oparciu o model Malthusa zakładający nieograniczoność pojemności środowiska dla tego gatunku, co przy niewielkiej liczebności drapieżników prowadzi do gwałtownego wzrostu populacji ofiar i przyczynia się do okresowości rozwiązań. Załóżmy więc, że wewnętrzna dynamika gatunku $\mathcal{P}_{o}$ rządzi się innymi prawami, w szczególności liczebność tej populacji jest ograniczona przez pojemność środowiska. Najprostszym znanym nam modelem opisującym taką dynamikę jest równanie logistyczne, zbudujmy zatem model drapieżnik – ofiara w oparciu o to równanie. Wobec tego, przy oznaczeniach wprowadzonych w wyjściowym modelu Lotki – Volterry, model z ograniczoną pojemnością środowiska dla ofiar możemy opisać następującym układem równań

$\begin{array}[]{ccl}\dot{V}&=&rV\left(1-\frac{V}{K}\right)-aVP,\\ \dot{P}&=&-sP+abVP,\end{array}$

(7.1)

gdzie $K$ oznacza pojemność środowiska, analogicznie jak w modelu logistycznym. Dynamika rozwiązań układu (7.1) zależy w istotny sposób od tego parametru. Zauważmy, że postawowe własności takie jak istnienie, jednoznaczność, nieujemność i przedłużalność rozwiązań możemy zbadać w analogiczny sposób jak dla układu (6.1). Wnioskujemy zatem, że dla dowolnego nieujemnego warunku początkowego $(V_{0},P_{0})$ istnieje jednoznaczne, nieujemne rozwiązanie układu (7.1) określone dla wszystkich $t\geq 0$ . Dość prosto możemy wykazać też ograniczoność rozwiązań. Sformułujmy odpowiednie stwierdzenie

Stwierdzenie 7.1

Rozwiązania układu (7.1) dla nieujemnego warunku początkowego $(V_{0},P_{0})$ , $V_{0}$ , $P_{0}\geq 0$ , pozostają ograniczone. Co więcej, $V(t)\leq\max\{ V_{0},K\}$ dla $t\geq 0$ .

Z nieujemności rozwiązań układu (7.1) wynika następujące szacowanie $\dot{V}\leq rV\left(1-\tfrac{V}{K}\right)$ . Rzeczywiście, jeśli $V_{0}=0$ , to $V(t)\equiv 0$ , jak dla równania logistycznego, a jeśli $P_{0}=0$ , to $P(t)\equiv 0$ i $V$ spełnia równanie logistyczne, stąd $V(t)\leq\max(V_{0},K)$ . Dla $V_{0}$ , $P_{0}>0$ mamy $V(t)$ , $P(t)>0$ dla dowolnego $t>0$ . Rozważmy równanie logistyczne $\dot{x}=rx\left(1-\tfrac{x}{K}\right)$ z warunkiem początkowym $x_{0}=V_{0}>0$ oraz różnicę $y(t)=x(t)-V(t)$ . Zauważmy, że $\dot{y}(0)=\dot{x}(0)-\dot{V}(0)=aV_{0}P_{0}>0$ . Zatem na pewnym odcinku $(0,t^{*})$ zachodzi $x(t)>V(t)$ . Jeśli istnieje $t>0$ , takie że $V(t)>x(t)$ , to istnieje pierwsza taka chwila $\bar{t}$ , że $V(\bar{t})=x(\bar{t})$ oraz $\dot{y}(\bar{t})\leq 0$ , ale $\dot{y}(\bar{t})=aV(\bar{t})x(\bar{t})>0$ , co jest sprzeczne z definicją punktu $\bar{t}$ . Wynika stąd, że zmienna $V$ nie przekracza rozwiązań równania logistycznego z warunkiem początkowym $V_{0}$ . Zatem $V(t)\leq V_{{\max}}=\max\{ V_{0},K\}$ . Załóżmy teraz, że zmienna $P$ rośnie nieograniczenie. Mamy wtedy asymptotycznie $\dot{V}(t)\sim C_{1}-C_{2}P(t)$ , gdzie stałe $C_{1}$ , $C_{2}>0$ zależą od parametrów i warunku początkowego. Stąd $\dot{V}(t)\to-\infty$ przy $t\to+\infty$ , co przeczy ograniczoności $V$ . Zatem zmienna $P$ także pozostaje ograniczona.

∎

Zajmiemy się teraz analizą portretu fazowego. Izokliny zerowe są następującymi prostymi

dla zmiennej $V$ : $V=0$ lub $P=\tfrac{r}{a}\left(1-\frac{V}{K}\right)$ ;
dla zmiennej $P$ : $P=0$ lub $V=\tfrac{s}{ab}$ .

W zależności od parametru $K$ — nietrywialne izokliny mogą się przeciąć w analizowanej przestrzeni fazowej, gdy $K$ jest relatywnie duże, albo nie, gdy $K$ jest małe. Zatem liczba stanów stacjonarnych zależy od $K$ : zawsze są 2 stany $(0,0)$ i $(K,0)$ , natomiast dla $K>\tfrac{s}{ab}$ istnieje dodatni stan stacjonarny $(\bar{V},\bar{P})$ , gdzie $\bar{V}=\tfrac{s}{ab}$ i $\bar{P}=\tfrac{r}{a}\left(1-\frac{s}{abK}\right)$ . Wyznaczmy wartości własne dla poszczególnych stanów stacjonarnych. Macierz Jacobiego układu (7.1) ma postać

$MJ(V,P)=\left(\begin{array}[]{cc}r\left(1-2\frac{V}{K}\right)-aP&-aV\\ abP&-s+abV\end{array}\right).$

Stąd wyznaczamy wartości własne

dla $(0,0)$ mamy $\lambda _{1}=r>0$ , $\lambda _{2}=-s<0$ ;
dla $(K,0)$ : $\lambda _{1}=-r<0$ , $\lambda _{2}=abK-s$ ;
dla $(\bar{V},\bar{P})$ wartości własne są rozwiązaniem równania charakterystycznego

$\lambda^{2}+\frac{r\bar{V}}{K}\lambda+a^{2}b\bar{V}\bar{P}=0,$

więc w zależności od parametrów albo $\lambda _{1}$ i $\lambda _{2}$ są rzeczywiste ujemne, albo są zespolone i mają ujemne części rzeczywiste.

Wobec tego stan stacjonarny $(0,0)$ jest siodłem bez względu na parametry układu, z kolei $(K,0)$ jest albo węzłem stabilnym, o ile $K<\tfrac{s}{ab}$ , czyli dodatni stan stacjonarny nie istnieje, albo siodłem, gdy $K>\tfrac{s}{ab}$ , a wtedy istnieje $(\bar{V},\bar{P})$ , który jest albo węzłem, albo ogniskiem stabilnym. Typ punktu krytycznego $(\bar{V},\bar{P})$ zależy także od $K$ , jeśli $K$ jest duże, tak że $4a^{2}b^{2}K^{2}-4asbK-rs>0$ , to obserwujemy ognisko i gasnące oscylacje, w przeciwnym przypadku — węzeł i zbieżność do stanu stacjonarnego jest od pewnego momentu monotoniczna.

$Szkice portretu fazowego układu~\eqref{LV-poj} z zaznaczonymi izoklinami zerowymi i kierunkami przebiegu poszczególnych zmiennych w przypadku gdy istnieje dodatni punkt stacjonarny (po lewej) oraz w przypadku gdy dodatni punkt stacjonarny nie istnieje (po prawej).$

Rys. 7.1. Szkice portretu fazowego układu (7.1) z zaznaczonymi izoklinami zerowymi i kierunkami przebiegu poszczególnych zmiennych w przypadku gdy istnieje dodatni punkt stacjonarny (po lewej) oraz w przypadku gdy dodatni punkt stacjonarny nie istnieje (po prawej).

$Portret fazowy układu~\eqref{LV-poj} wraz z przykładowymi krzywymi fazowymi w przypadku gdy istnieje dodatni punkt stacjonarny (wykresy u góry) oraz w przypadku gdy dodatni punkt stacjonarny nie istnieje (wykres u dołu).$

Rys. 7.2. Portret fazowy układu (7.1) wraz z przykładowymi krzywymi fazowymi w przypadku gdy istnieje dodatni punkt stacjonarny (wykresy u góry) oraz w przypadku gdy dodatni punkt stacjonarny nie istnieje (wykres u dołu).

Przystępując do narysowania portretu fazowego zauważmy, że trywialne izokliny zerowe stanowią trajektorie układu. Rzeczywiście: jeśli $(0,P_{0})$ , $P_{0}>0$ , jest warunkiem początkowym, to oczywiście $\dot{V}=0$ , czyli zmienna $V$ nie zmienia się, $V(t)=0$ dla $t>0$ , natomiast dla zmiennej $P$ mamy równanie $\dot{P}=-sP$ , czyli $P(t)=P_{0}\exp(-st)$ dla $t>0$ . Podobnie dla warunku początkowego $(V_{0},0)$ , $V_{0}>0$ , zachodzi $\dot{P}=0$ , czyli $P(t)=0$ dla $t>0$ oraz $\dot{V}=rV\left(1-\tfrac{V}{K}\right)$ , zatem $V(t)$ spełnia równanie logistyczne, co zauważyliśmy już wcześniej w dowodzie stwierdzenia 7.1.

Jeśli izokliny nie przecinają się wewnątrz przestrzeni fazowej, to przestrzeń tę możemy podzielić na trzy obszary, por. rys. 7.1 i rys. 7.2: $A$ pod izokliną dla $V$ , $B$ nad tą izokliną i na lewo od izokliny dla $P$ oraz $C$ — na prawo od izokliny dla $P$ . Zauważmy, że rozwiązania dla warunku początkowego z $C$ muszą przejść do $B$ . Gdyby tak się nie stało, to mielibyśmy obie zmienne monotoniczne, przy czym zmienna $V$ pozostawałaby malejąca i ograniczona, a zmienna $P$ rosłaby nieograniczenie. W takim przypadku z pierwszego równania asymptotycznie dostajemy $\dot{V}\sim c_{1}-c_{2}P$ , gdzie $c_{1}$ i $c_{2}$ są dodatnimi stałymi, czyli $\dot{V}\to-\infty$ , co przeczy ograniczoności $V$ . Zatem rozwiązanie przechodzi do $B$ — tu obie zmienne są malejące, zatem albo pozostają w $B$ i wobec tego mają granicę, czyli zbiegają do stanu stacjonarnego $(K,0)$ na brzegu obszaru, albo przechodzą do $A$ . W $A$ obie zmienne także są monotoniczne, tyle że $V$ rośnie. Przebieg pola wektorowego pokazuje, że rozwiązania nie mogą wyjść z obszaru $A$ , więc także są zbieżne i zbiegają do $(K,0)$ .

W drugim przypadku, gdy istnieje dodatni stan stacjonarny, przebieg pola wektorowego wskazuje, podobnie jak dla modelu Lotki – Volterry, że rozwiązania okrążają punkt krytyczny $(\bar{V},\bar{P})$ w przestrzeni fazowej, por. rys. 7.1 i rys. 7.2. Wartości własne dla tego punktu implikują, że jest on lokalnie stabilny, ale na takiej podstawie nie możemy wykluczyć istnienia orbit zamkniętych, czyli cykli granicznych. Można to próbować wykazać na różne sposoby, w tym przypadku dzięki prostocie modelu działa kilka standardowych metod. Najprostszy sposób w przypadku układów dwuwymiarowych polega na zastosowaniu kryterium Dulaca – Bendixsona. Dla danego układu o prawej stronie opisanej funkcją $G(V,P)=\left(G_{1}(V,P),G_{2}(V,P)\right)^{T}$ rozważamy pewną funkcję $B:(\mathbb{R}^{+})^{2}\to\mathbb{R}$ klasy $\mathbf{C}^{1}$ , przy czym

$\text{div}BG=\frac{\partial BG_{1}}{\partial V}+\frac{\partial BG_{1}}{\partial P}$

nie zmienia znaku w $(\mathbb{R}^{+})^{2}$ . Wtedy w $(\mathbb{R}^{+})^{2}$ nie ma cykli, w szczególności nie ma cykli granicznych. Typowo stosowana funkcja $B$ ma postać $B(x,y)=\tfrac{1}{xy}$ . W przypadku układu (7.1) dostajemy

$\text{div}BG(V,P)=\tfrac{\partial}{\partial V}\left(\frac{r}{P}\left(1-\frac{V}{K}\right)-a\right)+\tfrac{\partial}{\partial P}\left(-\frac{s}{V}+ab\right)=-\frac{r}{PK}<0,$

z czego wnioskujemy globalną stabilność stanu stacjonarnego $(\bar{V},\bar{P})$ , gdyż jeśli nie ma cyklu granicznego, to z twierdzenia Poincarégo – Bendixsona wszystkie rozwiązania zbiegają do stanu stacjonarnego, w tym przypadku do jedynego stanu stabilnego $(\bar{V},\bar{P})$ .

Bardziej uniwersalna metoda (działająca także w przypadku modeli o większej liczbie zmiennych) polega na znalezieniu funkcjonału Lapunowa. W tym przypadku można zastosować funkcjonał Lapunowa zaproponowany dla modelu Lotki – Volterry. W celu uproszczenia obliczeń przeprowadźmy najpierw zamianę zmiennych $x=V-\bar{V}$ , $y=P-\bar{P}$ , która przesuwa stan stacjonarny do punktu $(0,0)$ . Zauważmy, że po zamianie zmiennych powstaje nowa przestrzeń fazowa $\mathcal{DF}=(-\bar{V},+\infty)\times(-\bar{P},+\infty)$ , a układ (7.1) przyjmuje postać

$\begin{array}[]{cclcl}\dot{x}&=&r(x+\bar{V})\left(1-\frac{x+\bar{V}}{K}\right)-a(x+\bar{V})(y+\bar{P})&=&-(x+\bar{V})\left(\frac{rx}{K}+ay\right),\\ \dot{y}&=&-s(y+\bar{P})+ab(x+\bar{V})(y+\bar{P})&=&abx(y+\bar{P}).\end{array}$

(7.2)

Funkcjonał Lapunowa dla układu (7.2) definiujemy jako

$L(x,y)=\alpha _{1}\left(x-\bar{V}\ln\frac{x+\bar{V}}{\bar{V}}\right)+\alpha _{2}\left(y-\bar{P}\ln\frac{y+\bar{P}}{\bar{P}}\right).$

Należy teraz sprawdzić odpowiednie własności funkcjonału $L$ : w szczególności mamy $L(x,y)\geq 0$ w przestrzeni $\mathcal{DF}$ oraz $L(x,y)=0\ \Longleftrightarrow\ x=0\ \text{i}\ y=0$ (czyli $V=\bar{V}$ i $P=\bar{P}$ ). Policzmy pochodną $L$ wzdłuż trajektorii układu (7.2)

$\begin{split}\dot{L}(x(t),y(t))&=\alpha _{1}\frac{x}{x+\bar{V}}\left(-(x+\bar{V})\left(\frac{rx}{K}+ay\right)\right)+\alpha _{2}\frac{y}{y+\bar{P}}\left(abx(y+\bar{P})\right)=\\ &=-\alpha _{1}\frac{rx^{2}}{K}\leq 0\\ \end{split}$

dla $\alpha _{1}=b$ i $\alpha _{2}=1$ . Wynika stąd globalna stabilność dodatniego stanu stacjonarnego. Aby wykazać jego globalną asymptotyczną stabilność wystarczy udowodnić, że $L(V(t),P(t))$ jest ściśle malejąca, co pozostawiamy jako ćwiczenie.

$Przebieg przykładowych rozwiązań układu~\eqref{LV-poj}.$

Rys. 7.3. Przebieg przykładowych rozwiązań układu (7.1). Wykresy znajdujące się na górze przedstawiają rozwiązania układu (7.1) w przypadku gdy istnieje dodatni punkt stacjonarny i jest on odpowiednio ogniskiem stabilnym (po lewej) lub węzłem stabilnym (po prawej). Wykres na dole przedstawia rozwiązanie gdy nie istnieje dodatni punkt stacjonarny.

Zauważmy, że tak zdefiniowany układ ma stany stacjonarne asymptotycznie stabilne, przy czym dla $K<\tfrac{s}{ab}$ wszystkie rozwiązania (oprócz rozwiązań na brzegu przestrzeni fazowej) zbiegają do stanu stacjonarnego $(K,0)$ , natomiast dla $K>\tfrac{s}{ab}$ mamy zbieżność do dodatniego stanu stacjonarnego $(\bar{V},\bar{P})$ , a przy $K=\tfrac{s}{ab}$ następuje bifurkacja siodło – węzeł. Taki typ bifurkacji oznacza, że zmienia się charakter punktu krytycznego z siodła na węzeł, przy czym dla bifurkacyjnej wartości $K=\frac{s}{ab}$ parametru $K$ ze stanu stacjonarnego $(K,0)$ bifurkuje dodatni stan stacjonarny i dla takiego $K$ rozwiązania te sklejają się w jedno i choć nie przeprowadziliśmy analizy portretu fazowego dla tego przypadku, to doświadczenia związane z badaniem przebiegu orbit w pozostałych przypadkach pozwalają nam wnioskować, że także i dla bifurkacyjnej wartości $K$ mamy globalną stabilność stanu stacjonarnego $(K,0)$ . Zauważmy dalej, że niewielka zmiana prawej strony układu nie prowadzi w przypadku tego modelu do drastycznych zmian zachowania rozwiązań, zatem taki układ nie ma już niepożądanej własności niestabilności strukturalnej charakteryzującej układ (6.1).

7.2. Model z kryjówkami dla ofiar

Inny sposób na uzyskanie stabilności strukturalnej stanowi przyjęcie założenia, że część ofiar jest stale niedostępna dla drapieżników — w środowisku występuje pewna liczba kryjówek, w których może się schować określona liczba ofiar. Niech liczba ofiar, która może się schować przed drapieżnikiem wynosi $K$ . Jeśli zatem w chwili początkowej $t=0$ liczebność całkowita ofiar $V(t)$ nie przekracza $K$ , to wszystkie ofiary chowają się przed drapieżnikami i drapieżniki nie są w stanie nic upolować. W takiej sytuacji liczebność ofiar rośnie i w pewnym momencie przekracza $K$ — dopiero wtedy możemy mówić o właściwym układzie drapieżnik – ofiara i zastosować odpowiedni model. W modelu tym zakładamy, że liczba ofiar dostępnych dla drapieżników wynosi $V-K$ , $V>K$ , zatem tylko takie ofiary mogą zostać upolowane i wpłynąć na rozwój populacji drapieżników. Wyjściowy model Lotki – Volterry przyjmuje przy takim założeniu postać

$\begin{array}[]{ccl}\dot{V}&=&rV-a(V-K)P,\\ \dot{P}&=&-sP+ab(V-K)P,\end{array}$

(7.3)

przy czym parametry $r$ , $a$ , $s$ , $b$ mają takie same znaczenie jak dla układu (6.1), natomiast $K$ oznacza tu liczbę kryjówek (w odróżnieniu od typowego oznaczenia pojemności środowiska).

Lokalne istnienie i jednoznaczność rozwiązań układu (7.3) ponownie wynikają z wielomianowej postaci prawej strony układu. Ponieważ model ma sens dla $V\geq K$ , to ograniczymy nasze rozważania do przestrzeni fazowej $\mathcal{DF}=[K,+\infty)\times\mathbb{R}^{+}$ . Zauważmy, że

Stwierdzenie 7.2

Zbiór $\mathcal{DF}$ jest przestrzenią niezmienniczą dla układu (7.3).

Niech $(V_{0},P_{0})\in\mathcal{DF}$ będzie warunkiem początkowym. Jeśli $P_{0}=0$ , to $\dot{P}\equiv 0$ i rozwiązanie pozostaje na osi poziomej, co więcej, żadne rozwiązanie z wnętrza przestrzeni fazowej $\mathcal{DF}$ nie osiągnie tego brzegu ze względu na jednoznaczność rozwiązań. Niech $V_{0}=K$ . Wtedy $\dot{V}(0)=rK>0$ , zatem rozwiązanie rośnie w pewnym przedziale czasu. Jeśli istnieje chwila $\tilde{t}>0$ , w której rozwiązanie osiąga punkt na brzegu, $V(\tilde{t})=K$ , to oczywiście $\dot{V}(\tilde{t})=rK>0$ , zatem takie rozwiązanie nie może przekroczyć tego brzegu tego brzegu z prawej strony na lewą. Stąd $V(t)>K$ dla wszystkich $t>0$ .

∎

Z niezmienniczości zbioru $\mathcal{DF}$ wynika od razu globalne istnienie rozwiązań: dla zmiennej $V$ dostajemy oszacowanie liniowe $\dot{V}\leq rV$ , czyli $V(t)\leq V_{0}\text{e}^{{rt}}$ i dalej na każdym ograniczonym przedziale $[0,t^{*})$ szacujemy $\dot{P}\leq AP$ , gdzie $A=abV_{0}\text{e}^{{rt^{*}}}-s$ , podobnie jak w przypadku układu (6.1).

$Portret fazowy układu~\eqref{LV-kryj} wraz z przykładowymi krzywymi fazowymi.$

Rys. 7.4. Portret fazowy układu (7.3) w przypadku gdy dodatni stan stacjonarny jest stabilnym ogniskiem (po lewej) lub węzłem (po prawej) oraz przykładowe krzywe fazowe w obu przypadkach.

Asymptotykę układu (7.3) badamy znanymi już metodami. Przeprowadzenie analizy portretu fazowego z zastosowaniem np. kryterium Dulaca – Bedixsona pozwala wysnuć następujący wniosek

Wniosek 7.1

Wszystkie rozwiązania układu (7.3) z warunkiem początkowym $(V_{0},P_{0})$ , takim że $V_{0}\geq K$ , $P_{0}>0$ , zbiegają asymptotycznie do dodatniego stanu stacjonarnego $(\bar{V},\bar{P})$ , $\bar{V}=K+\tfrac{s}{ab}$ , $\bar{P}=\tfrac{rb}{s}\bar{V}$ .

Zauważmy też, że podobnie jak dla układu (7.1), dodatni stan stacjonarny w zależności od parametrów (np. liczby ofiar mogących schować się w kryjówkach) może być albo ogniskiem, albo węzłem stabilnym, por. rys. 7.4.

Podsumowując — zarówno w przypadku modelu z kryjówkami dla ofiar jak i z pojemnością środowiska dostaliśmy globalną stabilność jednego ze stanów stacjonarnych, przy czym w przypadku modelu z kryjówkami jest to dodatni stan stacjonarny, natomiast w przypadku modelu z ograniczoną pojemnością środowiska może to być albo dodatni stan stacjonarny, albo stan brzegowy opisujący wyginięcie (ekstynkcję) populacji drapieżników. Jeśli stabilny stan stacjonarny jest węzłem, to rozwiązania od pewnego momentu $\bar{t}$ pozostają monotoniczne, natomiast jeśli jest ogniskiem, to obserwujemy oscylacje i — ze względu na globalną stabilność — są to oscylacje gasnące.

Treść automatycznie generowana z plików źródłowych LaTeXa za pomocą oprogramowania wykorzystującego LaTeXML.

Modelowanie matematyczne w biologii i medycynie