Zagadnienia

6.1 Model Lotki – Volterry
6.2 Konstruktywna krytyka modelu Lotki – Volterry

6. Modele oddziaływań między dwiema populacjami I

Przejdziemy teraz do omówienia zagadnień związanych z modelowaniem ekosystemu, w którym występują dwie populacje. W zależności od oddziaływań między tymi populacjami, rozróżniamy trzy typy ekosystemów

układ drapieżnik – ofiara, w którym jeden z gatunków ekosystemu jest pożywieniem dla drugiego;
konkurencja, w którym gatunki rywalizują o zasoby środowiska (pożywienie, miejsca lęgowe itp.);
symbioza — zjawisko współżycia przynajmniej dwóch gatunków, które wpływają pozytywnie na siebie (mutualizm) lub wyraźne korzyści odnosi jedynie jedna ze stron, nie szkodząc drugiej (komensalizm).

W tym rozdziale zajmiemy się opisem średnich zagęszczeń obu populacji, nie uwzględnimy przestrzennego rozmieszczenia osobników w ekosystemie, co oczywiście kryje dodatkowe założenie, że osobniki są rozmieszczone jednorodnie i jest ich dostatecznie dużo, aby można było mówić o średnim zagęszczeniu.

6.1. Model Lotki – Volterry

Wprowadzimy najpierw pierwszy historycznie model opisujący oddziaływania dwóch populacji w ekosystemie. Dotyczy on układu drapieżnik – ofiara i został zaproponowany równolegle jako model populacyjny przez Vito Volterrę [17] oraz jako model łańcucha reakcji biochemicznych przez Alfreda Lotkę [10]. Volterra zaproponował ten model w celu wyjaśnienia pewnego, zdawałoby się paradoksalnego zjawiska, które zostało zaobserwowane po I wojnie światowej. Po ustaniu działań wojennych, kiedy ludzie na powrót zajęli się uprawianiem swoich zawodów, rybacy odkryli, że populacja ryb drapieżnych w Adriatyku zwiększyła się. Uznano to za paradoks, gdyż zdawałoby się, że wszystkie gatunki powinny ucierpieć w wyniku działań wojennych. Volterra na bazie swojego modelu wykazał, że wzrost liczebności drapieżników był całkiem naturalny, ponieważ w czasie wojny zaprzestano połowów i dzięki temu populacja drapieżników mogła wrócić do stanu naturalnego. Co więcej, model ten odzwierciedla znane w ekologii prawo zachowania średnich, które mówi, że w naturalnych siedliskach zmiany liczebności populacji w czasie zachodzą tak, że zachowana zostaje liczebność średnia.

Zaczniemy od przedstawienia modelu heurystycznego, na bazie którego wyprowadzamy równania modelu Lotki – Volterry. Zakładamy, że w ekosystemie występują dwa gatunki $\mathcal{P}_{d}$ i $\mathcal{P}_{o}$ , przy czym osobniki drugiego gatunku stanowią pożywienie osobników pierwszego gatunku, czyli drapieżników. Jeśli nie ma drapieżników, to gatunek $\mathcal{P}_{o}$ ma bardzo dobre warunki i może się rozwijać, podlegając prawu Malthusa ze współczynnikiem rozrodczości $r>0$ . Natomiast drapieżniki, jeśli nie ma ofiar, to nie mają co jeść, więc giną z głodu. Jeśli w środowisku są osobniki obu gatunków, to drapieżniki polują na ofiary. Zakładamy, że polowanie jest możliwe tylko w przypadku bezpośredniego kontaktu, osobniki poruszają się losowo, zatem liczba kontaktów jest proporcjonalna do liczebności obu gatunków. Zauważmy, że w takiej sytuacji dla pojedynczego drapieżnika mamy odpowiedź funkcjonalną Hollinga typu I. W związku z polowaniem ubywa osobników gatunku $\mathcal{P}_{o}$ , proporcjonalnie do liczby spotkań, a współczynnik proporcjonalności odzwierciedla skuteczność drapieżnika. Po upolowaniu ofiary drapieżnik zyskuje energię, której część przeznacza na rozmnażanie.

Niech $P(t)$ oznacza zagęszczenie drapieżników, a $V(t)$ — zagęszczenie ofiar. Na podstawie powyższego modelu heurystycznego zapiszemy układ równań różniczkowych dynamiki obu populacji

$\begin{array}[]{ccl}\dot{V}&=&rV-aVP,\\ \dot{P}&=&-sP+abVP,\end{array}$

(6.1)

gdzie dla uproszczenia zapisu pomijamy zmienną niezależną $t$ . Poszczególne składniki układu (6.1) mają następującą interpretację

$rV$ i $-sP$ opisują wewnętrzną dynamikę odpowiednio gatunku $\mathcal{P}_{o}$ i $\mathcal{P}_{d}$ , $r$ jest współczynnikiem rozrodczości ofiar, a $s$ — współczynnikiem śmiertelności drapieżników;
$aVP$ odzwierciedla liczbę losowych spotkań osobników obu gatunków, $a$ jest współczynnikiem skuteczności polowania, składnik ten interpretujemy także jako biomasę upolowanych ofiar;
współczynnik $b$ odzwierciedla część upolowanej biomasy, którą gatunek $\mathcal{P}_{d}$ przeznacza na reprodukcję.

Przejdziemy teraz do analizy zachowania rozwiązań układu (6.1). Ponieważ układ (6.1) jest autonomiczny, więc bez straty ogólności przyjmujemy, że $t_{0}=0$ i $V(0)=V_{0}$ , $P(0)=P_{0}$ . Sformułujemy najpierw następujące stwierdzenie.

Stwierdzenie 6.1

Rozwiązania układu (6.1) dla nieujemnego warunku początkowego $(V_{0},P_{0})$ , $V_{0}$ , $P_{0}\geq 0$ , istnieją dla wszystkich $t>0$ , są jednoznaczne i nieujemne.

Ponieważ prawa strona układu (6.1) jest klasy $\mathbf{C}^{1}$ (nawet analityczna) w $\mathbb{R}^{2}$ , to lokalne rozwiązania istnieją dla dowolnego warunku początkowego i są one jednoznaczne. Niech $[0,t^{*})$ będzie przedziałem określoności rozwiązania dla warunku początkowego $(V_{0},P_{0})$ . Pokażemy, że rozwiązanie to jest nieujemne. Niech $V_{0}=0$ . Wtedy $\dot{V}(0)=0$ i zmienna $V$ nie zmienia się. Mamy więc $V(t)\equiv 0$ i wtedy $\dot{P}=-sP$ , czyli $P(t)=P_{0}\text{e}^{{-st}}$ . Analogicznie dla $P_{0}=0$ mamy $P(t)\equiv 0$ i $V(t)=V_{0}\text{e}^{{rt}}$ . Niech teraz $V_{0}$ , $P_{0}>0$ . Załóżmy, że istnieje taka chwila $\bar{t}>0$ , że zmienna $V$ staje się ujemna. Wtedy $V(\bar{t})=0$ , co jest sprzeczne z jednoznacznością rozwiązań, czyli $V(t)>0$ dla wszystkich $t>0$ , dla których rozwiązanie istnieje. Podobnie $P(t)>0$ dla $t>0$ .

Pozostaje jeszcze udowodnić istnienie rozwiązań dla wszystkich $t>0$ . Załóżmy, że rozwiązanie istnieje na pewnym skończonym przedziale $[0,t^{*})$ Wiemy już, że rozwiązanie jest nieujemne, więc

$\dot{V}\leq rV.$

Po pomnożeniu powyższej nierówności przez $\text{e}^{{-rt}}$ dostajemy

$\frac{dV}{dt}\text{e}^{{-rt}}-rV\text{e}^{{-rt}}\leq 0\ \Longrightarrow\ \frac{d}{dt}\left(V\text{e}^{{-rt}}\right)\leq 0,$

czyli $V(t)\text{e}^{{-rt}}\leq V_{0}$ dla $t\geq 0$ , więc $V(t)\leq V_{0}\text{e}^{{rt}}$ . Stąd dla dowolnego skończonego $t^{*}$ zachodzi $V(t^{*})\leq V_{{\max}}=V_{0}\text{e}^{{rt}}$ oraz $\dot{V}(t^{*})\leq rV_{{\max}}$ . Teraz możemy oszacować pochodną zmiennej $P$ . Dla $t\leq t^{*}$ otrzymujemy $\dot{P}\leq abV_{{\max}}P$ , czyli $P(t)\leq P_{0}\text{e}^{{abV_{{\max}}t}}\leq P_{{\max}}=P_{0}\text{e}^{{abV_{{\max}}t^{*}}}$ oraz $\dot{P}\leq abV_{{\max}}P_{{\max}}$ . Szacując pochodne z dołu dostajemy $\dot{V}\geq-aV_{{\max}}P_{{\max}}$ oraz $\dot{P}\geq-sP_{{\max}}$ . Mamy więc funkcje ograniczone o ograniczonych pochodnych, z czego wynika, że istnieją granice $\lim _{{t\to t^{*}}}V(t)=V_{g}$ oraz $\lim _{{t\to t^{*}}}P(t)=P_{g}$ . Gdyby któraś z tych granic nie istniała, to możemy wybrać dwa podciągi zbieżne do różnych granic. Niech np. $\lim _{{n_{1}\to\infty}}V(t_{{n_{1}}})=V_{{g_{1}}}$ oraz $\lim _{{n_{2}\to\infty}}V(t_{{n_{2}}})=V_{{g_{2}}}$ dla $t_{{n_{1}}}$ , $t_{{n_{2}}}\to t^{*}$ i $V_{{g_{1}}}\ne V_{{g_{2}}}$ . Wtedy z twierdzenia o wartości średniej $\dot{V}(t_{{n_{{\xi}}}})=\tfrac{V(t_{{n_{1}}})-V(t_{{n_{2}}})}{t_{{n_{1}}}-t_{{n_{2}}}}\to\infty$ przy $n_{{\xi}}\to\infty$ , co jest sprzeczne z ograniczonością pochodnej. Biorąc warunek początkowy $t_{0}=t^{*}$ , $V(t_{0})=V_{g}$ , $P(t_{0})=P_{g}$ i korzystając z twierdzenia o istnieniu i jednoznaczności rozwiązań przedłużamy rozwiązanie dla $t>t^{*}$ . Wobec dowolności $t^{*}$ rozwiązanie istnieje dla wszystkich $t\geq 0$ .

∎

Zauważmy, że gładkość prawej strony i liniowe oszacowanie pochodnej zawsze gwarantuje przedłużalność rozwiązań.

Przystąpimy teraz do analizy asymptotyki rozwiązań układu (6.1). Wyznaczymy najpierw stany stacjonarne $(\bar{V},\bar{P})$ . Spełniają one układ równań

$\begin{array}[]{ccl}0&=&\bar{V}(r-a\bar{P}),\\ 0&=&\bar{P}(-s+ab\bar{V}),\end{array}$

zatem albo $(\bar{V},\bar{P})=(0,0)$ , albo $(\bar{V},\bar{P})=\left(\tfrac{s}{ab},\tfrac{r}{a}\right)$ . Lokalną stabilność spróbujemy zbadać stosując metodę linearyzacji. Macierz Jacobiego układu (6.1) ma postać

$\left(\begin{array}[]{cc}r-aP&-aV\\ abP&-s+abV\end{array}\right).$

Dla zerowego stanu stacjonarnego otrzymujemy

$\left(\begin{array}[]{cc}r&0\\ 0&-s\end{array}\right),$

widzimy więc, że jest to punkt siodłowy, gdyż wartości własne $\lambda _{1}=r$ , $\lambda _{2}=-s$ mają różne znaki. Z kolei dla dodatniego stanu stacjonarnego

$\left(\begin{array}[]{cc}0&-a\bar{V}\\ ab\bar{P}&0\end{array}\right),$

czyli wartości własne macierzy Jacobiego $\lambda _{{1,2}}=\pm i\sqrt{rs}$ są czysto urojone, zatem nie jest spełnione założenie twierdzenia o linearyzacji dotyczące hiperboliczności układu [13] i nie możemy w tym przypadku skorzystać z tej metody.

Spróbujmy naszkicować portret fazowy. Zaczniemy od wyznaczenia izoklin zerowych, czyli krzywych, na których $\dot{V}=0$ (izoklina zerowa zmiennej $V$ ) albo $\dot{P}=0$ (izoklina zerowa zmiennej $P$ ). Mamy

$\dot{V}=0\Rightarrow V=0\quad\text{albo}\quad P=\frac{r}{a}$

oraz

$\dot{P}=0\Rightarrow P=0\quad\text{albo}\quad V=\frac{s}{ab}.$

Izokliny zerowe dzielą przestrzeń fazową $\left(\mathbb{R}^{+}\right)^{2}$ na cztery podprzestrzenie, w których poszczególne zmienne mają ustalony kierunek przebiegu, por. rys. 6.1. Widzimy, że otrzymane w taki sposób pole wektorowe pokazuje, że rozwiązania okrążają dodatni stan stacjonarny $(\bar{V},\bar{P})$ .

$Szkic portretu fazowego układu~\eqref{LV}.$

Rys. 6.1. Szkic portretu fazowego układu (6.1) z zaznaczonymi izoklinami zerowymi i kierunkami przebiegu poszczególnych zmiennych.

Można pokazać, że wszystkie rozwiązania są okresowe, co wynika z własności całki pierwszej [13] układu (6.1). Znalezienie całki pierwszej i zbadanie jej własności pozostawiamy jako ćwiczenie.

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

Rys. 6.2. Portret fazowy układu (6.1) wraz z przykładowymi krzywymi fazowymi.

Na koniec zauważmy, że wszystkie rozwiązania oscylują wokół stanu stacjonarnego $(\bar{V},\bar{P})$ , por. rys. 6.2. Pokażemy, że każde rozwiązanie ma wartość średnią równą współrzędnym tego stanu. Niech $T$ oznacza okres rozwiązania $(V(t),P(t))$ .

$\frac{dV}{V}=r-aP\ \Longrightarrow\ \text{ln}\frac{V(T)}{V_{0}}=rT-a\int _{0}^{T}P(s)ds,$

a ponieważ rozwiązanie jest okresowe, to $V(T)=V_{0}$ , czyli

$\int _{0}^{T}P(s)ds=\frac{rT}{a}\ \Longrightarrow\ P_{{\text{śr}}}=\frac{1}{T}\int _{0}^{T}P(s)ds=\frac{r}{a}.$

Analogicznie pokazujemy, że $V_{{\text{śr}}}=\tfrac{s}{ab}$ . Rozwiązania oscylują wokół stanu stacjonarnego, przy czym oscylacje populacji drapieżników i ofiar są przesunięte w fazie, por. rys. 6.3.

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

Rys. 6.3. Przebieg przykładowych rozwiązań układu (6.1).

Na tej podstawie wyjaśnimy paradoks związany z populacją ryb drapieżnych w Adriatyku. Zauważmy, że jeśli odławiamy ryby, to zakładając jednakową intensywność odłowu dla obu populacji $\mathcal{P}_{o}$ i $\mathcal{P}_{d}$ , otrzymamy następującą modyfikację układu (6.1)

$\begin{array}[]{ccl}\dot{V}&=&(r-\gamma)V-aVP,\\ \dot{P}&=&-(s+\gamma)P+abVP,\end{array}$

(6.2)

gdzie $\gamma$ oznacza współczynnik odłowu. Zauważmy dalej, że aby populacja ofiar nie wyginęła, musimy dokonywać odłowów w sposób sensowny, tak by $\gamma<r$ . Jeśli łowimy zbyt dużo, $\gamma>r$ , to $\dot{V}<0$ dla wszystkich $t$ i populacja ginie. Natomiast gdy $\gamma<r$ , to układ (6.2) ma dokładnie taką samą postać jak (6.1), więc odławiane populacje mają średnie zagęszczenia równe

$V_{{\text{śr}}}^{{\text{odławiane}}}=\frac{s+\gamma}{ab}>V_{{\text{śr}}},\quad P_{{\text{śr}}}^{{\text{odławiane}}}=\frac{r-\gamma}{a}<P_{{\text{śr}}},$

zatem odławianie powoduje zmniejszenie średniej liczebności populacji drapieżników i zwiększenie średniej liczebności populacji ofiar. W związku z tym w czasie gdy nie było połowów, populacje wróciły do stanu naturalnego, więc liczebność drapieżników wzrosła.

6.2. Konstruktywna krytyka modelu Lotki – Volterry

Chociaż model Lotki – Volterry odzwierciedla obserwowane w naturze oscylacje w układach drapieżnik – ofiara, a także w sposób analityczny dowodzi ekologicznego prawa zachowania średnich w tych układach, to jednak nie jest pozbawiony wad. Przede wszystkim z matematycznego punktu widzenia jego główną wadę stanowi brak stabilności strukturalnej. Pojęcie stabilności strukturalnej wykracza poza ramy tego wykładu — w skrócie stabilność strukturalna oznacza, że niewielka zmiana prawej strony układu nie prowadzi do radykalnych zmian w dynamice rozwiązań. Widzimy, że w przypadku układu (6.1) nawet bardzo mała zmiana prawej strony może skutkować zmianą typu stanu stacjonarnego — stan typu centrum bardzo łatwo zaburzyć. Dodatkowo model ma też pewną własność, która podlega krytyce z ekologicznego punktu widzenia. Rozwiązania układu (6.1) oscylują w taki sposób, że wzrost populacji drapieżników poprzedza wzrost populacji ofiar, natomiast w większości rzeczywistych układów drapieżnik – ofiara takie oscylacje są przesunięte w stosunku do tych rozwiązań; najpierw obserwujemy wzrost populacji ofiar, a potem następujący po nim wzrost populacji drapieżników. Należy też zauważyć, że w modelu heurystycznym pominięto wiele istotnych dla układów drapieżnik – ofiara czynników. W związku z tym model Lotki – Volterry był na wiele różnych sposobów modyfikowany i poniżej omówimy pokrótce kilka z tych modyfikacji.

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

Modelowanie matematyczne w biologii i medycynie

Zagadnienia

6. Modele oddziaływań między dwiema populacjami I

6.1. Model Lotki – Volterry

Stwierdzenie 6.1

6.2. Konstruktywna krytyka modelu Lotki – Volterry