Model Lotki – Volterry

Wprowadzimy najpierw pierwszy historycznie model opisujący oddziaływania dwóch popu-lacji 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 biochemicz-nych 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 para-doks, 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ów-nania modelu Lotki – Volterry. Zakładamy, że w ekosystemie występują dwa gatunki P_d i 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 P_o ma bardzo dobre warunki i może się rozwijać, podlegając prawu Malthusa ze współczynnikiem rozrodczości r > 0. Natomiast dra-pież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

Modelowanie matematyczne w biologii i medycynie c U. Foryś, Uniwersytet Warszawski, 2011.

6.1. Model Lotki – Volterry 37

drapieżnika mamy odpowiedź funkcjonalną Hollinga typu I. W związku z polowaniem ubywa osobników gatunku 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 po-pulacji ˙ V = rV − aV P, ˙ P = −sP + abV P, (6.1)

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

— rV i −sP opisują wewnętrzną dynamikę odpowiednio gatunku P_o i P_d, r jest współczynni-kiem rozrodczości ofiar, a s — współczynniwspółczynni-kiem śmiertelności drapieżników;

— aV P odzwierciedla liczbę losowych spotkań osobników obu gatunków, a jest współczynni-kiem skuteczności polowania, składnik ten interpretujemy także jako biomasę upolowanych ofiar;

— współczynnik b odzwierciedla część upolowanej biomasy, którą gatunek 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 i V (0) = V₀, P (0) = P₀. Sformułujemy najpierw następujące stwierdzenie.

Stwierdzenie 6.1. Rozwiązania układu (6.1) dla nieujemnego warunku początkowego (V₀, P0),

V0, P₀ ­ 0, istnieją dla wszystkich t > 0, są jednoznaczne i nieujemne.

Dowód. Ponieważ prawa strona układu (6.1) jest klasy C¹ (nawet analityczna) w R², 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₀, P₀). Pokażemy, że rozwiązanie to jest nieujemne. Niech V₀ = 0. Wtedy ˙V (0) = 0 i zmienna V nie zmienia się.

Mamy więc V (t) ≡ 0 i wtedy ˙P = −sP , czyli P (t) = P₀e^−st. Analogicznie dla P₀ = 0 mamy

P (t) ≡ 0 i V (t) = V₀ert. Niech teraz V₀, P₀ > 0. Załóżmy, że istnieje taka chwila ¯t > 0, że

zmienna V staje się ujemna. Wtedy V (¯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ąza-nie istrozwiąza-nieje na pewnym skończonym przedziale [0, t^∗) Wiemy już, że rozwiązanie jest nieujemne, więc

˙

V ¬ rV.

Po pomnożeniu powyższej nierówności przez e^−rt dostajemy

dV dt^e −rt− rV e^−rt¬ 0 =⇒ ^d dt  V e^−rt¬ 0,

czyli V (t)e^−rt ¬ V₀ dla t ­ 0, więc V (t) ¬ V₀e^rt. Stąd dla dowolnego skończonego t^∗ zachodzi

V (t^∗) ¬ V_max = V₀e^rt oraz ˙V (t^∗) ¬ rV_max. Teraz możemy oszacować pochodną zmiennej P . Dla t ¬ t^∗ otrzymujemy ˙P ¬ abV_maxP , czyli P (t) ¬ P₀eabVmaxt ¬ P_max = P₀eabVmaxt^∗ oraz

˙

P ¬ abVmaxPmax. Szacując pochodne z dołu dostajemy ˙V ­ −aVmaxPmax oraz ˙P ­ −sPmax. Mamy więc funkcje ograniczone o ograniczonych pochodnych, z czego wynika, że istnieją gra-nice lim_t→t∗V (t) = V_g oraz lim_t→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_n1→∞V (tn₁) = V_g1 oraz

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

lim_n2→∞V (tn₂) = V_g2 dla t_n1, t_n2 → t^∗ i V_g1 6= V_g2. Wtedy z twierdzenia o wartości średniej ˙

V (t_nξ) = ^{V (t}n1^{)−V (t}n2⁾

t_n1−t_n2 → ∞ przy n_ξ → ∞, co jest sprzeczne z ograniczonością pochodnej.

Biorąc warunek początkowy t₀ = t^∗, V (t₀) = V_g, P (t₀) = 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 ­ 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 ( ¯V , ¯P ). Spełniają one układ równań

0 = V (r − a ¯^¯ P ),

0 = P (−s + ab ¯^¯ V ),

zatem albo ( ¯V , ¯P ) = (0, 0), albo ( ¯V , ¯P ) = _ab^s,^r_a

. Lokalną stabilność spróbujemy zbadać stosu-jąc metodę linearyzacji. Macierz Jacobiego układu (6.1) ma postać

   r − aP −aV abP −s + abV   .

Dla zerowego stanu stacjonarnego otrzymujemy

   r 0 0 −s   ,

widzimy więc, że jest to punkt siodłowy, gdyż wartości własne λ₁ = r, λ₂ = −s mają różne znaki. Z kolei dla dodatniego stanu stacjonarnego

   0 −a ¯V ab ¯P 0   ,

czyli wartości własne macierzy Jacobiego λ_1,2 = ±i√

rs są czysto urojone, zatem nie jest

speł-nione założenie twierdzenia o linearyzacji dotyczącehiperboliczności układu [14] 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 ˙V = 0 (izoklina zerowa zmiennej V ) albo ˙P = 0 (izoklina zerowa zmiennej P ). Mamy ˙ V = 0 ⇒ V = 0 albo P = ^r a oraz ˙ P = 0 ⇒ P = 0 albo V = ^s ab^.

Izokliny zerowe dzielą przestrzeń fazową R⁺
2

na cztery podprzestrzenie, w których poszcze-gó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 ( ¯V , ¯P ).

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

6.1. Model Lotki – Volterry 39

Rysunek 6.1. Szkic portretu fazowego układu (6.1) z zaznaczonymi izoklinami zerowymi i kie-runkami przebiegu poszczególnych zmiennych.

Rysunek 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 ( ¯V , ¯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)).

dV V ^{= r − aP =⇒ ln} V (T ) V₀ ^{= rT − a} Z T 0 P (s)ds,

a ponieważ rozwiązanie jest okresowe, to V (T ) = V₀, czyli

Z T 0 P (s)ds = ^rT a ^{=⇒ Pśr=} 1 T Z T 0 P (s)ds = ^r a^.

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

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

Rysunek 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 P_o i P_d, otrzymamy następującą modyfikację układu (6.1)

˙

V = (r − γ)V − aV P, ˙

P = −(s + γ)P + abV P,