Ekologia i dynamika populacji — cykle drapieżnik-ofiara, sieci pokarmowe, kaskady troficzne i cykl węglowy

Matematyka żywych systemów

Modele ekologiczne zajmują niezwykłą pozycję w nauce: są wystarczająco proste, by zapisać je w jednej linijce, a mimo to generują zachowanie — oscylacje, chaos, wymieranie, koegzystencję — które co roku zaskakuje ekologów. Równania drapieżnik-ofiara Lotki-Volterry pochodzą z lat 1925 i 1926. Sto lat później nadal pozostają fundamentem ekologii populacyjnej, rozszerzane i udoskonalane, ale nigdy niezastąpione.

To, co łączy symulacje w tej kolekcji, to wspólny język matematyczny: sprzężone nieliniowe równania różniczkowe zwyczajne (ODE) na zmiennych populacyjnych. Wilki i łosie, plankton i kryl, drzewa i atmosferyczne CO₂ — wszystkie są opisane tą samą formą równania, różniącą się jedynie liczbą oddziałujących zmiennych i strukturą wyrazów interakcji.

Symulacja 1: Dynamika drapieżnik-ofiara

Lotka-Volterra — oscylacje bez zegara

Alfred Lotka (1925) i Vito Volterra (1926) wyprowadzili niezależnie tę samą parę równań. Volterra próbował wyjaśnić, dlaczego odsetek ryb drapieżnych na Morzu Adriatyckim rósł w czasie I wojny światowej — gdy zawieszono połowy, wzrosła zarówno populacja drapieżników, jak i ofiar, ale nieproporcjonalnie. Równania wyjaśniły dlaczego: oba gatunki są uwięzione w cyklu rządzonym przez wielkość zachowaną, analogiczną do energii w mechanice klasycznej. Przerwij zewnętrzne zaburzenie (połowy), a cykl wznawia się.

Zmienne stanu:
  N = populacja ofiar (króliki, ryby, jelenie, ...)
  P = populacja drapieżników (lisy, rekiny, wilki, ...)

Równania różniczkowe:
  dN/dt = α·N  - β·N·P     (ofiara: rośnie, jest zjadana)
  dP/dt = δ·N·P - γ·P      (drapieżnik: zyskuje z ofiar, umiera)

Parametry:
  α = tempo rozrodu ofiar         [dzień⁻¹]
  β = tempo drapieżnictwa         [dzień⁻¹ na drapieżnika]
  γ = tempo śmiertelności drapieżnika [dzień⁻¹]
  δ = efektywność (ofiara→drapieżnik) [bezwymiarowa konwersja]

Punkty stałe:
  Trywialny: (N, P) = (0, 0)   — wymarcie (niestabilny)
  Koegzystencja: (N*, P*) = (γ/δ, α/β)   — centrum neutralne

Wielkość zachowana (funkcja Lapunowa):
  V(N,P) = δN - γ·ln N + βP - α·ln P = const wzdłuż trajektorii
  → zamknięte orbity na płaszczyźnie fazowej (neutralnie stabilne cykle)

Podsumowanie: Układ Lotki-Volterry jest „strukturalnie niestabilny" —
każda mała modyfikacja (zagęszczenie ofiar, nasycające drapieżnictwo)
zamienia neutralne cykle w spirale zbieżne lub rozbieżne. Prawdziwe
układy potrzebują odpowiedzi funkcjonalnych typu Hollinga II lub III
dla stabilności.

Symulacja Lotki-Volterry wyświetla jednocześnie szereg czasowy (populacja vs. czas) i portret fazowy (populacja ofiar na osi x, drapieżnika na osi y). Widać zamknięte orbity w przestrzeni fazowej — matematyczny podpis oscylacji zachowawczych. Przesuń warunki początkowe wokół punktu stałego i obserwuj, jak zmienia się amplituda, podczas gdy okres ledwo się przesuwa: cecha charakterystyczna równań LV, którą zniszczyłaby jakakolwiek modyfikacja z zależnością gęstościową w którymkolwiek z gatunków.

Symulacja 2: Dynamika sieci pokarmowej

Sieć pokarmowa — wielogatunkowe sieci interakcji

Rzeczywiste ekosystemy zawierają nie dwa gatunki, lecz setki, połączone złożoną siecią relacji kto-kogo-zjada. Sieć pokarmowa to graf skierowany, w którym każdy węzeł to gatunek (lub gildia troficzna), a każda krawędź reprezentuje przepływ energii. Struktura sieci — jej łączność, poziomy troficzne oraz rozkład siły interakcji — decyduje o tym, czy ekosystem jest stabilny, czy kruchy.

Dla n gatunków z wektorem populacji x = (x₁, x₂, ..., xₙ):

  dxᵢ/dt = xᵢ (rᵢ + Σⱼ aᵢⱼ·xⱼ)

  rᵢ  = wewnętrzne tempo wzrostu (+ dla producentów, - dla drapieżników szczytowych)
  aᵢⱼ = macierz współczynników interakcji:
         aᵢⱼ > 0: gatunek j przynosi korzyść gatunkowi i (ofiara zjadana przez i)
         aᵢⱼ < 0: gatunek j szkodzi gatunkowi i (konkurent, drapieżnik)
         aᵢⱼ = 0: brak bezpośredniej interakcji

Kryterium stabilności (May 1972):
  Losowy ekosystem z n gatunkami, średnią siłą interakcji s,
  i łącznością C jest stabilny, jeśli:
         s · √(n·C) < 1
  (stabilność maleje wraz ze wzrostem n, C lub s)

Ten słynny wynik sugerował, że różnorodne, silnie połączone
ekosystemy powinny być niestabilne — sprzecznie z intuicją terenową
ekologów. Rozwiązanie: prawdziwe sieci pokarmowe NIE są losowe;
mają specyficzne struktury (słabe powiązania, wszystkożerność, pętle),
które nadają im stabilność.

Poziomy troficzne:
  L1: Producenci pierwotni     (rośliny, fitoplankton) — rᵢ > 0
  L2: Konsumenci pierwotni     (roślinożercy)
  L3: Konsumenci wtórni        (drapieżniki roślinożerców)
  L4: Drapieżniki szczytowe    — rᵢ < 0, całkowicie zależne od ofiar
  Detrytofagi: rozkładają martwą materię → recykling do L1

Symulacja sieci pokarmowej buduje wielogatunkowy ekosystem z maksymalnie 8 gatunkami na 3–4 poziomach troficznych. Usuń gatunek jednym kliknięciem i obserwuj kaskadę: usunięcie drapieżnika szczytowego może spowodować „uwolnienie mezodrapieżnika", gdzie drapieżniki średniego szczebla eksplodują liczebnie. Dodaj dopływ składników odżywczych do producenta i obserwuj, jak wzbogacenie destabilizuje sieć. Symulacja oblicza wartości własne macierzy społeczności w czasie rzeczywistym i wyświetla wskaźnik stabilności.

Symulacja 3: Kaskady troficzne

Kaskada troficzna — kontrola odgórna i gatunki kluczowe

Kaskada troficzna zachodzi, gdy drapieżnik pośrednio przynosi korzyść roślinom, tłumiąc roślinożerców. Ponowne wprowadzenie wilków w Yellowstone to najsłynniejszy przykład lądowy; wydry morskie, „wilki lasów wodorostów", utrzymujące jeżowce w ryzach, są morskim odpowiednikiem. Zrozumienie kaskad troficznych ma kluczowe znaczenie dla restytucji ekologicznej i biologii konserwatorskiej.

Trzy gatunki: Roślina (V), Roślinożerca (H), Mięsożerca (C)

dV/dt = r·V·(1 - V/K)  - f(V)·H
dH/dt = e₁·f(V)·H      - g(H)·C  - d_H·H
dC/dt = e₂·g(H)·C      - d_C·C

Odpowiedzi funkcjonalne:
  f(V) = a·V / (1 + a·h·V)   — Holling typu II (nasycające drapieżnictwo)
  g(H) = b·H / (1 + b·k·H)   — Holling typu II dla mięsożerców
  a, b = tempo ataku
  h, k = czas obsługi (max tempo pobierania pokarmu = 1/h)

Mechanizm kaskady:
  Wzrost C → spadek H → wzrost V   (dodatni efekt pośredni)
  To łańcuch ujemny-dodatni-ujemny → efekt netto dodatni (kaskada)

Kryterium gatunku kluczowego:
  Gatunek jest „kluczowy", jeśli jego usunięcie ma nieproporcjonalnie
  duży wpływ w stosunku do jego biomasy (Paine 1969).
  Gatunki kluczowe to zazwyczaj drapieżniki szczytowe o niskiej biomasie,
  ale wysokiej sile interakcji per capita.

Tłumienie behawioralne („krajobraz strachu"):
  Drapieżniki nie muszą zjadać ofiar, by je kontrolować.
  Strach przed drapieżnictwem ogranicza miejsca żerowania roślinożerców,
  redukując wpływ wypasu w obszarach uczęszczanych przez drapieżniki.
  Ten „krajobraz strachu" podwaja realizowany efekt odgórny.

Symulacja 4: Zmierzchowa strefa oceanu

Strefa afotyczna — życie bez światła

Poniżej 200 metrów pochłaniany jest ostatni foton światła słonecznego. Fotosynteza nie jest możliwa. Mimo to głęboki ocean nie jest pozbawiony życia: „pompa biologiczna" — opadanie materii organicznej z oświetlonej strefy fotycznej — dostarcza strumień energii, który podtrzymuje jeden z największych i najmniej zbadanych ekosystemów Ziemi. Śnieg morski, grudki kałowe, martwe organizmy i migrujący zooplankton przenoszą węgiel z powierzchni do głębin, sekwestrując CO₂ w skalach czasowych wieków.

Strumień eksportu węgla (prawo potęgowe Martina 1987):
  F(z) = F₀ · (z / z₀)^(-b)
  F₀ = strumień eksportu przy głębokości referencyjnej z₀ (zwykle 100 m)
  b  ≈ 0,86 (średnia globalna, waha się 0,4–1,4)
  z  = głębokość [m]

Przykład:
  F(100m)  = 10 g C/m²/rok  → eksport powierzchniowy
  F(500m)  = 10 × (5)^-0,86 ≈ 2,3 g C/m²/rok
  F(4000m) = 10 × (40)^-0,86 ≈ 0,5 g C/m²/rok  (tylko 5% dociera do dna)

Organizmy na głębokości:
  Strefa mezopelagiczna (200–1000 m):
    Myktofidy, kałamarnice, widłonogi — dobowa migracja pionowa
    Codziennie pokonują 200-300 m, by żerować na powierzchni i wracać o świcie
    Ta migracja transportuje ~1–5 Gt C rocznie globalnie

  Strefa batypelagiczna (1000–4000 m):
    Żabnice, viperfish — chemoreceptory, bioluminescencja
    Brak światła słonecznego → brak fotosyntezy → cała energia ze śniegu morskiego

Koszt energetyczny dobowej migracji pionowej:
  Widłonóg migrujący 200 m/noc wydatkuje ~15% dziennego budżetu energii
  na pływanie, ale zyskuje ~30% więcej dostępu do pokarmu — netto dodatni

Symulacja 5: Globalny cykl węglowy

Cykl węglowy — najważniejszy przepływ ekologiczny Ziemi

Węgiel krąży przez pięć głównych rezerwuarów: atmosferę, biosferę lądową (rośliny i gleby), powierzchniową warstwę oceanu, głęboki ocean i litosferę (paliwa kopalne i skały węglanowe). Naturalne strumienie między tymi rezerwuarami były w przybliżeniu zrównoważone przed erą przemysłową. Od 1750 roku spalanie paliw kopalnych i zmiana użytkowania gruntów dodały ~640 Gt węgla do atmosfery — około połowa została pochłonięta przez ocean i biosferę lądową; reszta się zgromadziła, podnosząc atmosferyczne CO₂ z 280 do ponad 420 ppm.

Rezerwuary (Gt C, przybliżone wartości z 2024):
  Atmosfera:            850 Gt C     (CO₂: ~420 ppm)
  Rośliny lądowe:       550 Gt C
  Gleba + detrytus:    2500 Gt C
  Powierzchnia oceanu:   900 Gt C
  Głęboki ocean:      38 000 Gt C
  Paliwa kopalne:     ~1200 Gt C pozostałych do wydobycia

Roczne strumienie (naturalne, przybliżone):
  GPP (produkcja pierwotna brutto): ~120 Gt C/rok   (fotosynteza)
  Oddychanie ekosystemu:            ~118 Gt C/rok   (z powrotem do atm.)
  Pochłanianie oceaniczne (powietrze-morze): ~2,3 Gt C/rok  (netto do oceanu)
  Wietrzenie do oceanu:               ~0,4 Gt C/rok

Zaburzenie antropogeniczne (2023):
  Paliwa kopalne + przemysł:     +9,7 Gt C/rok
  Zmiana użytkowania gruntów:    +1,1 Gt C/rok
  Emisje całkowite:              +10,8 Gt C/rok
  Akumulacja atmosferyczna:      +5,1 Gt C/rok  (rezydualna)
  Pochłanianie oceaniczne:       +2,9 Gt C/rok
  Pochłanianie lądowe:           +2,8 Gt C/rok

Frakcja przenoszona w powietrzu CO₂ ≈ 47–50% (reszta pochłonięta przez ocean + ląd)
Zakwaszenie oceanu:  pH ≈ 8,2 (1850) → 8,1 (obecnie) → ~7,9 (2100 przy BAU)

Symulacja 6: Stada i zachowanie zbiorowe

Stado ptaków — emergentna struktura z lokalnych reguł

Mrowienie się szpaków — tysiące ptaków poruszających się w płynnych, skoordynowanych falach — nie ma dyrygenta. Każdy ptak przestrzega trzech lokalnych reguł (separacja, wyrównanie, kohezja), które razem tworzą skoordynowany ruch grupowy odporny na zaburzenia. Z perspektywy ekologicznej zachowanie stadne wyewoluowało głównie w celu unikania drapieżników: stado 10 000 ptaków jest 10 000 razy trudniejsze do zasadzki niż pojedynczy ptak, a efekt dezorientacji sprawia, że celowanie w konkretną jednostkę jest niemal niemożliwe.

Standardowe Boids (Reynolds 1987):
  f_sep = -Σⱼ (pⱼ - pᵢ) / |pⱼ - pᵢ|²   [separacja]
  f_ali = (Σⱼ vⱼ / |N|) - vᵢ            [wyrównanie]
  f_coh = (Σⱼ pⱼ / |N|) - pᵢ - vᵢ      [kohezja]

Unikanie drapieżnika (4. reguła):
  f_pred = -k_p · (pᵢ - p_pred) / |pᵢ - p_pred|³   jeśli |pᵢ - p_pred| < R_alert
  Odpychanie odwrotnie proporcjonalne do sześcianu daje ostre, nagłe manewry ucieczki

Zbiorowe zachowania antydrapieżnicze:
  Wakuola: stado dzieli się, pozostawia puste miejsce wokół drapieżnika
  Błyskawiczna ekspansja: wszystkie ptaki jednocześnie uciekają na zewnątrz
  Skoordynowany zwrot: fala zwrotu propaguje się z ~15 ms/ptak (dane szpaków)
  Efekt dezorientacji: czas namierzenia drapieżnika ∝ n (rozmiar stada)

Propagacja informacji:
  Reakcja „przestrachu" u jednego ptaka propaguje się do sąsiadów w ~13 ms.
  Przy 10 najbliższych interakcjach sąsiedzkich, stado 1000 ptaków aktualizuje się w
  ~125 ms — wystarczająco szybko, by wydawać się natychmiastowe dla drapieżnika.

Samolubne stado Hamiltona (1971):
  Każde zwierzę minimalizuje swoją „domenę zagrożenia" (obszar bliższy jemu
  niż jakiemukolwiek innemu zwierzęciu). To skłania zwierzęta do przemieszczania się
  ku centrum grupy — tworząc agregację bez siły kohezji.

Symulacja 7: Wzrost upraw i agronomia

Wzrost upraw — od nasiona do zbiorów poprzez dni stopniowzrostu

Rolnictwo to stosowana ekologia: zarządzanie dynamiką populacji roślin uprawnych w celu maksymalizacji produkcji żywności. Modele wzrostu upraw — od prostych liczników dni stopniowzrostu po pełne modele procesowe, takie jak APSIM i DSSAT — przewidują etapy fenologiczne (kiełkowanie, kwitnienie, dojrzałość), akumulację biomasy i plon na podstawie danych pogodowych. Są używane wszędzie: do optymalizacji terminów siewu, szacowania plonów w scenariuszach klimatycznych i planowania nawadniania.

Zgromadzone dni stopniowzrostu (GDD):
  GDD_dzień = max(T_baza, (T_max + T_min)/2) - T_baza
  GDD_skumulowane += GDD_dzień jeśli > 0

Progi fenologiczne (pszenica, przybliżone):
  Kiełkowanie:    ~150 GDD₅  (baza 5°C)
  Krzewienie:     ~450 GDD₅
  Kłoszenie:     ~1300 GDD₅
  Dojrzałość fizjologiczna: ~2000 GDD₅

Akumulacja biomasy (Monteith 1977):
  ΔB = ε_c · PAR_przechwycone · FPAR
  ε_c = efektywność wykorzystania promieniowania [g/MJ]  (pszenica ≈ 1,2 g/MJ)
  FPAR = ułamek PAR przechwycony przez koronę = 1 - e^(-k·LAI)
  k    = współczynnik ekstynkcji ≈ 0,5
  LAI  = wskaźnik powierzchni liści [m² liścia / m² gruntu]

Plon:
  Y = B_całkowita · HI
  HI = indeks zbiorczy (ziarno / całkowita biomasa) ≈ 0,45 pszenica, 0,50 kukurydza

Czynnik stresu wodnego:
  ks = min(1, ETa / ETc)
  Mnoży ΔB: ograniczona woda → ograniczony wzrost
  ETa = rzeczywista ewapotranspiracja, ETc = potencjalna

Struktura wszystkich tych modeli

Patrząc na wszystkie siedem symulacji, wyłania się wzorzec. Każdy model, od dwugatunkowej Lotki-Volterry po pięciorezerwuarowy cykl węglowy, ma tę samą formę: wektor zmiennych stanu x i prawą stronę F(x, parametry):

dx/dt = F(x, p)

Modele różnią się tym, co wchodzi do x (populacje gatunków, rozmiary rezerwuarów węgla, biomasa upraw), jaką strukturę ma F (wyrazy bilinearne dla drapieżnictwa, nasycające wyrazy Hollinga-II dla sytości drapieżnika, wyrazy logistyczne dla konkurencji zależnej od gęstości) oraz co oznaczają parametry. Ale zadanie matematyczne — całkowanie tych ODE w czasie, znajdowanie punktów stałych, ocena stabilności przez wartości własne jakobianu, badanie bifurkacji przy zmianie parametrów — jest zawsze takie samo.

Pełny katalog symulacji ekologicznych — przeglądaj wszystkie symulacje ekologiczne — obejmuje również modele afotycznej strefy globalnego oceanu, model erozji gleby RUSLE oraz model wzrostu roślin GDD. Spotlight #17 o ekologii i biologii obejmuje wcześniejsze symulacje w kolekcji. Dla leżącej u podstaw matematyki Learning #16 o ODE w biologii przeprowadza przez każdy typ równania wraz z wyprowadzeniami.