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 – 1. Modelowanie pojedynczej populacji I – MIM UW

Zagadnienia

1. Modelowanie pojedynczej populacji I

W pierwszej części wykładu przedstawimy podstawowe modele opisujące ekosystem, w którym występuje tylko jedna populacja i opiszemy dynamikę takiej populacji w zależności od przyjętego modelu heurystycznego i formalizmu matematycznego.

1.1. Model Malthusa

W najprostszym przypadku zakładamy, że opisywana populacja ma w danym siedlisku bardzo dobre warunki rozwoju, wyrażające się głównie w taki sposób, że każdy osobnik ma nieograniczony dostęp do pożywienia i miejsc lęgowych oraz że obserwujemy jedynie proces rozmnażania. Aby opis był możliwie najprostszy, osobniki w danej populacji powinny być jednakowe i podlegać jednakowym prawom. Co więcej, powinny być równomiernie rozłożone w przestrzeni, żeby nie wystąpiła konieczność analizowania zależności od przestrzeni, a tylko od czasu. Przyjmujemy zatem następujący opis heurystyczny dla populacji \mathcal{P}

  • populacja jest jednorodna, składa się z genetycznie identycznych osobników rozmnażających się partenogenetycznie;

  • osobnik rodzi się w pełni ukształtowany, zdolny do rozrodu i może rozmnażać się w dowolnym wieku;

  • momenty rozmnażania są w dowolnym przedziale czasu rozłożone jednostajnie;

  • każdy osobnik wydaje na świat potomstwo co \tau jednostek czasu, \tau jest ustalone i jednakowe dla wszystkich osobników;

  • każdorazowo jeden rodzic ma \lambda osobników potomnych.

Na podstawie powyższych informacji spróbujemy zaproponować równanie opisujące średnią liczebność populacji N(t) w chwili t. Dokładniej, załóżmy że znamy liczebność w pewnej ustalonej chwili t i chcemy obliczyć N(t+\Delta t) po upływie czasu \Delta t. Zauważmy, że w przedziale czasu (t,t+\Delta t) jest

\frac{\Delta t}{\tau}

momentów rozmnażania. Każdy rodzic ma w tym przedziale czasu średnio

\frac{\Delta t}{\tau}\lambda

potomków. Liczba osobników zdolnych do rozmnażania w chwili t wynosi N(t), zatem przyrost liczebności w krótkim przedziale czasu (t,t+\Delta t) możemy przybliżyć równaniem

N(t+\Delta t)-N(t)=\frac{\Delta t}{\tau}\lambda N(t). (1.1)

Dzieląc stronami przez \Delta t i przechodząc do granicy otrzymujemy

\lim\limits _{{\Delta t\to 0}}\frac{N(t+\Delta t)-N(t)}{\Delta t}=\frac{d}{dt}N(t)=\frac{\lambda}{\tau}N(t).

Przyjmijmy oznaczenia \tfrac{d}{dt}N(t)=\dot{N}(t), \tfrac{\lambda}{\tau}=r i przepiszmy powyższe równanie w postaci

\dot{N}(t)=rN(t), (1.2)

gdzie r nazwiemy współczynnikiem rozrodczości populacji \mathcal{P}. Zbudowaliśmy w ten sposób model ciągły (tzn. zmienna czasowa t zmienia się w sposób ciągły i model heurystyczny jest opisany za pomocą równania różniczkowego). Jeśli natomiast założymy, że \Delta t=1, gdzie 1 oznacza ustaloną jednostkę czasu i zastosujemy oznaczenie N(t)=N_{t}, to równanie (1.1) możemy zapisać jako

N_{{t+1}}=(1+r)N_{t}, (1.3)

otrzymujemy więc model dyskretny (zmienna czasowa zmienia się w sposób dyskretny, najczęściej t\in\mathbb{N}, a liczebność populacji opisana jest za pomocą równania rekurencyjnego).

Obydwa równania (1.2) i (1.3) nazywamy modelem Malthusa, odpowiednio w wersji ciągłej i dyskretnej. Thomas Malthus, angielski ekonomista i demograf, pod koniec XVIII wieku zwrócił uwagę na zbyt szybki przyrost liczebności populacji ludzkiej.

W pracy zatytułowanej ,,An Essay on the Principle of Population” sformułował swoje słynne prawo, które mówi, że liczba ludności wzrasta w tempie geometrycznym, natomiast zasoby żywności w tempie arytmetycznym, co w oczywisty sposób musi doprowadzić do katastrofy. Zauważmy, że równanie dyskretne (1.3) w bezpośredni sposób odzwierciedla fragment prawa Malthusa odnoszący się do liczebności populacji ludzkiej, natomiast równanie (1.2) jest jego ciągłym odpowiednikiem.

Każde z równań (1.2) i (1.3) możemy w prosty sposób rozwiązać, korzystając z ich liniowości. Przypomnijmy, że dla układów liniowych równań różniczkowych poszukujemy rozwiązań bazowych w postaci wykładniczej

N(t)=N_{0}\text{e}^{{\lambda t}},

a dla równań różnicowych — w postaci potęgowej

N_{t}=N_{0}\lambda^{t},

przy czym stałą \lambda nazywamy wartością własną. Oczywiście w przypadku jednej zmiennej równanie (1.2) rozwiązujemy wprost jako równanie o zmiennych rozdzielonych

\int _{{N_{0}}}^{{N(t)}}\frac{dN}{N}=r\int _{0}^{t}ds\Rightarrow\text{ln}\frac{N(t)}{N_{0}}=rt,

czyli N(t)=N_{0}\text{e}^{{rt}}. Natomiast równanie (1.3) opisuje ciąg geometryczny o ilorazie \tfrac{N_{{t+1}}}{N_{t}}, którego poszczególne wyrazy obliczamy ze wzoru ogólnego

N_{t}=N_{0}(1+r)^{t},

co można łatwo udowodnić stosując zasadę indukcji matematycznej.

Jakie są własności rozwiązań obu wersji modelu Malthusa? Widzimy, że zarówno funkcja N_{0}\text{e}^{{rt}} jak i ciąg geometryczny N_{0}(1+r)^{t}, r>0, są rosnące nieograniczenie, por. rys. 1.1. Co więcej, dla ustalonego N_{0}, jeśli współczynnik rozrodczości w modelu ciągłym oznaczymy przez r_{c}, a w dyskretnym — r_{d}, to możemy dobrać r_{d} do r_{c} (albo odwrotnie), tak by rozwiązania dla t\in\mathbb{N} pokrywały się. Mianowicie, chcemy mieć

N_{0}\text{e}^{{r_{c}n}}=N_{0}(1+r_{d})^{n},\quad n\in\mathbb{N},

więc \text{e}^{{r_{c}}}=1+r_{d}, czyli r_{d}=\text{e}^{{r_{c}}}-1, przy czym oczywiście r_{d}>0\Leftrightarrow r_{c}>0.

Porównywanie rozwiązań dyskretnej i ciągłej wersji modelu Malthusa.
Rys. 1.1. Porównanie rozwiązań wersji ciągłej (1.2) oraz dyskretnej (1.3) modelu Malthusa dla różnych wartości parametru r_{c} przy założeniu, że r_{d}=\text{e}^{{r_{c}}}-1.

Z dynamiki modeli (1.2) i (1.3) wynika jasno, że mogą one mieć tylko ograniczony zakres stosowalności. Jedynie hipotetyczna populacja może rozrastać się nieograniczenie, natomiast w rzeczywistości zawsze działają różne mechanizmy ograniczające taki przyrost, w szczególności pojemność siedliska, w którym występuje populacja \mathcal{P} jest zawsze ograniczona i zbyt duża liczba osobników w tym siedlisku się nie zmieści. Pokreślić należy, już sam Malthus zwrócił uwagę na fakt, że ilość żywności przyrasta znacznie wolniej niż liczba ludności, co w konsekwencji musi prowadzić do konkurencji o pożywienie. Zawsze także obok procesu rozrodczości występuje proces śmiertelności.

1.2. Procesy rozrodczości i śmiertelności

Ponieważ uświadomiliśmy sobie, że w naszym modelu dynamiki populacji musimy uwzględnić nie tylko rozrodczość, ale także i inne procesy, w szczególności śmiertelność osobników danej populacji, a chcemy to ponownie zrobić w jak najprostszy sposób, to odpowiedni model heurystyczny budujemy analogicznie jak w przypadku procesu rozrodczości, przyjmując takie same założenia o jednorodności populacji. Wobec tego zamiast równania (1.1) dostaniemy w rezultacie

N(t+\Delta t)-N(t)=r\Delta tN(t)-s\Delta tN(t),

gdzie s nazywamy współczynnikiem śmiertelności. Odzwierciedla on procent osobników umierających w jednostce czasu. Zauważmy, że \tfrac{1}{s} możemy interpretować jako średnią długość życia pojedynczego osobnika.

Ostatecznie z analitycznego punktu widzenia, po uwzględnieniu procesu śmiertelności w modelu Malthusa dostajemy taką samą strukturę matematyczną jak poprzednio. Oznaczmy przez r_{n}=r-s współczynnik rozrodczości ,,netto”, czyli różnicę między współczynnikiem rozrodczości a śmiertelności. Często r_{n} nazywamy współczynnikiem reprodukcji lub współczynnikiem przyrostu naturalnego. Wtedy mamy

\dot{N}(t)=r_{{n}}N(t)\Rightarrow N(t)=N_{0}\text{e}^{{r_{{n}}t}} (1.4)

w modelu ciągłym oraz

N_{{t+1}}=(1+r_{{n}})N_{t}\Rightarrow N_{t}=N_{0}(1+r_{{n}})^{t} (1.5)

w modelu dyskretnym.

Widzimy, że jeśli przyrost naturalny jest dodatni, r_{{n}}>0, to liczebność populacji rośnie, choć oczywiście nieco wolniej niż dla s=0. Natomiast gdy r_{{n}}<0, to charakter rozwiązań obu modeli zmienia się. W modelu ciągłym mamy

\dot{N}(t)<0\quad\text{dla}\quad N(t)>0,

zatem liczebność populacji spada. Rozwiązaniem jest funkcja wykładnicza o ujemnym wykładniku, czyli

\lim\limits _{{t\to\infty}}N(t)=0,

wobec tego asymptotycznie populacja wymiera. Analogicznie w modelu dyskretnym mamy ciąg geometryczny o ilorazie mniejszym niż 1, czyli także N_{t} maleje do zera przy t\to\infty, por. rys. 1.2.

Porównywanie rozwiązań dyskretnej i ciągłej wersji modelu Malthusa z uwzględnieniem śmiertelności.
Rys. 1.2. Porównanie rozwiązań wersji ciągłej (1.4) oraz dyskretnej (1.5) modelu Malthusa z uwzględnieniem śmiertelności dla różnych wartości współczynnika r_{{nc}} przy założeniu, że r_{{nd}}=\text{e}^{{r_{{nc}}}}-1, gdzie r_{{nc}} i r_{{nd}} to współczynniki rozrodczości ,,netto” odpowiednio dla modelu ciągłego i dyskretnego.

1.3. Migracje

Kolejnym krokiem przybliżającym nasz model do rzeczywistości może być uwzględnienie emigracji osobników w związku z ograniczoną pojemnością siedliska. W najprostszych modelach zakłada się dwa typy migracji

  1. migracja stała w czasie;

  2. migracja zależna od zagęszczenia.

Zauważmy, że przy drugim typie migracji w równaniu (1.4) ponownie musimy uwzględnić składnik \alpha N(t), gdzie \alpha oznacza frakcję osobników migrujących przypadających na jednostkę czasu, zatem

\dot{N}(t)=r_{n}N(t)-\alpha N(t)

i taki sposób migracji nie wpływa na ogólną postać rozwiązań modelu — rozwiązaniem jest zawsze funkcja wykładnicza, przy czym jeśli emigracja nie jest zbyt duża, czyli r_{n}>\alpha, to populacja w dalszym ciągu rośnie nieograniczenie, jeśli przyrost naturalny jest dodatni, natomiast zbyt duża emigracja \alpha>r_{n} spowoduje, że siedlisko zostanie opuszczone, ale matematycznie realizuje się to w nieskończonym czasie. Wobec tego dynamika nie różni się od modelu rozrodczości/śmiertelności omawianego poprzednio. Tak samo przedstawia się dynamika wersji dyskretnej.

Zastanówmy się wobec tego jak zmieni się model, jeśli założymy, że emigracja jest stała w czasie i nie zależy od zagęszczenia populacji. Wtedy współczynnik emigracji m odzwierciedla liczbę osobników, które emigrują (zatem ubywa osobników w siedlisku) w jednostce czasu (na jednostkę powierzchni, ponieważ opisujemy zagęszczenie). Zamiast równania (1.4) otrzymamy

\dot{N}(t)=r_{n}N(t)-m,\quad m>0. (1.6)

Zwróćmy uwagę, że z matematycznego punktu widzenia możemy rozważać różne układy znaków parametrów r_{n} i \alpha w równaniu

\dot{N}(t)=r_{n}N(t)+\alpha, (1.7)

odzwierciedlające odpowiednio

  1. r_{n}>0 i \alpha>0 — dodatni przyrost naturalny i imigrację;

  2. r_{n}>0 i \alpha<0 — dodatni przyrost naturalny i emigrację;

  3. r_{n}<0 i \alpha<0 — ujemny przyrost naturalny i emigrację;

  4. r_{n}<0 i \alpha>0 — ujemny przyrost naturalny i imigrację.

Model (1.7) możemy także interpretować w inny sposób. Współczynnik \alpha>0 oznacza wprowadzanie (introdukcję) nowych osobników do siedliska, natomiast \alpha<0 oznacza odławianie osobników. Wyobraźmy sobie zatem staw z rybami i spróbujmy zinterpretować każdy przypadek. Najłatwiejsza wydaje się interpretacja drugiego przypadku — populacja ryb zwiększa się, więc je odławiamy. Przypadek pierwszy jest z tego punktu widzenia najbardziej problematyczny, ale zakładając niewielką wartość r_{n} możemy pewien czas możemy nie zauważać istotnego przyrostu, więc introdukcja w celu podtrzymania gatunku może wydawać się konieczna. Podobnie w trzecim przypadku — przy małym ujemnym przyroście naturalnym możemy przez pewien czas nie zauważać zmniejszania się populacji z tego powodu i eksploatować tę populację przez odłów (czasem zdarza się tak, że zauważamy ujemny przyrost naturalny, ale z różnych powodów w dalszym ciągu eksploatujemy populację). Przypadek czwarty oznacza oczywiście pozytywną ingerencję człowieka — ponieważ przyrost naturalny jest ujemny, więc zarybiamy staw, by utrzymać populację w siedlisku.

Porównywanie rozwiązań ciągłej wersji modelu Malthusa z uwzględnieniem śmiertelności i migracji dla różnych wartości parametrów $\alpha$.
Rys. 1.3. Przykładowe rozwiązania równania (1.7) dla różnych układów znaków parametrów r_{n} oraz \alpha.

Z punktu widzenia interesującego nas teraz zagadnienia ograniczenia pojemności siedliska najważniejszy jest przypadek drugi, zatem przeprowadzimy jego analizę, a pozostałe trzy należy potraktować jako ćwiczenie. Zajmiemy się więc teraz równaniem (1.6) z dodatnimi parametrami r_{n}>0 i m>0. Narysujmy najpierw zależność (N,\dot{N}) w przestrzeni (\mathbb{R}^{+})^{2}, gdzie \mathbb{R} oznacza liczby rzeczywiste nieujemne, por. rys. 1.4.

Zależność $(N,\dot{N})$ w przestrzeni  $(\mathbb{R}^{+})^{2}$, gdzie $\mathbb{R}^{+}$ oznacza liczby rzeczywiste nieujemne (po lewej) oraz odpowiadający tej zależności portret fazowy (po prawej).
Rys. 1.4. Zależność (N,\dot{N}) dla równania (1.6) w przestrzeni (\mathbb{R}^{+})^{2}, gdzie \mathbb{R}^{+} oznacza liczby rzeczywiste nieujemne (po lewej) oraz odpowiadający tej zależności portret fazowy (po prawej).

Zauważmy, że pochodna \dot{N}(t) jest dla N<\tfrac{m}{r_{n}} ujemna, co oznacza, że liczebność maleje, natomiast dla N>\tfrac{m}{r_{n}} mamy \dot{N}(t)>0, czyli liczebność rośnie. Jakie wnioski możemy wysnuć z tej prostej analizy? Musimy zacząć od przypomnienia sobie podstawowych własności rozwiązań równań różniczkowych zwyczajnych (RRZ). W naszym przypadku prawa strona równania (1.6) jest funkcją klasy \mathbf{C}^{1} (nawet analityczną), więc rozwiązania tego równania są jednoznaczne. Co więcej, ponieważ równanie (1.6) jest liniowe, więc dla dowolnego N_{0}\geq 0 (tylko taki warunek początkowy ma sens biologiczny, z analitycznego punktu widzenia możemy rozpatrywać dowolne N_{0}\in\mathbb{R}) rozwiązanie istnieje dla wszystkich t\in\mathbb{R}, choć zwykle interesuje nas przewidywanie dynamiki populacji w przyszłości, więc ograniczamy się do t\geq 0. Zauważmy dalej, że \bar{N}=\tfrac{m}{r_{n}} jest rozwiązaniem stacjonarnym, czyli jeśli N_{0}=\tfrac{m}{r_{n}}, to N(t)=\tfrac{m}{r_{n}} dla dowolnego t, co oznacza, że rozwiązanie nie zmienia się w czasie. Jest to jedyne rozwiązanie stacjonarne. Stąd jeśli N_{0}<\bar{N}, to liczebność populacji maleje i nie ma żadnej ,,bariery”, która mogłaby ten spadek zahamować. Z teorii RRZ wynika, że taką barierę może stanowić tylko rozwiązanie stacjonarne. Faktycznie — załóżmy, że N(t)\to N_{g} przy t\to\infty. Wtedy \dot{N}(t)\to r_{n}N_{g}-m. Ponieważ N_{g}<\bar{N}=\tfrac{m}{r_{n}}, to od pewnego momentu \bar{t} zachodzi \dot{N}(t)<-a<0 z czego wynika, że N(t)\leq N(\bar{t})-a(t-\bar{t}), a to implikuje, że N(t)\to-\infty. Wobec tego nie może zachodzić nierówność N_{g}<\bar{N}. Wnioskujemy więc, że N(t)\to-\infty, zatem istnieje taki moment t_{k}>0, dla którego N(t_{k})=0. Ze względu na liniowość równania (1.6) możemy je łatwo rozwiązać (jak poprzednio jest to równanie o zmiennych rozdzielonych lub też możemy skorzystać z metody uzmienniania stałej)

N(t)=\frac{m}{r_{n}}+\left(N_{0}-\frac{m}{r_{n}}\right)\text{e}^{{r_{n}t}}.

Analityczna postać rozwiązania tylko potwierdza wcześniejsze wnioski płynące z zależności (N,\dot{N}). Zauważmy dalej, że obliczając drugą pochodną dostajemy

\ddot{N}(t)=r_{n}\dot{N},

zatem jeśli funkcja N(t) jest rosnąca, to rośnie w sposób wypukły (jest to w zasadzie wzrost wykładniczy), natomiast gdy jest malejąca — jest także wklęsła, co oczywiście oznacza, że w skończonym czasie t_{k} przekracza oś poziomą (co już wcześniej zauważyliśmy). Znając postać rozwiązania możemy wyznaczyć t_{k} jako taki czas, dla którego liczebność populacji zeruje się. Skoro N(t_{k})=0, to

0=\frac{m}{r_{n}}+\left(N_{0}-\frac{m}{r_{n}}\right)\text{e}^{{r_{n}t_{k}}}\Longrightarrow t_{k}=\frac{1}{r_{n}}\text{ln}\frac{m}{m-r_{n}N_{0}},\quad\text{o ile }\  N_{0}<\frac{m}{r_{n}}.

Przykładowe rozwiązania tego modelu zostały przedstawione na prawym górnym wykresie na rys. 1.3. Widzimy, że w zależności od wielkości emigracji populacja albo rośnie nieograniczenie, jak w przypadku bazowego modelu Malthusa, albo wymiera w skończonym czasie, jeśli zbyt dużo osobników emigruje. Pozostałe przypadki migracji analizujemy w analogiczny sposób, pozostawiamy tę analizę i płynące z niej wnioski jako ćwiczenie.

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.