Симуляція нейрона Годжкіна-Хакслі у браузері

1. Чотири рівняння Годжкіна-Хакслі

Мембрана поводиться як конденсатор паралельно з трьома змінними резисторами (йонні канали). Баланс струмів дає головне рівняння для напруги V, плюс по одному рівнянню кінетики першого порядку для кожної ворітної змінної:

Стандартні параметри для аксона кальмара (Hodgkin & Huxley, 1952):

• C_m = 1 мкФ/см² — ємність мембрани
• g_Na = 120 мСм/см², E_Na = 50 мВ — провідність та рівноважний потенціал натрію
• g_K = 36 мСм/см², E_K = −77 мВ — провідність та рівноважний потенціал калію
• g_L = 0.3 мСм/см², E_L = −54.4 мВ — провідність та потенціал витоку

m³h моделює швидку активацію Na⁺ (m, три незалежних ворота) у поєднанні з повільнішою інактивацією Na⁺ (h, одні ворота); n⁴ моделює активацію K⁺ із затримкою (чотири ворота). Ці конкретні показники степеня були підібрані Годжкіним і Хакслі емпірично, і вони й досі напрочуд добре відтворюють реальні дані аксона.

2. Функції швидкості α/β

Швидкість відкриття α та закриття β кожних воріт — це функції, залежні від напруги, підібрані емпірично з експериментів voltage-clamp:

function alphaM(V) { return 0.1 * (V + 40) / (1 - Math.exp(-(V + 40) / 10)); }
function betaM(V)  { return 4.0 * Math.exp(-(V + 65) / 18); }

function alphaH(V) { return 0.07 * Math.exp(-(V + 65) / 20); }
function betaH(V)  { return 1.0 / (1 + Math.exp(-(V + 35) / 10)); }

function alphaN(V) { return 0.01 * (V + 55) / (1 - Math.exp(-(V + 55) / 10)); }
function betaN(V)  { return 0.125 * Math.exp(-(V + 65) / 80); }

3. Вектор стану та функція похідних

Об'єднайте [V, m, h, n] в один масив і напишіть одну функцію, що повертає всі чотири похідні — саме її буде викликати інтегратор RK4 чотири рази на крок:

const params = {
  Cm: 1.0, gNa: 120, gK: 36, gL: 0.3,
  ENa: 50, EK: -77, EL: -54.4
};

function derivatives(state, Iext) {
  const [V, m, h, n] = state;
  const { Cm, gNa, gK, gL, ENa, EK, EL } = params;

  const iNa = gNa * m**3 * h * (V - ENa);
  const iK  = gK  * n**4 * (V - EK);
  const iL  = gL  * (V - EL);

  const dV = (Iext - iNa - iK - iL) / Cm;
  const dm = alphaM(V)*(1-m) - betaM(V)*m;
  const dh = alphaH(V)*(1-h) - betaH(V)*h;
  const dn = alphaN(V)*(1-n) - betaN(V)*n;

  return [dV, dm, dh, dn];
}

Початковий стан спокою при V = −65 мВ: ворота встановлюються у свій стаціонарний стан, m₀ ≈ 0.05, h₀ ≈ 0.6, n₀ ≈ 0.32 — обчислені один раз з x_∞(−65) = α(−65)/(α(−65)+β(−65)) для кожних воріт.

4. Інтегрування RK4

Ворота активації Na⁺ (m) відкриваються менш ніж за мілісекунду — явному методу Ейлера потрібен непрактично малий крок часу, щоб лишатися стабільним. Метод Рунге-Кутта четвертого порядку (RK4) дає точність четвертого порядку на кожному кроці, дозволяючи комфортно використовувати dt ≈ 0.01 мс:

function rk4Step(state, Iext, dt) {
  const add = (a, b, scale) => a.map((v, i) => v + b[i] * scale);

  const k1 = derivatives(state, Iext);
  const k2 = derivatives(add(state, k1, dt/2), Iext);
  const k3 = derivatives(add(state, k2, dt/2), Iext);
  const k4 = derivatives(add(state, k3, dt), Iext);

  return state.map((v, i) =>
    v + (dt/6) * (k1[i] + 2*k2[i] + 2*k3[i] + k4[i])
  );
}

5. Стимуляція зовнішнім струмом

Цікава поведінка з'являється, щойно ви подаєте струм. Нижче порогу (приблизно 6–7 мкА/см² для короткого імпульсу з цими параметрами) збурення просто загасає. Вище порогу петля позитивного зворотного зв'язку між деполяризацією та активацією Na⁺ (m зростає → більше вхідного струму → V зростає швидше → m зростає ще більше) запускає повний, стереотипний спайк — закон «все або нічого»:

function stimulusCurrent(t) {
  // 1 мс імпульс 10 мкА/см², що починається при t = 5 мс
  return (t >= 5 && t < 6) ? 10 : 0;
}

let state = [-65, 0.05, 0.6, 0.32]; // [V, m, h, n]
let t = 0;
const dt = 0.01; // мс — достатньо малий для стабільної кінетики Na+

function tick() {
  const Iext = stimulusCurrent(t);
  state = rk4Step(state, Iext, dt);
  t += dt;
}

6. Відображення на канвасі

Додавайте V(t) у циклічний буфер і малюйте його як лінійний графік, що прокручується — точно як осцилограма. Накладіть m, h, n, масштабовані до [0,1], на другій осі, щоб побачити, як ворота відкриваються й закриваються синхронно зі спайком:

const ctx = canvas.getContext('2d');
const history = []; // зразки { V, m, h, n }

function render() {
  ctx.clearRect(0, 0, canvas.width, canvas.height);
  const w = canvas.width, h = canvas.height;

  // Графік напруги: перетворення V у [-90, 40] мВ на Y канваса
  ctx.beginPath();
  history.forEach((s, i) => {
    const x = (i / history.length) * w;
    const y = h - ((s.V + 90) / 130) * h;
    i === 0 ? ctx.moveTo(x, y) : ctx.lineTo(x, y);
  });
  ctx.strokeStyle = '#60a5fa';
  ctx.stroke();
}

7. Рефрактерний період і повторне спрацювання

Одразу після спайку h (інактивація Na⁺) близька до нуля, а n (активація K⁺) підвищена — нейрон тимчасово не може спрацювати знову навіть за сильного стимулу. Це абсолютний рефрактерний період. Коли h відновлюється, а n релаксує назад протягом декількох мілісекунд, нейрон входить у відносний рефрактерний період, коли лише сильніший за звичайний стимул може запустити наступний спайк.

Подайте на модель тривалий (а не імпульсний) струм вище другого, вищого порогу — і отримаєте періодичне повторне спрацювання; частота зростає з Iext до точки насичення, визначеної тим, наскільки швидко можуть відновитися ворітні змінні — спрощена версія кодування частоти спрацювання, яку реальні нейрони використовують для представлення інтенсивності стимулу.