Чому ОДР — це мова біології
Диференціальне рівняння запитує: як швидко змінюється ця величина? Біологія повна швидкостей — швидкостей поділу клітин, швидкостей хижацтва, швидкостей відкриття іонних каналів, швидкостей ферментативного каталізу. Щоразу, коли швидкість залежить від поточного стану системи, ви маєте справу з диференціальним рівнянням. Коли кілька швидкостей залежать одна від одної, ви отримуєте систему пов'язаних ОДР — і саме тут виникає цікава динаміка.
На відміну від фізики, біологічні моделі рідко виводяться виключно з перших принципів. Вони поєднують обмеження з перших принципів (збереження маси, термодинаміка) з емпіричними законами швидкості (кінетика Хілла, Міхаеліса-Ментен, логістичне зростання), які обґрунтовуються відповідністю даним. Результат — прагматична математична біологія: моделі, які корисно прогнозують, не будучи точними.
Частина 1: динаміка популяцій
Експоненційне зростання та логістичне рівняння
Найпростіша модель популяції починається з експоненційного зростання: dN/dt = rN. Кожна особина розмножується зі швидкістю r, тож швидкість зростання пропорційна поточному розміру. Це дає N(t) = N₀ e^(rt) — необмежене зростання. Реальні популяції обмежені ресурсами, тому Ферхюльст (1838) додав член насичення: dN/dt = rN (1 − N/K), де K — ємність середовища. Логістичне рівняння має S-подібний (сигмоїдний) розв'язок, що виходить на плато при K.
Логістичне зростання
dN/dt = rN(1 − N/K)
r = внутрішня швидкість зростання (рік⁻¹)
K = ємність середовища (особин)
Розв'язок: N(t) = K / (1 + ((K − N₀)/N₀) · e^(−rt))
Нерухомі точки: N* = 0 (нестійка), N* = K (стійка)
Максимальна швидкість зростання при N = K/2 (точка перегину)
Час подвоєння при низькій густині ≈ ln(2) / r
Лотка-Вольтерра: система «хижак–жертва» для двох видів
Додайте другий вид — хижака — і рівняння перетворюються на 2D динамічну систему. Модель Лотка-Вольтерра складається з двох пов'язаних нелінійних ОДР. Жертва зростає експоненційно за відсутності хижаків; хижаки вимирають експоненційно за відсутності жертви. Члени хижацтва пов'язують їх: βNP зменшує популяцію жертви, δNP збільшує популяцію хижаків.
Система Лотка-Вольтерра
dN/dt = αN − βNP (жертва)
dP/dt = δNP − γP (хижак)
Нуль-ізокліни: N* = γ/δ (вертикальна), P* = α/β (горизонтальна)
Нерухома точка (N*, P*) — це центр, нейтральна стійкість
Траєкторії у фазовому просторі — замкнені криві (величина H зберігається)
H = δN − γ ln N + βP − α ln P = const.
Реальні системи: затухаючі спіралі до рівноваги (ємність середовища + шум)
Симуляція «хижак–жертва»
Агентна модель Лотка-Вольтерра: кролики пасуться й розмножуються, лисиці полюють і голодують. Спостерігайте класичний контрзгодинниковий цикл у фазовому просторі. Змінюйте швидкості народжуваності й смертності, щоб дослідити різні динамічні режими — стійкі цикли, вимирання, співіснування.
Харчова мережа (6 видів)
Розширена система ОДР із шістьма пов'язаними видами: трава, кущі, кролик, олень, лисиця, вовк. Інтегрування RK4. Мережева візуалізація, де розмір вузла відображає популяцію. Запустіть трофічний каскад, видаливши вид, і спостерігайте реакцію всієї системи.
Аналіз фазової площини — ключовий інструмент. Замість того щоб будувати N(t) і P(t) окремо, побудуйте P відносно N. На цій фазовій площині траєкторія розкриває структуру атрактора: чи є рівновага стійкою спіраллю (популяції сходяться), центром (нейтральні коливання), чи сідлом (одна популяція колапсує)? Форма траєкторії розповідає все про довгострокову поведінку без аналітичного розв'язання ОДР.
Частина 2: нейронаука — модель Ходжкіна-Хакслі
У 1952 році Алан Ходжкін та Ендрю Хакслі опублікували першу кількісну модель потенціалу дії — електричного імпульсу, яким нейрони обмінюються сигналами. Вони вводили мікроелектроди в гігантський аксон кальмара, фіксували напругу мембрани на різних рівнях та вимірювали струми. З цих експериментів вони вивели криві провідності для натрієвих і калієвих каналів, а потім записали чотири пов'язані ОДР, які з дивовижною точністю відтворювали потенціал дії. У 1963 році вони отримали Нобелівську премію з фізіології або медицини.
Модель Ходжкіна-Хакслі — шедевр механістичного біологічного моделювання. Мембрана розглядається як електричний контур: ємність мембрани C_m накопичує заряд, а провідності g_Na, g_K і g_L (витік) переносять струм. Напруго-залежні провідності описуються керуючими змінними — m, h (натрій) та n (калій), кожна з яких підпорядковується власному ОДР першого порядку.
Система Ходжкіна-Хакслі (4 пов'язані ОДР)
C_m · dV/dt = I_зовн − 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 = напруга мембрани (мВ), m = активація Na, h = інактивація Na
n = активація K; усі керуючі змінні ∈ [0, 1]
g_Na = 120, g_K = 36, g_L = 0,3 мС/см² (значення Ходжкіна-Хакслі, 1952)
E_Na = +50 мВ, E_K = −77 мВ, E_L = −54,4 мВ (потенціали реверсії)
Нейрон Ходжкіна-Хакслі
Повна модель нейрона на основі провідності: V(t) зі струмами Na⁺, K⁺, витоку. Керуючі змінні m, h, n у реальному часі. Готові сценарії: одиночний спайк, повторюваний розряд, підпороговий режим, високочастотна серія. Інтегрування Ейлера при dt = 0,01 мс. Частота спайків і виявлення порогу в реальному часі.
Потенціал дії серця
Потенціал дії шлуночка: п'ять фаз (спокій, наростання, рання швидка реполяризація, плато, кінцева реполяризація). Провідності каналів Na⁺/Ca²⁺/K⁺, рефрактерний період, автоматизм пейсмейкера. Порівняйте форми хвилі шлуночка, передсердя й синусового вузла.
Сила моделі Ходжкіна-Хакслі — у її загальності. Та сама формальна структура — ємнісна мембрана, напруго-керовані провідності, кінетика керуючих змінних першого порядку — описує нейрони в різних типах, серцеві клітини, бета-клітини підшлункової залози та волокна скелетних м'язів. Сучасна обчислювальна нейронаука використовує багатокомпартментні розширення з десятками типів іонних каналів, але математичний каркас завжди залишається моделлю Ходжкіна-Хакслі.
Частина 3: ферментативна кінетика — Міхаеліс-Ментен
Ферменти — це біологічні каталізатори. Вони зв'язують субстрат, перетворюють його на продукт і вивільняють продукт — усе це, не витрачаючись самі. Швидкість ферментативної реакції математично описали Леонор Міхаеліс і Мод Ментен у 1913 році. Їхня модель передбачає двоетапний механізм: зв'язування (оборотне) і каталіз (необоротний), з швидкою рівновагою між комплексом фермент-субстрат і вільним ферментом.
Кінетика Міхаеліса-Ментен
Реакція: E + S ⇌ ES → E + P
Швидкість реакції: v = V_max · [S] / (K_m + [S])
V_max = k_cat · [E_загальний] (максимальна швидкість)
K_m = (k₋₁ + k_cat) / k₁ (константа Міхаеліса)
При [S] = K_m: v = V_max / 2 (напівмаксимальна швидкість)
При [S] ≪ K_m: v ≈ (V_max / K_m) · [S] (перший порядок за S)
При [S] ≫ K_m: v ≈ V_max (нульовий порядок, насичення)
Лайнуівер-Берк: 1/v = (K_m / V_max) · (1/[S]) + 1/V_max
Ферментативна кінетика
Крива швидкості Міхаеліса-Ментен v відносно [S], подвійний зворотний графік Лайнуівера-Берка, конкурентне / неконкурентне / антиконкурентне інгібування, крива виснаження субстрату. Повзунки V_max і K_m. Подивіться, як інгібітори по-різному зсувають видимі K_m і V_max.
Кінетика хімічних реакцій
Константа швидкості Арреніуса k = A·e^(−Ea/RT), послідовні реакції A→B→C, вплив каталізатора на енергію активації. Кінетика першого й другого порядку. Спостерігайте, як криві концентрації від часу змінюються при налаштуванні енергії активації й температури.
Кінетика Хілла узагальнює Міхаеліса-Ментен для кооперативних систем. Багато біологічних реакцій мають не гіперболічну, а сигмоїдну форму — рівняння Хілла v = V_max · [S]ⁿ / (K_d + [S]ⁿ), де n > 1, описує кооперативне зв'язування. Гемоглобін зв'язує кисень кооперативно (n ≈ 2,8): щойно одна молекула O₂ зв'язана, наступне зв'язування полегшується. Це створює різкий перемикальний відгук, який робить гемоглобін ефективним переносником кисню. Та сама кооперативна логіка з'являється у зв'язуванні факторів транскрипції, керуванні іонними каналами й репресії циркадного годинника.
Частина 4: епідеміологія — модель SIR
Модель SIR ділить популяцію на три відсіки: сприйнятливі (можуть захворіти), інфіковані (наразі заразні) та одужалі (імунні). Кермак і Мак-Кендрик опублікували цю модель у 1927 році, натхненні пандемією грипу 1918 року. Попри свою простоту, модель відображає основну динаміку спалахів інфекційних хвороб — включно з критичною концепцією колективного імунітету.
Відсікова модель SIR
dS/dt = −β · S · I / N
dI/dt = β · S · I / N − γ · I
dR/dt = γ · I
N = S + I + R = const. (закрита популяція, без народжень/смертей)
β = швидкість передачі (контакти × ймовірність на контакт)
γ = швидкість одужання (= 1 / період заразності)
R₀ = β / γ (базове репродуктивне число)
Епідемічний поріг: R₀ > 1 (спалах відбувається)
Поріг колективного імунітету: частка імунних p_c = 1 − 1/R₀
Рівняння кінцевого розміру: S_∞ = N · e^(−R₀(1 − S_∞/N))
Епідемічна модель SIR
Живі криві S/I/R і 2D сітка агентів. Повзунки для β (передача) і γ (одужання). Спостерігайте епідемічний поріг: нижче R₀ = 1 інфекція вимирає; вище нього — хвиля епідемії прокочується популяцією. Увімкніть вакцинацію, щоб побачити колективний імунітет у дії.
Циркадний осцилятор
Осцилятор Гудвіна, що виробляє 24-годинні цикли мелатоніну, кортизолу й температури. Класичний приклад біологічного контуру зворотного зв'язку ОДР. Демонструє стійкі коливання від контуру негативного зворотного зв'язку із затримкою в часі — принцип, що живить усі циркадні годинники.
Частина 5: розв'язувачі ОДР — правильна чисельна реалізація
Біологічні ОДР рідко розв'язуються аналітично. Ми розв'язуємо їх чисельно: починаючи з початкових умов, крокуємо вперед у часі, використовуючи наближене правило. Найпростіший — метод Ейлера: x(t + dt) = x(t) + dt · f(x, t). Його легко реалізувати, але він швидко накопичує похибку і може бути нестійким для жорстких систем.
Метод Рунге-Кутта 4-го порядку (RK4) — робочий кінь біологічного моделювання. Він обчислює праву частину чотири рази за крок і поєднує результати з оптимальними вагами — даючи точність 4-го порядку за помірної вартості. Для жорстких систем (з різко різними часовими масштабами, як у ферментативних каскадах) потрібні неявні методи, такі як VODE, або формули чисельного диференціювання (NDF), щоб уникнути чисельного вибуху при розумних розмірах кроку.
Крок інтегрування 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₄)
Локальна похибка усічення: O(dt⁵); глобальна похибка: O(dt⁴)
Похибка методу Ейлера: O(dt²) за крок, O(dt) глобально
Жорсткі системи: використовуйте dt ≪ 1/λ_max (найбільше власне значення)
Жорсткість — прихований ворог біологічних симуляцій. Натрієвий канал у моделі Ходжкіна-Хакслі активується за ~0,2 мс, але симуляція нейрона може тривати 1000 мс. Відношення масштабів часу — близько 5000. З методом Ейлера dt має бути меншим за найшвидший масштаб часу (0,01 мс), що вимагає 100 000 кроків. RK4 при dt = 0,01 мс стійкий і точний. Але для ферментативних каскадів із масштабами часу від мілісекунд до годин необхідні адаптивні неявні розв'язувачі — інакше симуляції або нестійкі, або неприпустимо повільні.
Складаємо все разом: як читати фазову діаграму
Фазова площина — фундаментальний інструмент візуалізації для 2D систем ОДР. Побудуйте одну змінну відносно іншої; траєкторія системи прокреслює криву. Рівноваги з'являються як нерухомі точки; їхня стійкість визначається власними значеннями матриці Якобі в цій точці. Від'ємні дійсні частини = стійкість (спіраль чи вузол); додатні дійсні частини = нестійкість; нульові дійсні частини = центр (нейтральність).
З більш ніж двома змінними ми не можемо безпосередньо візуалізувати повний фазовий простір; натомість використовуємо перерізи Пуанкаре (площину крізь траєкторію), біфуркаційні діаграми (рівновага відносно параметра) та чисельне продовження. Ці інструменти показують, як якісна поведінка змінюється, коли параметр перетинає точку біфуркації: стійка рівновага стає нестійкою й породжує граничний цикл (біфуркація Хопфа), або дві рівноваги стикаються й взаємознищуються (сідло-вузлова біфуркація). Циркадний годинник, серцебиття та потенціал дії нейрона найкраще розуміти саме через цю призму.