Hodgkin-Huxley: Action Potential

// Hodgkin-Huxley action potential model (1952).
//
// State (4 variables):
//   V    membrane voltage (mV)
//   m    sodium activation gate         in [0,1]
//   h    sodium inactivation gate        in [0,1]
//   n    potassium activation gate       in [0,1]
//
// Currents (uA/cm^2):
//   I_Na = gNa * m^3 * h * (V - E_Na)
//   I_K  = gK  * n^4     * (V - E_K)
//   I_L  = gL          * (V - E_L)
//   C * dV/dt = I_inj - I_Na - I_K - I_L
//
// Squid giant axon parameters (Hodgkin & Huxley 1952), modern sign
// convention where E_Na > 0. Threshold sits around V = -55 mV.
//
// Integrated with classical RK4 at dt = 0.01 ms; many substeps per frame
// for stability through the spike upstroke.

// ---------- HH constants ----------
const C_M  = 1.0;     // membrane capacitance (uF/cm^2)
const G_NA = 120.0;   // mS/cm^2
const G_K  = 36.0;
const G_L  = 0.3;
const E_NA = 50.0;    // mV
const E_K  = -77.0;
const E_L  = -54.387;
const V_REST = -65.0;

// Integrator
const DT_SIM = 0.01;       // ms per RK4 step
const STEPS_PER_FRAME = 28; // sim ms per visible frame = STEPS*DT = 0.28 ms

// Trace buffer (each entry = 1 RK4 step, so 1 sample per 0.01 ms)
const TRACE_MS = 250;                          // visible window
const TRACE_LEN = Math.ceil(TRACE_MS / DT_SIM); // 25000 samples
let trace;        // Float32Array of V samples, ring buffer
let mTrace, hTrace, nTrace; // Float32Array gating traces, downsampled
const GATE_TRACE_LEN = 600; // smaller ring buffer for the gate plots
let gateHead = 0;
let gateCount = 0;
let gateAccum = 0;
let head = 0;
let count = 0;

// State
let V, m, h, n;
let simTimeMs = 0;
let spikeCount = 0;
let lastV = V_REST;
let aboveZero = false;

// Injected current (uA/cm^2): baseline + sustained from mouseY + transient pulse.
let iPulseRemain = 0;     // ms remaining of click pulse
let iPulseAmp = 0;        // amplitude of active click pulse
let iSustained = 0;       // continuous current from mouseY
let mouseYNorm = 0.5;     // 0=top, 1=bottom (normalized fraction)

// Layout
let W, H, dpr = 1;
let input = null;

// Bands within the canvas (top membrane, middle voltage trace, bottom gates)
let bandMembrane, bandTrace, bandGates;

// Colors
const COL_BG       = '#06080d';
const COL_BG2      = '#0d1117';
const COL_AXIS     = 'rgba(120,140,170,0.35)';
const COL_AXIS_DIM = 'rgba(120,140,170,0.15)';
const COL_TEXT     = 'rgba(190,200,220,0.85)';
const COL_DIM      = 'rgba(140,155,180,0.6)';
const COL_V        = '#7ed6ff';      // voltage trace
const COL_V_GLOW   = 'rgba(126,214,255,0.35)';
const COL_NA       = '#ff8a5b';      // sodium
const COL_K        = '#7afeb1';      // potassium
const COL_M        = '#ff8a5b';
const COL_H        = '#ffd06a';
const COL_N        = '#7afeb1';
const COL_THRESH   = 'rgba(255,180,90,0.4)';
const COL_REST     = 'rgba(120,170,220,0.28)';
const COL_MEMBRANE = '#3a4763';
const COL_LIPID    = '#1a2031';

// Animated channel pictograms — these are bands of "channels" sliding up the
// membrane to convey ion flux direction. We don't simulate individual ions —
// we visually link channel openness + flux direction to gating variables.
const CHANNEL_ROWS = 5; // rows of Na+ and K+ channels stacked vertically
let channelAnim = 0;     // phase for ion-flow ticks

// ---------- HH gating functions ----------

function alpha_m(V) {
  const x = -(V + 40);
  // Limit -> 1.0 as V -> -40 to avoid 0/0
  if (Math.abs(x) < 1e-6) return 1.0;
  return 0.1 * x / (Math.exp(x / 10) - 1);
}
function beta_m(V)  { return 4 * Math.exp(-(V + 65) / 18); }
function alpha_h(V) { return 0.07 * Math.exp(-(V + 65) / 20); }
function beta_h(V)  { return 1 / (Math.exp(-(V + 35) / 10) + 1); }
function alpha_n(V) {
  const x = -(V + 55);
  if (Math.abs(x) < 1e-6) return 0.1;
  return 0.01 * x / (Math.exp(x / 10) - 1);
}
function beta_n(V)  { return 0.125 * Math.exp(-(V + 65) / 80); }

// Steady-state at rest (used for init)
function steadyState(V) {
  const am = alpha_m(V), bm = beta_m(V);
  const ah = alpha_h(V), bh = beta_h(V);
  const an = alpha_n(V), bn = beta_n(V);
  return {
    m: am / (am + bm),
    h: ah / (ah + bh),
    n: an / (an + bn),
  };
}

// ---------- HH derivatives ----------

function derivs(V, m, h, n, I) {
  const INa = G_NA * m * m * m * h * (V - E_NA);
  const IK  = G_K  * n * n * n * n * (V - E_K);
  const IL  = G_L  *                 (V - E_L);
  const dV  = (I - INa - IK - IL) / C_M;

  const am = alpha_m(V), bm = beta_m(V);
  const ah = alpha_h(V), bh = beta_h(V);
  const an = alpha_n(V), bn = beta_n(V);

  return [
    dV,
    am * (1 - m) - bm * m,
    ah * (1 - h) - bh * h,
    an * (1 - n) - bn * n,
  ];
}

// Classical RK4 step.
function rk4Step(I) {
  const k1 = derivs(V, m, h, n, I);
  const k2 = derivs(
    V + 0.5 * DT_SIM * k1[0],
    m + 0.5 * DT_SIM * k1[1],
    h + 0.5 * DT_SIM * k1[2],
    n + 0.5 * DT_SIM * k1[3],
    I
  );
  const k3 = derivs(
    V + 0.5 * DT_SIM * k2[0],
    m + 0.5 * DT_SIM * k2[1],
    h + 0.5 * DT_SIM * k2[2],
    n + 0.5 * DT_SIM * k2[3],
    I
  );
  const k4 = derivs(
    V + DT_SIM * k3[0],
    m + DT_SIM * k3[1],
    h + DT_SIM * k3[2],
    n + DT_SIM * k3[3],
    I
  );
  V = V + (DT_SIM / 6) * (k1[0] + 2 * k2[0] + 2 * k3[0] + k4[0]);
  m = m + (DT_SIM / 6) * (k1[1] + 2 * k2[1] + 2 * k3[1] + k4[1]);
  h = h + (DT_SIM / 6) * (k1[2] + 2 * k2[2] + 2 * k3[2] + k4[2]);
  n = n + (DT_SIM / 6) * (k1[3] + 2 * k2[3] + 2 * k3[3] + k4[3]);
  // Clamp gates to physical [0,1]
  if (m < 0) m = 0; else if (m > 1) m = 1;
  if (h < 0) h = 0; else if (h > 1) h = 1;
  if (n < 0) n = 0; else if (n > 1) n = 1;
  // Clamp voltage to a generous range (depolarization block can push high)
  if (V > 80)  V = 80;
  if (V < -100) V = -100;
}

// ---------- layout ----------

function computeLayout() {
  // Vertical bands. Membrane on top (~32%), V trace in middle (~38%), gates bottom (~30%).
  const topH    = Math.max(120, Math.floor(H * 0.32));
  const traceH  = Math.max(120, Math.floor(H * 0.38));
  const gatesH  = H - topH - traceH;
  bandMembrane = { x: 0, y: 0,           w: W, h: topH };
  bandTrace    = { x: 0, y: topH,        w: W, h: traceH };
  bandGates    = { x: 0, y: topH + traceH, w: W, h: gatesH };
}

// ---------- init / tick ----------

function init({ canvas, ctx, width, height, input: inp }) {
  W = width; H = height;
  input = inp;
  computeLayout();

  // Initial state at rest
  V = V_REST;
  const ss = steadyState(V_REST);
  m = ss.m; h = ss.h; n = ss.n;

  trace = new Float32Array(TRACE_LEN);
  for (let i = 0; i < TRACE_LEN; i++) trace[i] = V_REST;
  head = 0;
  count = TRACE_LEN;

  mTrace = new Float32Array(GATE_TRACE_LEN);
  hTrace = new Float32Array(GATE_TRACE_LEN);
  nTrace = new Float32Array(GATE_TRACE_LEN);
  for (let i = 0; i < GATE_TRACE_LEN; i++) {
    mTrace[i] = m; hTrace[i] = h; nTrace[i] = n;
  }
  gateHead = 0; gateCount = GATE_TRACE_LEN; gateAccum = 0;

  simTimeMs = 0;
  spikeCount = 0;
  iPulseRemain = 0;
  iPulseAmp = 0;
  iSustained = 0;
  mouseYNorm = 0.5;
  channelAnim = 0;

  ctx.fillStyle = COL_BG;
  ctx.fillRect(0, 0, W, H);
}

function handleInput() {
  if (!input) return;
  // mouseY scrubs sustained current. Top of screen = 0, bottom = ~14 uA/cm^2.
  // Use the full canvas height so the user can scrub regardless of where the
  // cursor falls within the layout bands.
  if (typeof input.mouseY === 'number' && input.mouseY >= 0 && input.mouseY <= H) {
    mouseYNorm = input.mouseY / H;
    // Map: top half (0..0.5) -> 0..6 uA, bottom half (0.5..1) -> 6..16 uA
    // This gives a "no-spike / firing / block" feel.
    iSustained = mouseYNorm * 16;
  }
  // Clicks inject short pulse (1.5 ms, ~12 uA/cm^2). Multiple clicks stack.
  // consumeClicks() returns a Number (count); truthy when > 0.
  const clicks = (typeof input.consumeClicks === 'function') ? input.consumeClicks() : 0;
  if (clicks > 0) {
    iPulseRemain = 1.5;       // ms
    iPulseAmp = 12;           // uA/cm^2 — enough to fire from rest with no background
  }
}

function tick(args) {
  const ctx = args.ctx;
  if (args.width !== W || args.height !== H) {
    W = args.width; H = args.height;
    computeLayout();
  }

  handleInput();

  // ---- Integrate ----
  for (let s = 0; s < STEPS_PER_FRAME; s++) {
    let I = iSustained;
    if (iPulseRemain > 0) {
      I += iPulseAmp;
      iPulseRemain -= DT_SIM;
      if (iPulseRemain < 0) iPulseRemain = 0;
    }
    lastV = V;
    rk4Step(I);
    simTimeMs += DT_SIM;
    // Spike detection: crossing 0 mV from below
    if (!aboveZero && V > 0) { aboveZero = true; spikeCount++; }
    if (aboveZero && V < -20) aboveZero = false;
    // Append to voltage ring
    trace[head] = V;
    head = (head + 1) % TRACE_LEN;
    if (count < TRACE_LEN) count++;
    // Gate trace downsample: 1 sample per ~5 steps (~0.05 ms apart)
    gateAccum++;
    if (gateAccum >= 5) {
      gateAccum = 0;
      mTrace[gateHead] = m;
      hTrace[gateHead] = h;
      nTrace[gateHead] = n;
      gateHead = (gateHead + 1) % GATE_TRACE_LEN;
      if (gateCount < GATE_TRACE_LEN) gateCount++;
    }
  }
  channelAnim += 0.04;

  // ---- Draw ----
  drawBackground(ctx);
  drawMembrane(ctx);
  drawVoltageTrace(ctx);
  drawGates(ctx);
  drawHUD(ctx);
}

// ---------- drawing ----------

function drawBackground(ctx) {
  // Vertical gradient: a hair lighter near the membrane band, darker near the bottom.
  const g = ctx.createLinearGradient(0, 0, 0, H);
  g.addColorStop(0,    '#0a0e16');
  g.addColorStop(0.45, '#070a11');
  g.addColorStop(1,    '#05070c');
  ctx.fillStyle = g;
  ctx.fillRect(0, 0, W, H);

  // Subtle separators between the three bands
  ctx.strokeStyle = 'rgba(80,100,140,0.18)';
  ctx.lineWidth = 1;
  ctx.beginPath();
  ctx.moveTo(0, bandTrace.y + 0.5); ctx.lineTo(W, bandTrace.y + 0.5);
  ctx.moveTo(0, bandGates.y + 0.5); ctx.lineTo(W, bandGates.y + 0.5);
  ctx.stroke();
}

// ---- Top band: membrane schematic ----
function drawMembrane(ctx) {
  const b = bandMembrane;
  ctx.save();
  // Clip to band so glows don't bleed
  ctx.beginPath();
  ctx.rect(b.x, b.y, b.w, b.h);
  ctx.clip();

  const cy = b.y + b.h * 0.5;
  // Membrane band thickness scales with min dimension
  const memT = Math.max(38, Math.min(64, b.h * 0.42));
  const memTop = cy - memT * 0.5;
  const memBot = cy + memT * 0.5;

  // Extracellular side (top) — pale blue tint
  const gExt = ctx.createLinearGradient(0, b.y, 0, memTop);
  gExt.addColorStop(0, 'rgba(80,140,200,0.10)');
  gExt.addColorStop(1, 'rgba(80,140,200,0.02)');
  ctx.fillStyle = gExt;
  ctx.fillRect(b.x, b.y, b.w, memTop - b.y);

  // Intracellular side (bottom) — warm rose tint
  const gInt = ctx.createLinearGradient(0, memBot, 0, b.y + b.h);
  gInt.addColorStop(0, 'rgba(200,120,140,0.10)');
  gInt.addColorStop(1, 'rgba(200,120,140,0.02)');
  ctx.fillStyle = gInt;
  ctx.fillRect(b.x, memBot, b.w, b.y + b.h - memBot);

  // Voltage gradient across membrane: tint the bilayer based on V.
  // Map V in [-90, +50] to a color blend.
  const vNorm = Math.max(0, Math.min(1, (V + 90) / 140));
  const bilayerR = Math.floor(28 + vNorm * 70);
  const bilayerG = Math.floor(32 + vNorm * 30);
  const bilayerB = Math.floor(50 - vNorm * 10);

  // Lipid bilayer: two darker bands sandwiching a slightly lighter core
  ctx.fillStyle = `rgb(${bilayerR}, ${bilayerG}, ${bilayerB})`;
  ctx.fillRect(b.x, memTop, b.w, memT);

  // Inner highlight line
  ctx.strokeStyle = 'rgba(140,160,200,0.20)';
  ctx.lineWidth = 1;
  ctx.beginPath();
  ctx.moveTo(b.x, memTop + 0.5); ctx.lineTo(b.x + b.w, memTop + 0.5);
  ctx.moveTo(b.x, memBot - 0.5); ctx.lineTo(b.x + b.w, memBot - 0.5);
  ctx.stroke();

  // --- Phospholipid head/tail beads sprinkled along the bilayer ---
  // Cheap: a few dotted lines.
  const beadSpace = 12;
  ctx.fillStyle = 'rgba(190,210,240,0.35)';
  for (let x = b.x + 6; x < b.x + b.w; x += beadSpace) {
    ctx.fillRect(x, memTop + 2, 2, 2);
    ctx.fillRect(x, memBot - 4, 2, 2);
  }

  // --- Na+ and K+ channels ---
  // Lay out channels horizontally. Na channels on the left half, K on the right.
  // Each is rendered as a "pore" through the bilayer; openness = m^3*h for Na, n^4 for K.
  // Ion arrows show flux direction relative to the driving force.

  const halfW = b.w * 0.5;
  const naCount = Math.max(3, Math.floor(halfW / 95));
  const kCount  = Math.max(3, Math.floor(halfW / 95));

  const naOpen = clamp01(Math.pow(m, 3) * h);
  const kOpen  = clamp01(Math.pow(n, 4));

  // Driving force determines flux direction & magnitude (sign).
  const dfNa = (V - E_NA);  // negative when below E_Na -> Na flows INWARD (down arrow)
  const dfK  = (V - E_K);   // positive when above E_K -> K flows OUTWARD (up arrow)

  // Na channels
  for (let i = 0; i < naCount; i++) {
    const cx = b.x + (i + 0.5) * (halfW / naCount);
    drawChannel(ctx, cx, cy, memT, naOpen, dfNa, COL_NA, 'Na', i);
  }
  // K channels
  for (let i = 0; i < kCount; i++) {
    const cx = b.x + halfW + (i + 0.5) * (halfW / kCount);
    drawChannel(ctx, cx, cy, memT, kOpen, dfK, COL_K, 'K', i + 100);
  }

  // Labels: extracellular / intracellular
  ctx.font = '10px ui-monospace, monospace';
  ctx.fillStyle = COL_DIM;
  ctx.textAlign = 'left';
  ctx.fillText('extracellular', b.x + 10, b.y + 14);
  ctx.fillText('intracellular', b.x + 10, b.y + b.h - 6);

  // Section labels
  ctx.fillStyle = 'rgba(255,138,91,0.85)';
  ctx.font = 'bold 11px ui-monospace, monospace';
  ctx.fillText('Na+', b.x + 10, cy + 4);
  ctx.fillStyle = 'rgba(122,254,177,0.85)';
  ctx.fillText('K+',  b.x + halfW + 10, cy + 4);

  // Voltage indicator on the right edge of the band: little vertical thermometer.
  drawVoltageBar(ctx, b.x + b.w - 26, b.y + 12, 14, b.h - 24, V);

  ctx.restore();
}

function drawChannel(ctx, cx, cy, memT, openness, driveForce, color, kind, seed) {
  const poreW = 14;
  const poreH = memT - 4;
  const top = cy - poreH * 0.5;
  // Pore body
  ctx.fillStyle = '#161c2a';
  ctx.fillRect(cx - poreW * 0.5, top, poreW, poreH);
  ctx.strokeStyle = 'rgba(150,170,200,0.45)';
  ctx.lineWidth = 1;
  ctx.strokeRect(cx - poreW * 0.5 + 0.5, top + 0.5, poreW - 1, poreH - 1);

  // Inner channel: openness controls how much "open" the gate is.
  // Render an open region from the bottom (intracellular) upward.
  const innerW = poreW - 4;
  const innerH = poreH - 4;
  const openH = innerH * openness;
  const innerX = cx - innerW * 0.5;
  const innerY = top + 2;

  // Closed remainder (dark)
  ctx.fillStyle = '#0b0e16';
  ctx.fillRect(innerX, innerY, innerW, innerH - openH);
  // Open region (colored, glowing)
  if (openH > 0.5) {
    const g = ctx.createLinearGradient(innerX, innerY + innerH - openH, innerX, innerY + innerH);
    g.addColorStop(0, hexA(color, 0.95));
    g.addColorStop(1, hexA(color, 0.55));
    ctx.fillStyle = g;
    ctx.fillRect(innerX, innerY + innerH - openH, innerW, openH);
    // Glow halo
    ctx.fillStyle = hexA(color, 0.12 * Math.min(1, openness * 1.5));
    ctx.fillRect(innerX - 3, innerY + innerH - openH - 3, innerW + 6, openH + 6);
  }

  // Ion flux arrows: only render if channel is open enough AND there is drive.
  // Na: dfNa < 0 -> inward (downward arrow); K: dfK > 0 -> outward (upward arrow).
  const fluxStrength = clamp01(openness * Math.min(1, Math.abs(driveForce) / 80));
  if (fluxStrength > 0.05) {
    const dir = driveForce < 0 ? +1 : -1; // +1 = down (into cell), -1 = up (out of cell)
    const numIons = 3;
    const phase = (channelAnim + seed * 0.37) % 1;
    for (let k = 0; k < numIons; k++) {
      const t = (phase + k / numIons) % 1;
      // Travel from above pore (dir=+1) to below pore, or vice versa.
      let iy;
      if (dir > 0) {
        iy = (cy - poreH * 0.5 - 12) + t * (poreH + 20);
      } else {
        iy = (cy + poreH * 0.5 + 12) - t * (poreH + 20);
      }
      const a = 0.85 * fluxStrength * Math.sin(Math.PI * t); // fade in/out
      ctx.fillStyle = hexA(color, a);
      ctx.beginPath();
      ctx.arc(cx, iy, 2.4, 0, Math.PI * 2);
      ctx.fill();
    }
  }

  // Tiny label below the pore: openness percent (only on larger canvases)
  if (memT >= 50) {
    ctx.fillStyle = hexA(color, 0.7);
    ctx.font = '9px ui-monospace, monospace';
    ctx.textAlign = 'center';
    ctx.fillText(`${Math.round(openness * 100)}%`, cx, cy + poreH * 0.5 + 12);
    ctx.textAlign = 'left';
  }
}

function drawVoltageBar(ctx, x, y, w, h, V) {
  // Map V from [-90, +60] to height fraction (bottom = -90, top = +60)
  const minV = -90, maxV = 60;
  const frac = (V - minV) / (maxV - minV);
  const fillH = Math.max(0, Math.min(1, frac)) * h;

  // Frame
  ctx.fillStyle = 'rgba(20,25,40,0.85)';
  ctx.fillRect(x, y, w, h);
  ctx.strokeStyle = 'rgba(140,160,200,0.45)';
  ctx.lineWidth = 1;
  ctx.strokeRect(x + 0.5, y + 0.5, w - 1, h - 1);

  // Fill from bottom upward
  const fy = y + h - fillH;
  const g = ctx.createLinearGradient(x, fy, x, y + h);
  // Hue depends on whether we're depolarized or hyperpolarized vs rest
  let topCol, botCol;
  if (V > -55) { topCol = 'rgba(255,138,91,0.85)'; botCol = 'rgba(255,80,40,0.55)'; }
  else         { topCol = 'rgba(126,214,255,0.8)'; botCol = 'rgba(80,160,220,0.5)'; }
  g.addColorStop(0, topCol);
  g.addColorStop(1, botCol);
  ctx.fillStyle = g;
  ctx.fillRect(x + 1, fy, w - 2, fillH);

  // Tick at rest level
  const restY = y + h - ((V_REST - minV) / (maxV - minV)) * h;
  ctx.strokeStyle = COL_REST;
  ctx.beginPath();
  ctx.moveTo(x - 2, restY); ctx.lineTo(x + w + 2, restY);
  ctx.stroke();
  // Tick at threshold
  const thrY = y + h - ((-55 - minV) / (maxV - minV)) * h;
  ctx.strokeStyle = COL_THRESH;
  ctx.setLineDash([2, 3]);
  ctx.beginPath();
  ctx.moveTo(x - 2, thrY); ctx.lineTo(x + w + 2, thrY);
  ctx.stroke();
  ctx.setLineDash([]);

  ctx.fillStyle = COL_DIM;
  ctx.font = '9px ui-monospace, monospace';
  ctx.textAlign = 'right';
  ctx.fillText(`${V.toFixed(0)} mV`, x - 4, restY + 3);
  ctx.textAlign = 'left';
}

// ---- Middle band: voltage trace ----
function drawVoltageTrace(ctx) {
  const b = bandTrace;
  const padL = 44, padR = 16, padT = 16, padB = 22;
  const plotX = b.x + padL;
  const plotY = b.y + padT;
  const plotW = b.w - padL - padR;
  const plotH = b.h - padT - padB;

  // Plot frame
  ctx.fillStyle = 'rgba(8,11,18,0.65)';
  ctx.fillRect(plotX, plotY, plotW, plotH);

  // V range: -90 to +60 mV
  const vmin = -90, vmax = 60;
  const vToY = (v) => plotY + (1 - (v - vmin) / (vmax - vmin)) * plotH;

  // Horizontal reference lines
  ctx.font = '9px ui-monospace, monospace';
  ctx.textAlign = 'right';
  const refs = [
    { v: 50,    label: '+50',         col: 'rgba(120,140,170,0.15)' },
    { v: 0,     label: '  0',         col: 'rgba(120,140,170,0.22)' },
    { v: -55,   label: '-55 thresh',  col: COL_THRESH },
    { v: -65,   label: '-65 rest',    col: COL_REST },
    { v: -77,   label: '-77 EK',      col: 'rgba(122,254,177,0.18)' },
  ];
  for (const r of refs) {
    const y = vToY(r.v);
    ctx.strokeStyle = r.col;
    ctx.setLineDash(r.label.includes('thresh') ? [3, 3] : []);
    ctx.beginPath();
    ctx.moveTo(plotX, y); ctx.lineTo(plotX + plotW, y);
    ctx.stroke();
    ctx.setLineDash([]);
    ctx.fillStyle = r.label.includes('thresh') ? 'rgba(255,180,90,0.7)' :
                    r.label.includes('rest')   ? 'rgba(120,170,220,0.75)' :
                    'rgba(140,155,180,0.55)';
    ctx.fillText(r.label, plotX - 4, y + 3);
  }
  ctx.textAlign = 'left';

  // x-axis (time) ticks every 50 ms
  const tMaxMs = TRACE_MS;
  ctx.strokeStyle = COL_AXIS_DIM;
  ctx.fillStyle = COL_DIM;
  for (let t = 0; t <= tMaxMs; t += 50) {
    const x = plotX + (t / tMaxMs) * plotW;
    ctx.beginPath();
    ctx.moveTo(x, plotY + plotH - 3);
    ctx.lineTo(x, plotY + plotH);
    ctx.stroke();
    if (t !== 0 && t !== tMaxMs) {
      const label = `${-tMaxMs + t}`;
      ctx.textAlign = 'center';
      ctx.fillText(label, x, plotY + plotH + 11);
    }
  }
  ctx.textAlign = 'left';
  ctx.fillText('-250 ms', plotX, plotY + plotH + 11);
  ctx.textAlign = 'right';
  ctx.fillText('now', plotX + plotW, plotY + plotH + 11);
  ctx.textAlign = 'left';

  // Draw the trace. trace[head-1] is the newest sample. We map left edge =
  // oldest, right edge = newest.
  const startIdx = (head - count + TRACE_LEN) % TRACE_LEN;
  // Downsample: only one point per pixel.
  const samplesPerPx = count / plotW;

  // Underlay: area fill below 0 mV line for visual "spike" emphasis
  ctx.beginPath();
  ctx.moveTo(plotX, vToY(V_REST));
  let firstY = null;
  for (let px = 0; px < plotW; px++) {
    const s0 = Math.floor(px * samplesPerPx);
    let vmaxLocal = -1e9;
    let vminLocal = 1e9;
    const sEnd = Math.min(count, Math.floor((px + 1) * samplesPerPx));
    for (let i = s0; i < sEnd; i++) {
      const v = trace[(startIdx + i) % TRACE_LEN];
      if (v > vmaxLocal) vmaxLocal = v;
      if (v < vminLocal) vminLocal = v;
    }
    if (vmaxLocal === -1e9) {
      const v = trace[(startIdx + s0) % TRACE_LEN];
      vmaxLocal = v; vminLocal = v;
    }
    const y = vToY(vmaxLocal);
    if (firstY === null) { firstY = y; ctx.moveTo(plotX, y); }
    ctx.lineTo(plotX + px, y);
  }
  // Close path to baseline for fill
  ctx.lineTo(plotX + plotW, plotY + plotH);
  ctx.lineTo(plotX, plotY + plotH);
  ctx.closePath();
  const fillG = ctx.createLinearGradient(0, plotY, 0, plotY + plotH);
  fillG.addColorStop(0, 'rgba(126,214,255,0.20)');
  fillG.addColorStop(0.5, 'rgba(126,214,255,0.06)');
  fillG.addColorStop(1, 'rgba(126,214,255,0)');
  ctx.fillStyle = fillG;
  ctx.fill();

  // Main trace stroke (min/max per pixel as a thin vertical bar, then top line)
  ctx.strokeStyle = COL_V;
  ctx.lineWidth = 1.4;
  ctx.beginPath();
  for (let px = 0; px < plotW; px++) {
    const s0 = Math.floor(px * samplesPerPx);
    const sEnd = Math.min(count, Math.floor((px + 1) * samplesPerPx));
    let vmaxLocal = -1e9;
    let vminLocal = 1e9;
    for (let i = s0; i < sEnd; i++) {
      const v = trace[(startIdx + i) % TRACE_LEN];
      if (v > vmaxLocal) vmaxLocal = v;
      if (v < vminLocal) vminLocal = v;
    }
    if (vmaxLocal === -1e9) {
      const v = trace[(startIdx + s0) % TRACE_LEN];
      vmaxLocal = v; vminLocal = v;
    }
    if (px === 0) ctx.moveTo(plotX + px, vToY(vmaxLocal));
    else ctx.lineTo(plotX + px, vToY(vmaxLocal));
    // For pixels containing significant variation, also draw min so we don't
    // alias spikes into hairlines.
    if (vmaxLocal - vminLocal > 4) {
      ctx.stroke();
      ctx.beginPath();
      ctx.moveTo(plotX + px, vToY(vmaxLocal));
      ctx.lineTo(plotX + px, vToY(vminLocal));
      ctx.stroke();
      ctx.beginPath();
      ctx.moveTo(plotX + px, vToY(vmaxLocal));
    }
  }
  ctx.stroke();

  // Subtle glow line (one px wider, lower alpha)
  ctx.strokeStyle = COL_V_GLOW;
  ctx.lineWidth = 3;
  ctx.beginPath();
  for (let px = 0; px < plotW; px += 2) {
    const s0 = Math.floor(px * samplesPerPx);
    const v = trace[(startIdx + s0) % TRACE_LEN];
    if (px === 0) ctx.moveTo(plotX + px, vToY(v));
    else ctx.lineTo(plotX + px, vToY(v));
  }
  ctx.stroke();

  // Current voltage marker (right edge)
  const yNow = vToY(V);
  ctx.fillStyle = '#ffffff';
  ctx.beginPath();
  ctx.arc(plotX + plotW - 1, yNow, 3, 0, Math.PI * 2);
  ctx.fill();

  // Title
  ctx.fillStyle = COL_TEXT;
  ctx.font = '11px ui-monospace, monospace';
  ctx.fillText('V(t)  membrane voltage (mV)', plotX, plotY - 4);

  // Pulse indicator: a small triangle on the right edge if a pulse is active
  if (iPulseRemain > 0) {
    ctx.fillStyle = 'rgba(255,200,90,0.9)';
    ctx.beginPath();
    ctx.moveTo(plotX + plotW + 4, yNow - 4);
    ctx.lineTo(plotX + plotW + 4, yNow + 4);
    ctx.lineTo(plotX + plotW + 10, yNow);
    ctx.closePath();
    ctx.fill();
  }
}

// ---- Bottom band: m, h, n traces ----
function drawGates(ctx) {
  const b = bandGates;
  if (b.h < 40) return;
  const padL = 44, padR = 16, padT = 8, padB = 14;
  const plotX = b.x + padL;
  const plotY = b.y + padT;
  const plotW = b.w - padL - padR;
  const plotH = b.h - padT - padB;

  ctx.fillStyle = 'rgba(8,11,18,0.65)';
  ctx.fillRect(plotX, plotY, plotW, plotH);

  // 0 and 1 reference lines
  ctx.strokeStyle = COL_AXIS_DIM;
  ctx.beginPath();
  ctx.moveTo(plotX, plotY + 0.5); ctx.lineTo(plotX + plotW, plotY + 0.5);
  ctx.moveTo(plotX, plotY + plotH - 0.5); ctx.lineTo(plotX + plotW, plotY + plotH - 0.5);
  ctx.moveTo(plotX, plotY + plotH * 0.5); ctx.lineTo(plotX + plotW, plotY + plotH * 0.5);
  ctx.stroke();

  ctx.fillStyle = COL_DIM;
  ctx.font = '9px ui-monospace, monospace';
  ctx.textAlign = 'right';
  ctx.fillText('1.0', plotX - 4, plotY + 7);
  ctx.fillText('0.5', plotX - 4, plotY + plotH * 0.5 + 3);
  ctx.fillText('0.0', plotX - 4, plotY + plotH - 2);
  ctx.textAlign = 'left';

  // Draw three traces
  const startIdx = (gateHead - gateCount + GATE_TRACE_LEN) % GATE_TRACE_LEN;
  const valAtPx = (arr, px) => {
    const i = Math.floor((px / plotW) * gateCount);
    return arr[(startIdx + Math.min(gateCount - 1, i)) % GATE_TRACE_LEN];
  };
  function plot(arr, color, label, labelX) {
    ctx.strokeStyle = color;
    ctx.lineWidth = 1.3;
    ctx.beginPath();
    for (let px = 0; px < plotW; px++) {
      const v = valAtPx(arr, px);
      const y = plotY + (1 - v) * plotH;
      if (px === 0) ctx.moveTo(plotX + px, y);
      else ctx.lineTo(plotX + px, y);
    }
    ctx.stroke();
    // Label at right edge (current value)
    const vNow = arr[(gateHead - 1 + GATE_TRACE_LEN) % GATE_TRACE_LEN];
    const yNow = plotY + (1 - vNow) * plotH;
    ctx.fillStyle = color;
    ctx.font = '10px ui-monospace, monospace';
    ctx.fillText(label, labelX, yNow + 4);
  }
  plot(mTrace, COL_M, `m=${m.toFixed(2)}`, plotX + plotW + 4);
  plot(hTrace, COL_H, `h=${h.toFixed(2)}`, plotX + plotW + 4);
  plot(nTrace, COL_N, `n=${n.toFixed(2)}`, plotX + plotW + 4);

  // Title
  ctx.fillStyle = COL_TEXT;
  ctx.font = '11px ui-monospace, monospace';
  ctx.fillText('gating variables  m (Na act)  h (Na inact)  n (K act)', plotX, plotY - 1);

  // Tiny legend dots (anchored to the right of the plot so they never overflow
  // the title text on a narrow canvas).
  if (plotW > 240) {
    const baseX = plotX + plotW - 96;
    ctx.fillStyle = COL_M; ctx.fillRect(baseX,      plotY - 9, 8, 2);
    ctx.fillStyle = COL_H; ctx.fillRect(baseX + 32, plotY - 9, 8, 2);
    ctx.fillStyle = COL_N; ctx.fillRect(baseX + 64, plotY - 9, 8, 2);
  }
}

// ---- HUD overlay ----
function drawHUD(ctx) {
  // Top-left: live readout
  ctx.font = '10px ui-monospace, monospace';
  ctx.fillStyle = 'rgba(190,200,220,0.85)';
  const iTotal = iSustained + (iPulseRemain > 0 ? iPulseAmp : 0);
  ctx.fillText(`I_inj = ${iTotal.toFixed(1)} uA/cm^2`, 10, H - 30);
  ctx.fillStyle = COL_DIM;
  ctx.fillText(`spikes: ${spikeCount}`, 10, H - 18);

  // Bottom-right: regime label based on iSustained
  let regime = 'silent';
  let regimeCol = 'rgba(140,160,200,0.7)';
  if (iSustained > 9.5)      { regime = 'depolarization block'; regimeCol = 'rgba(255,138,91,0.85)'; }
  else if (iSustained > 6.2) { regime = 'sustained firing';     regimeCol = 'rgba(255,205,90,0.9)'; }
  else if (iSustained > 2.3) { regime = 'subthreshold';         regimeCol = 'rgba(180,200,230,0.85)'; }
  ctx.textAlign = 'right';
  ctx.fillStyle = regimeCol;
  ctx.fillText(regime, W - 10, H - 8);
  ctx.textAlign = 'left';
}

// ---------- utilities ----------

function clamp01(x) { return x < 0 ? 0 : x > 1 ? 1 : x; }

// Build "color with alpha" from a hex string like "#7afeb1"
function hexA(hex, a) {
  const h = hex.replace('#', '');
  const r = parseInt(h.slice(0, 2), 16);
  const g = parseInt(h.slice(2, 4), 16);
  const b = parseInt(h.slice(4, 6), 16);
  return `rgba(${r},${g},${b},${a})`;
}