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 – 4. Modele pojedynczej populacji z uwzględnieniem wieku I – MIM UW

Zagadnienia

4. Modele pojedynczej populacji z uwzględnieniem wieku I

Przejdziemy teraz do opisu sytuacji, gdy populacja nie jest jednorodna — rozróżniamy osobniki w różnym wieku.

4.1. Ciąg Fibonacciego

W najprostszym przypadku w populacji wyodrębniamy dwie grupy wiekowe:

  • osobników niedojrzałych (w wieku przed reprodukcyjnym);

  • osobników dojrzałych (w wieku reprodukcyjnym).

Co ciekawe, model opisujący taką sytuację jest najstarszym znanym modelem populacyjnym i znamy go pod nazwą ciągu Fibonacciego.

W I połowie XIII wieku Leonardo z Pizy (Fibonacci, czyli ,,filus Bonacci” — syn Bonacciego) zastosował ten ciąg do opisu następującego zagadnienia populacyjnego.

,,Pewien człowiek wziął parę królików i umieścił je w miejscu otoczonym ze wszystkich stron murem. Ile par królików urodzi się z tej pary w ciągu roku, jeśli założymy, że z każdej pary po miesiącu rodzi się nowa para, która staje się płodna po upływie kolejnego miesiąca?”

Liber abaci rozdział III.

Ilustracja zagadnienia populacyjnego opisanego przez Leonarda z Pizy.
Rys. 4.1. Ilustracja zagadnienia populacyjnego opisanego przez Leonarda z Pizy.

Jako rozwiązanie tego zagadnienia Fibonacci zaproponował ciąg:

1,\; 2,\; 3,\; 5,\; 8,\; 13,\; 21,\; 34,\; 55,\; 89,\; 144\ldots,

znany dziś właśnie jako ciąg Fibonacciego. Ogólny wyraz tego ciągu opisujemy za pomocą formuły rekurencyjnej

F_{{n+1}}=F_{n}+F_{{n-1}},\quad F_{1}=1,\  F_{2}=1. (4.1)

W Liber abaci Fibonacci pominął pierwszy wyraz ciągu — zauważmy, że powinny być dwa wyrazy równe 1, gdyż pierwsza para też staje się płodna dopiero po upływie miesiąca, por. rys. 4.1. Analogiczny schemat opisuje liczbę przodków pojedynczego trutnia w roju pszczół, por. rys. 4.2.

Ilustracja przedstawiająca schemat opisujący liczbę przodków pojedynczego trutnia w roju pszczół.
Rys. 4.2. Ilustracja przedstawiająca schemat opisujący liczbę przodków pojedynczego trutnia w roju pszczół.

Nazwijmy kolejne wyrazy ciągu Fibonacciego liczbami Fibonacciego. Okazuje się, że liczby Fibonacciego opisują różne wielkości przyrodnicze, jak np. liczbę płatków kwiatów (stokrotki stanowią sztandarowy przykład, typowo liczba płatków równa jest 34, 55 albo 89!), liczbę pędów roślin w kolejnych fazach wzrostu, czy też liczbę spiral w różnych konstrukcjach spiralnych (spirale prawoskrętne i lewoskrętne), takich jak kwiatostany słonecznika, owoc ananasa czy szyszki, por. rys. 4.3.

Zdjęcia przedstawiające wielkości przyrodnicze, które są opisywane przez liczby Fibonacciego. Od lewej: liczba płatków stokrotki; liczba pędów roślin; liczba spiral w kwiatostanie słonecznika.
Rys. 4.3. Zdjęcia przedstawiające wielkości przyrodnicze, które są opisywane przez liczby Fibonacciego. Od lewej: liczba płatków stokrotki; liczba pędów roślin; liczba spiral w kwiatostanie słonecznika.

Równanie (4.1) jest liniowym równaniem dyskretnym, przy czym możemy je traktować jako równanie z opóźnionym argumentem — typowo dyskretny układ dynamiczny bez opóźnienia zapisujemy w postaci

X_{{n+1}}=G(X_{n}), (4.2)

gdzie X_{n} opisuje liczebność populacji (czy inną wielkość biofizyczną) w chwili n, a G odzwierciedla prawa rządzące dynamiką tej populacji. W ogólnym przypadku X_{n}\in\mathbb{R}^{n}. Oczywiście postać funkcji G zależy od modelu heurystycznego (np. dla modelu Malthusa (1.3) funkcja G jest liniowa) i w istotny sposób wpływa na przebieg rozwiązań, czyli dynamikę ciągu \left(X_{n}\right)_{{n=1}}^{{\infty}}. Dyskretny układ z opóźnieniem zawsze możemy sprowadzić do układu postaci (4.2) wprowadzając pomocnicze zmienne. Przyjrzyjmy się tej procedurze na przykładzie ciągu Fibonacciego (4.1). Niech Y_{n} oznacza wektor kolejnych liczebności Y_{n}=(X_{n},X_{{n-1}})^{{\text{T}}}. Wtedy

Y_{{n+1}}=\left(\begin{array}[]{cc}1&1\\
1&0\end{array}\right)Y_{n},

zatem jedno równanie z opóźnieniem rzędu jednej jednostki czasu zamienia się na układ dwóch równań postaci (4.2). Liczba równań w układzie wyjściowym zależy od wielkości opóźnień w układzie początkowym.

Jak już przypominaliśmy przy rozwiązywaniu równania (1.3), rozwiązań układu (4.2) dla liniowej funkcji G szukamy w postaci X_{n}=\lambda^{n}X_{0}, gdzie \lambda jest wartością własną macierzy przekształcenia G. Podstawiając taką postać rozwiązania do zależności rekurencyjnej (4.1) otrzymujemy

\lambda^{{n+2}}X_{0}=\lambda^{{n+1}}X_{0}+\lambda^{n}X_{0},

zatem \lambda spełnia równanie

\lambda^{2}-\lambda-1=0, (4.3)

czyli równanie charakterystyczne dla macierzy \left(\begin{array}[]{cc}1&1\\
1&0\end{array}\right). Rozwiązując równanie (4.3) dostajemy

\lambda _{1}=\frac{1+\sqrt{5}}{2},\quad\lambda _{2}=\frac{1-\sqrt{5}}{2}.

Wobec tego wyraz ogólny ciągu Fibonacciego obliczamy jako kombinację liniową

F_{n}=a\lambda _{1}^{n}+b\lambda _{2}^{n}, (4.4)

gdzie a i b musimy wyznaczyć korzystając z warunku początkowego F_{1}=F_{2}=1. Dostajemy zatem układ równań

\begin{array}[]{lclcc}a(1+\sqrt{5})&+&b(1-\sqrt{5})&=&2,\\
a{(1+\sqrt{5})^{2}}&+&b{(1-\sqrt{5})^{2}}&=&4,\\
\end{array}

skąd obliczamy a=\tfrac{\sqrt{5}}{5}=-b i ostatecznie ogólny wyraz ciągu Fibonacciego

F_{n}=\frac{\sqrt{5}}{5}\cdot\frac{(1+\sqrt{5})^{n}-(1-\sqrt{5})^{n}}{2^{n}}. (4.5)

Zauważmy, że ciąg \left(F_{n}\right)_{{n=1}}^{{\infty}} jest rosnący i F_{n}\to\infty przy n\to\infty. Co więcej, tę samą własność mają wszystkie ciągi określone formułą rekurencyjną (4.1) bez względu na warunek początkowy, o ile tylko F_{1}, F_{2}>0. Jeśli wzór (4.5) przepiszemy w postaci

F_{n}=\frac{\sqrt{5}}{5}\left(\frac{1+\sqrt{5}}{2}\right)^{n}\left(1-\left(\frac{1-\sqrt{5}}{1+\sqrt{5}}\right)^{n}\right),

to możemy zauważyć, że wraz z rosnącym n dynamika ciągu \left(F_{n}\right)_{{n=1}}^{{\infty}} jest coraz bliższa dynamice ciągu \left(\frac{\sqrt{5}}{5}\left(\frac{1+\sqrt{5}}{2}\right)^{n}\right)_{{n=1}}^{{\infty}} determinowanej przez dominującą wartość własną \lambda _{1}. Wnioskujemy więc, że wraz z rosnącym n dynamika ta coraz mniej różni się od dynamiki dyskretnego modelu Malthusa (1.3). Podobne zachowanie zaobserwujemy dla większości wartości początkowych F_{1}, F_{2}>0.

Widzimy więc, że dynamika ciągu Fibonacciego jest raczej ,,uboga”. Spójrzmy jednak na ciąg \tfrac{F_{n}}{F_{{n-1}}}. Zdefiniujmy

x_{n}=\frac{F_{n}}{F_{{n-1}}}=\frac{F_{{n-1}}+F_{{n-2}}}{F_{{n-1}}}=1+\frac{1}{\frac{F_{{n-1}}}{F_{{n-2}}}}=1+\frac{1}{x_{{n-1}}}.

Możemy sprawdzić, że ciąg (x_{n})_{{n=1}}^{{\infty}} oscyluje i dzieli się na dwa podciągi (x_{{2k}})_{{k=1}}^{{\infty}}, (x_{{2k+1}})_{{k=1}}^{{\infty}}, monotoniczne, zbieżne do tej samej liczby \tfrac{1+\sqrt{5}}{2}, z czego wynika, że cały ciąg zbiega do tej liczby.

Początkowe wyrazy ciągu $x_{n}=\frac{F_{n}}{F_{{n-1}}}$.
Rys. 4.4. Początkowe wyrazy ciągu x_{n}=\frac{F_{n}}{F_{{n-1}}}.

Liczba ta znana jest od starożytności, oznaczamy ją \Phi i definiuje ona tzw. złotą proporcję lub złoty podział. Nazwa ta wiąże się z podziałem odcinka na dwie części, z których dłuższa do krótszej ma się tak jak cały odcinek do dłuższej. W starożytności mówiono także o ,,boskiej proporcji” i wykorzystywano ją w wielu różnych konstrukcjach, także przy budowie Partenonu. Również w czasach nowożytnych wielu ludzi doszukuje się liczb Fibonacciego czy złotej proporcji w różnych zjawiskach, najczęściej jednak nie ma to prawie żadnych podstaw naukowych.

4.2. Macierze Lesliego

W ogólnym przypadku możemy wyodrębnić więcej grup wiekowych. Niech populacja \mathcal{P} będzie podzielona na k grup wiekowych. W tej sytuacji liczebność całej populacji jest opisana za pomocą wektora \mathbf{N}_{t}=\left(N^{0}_{t},N^{2}_{t},\ldots,N^{k}_{t}\right)^{{\text{T}}}, gdzie N^{i}_{t} oznacza liczebność grupy wiekowej i w chwili t, dla i=0 mamy osobniki nowonarodzone, a dla i=k — osobniki najstarsze, które nie mają szansy dożyć do następnej chwili. Podobnie jak w przypadku modelu Malthusa chcemy na podstawie znajomości wektora liczebności \mathbf{N}_{t} określić \mathbf{N}_{{t+1}}. Jak zwykle musimy zacząć od modelu heurystycznego, czyli omówienia procesów, które są według nas istotne przy opisie dynamiki populacji \mathcal{P}. Ponownie podobnie jak w modelu Malthusa zakładamy, że osobniki są jednorodne w obrębie każdej grupy wiekowej, w każdej grupie wiekowej mamy do czynienia tylko z dwoma procesami — rozrodczości i starzenia się/śmiertelności, przy czym współczynniki rozrodczości i starzenia się są charakterystyczne dla danej grupy wiekowej. Analogicznie jak w przypadku modelu Malthusa w obrębie danej grupy wiekowej oba procesy przebiegają tak samo dla wszystkich osobników i są rozłożone równomiernie w czasie. Zakładamy też, że jednostka czasu jest równa jednostce zmiany wieku (starzenia się), czyli że po upływie np. jednego roku osobnik starzeje się o rok, co więcej po upływie jednej jednostki czasu przechodzi do następnej grupy wiekowej albo umiera. Dostajemy więc:

N_{{t+1}}^{{i+1}}=\gamma _{i}N_{t}^{i},\  i\in\mathbb{N},\quad N_{{t+1}}^{0}=\sum\limits _{{i=0}}^{k}\alpha _{i}N_{t}^{i},

gdzie \gamma _{i} odzwierciedla przeżywalność osobników z grupy wiekowej i (zatem 1-\gamma _{i} to śmiertelność), a \alpha _{i} — rozrodczość danej grupy wiekowej. W związku z tym formalnie model Lesliego można przedstawić jako trójkę (\mathbf{N}_{0},\alpha,\eta), gdzie wektory \mathbf{N}_{0}=(N_{0}^{0},\ldots,N_{0}^{k})^{T}, \alpha=(\alpha _{0},\ldots,\alpha _{k}), \eta=(\eta _{0},\ldots,\eta _{{k-1}}) oznaczają odpowiednio liczebności początkowe, współczynniki rozrodczości oraz śmiertelności, \eta _{i}=1-\gamma _{i}.

Zapisując ten model w postaci (4.2) otrzymujemy

\mathbf{N}_{{t+1}}=\mathbf{M}\mathbf{N}_{t}, (4.6)

gdzie

\mathbf{M}=\left(\begin{array}[]{ccccc}\alpha _{0}&\alpha _{1}&\alpha _{2}&\ldots&\alpha _{k}\\
\gamma _{0}&0&0&\ldots&0\\
0&\gamma _{1}&0&\ldots&0\\
\vdots&0&\ddots&\ddots&\vdots\\
0&\ldots&\ldots&\gamma _{{k-1}}&0\end{array}\right),

przy czym macierz \mathbf{M} nazywamy macierzą Lesliego, a model (4.6) modelem (macierzowym) Lesliego (patrz [8, 9, 2]). Jeśli mamy zadany początkowy rozkład wieku \mathbf{N}_{0}, to rozkład w dowolnej chwili t\in\mathbb{N} możemy obliczyć ze wzoru

\mathbf{N}_{t}=\mathbf{M}^{t}\mathbf{N}_{0}.

Dla dowolnych t, p\in\mathbb{N} zachodzi też \mathbf{N}_{{t+p}}=\mathbf{M}^{p}\mathbf{N}_{t}.

Definicja 4.1

Stanem równowagi modelu Lesliego nazwiemy parę (\bar{\mathbf{N}},r) spełniającą zależność

r\bar{\mathbf{N}}=\mathbf{M}\bar{\mathbf{N}},

gdzie \bar{\mathbf{N}} jest wektorem liczebności, który nazywamy ustaloną strukturą wieku, a r\in\mathbb{R} oznacza stały współczynnik wzrostu.

W przypadku macierzy Lesliego liniowa niezależność wierszy i kolumn pozwala zapisać \mathbf{M} w postaci

\mathbf{M}=\mathbf{V}\Lambda\mathbf{V}^{{-1}},

gdzie \Lambda to diagonalna macierz wartości własnych, a \mathbf{V} macierz kolumnowych wektorów własnych V_{0},\ldots,V_{{k}}. Wynika stąd, że

\mathbf{M}^{t}=(\mathbf{V}\Lambda\mathbf{V}^{{-1}})^{t}=\mathbf{V}\Lambda^{t}\mathbf{V}^{{-1}}\ \Longrightarrow\ \mathbf{N}_{t}=\mathbf{V}\Lambda^{t}\mathbf{V}^{{-1}}\mathbf{N}_{0}.

Oznaczmy W=\mathbf{V}^{{-1}}\mathbf{N}_{0}, W=(w_{0},\ldots,w_{{k}}). Wtedy

\mathbf{N}_{t}=\sum\limits _{{i=0}}^{{k}}V_{i}\lambda _{i}^{t}w_{i}. (4.7)

Własności rozwiązań zależą zarówno od postaci macierzy \mathbf{M} jak i początkowego rozkładu wieku. Twierdzenie Frobeniusa – Perrona gwarantuje istnienie dominującej rzeczywistej wartości własnej \lambda _{0}\geq|\lambda _{i}|, i=1\ldots,k. Dla tej wartości własnej zarówno prawo jak i lewostronny wektor własny jest rzeczywisty i ma nieujemne współrzędne. Jeśli macierz \mathbf{M} spełnia pewne dodatkowe założenia, np. gdy wskaźniki i, dla których \alpha _{i}>0 nie mają większego wspólnego dzielnika niż 1, to \lambda _{0}>|\lambda _{i}|, i=1,\ldots,k i odpowiadający jej wektor własny \mathbf{M}V_{0}=\lambda _{0}V_{0} ma współrzędne nieujemne i dla większości warunków początkowych asymptotyczne zachowanie rozwiązań jest determinowane przez aV\lambda _{0}^{t}, gdzie a=\text{const} jest stałą zależną od warunku początkowego, ponieważ ze wzoru (4.7) dla dominującej wartości własnej \lambda _{0} wynika rozkład asymptotyczny

\frac{1}{\lambda _{0}^{t}}\mathbf{N}_{t}=V_{0}w_{0}+\sum\limits _{{i=1}}^{{k}}V_{i}\left(\frac{\lambda _{i}}{\lambda _{0}}\right)^{t}w_{i}\to V_{0}w_{0}\ \ \text{przy}\ \  t\to\infty.

Wobec tego dla większości warunków początkowych populacja asymptotycznie osiąga rozkład wieku V_{0} przy współczynniku rozrodczości \lambda _{0}, gdyż

\mathbf{N}_{t}\approx w_{0}V_{0}\lambda _{0}^{t}\quad\text{przy}\quad t\to\infty

i jeśli

  1. \lambda _{0}>1, to liczebności wszystkich grup wiekowych rosną nieograniczenie, mówimy że populacja jest rozwojowa;

  2. \lambda _{0}<1, to populacja jest wymierająca;

  3. \lambda _{0}=1, to liczebność całej populacji i poszczególnych grup wiekowych stabilizuje się — populacja jest stacjonarna, por rys. 4.5.

Przykładowe rozwiązania modelu~\eqref{Leslie} w przypadku trzech grup wiekowych w zależności od dominującej wartości własnej $\lambda _{0}$. Wstawki pokazują ewolucję procentowego udziału poszczególnych grup wiekowych w całości populacji.
Rys. 4.5. Przykładowe rozwiązania modelu (4.6) w przypadku trzech grup wiekowych w zależności od dominującej wartości własnej \lambda _{0}. Wstawki pokazują ewolucję procentowego udziału poszczególnych grup wiekowych w całości populacji.

Wektor V_{0} odzwierciedla rozkład liczebności populacji w stanie równowagi, przy czym jeśli ||V_{0}||=1, to współrzędne tego wektora odzwierciedlają procentowy wkład poszczególnych grup wiekowych w populację. W takiej sytuacji dynamika modelu Lesliego jest podobna do dynamiki modelu Malthusa, co więcej, jeśli rozpatrzymy całkowitą liczebność populacji N_{t}=\sum\limits _{{i=0}}^{k}N_{t}^{i}, to N_{t}\sim\lambda _{0}^{t}N_{0}.

Możliwe jest jednak także inne zachowanie rozwiązań, w szczególności cykliczne zmiany struktury wieku. Wiąże się to między innymi z dynamiką determinowaną przez początkowe rozkłady wieku należące do podprzestrzeni generowanych przez wektory własne odpowiadające wartościom własnym \lambda _{1},\ldots,\lambda _{k} w sytuacji, gdy wśród tych wartości własnych nie ma już kolejnej dominującej. Oczywiście zależy to w istotny sposób od postaci macierzy Lesliego. Rozpatrzmy przykład macierzy Lesliego odzwierciedlającej sytuację, gdy rozmnażają się tylko osobniki najstarsze. Rozpatrzmy przykład populacji z wydzielonymi dwoma grupami wiekowymi — rozróżniamy osobniki dojrzałe, zdolne do reprodukcji i osobniki niedojrzałe. Mamy więc

\mathbf{M}=\left(\begin{array}[]{cc}0&m\\
s&0\end{array}\right). (4.8)

Obliczając kolejne potęgi macierzy (4.6) dostajemy

\mathbf{M}^{{2t-1}}=\left(\begin{array}[]{cc}0&m^{t}s^{{t-1}}\\
s^{t}m^{{t-1}}&0\end{array}\right)\quad\text{oraz}\quad\mathbf{M}^{{2t}}=\left(\begin{array}[]{cc}m^{t}s^{{t}}&0\\
0&s^{t}m^{{t}}\end{array}\right),

co łatwo wykazać indukcyjnie. Zauważmy, że jeśli ms=1 i s<1, to \mathbf{N}_{{2t}}=\mathbf{N}_{0} oraz \mathbf{N}_{{2t-1}}=\left(\begin{array}[]{c}mN^{2}_{0}\\
sN^{1}_{0}\end{array}\right), czyli mamy zachowanie cykliczne, por. rys. 4.6.

Cykliczne zachowanie rozwiązań modelu~\eqref{Lieslie} gdy $ms=1$ i $s<1$.
Rys. 4.6. Cykliczne zachowanie rozwiązań modelu (4.8) gdy ms=1 i s<1.

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.