Notice: Undefined index: mode in /home/misc/mst/public_html/common.php on line 63 Notice: Undefined index: mode in /home/misc/mst/public_html/common.php on line 65 Notice: Undefined index: mode in /home/misc/mst/public_html/common.php on line 67 Notice: Undefined index: mode in /home/misc/mst/public_html/common.php on line 69 Notice: Undefined variable: base in /home/misc/mst/public_html/lecture.php on line 36 Modele matematyczne biologii i medycyny – 6. Modele oddziaływań między dwiema populacjami I – MIM UW

Zagadnienia

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

  1. układ drapieżnik – ofiara, w którym jeden z gatunków ekosystemu jest pożywieniem dla drugiego;

  2. konkurencja, w którym gatunki rywalizują o zasoby środowiska (pożywienie, miejsca lęgowe itp.);

  3. 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.

Projekt współfinansowany przez Unię Europejską w ramach Europejskiego Funduszu Społecznego.

Projekt współfinansowany przez Ministerstwo Nauki i Szkolnictwa Wyższego i przez Uniwersytet Warszawski.