Dlaczego równania różniczkowe to język biologii
Równanie różniczkowe pyta: jak szybko zmienia się dana wielkość? Biologia pełna jest tempa zmian — tempa podziałów komórkowych, tempa drapieżnictwa, tempa otwierania kanałów jonowych, tempa katalizy enzymatycznej. Ilekroć masz tempo zależne od bieżącego stanu układu, masz równanie różniczkowe. Gdy wiele takich temp zależy od siebie nawzajem, masz układ sprzężonych równań różniczkowych — i właśnie wtedy pojawia się interesująca dynamika.
W przeciwieństwie do fizyki, modele biologiczne rzadko wywodzą się wyłącznie z pierwszych zasad. Łączą one ograniczenia wynikające z pierwszych zasad (zachowanie masy, termodynamika) z empirycznymi prawami tempa (kinetyka Hilla, Michaelisa-Menten, wzrost logistyczny), których zasadność potwierdza dopasowanie do danych. Efektem jest pragmatyczna biologia matematyczna: modele użytecznie predykcyjne, choć niekoniecznie ścisłe.
Część 1: Dynamika populacji
Wzrost wykładniczy i równanie logistyczne
Najprostszy model populacji zaczyna się od wzrostu wykładniczego: dN/dt = rN. Każdy osobnik rozmnaża się z tempem r, więc tempo wzrostu jest proporcjonalne do bieżącej liczebności. Daje to N(t) = N₀ e^(rt) — wzrost nieograniczony. Rzeczywiste populacje są ograniczone zasobami, więc Verhulst (1838) dodał człon nasycenia: dN/dt = rN (1 − N/K), gdzie K to pojemność środowiska. Równanie logistyczne ma rozwiązanie w kształcie litery S (sigmoidalne), które wypłaszcza się przy wartości K.
Wzrost logistyczny
dN/dt = rN(1 − N/K)
r = wewnętrzne tempo wzrostu (rok⁻¹)
K = pojemność środowiska (liczba osobników)
Rozwiązanie: N(t) = K / (1 + ((K − N₀)/N₀) · e^(−rt))
Punkty stałe: N* = 0 (niestabilny), N* = K (stabilny)
Maksymalne tempo wzrostu przy N = K/2 (punkt przegięcia)
Czas podwojenia przy niskiej gęstości ≈ ln(2) / r
Lotka-Volterra: dwugatunkowy układ drapieżnik–ofiara
Dodaj drugi gatunek — drapieżnika — a równania staną się dwuwymiarowym układem dynamicznym. Model Lotki-Volterry składa się z dwóch sprzężonych nieliniowych równań różniczkowych. Ofiara rośnie wykładniczo pod nieobecność drapieżników; drapieżniki giną wykładniczo pod nieobecność ofiar. Sprzęga je człon drapieżnictwa: βNP zmniejsza liczebność ofiar, δNP zwiększa liczebność drapieżników.
Układ Lotki-Volterry
dN/dt = αN − βNP (ofiara)
dP/dt = δNP − γP (drapieżnik)
Izokliny zerowe: N* = γ/δ (pionowa), P* = α/β (pozioma)
Punkt stały (N*, P*) to centrum — stabilność neutralna
Trajektorie w przestrzeni fazowej są krzywymi zamkniętymi (wielkość zachowana H)
H = δN − γ ln N + βP − α ln P = const.
Układy rzeczywiste: tłumione spirale zbiegające do stanu równowagi (pojemność środowiska + szum)
Symulacja drapieżnik–ofiara
Model Lotki-Volterry oparty na agentach: króliki żerują i się rozmnażają, lisy polują i głodują. Obserwuj klasyczny przeciwny do ruchu wskazówek zegara cykl w przestrzeni fazowej. Dostosuj tempo narodzin i śmierci, aby zbadać różne reżimy dynamiczne — stabilne cykle, wymieranie, koegzystencję.
Sieć pokarmowa (6 gatunków)
Rozszerzony układ równań różniczkowych z sześcioma sprzężonymi gatunkami: trawa, krzewy, królik, jeleń, lis, wilk. Integracja RK4. Wizualizacja sieciowa, w której rozmiar węzła śledzi liczebność populacji. Wywołuj kaskady troficzne, usuwając gatunki, i obserwuj reakcję całego systemu.
Analiza płaszczyzny fazowej to kluczowe narzędzie. Zamiast wykreślać N(t) i P(t) osobno, wykreśl P względem N. Na tej płaszczyźnie fazowej trajektoria ujawnia strukturę atraktora: czy stan równowagi to stabilna spirala (populacje zbiegają do siebie), centrum (oscylacja neutralna), czy siodło (jedna populacja się załamuje)? Kształt trajektorii mówi wszystko o długoterminowym zachowaniu bez konieczności analitycznego rozwiązywania równania różniczkowego.
Część 2: Neuronauka — model Hodgkina-Huxleya
W 1952 roku Alan Hodgkin i Andrew Huxley opublikowali pierwszy ilościowy model potencjału czynnościowego — impulsu elektrycznego, którym komunikują się neurony. Wprowadzili mikroelektrody do olbrzymiego aksonu kalmara, zaciskali napięcie błony na różnych poziomach i mierzyli prądy. Na podstawie tych eksperymentów wyodrębnili krzywe przewodnictwa dla kanałów sodowych i potasowych, a następnie zapisali cztery sprzężone równania różniczkowe, które z niezwykłą wiernością odtwarzały potencjał czynnościowy. W 1963 roku otrzymali Nagrodę Nobla w dziedzinie fizjologii lub medycyny.
Model Hodgkina-Huxleya to arcydzieło mechanistycznego modelowania biologicznego. Błonę traktuje się jako obwód elektryczny: pojemność błonowa C_m magazynuje ładunek, a przewodnictwa g_Na, g_K i g_L (upływ) przenoszą prąd. Przewodnictwa zależne od napięcia opisują zmienne bramkujące — m, h (sód) i n (potas) — z których każda podlega własnemu równaniu różniczkowemu pierwszego rzędu.
Układ Hodgkina-Huxleya (4 sprzężone równania różniczkowe)
C_m · dV/dt = I_ext − g_Na·m³h·(V − E_Na) − g_K·n⁴·(V − E_K) − g_L·(V − E_L)
dm/dt = α_m(V)·(1−m) − β_m(V)·m
dh/dt = α_h(V)·(1−h) − β_h(V)·h
dn/dt = α_n(V)·(1−n) − β_n(V)·n
V = napięcie błonowe (mV), m = aktywacja Na, h = inaktywacja Na
n = aktywacja K; wszystkie zmienne bramkujące ∈ [0, 1]
g_Na = 120, g_K = 36, g_L = 0,3 mS/cm² (wartości Hodgkina-Huxleya z 1952 r.)
E_Na = +50 mV, E_K = −77 mV, E_L = −54,4 mV (potencjały równowagi)
Neuron Hodgkina-Huxleya
Pełny model neuronu oparty na przewodnictwie: V(t) z prądami Na⁺, K⁺ i upływu. Zmienne bramkujące m, h, n wykreślane w czasie rzeczywistym. Presety: pojedynczy spike, powtarzalne wystrzały, poniżej progu, wysokoczęstotliwościowy wybuch. Integracja Eulera przy dt = 0,01 ms. Tempo wystrzałów i wykrywanie progu na żywo.
Potencjał czynnościowy serca
Potencjał czynnościowy komorowy: pięć faz (spoczynek, narastanie, wczesna szybka repolaryzacja, plateau, końcowa repolaryzacja). Przewodnictwa kanałów Na⁺/Ca²⁺/K⁺, okres refrakcji, automatyzm rozrusznika. Porównaj przebiegi komorowe, przedsionkowe i węzła SA.
Siła modelu Hodgkina-Huxleya tkwi w jego ogólności. Ta sama struktura formalna — pojemnościowa błona, przewodnictwa bramkowane napięciem, kinetyka bramkowania pierwszego rzędu — opisuje neurony różnych gromad zwierząt, komórki serca, komórki beta trzustki i włókna mięśni szkieletowych. Współczesna neuronauka obliczeniowa wykorzystuje rozszerzenia wielokompartmentowe z dziesiątkami typów kanałów jonowych, ale szkielet matematyczny to zawsze Hodgkin-Huxley.
Część 3: Kinetyka enzymatyczna — Michaelis-Menten
Enzymy to katalizatory biologiczne. Wiążą substrat, przekształcają go w produkt i uwalniają produkt — a przy tym same nie ulegają zużyciu. Tempo reakcji katalizowanej enzymatycznie opisali matematycznie Leonor Michaelis i Maud Menten w 1913 roku. Ich model zakłada dwuetapowy mechanizm: wiązanie (odwracalne) i katalizę (nieodwracalną), przy szybkiej równowadze między kompleksem enzym-substrat a wolnym enzymem.
Kinetyka Michaelisa-Menten
Reakcja: E + S ⇌ ES → E + P
Tempo reakcji: v = V_max · [S] / (K_m + [S])
V_max = k_cat · [E_total] (prędkość maksymalna)
K_m = (k₋₁ + k_cat) / k₁ (stała Michaelisa)
Przy [S] = K_m: v = V_max / 2 (tempo połówkowe)
Przy [S] ≪ K_m: v ≈ (V_max / K_m) · [S] (pierwszego rzędu względem S)
Przy [S] ≫ K_m: v ≈ V_max (rzędu zerowego, nasycenie)
Lineweaver-Burk: 1/v = (K_m / V_max) · (1/[S]) + 1/V_max
Kinetyka enzymatyczna
Krzywa tempa Michaelisa-Menten v w funkcji [S], wykres podwójnie odwrotny Lineweavera-Burka, inhibicja kompetycyjna / niekompetycyjna / antykompetycyjna, krzywa wyczerpywania substratu. Suwaki V_max i K_m. Zobacz, jak inhibitory różnie przesuwają pozorne K_m i V_max.
Kinetyka reakcji chemicznych
Stała szybkości Arrheniusa k = A·e^(−Ea/RT), reakcje sekwencyjne A→B→C, wpływ katalizatora na energię aktywacji. Kinetyka pierwszego i drugiego rzędu. Obserwuj, jak krzywe stężenia w funkcji czasu zmieniają się przy regulacji energii aktywacji i temperatury.
Kinetyka Hilla uogólnia model Michaelisa-Menten na układy kooperatywne. Wiele odpowiedzi biologicznych ma kształt nie hiperboliczny, lecz sigmoidalny — opisuje je równanie Hilla v = V_max · [S]ⁿ / (K_d + [S]ⁿ), gdzie n > 1 oznacza wiązanie kooperatywne. Hemoglobina wiąże tlen kooperatywnie (n ≈ 2,8): po związaniu jednej cząsteczki O₂ kolejne wiązanie jest łatwiejsze. Powstaje w ten sposób stroma, przełącznikowa odpowiedź, która czyni hemoglobinę efektywnym nośnikiem tlenu. Ta sama kooperatywna logika pojawia się w wiązaniu czynników transkrypcyjnych, bramkowaniu kanałów jonowych i represji zegara dobowego.
Część 4: Epidemiologia — model SIR
Model SIR dzieli populację na trzy przedziały: podatnych (mogą się zarazić), zakażonych (aktualnie zakaźnych) i ozdrowiałych (odpornych). Kermack i McKendrick opublikowali te ramy w 1927 roku, zainspirowani pandemią grypy z 1918 roku. Mimo swojej prostoty model uchwytuje najważniejszą dynamikę wybuchów chorób zakaźnych — w tym kluczowe pojęcie odporności zbiorowej.
Model przedziałowy SIR
dS/dt = −β · S · I / N
dI/dt = β · S · I / N − γ · I
dR/dt = γ · I
N = S + I + R = const. (populacja zamknięta, bez narodzin/zgonów)
β = tempo transmisji (kontakty × prawdopodobieństwo na kontakt)
γ = tempo powrotu do zdrowia (= 1 / okres zakaźności)
R₀ = β / γ (podstawowa liczba reprodukcji)
Próg epidemiczny: R₀ > 1 (wybuch epidemii)
Próg odporności zbiorowej: frakcja odpornych p_c = 1 − 1/R₀
Równanie rozmiaru końcowego: S_∞ = N · e^(−R₀(1 − S_∞/N))
Model epidemii SIR
Krzywe S/I/R na żywo i siatka agentów 2D. Suwaki dla β (transmisja) i γ (powrót do zdrowia). Obserwuj próg epidemiczny: poniżej R₀ = 1 infekcja wygasa; powyżej — przez populację przetacza się fala epidemii. Włącz szczepienia, by zobaczyć odporność zbiorową w akcji.
Oscylator rytmu dobowego
Oscylator Goodwina generujący 24-godzinne cykle melatoniny, kortyzolu i temperatury. Klasyczny przykład biologicznej pętli sprzężenia zwrotnego opisanej równaniem różniczkowym. Pokazuje trwałą oscylację powstającą z ujemnego sprzężenia zwrotnego z opóźnieniem czasowym — zasadę napędzającą wszystkie zegary dobowe.
Część 5: Solvery równań różniczkowych — poprawna numeryka
Biologiczne równania różniczkowe rzadko dają się rozwiązać analitycznie. Rozwiązujemy je numerycznie: zaczynając od warunków początkowych, krok po kroku posuwamy się w czasie za pomocą reguły przybliżającej. Najprostsza to metoda Eulera: x(t + dt) = x(t) + dt · f(x, t). Łatwo ją zaimplementować, ale szybko kumuluje błąd i może być niestabilna dla układów sztywnych.
Metoda Rungego-Kutty czwartego rzędu (RK4) to koń roboczy symulacji biologicznych. Oblicza prawą stronę równania cztery razy na krok i łączy wyniki z optymalnymi wagami — dając dokładność czwartego rzędu przy umiarkowanym koszcie. Dla układów sztywnych (drastycznie różne skale czasowe, jak w kaskadach enzymatycznych) potrzebne są metody niejawne, takie jak VODE lub numeryczne formuły różniczkowania (NDF), aby uniknąć numerycznego wybuchu przy rozsądnych krokach czasowych.
Krok całkowania RK4
k₁ = f(xₙ, tₙ)
k₂ = f(xₙ + dt/2 · k₁, tₙ + dt/2)
k₃ = f(xₙ + dt/2 · k₂, tₙ + dt/2)
k₄ = f(xₙ + dt · k₃, tₙ + dt)
x_{n+1} = xₙ + (dt/6)·(k₁ + 2k₂ + 2k₃ + k₄)
Lokalny błąd obcięcia: O(dt⁵); błąd globalny: O(dt⁴)
Błąd metody Eulera: O(dt²) na krok, O(dt) globalnie
Układy sztywne: użyj dt ≪ 1/λ_max (największa wartość własna)
Sztywność to ukryty wróg symulacji biologicznych. Kanał sodowy Hodgkina-Huxleya aktywuje się w ~0,2 ms, ale symulacja neuronu może trwać 1000 ms. Stosunek skal czasowych wynosi ~5000. Przy metodzie Eulera dt musi być mniejsze od najszybszej skali czasowej (0,01 ms), co wymaga 100 000 kroków. RK4 przy dt = 0,01 ms jest stabilna i dokładna. Ale dla kaskad enzymatycznych o skalach czasowych sięgających od milisekund po godziny niezbędne stają się adaptacyjne solvery niejawne — w przeciwnym razie symulacje są albo niestabilne, albo nieakceptowalnie wolne.
Wszystko razem: czytanie diagramu fazowego
Płaszczyzna fazowa to podstawowe narzędzie wizualizacji dla dwuwymiarowych układów równań różniczkowych. Wykreśl jedną zmienną względem drugiej; trajektoria układu tworzy krzywą. Stany równowagi pojawiają się jako punkty stałe; ich stabilność wyznaczają wartości własne macierzy Jacobiego w punkcie stałym. Ujemne części rzeczywiste = stabilny (spirala lub węzeł); dodatnie części rzeczywiste = niestabilny; zerowe części rzeczywiste = centrum (neutralny).
Przy więcej niż dwóch zmiennych nie można bezpośrednio zwizualizować pełnej przestrzeni fazowej; zamiast tego używa się przekrojów Poincarégo (płaszczyzny przecinającej trajektorię), diagramów bifurkacyjnych (stan równowagi w funkcji parametru) oraz kontynuacji numerycznej. Narzędzia te pokazują, jak jakościowe zachowanie zmienia się, gdy parametr przekracza punkt bifurkacji: stabilny stan równowagi staje się niestabilny i rodzi cykl graniczny (bifurkacja Hopfa), albo dwa stany równowagi zderzają się i anihilują (bifurkacja siodło-węzeł). Zegar dobowy, bicie serca i potencjał czynnościowy neuronu najlepiej rozumieć właśnie przez ten pryzmat.