Hodgkin-Huxley Neuron Model Simulation

The Hodgkin-Huxley model describes how neurons generate action potentials through ion channel dynamics. It models:

• Membrane potential () changes over time.
• Sodium () and potassium () currents contributing to action potentials.
• Gating variables () that control ion channel opening.

Hodgkin-Huxley Equations

The membrane voltage follows:

where:

• is the membrane capacitance.
• is an external stimulus current.
• are sodium, potassium, and leak currents.

Each current follows Ohm’s law:

where:

• is conductance (voltage- and time-dependent).
• is the equilibrium potential for ion .

Ionic Currents

Sodium Current

Potassium Current

Leak Current

Gating Variables Dynamics

The activation () follows:

where and are voltage-dependent rates.

Numerical Solution

• We use Euler’s method for integration.
• External current is injected at multiple times to trigger multiple action potentials.
• Results plotted using Observable Plot.

TypeScript Code

This code:

1. Numerically solves the Hodgkin-Huxley model.
2. Generates multiple action potentials in response to periodic external current pulses.
3. Plots both the membrane potential and the injected current.

import {range} from "d3-array";
import * as Plot from "@observablehq/plot";

// Hodgkin-Huxley parameters
const params = {
  C_m: 1.0,
  g_Na: 120.0, E_Na: 50.0,
  g_K: 36.0, E_K: -77.0,
  g_L: 0.3, E_L: -54.387
};

// External current (stimulus)
function I_ext(t: number, amp = 10) {
  return t >= 10 && t <= 40 ? amp : 0;
}

// Alpha-beta functions
const alpha_m = (V: number) => 0.1 * (V + 40) / (1 - Math.exp(-(V + 40) / 10));
const beta_m  = (V: number) => 4 * Math.exp(-(V + 65) / 18);
const alpha_h = (V: number) => 0.07 * Math.exp(-(V + 65) / 20);
const beta_h  = (V: number) => 1 / (1 + Math.exp(-(V + 35) / 10));
const alpha_n = (V: number) => 0.01 * (V + 55) / (1 - Math.exp(-(V + 55) / 10));
const beta_n  = (V: number) => 0.125 * Math.exp(-(V + 65) / 80);

// Simulation step using Euler method
function simulateHH(dt = 0.01, T = 50, amp = 10) {
  const steps = Math.ceil(T / dt);
  let V = -65, m = 0.05, h = 0.6, n = 0.32;

  const results = range(steps).map(i => {
    const t = i * dt;
    
    const INa = params.g_Na * m**3 * h * (V - params.E_Na);
    const IK  = params.g_K * n**4 * (V - params.E_K);
    const IL  = params.g_L * (V - params.E_L);
    const dV = (I_ext(t, amp) - INa - IK - IL) / params.C_m;

    const dm = alpha_m(V)*(1-m) - beta_m(V)*m;
    const dh = alpha_h(V)*(1-h) - beta_h(V)*h;
    const dn = alpha_n(V)*(1-n) - beta_n(V)*n;

    V += dV * dt;
    m += dm * dt;
    h += dh * dt;
    n += dn * dt;

    return {t, V, m, h, n, INa, IK, gNa: params.g_Na*m**3*h, gK: params.g_K*n**4};
  });

  return results;
}

const data = simulateHH();

// 1. Membrane Potential vs Time
display(Plot.plot({
  marks: [
    Plot.line(data, {x: "t", y: "V"}),
  ],
  x: {label: "Time (ms)"},
  y: {label: "Voltage (mV)"},
  width: 700,
  height: 300,
  caption: "Membrane Potential vs. Time"
}))

// Ion Currents Plot
display(Plot.plot({
  marks: [
    Plot.line(data, {x: "t", y: "INa", stroke: "#f8766d", label: "Na⁺ Current"}),
    Plot.line(data, {x: "t", y: "IK", stroke: "#1f77b4", label: "K⁺ Current"}),
  ],
  x: {label: "Time (ms)"},
  y: {label: "Current (µA/cm²)"},
  width: 700,
  height: 300,
  caption: "Na⁺ and K⁺ Currents vs. Time"
}))

// Conductances
display(Plot.plot({
  marks: [
    Plot.line(data, {x: "t", y: "gNa", stroke: "#ff7f0e", label: "Na⁺ Conductance"}),
    Plot.line(data, {x: "t", y: "gK", stroke: "#2ca02c", label: "K⁺ Conductance"}),
  ],
  x: {label: "Time (ms)"},
  y: {label: "Conductance (mS/cm²)"},
  width: 700,
  height: 300,
  caption: "Ion Conductances vs. Time"
}))

// Phase-space plot
display(Plot.plot({
  marks: [Plot.line(data, {x: "m", y: "V", stroke: "#9467bd"})],
  x: {label: "m gating variable"},
  y: {label: "Membrane potential (mV)"},
  width: 500,
  height: 400,
  caption: "Phase-Space Plot (Voltage vs. m-gating variable)"
}))

// Stimulus Intensity Comparison
const amps = [6, 10, 20];
const intensityData = amps.flatMap(amp => simulateHH(0.05, 50, amp).map(d => ({...d, amp})));

display(Plot.plot({
  marks: [
    Plot.line(intensityData, {x: "t", y: "V", stroke: "amp", z: "amp"}),
    Plot.text(intensityData.filter(d => d.t===45), {x: "t", y: "V", text: d => `${d.amp} µA/cm²`, dx: 5})
  ],
  x: {label: "Time (ms)"},
  y: {label: "Voltage (mV)"},
  color: {legend: true},
  width: 700,
  height: 400,
  caption: "Neuron Response to Different Stimulus Intensities"
}))