LIF Population Raster

// Leaky integrate-and-fire population raster.
//
// N=100 neurons, each obeying  tau * dV/dt = -V + R*I_i  with heterogeneous
// drives I_i (Gaussian around a mean that crosses threshold for a fraction
// of the population). Sparse all-to-all synaptic coupling: when neuron j
// spikes, every other neuron's V receives a kick of J/N. V reaches V_th=1,
// spike registers, V reset to 0, refractory period applied.
//
// Coupling J knob (mouseY): low → asynchronous independent firing
// (Poisson-looking raster). Moderate → episodic synchronous bursts on top
// of the background. High → near-perfect periodic synchrony, vertical
// bands across the raster.
//
// Click: synchronous kick to a random 20% of neurons. Watch it propagate
// (or fail to propagate, at low J).
//
// Render: rolling raster (~1.5 s window, time left→right, neuron index
// top→bottom) above a smoothed instantaneous firing-rate trace.

const N = 100;
const TAU_MS = 20;        // membrane time constant
const V_TH = 1.0;         // threshold
const V_RESET = 0.0;
const REFRAC_MS = 3;
const DT_MS = 0.5;        // integration step (ms of sim time)
const STEPS_PER_FRAME = 4; // sim ms advanced per rendered frame = 2 ms
const WINDOW_MS = 1500;   // visible window
const RATE_BIN_MS = 2;    // bin for population rate

// state
let V;            // Float32Array(N)
let refrac;       // Int16Array(N) — ms remaining
let I;            // Float32Array(N) — heterogeneous drive
// Spike ring buffer: parallel typed arrays (avoids per-spike object allocation
// and the O(N) cost of Array.shift() to drop old spikes).
const SPIKE_CAP = 4096;
let spikeT;       // Float32Array — spike time (ms)
let spikeN;       // Int16Array — neuron index
let spikeHead = 0, spikeTail = 0, spikeCount = 0;
let rateHist;     // Float32Array of binned spike counts (per RATE_BIN_MS)
let rateSmoothed; // Float32Array reused for display smoothing
let rateHead;     // index into rateHist (newest)
let nowMs;        // sim time elapsed in ms
let coupling;     // J in [0, ~3.5]
let pulseQueued;  // ms remaining of a "kick"
let pulseAffected; // Uint8Array(N) — who got kicked this pulse
let lastClickAtMs;

function gauss() {
  // Box-Muller
  let u = 0, v = 0;
  while (u === 0) u = Math.random();
  while (v === 0) v = Math.random();
  return Math.sqrt(-2 * Math.log(u)) * Math.cos(2 * Math.PI * v);
}

function init() {
  V = new Float32Array(N);
  refrac = new Int16Array(N);
  I = new Float32Array(N);
  // mean drive of 1.05 (just above threshold for an isolated neuron) with
  // spread 0.15 → about 60% of neurons fire on their own, the rest are
  // subthreshold and only fire when pushed by network activity.
  for (let i = 0; i < N; i++) {
    I[i] = 1.05 + 0.15 * gauss();
    V[i] = Math.random() * V_TH * 0.8;  // desync initial phases
  }
  spikeT = new Float32Array(SPIKE_CAP);
  spikeN = new Int16Array(SPIKE_CAP);
  spikeHead = 0; spikeTail = 0; spikeCount = 0;
  const nBins = Math.ceil(WINDOW_MS / RATE_BIN_MS);
  rateHist = new Float32Array(nBins);
  rateSmoothed = new Float32Array(nBins);
  rateHead = 0;
  nowMs = 0;
  coupling = 0.6; // start in the interesting "moderate" regime
  pulseQueued = 0;
  pulseAffected = new Uint8Array(N);
  lastClickAtMs = -1e9;
}

function tick({ ctx, width: W, height: H, input, dt }) {
  // ----- input -----
  // mouseY scrubs coupling. Map top → high J (strong synchrony), bottom →
  // low J. Felt better this way: the dramatic synchronized state is "up".
  if (input && typeof input.mouseY === "number" && input.mouseY >= 0) {
    const t = 1 - Math.max(0, Math.min(1, input.mouseY / H));
    coupling = t * 3.5; // 0..3.5
  }
  if (input && typeof input.consumeClicks === "function") {
    const clicks = input.consumeClicks();
    if (clicks > 0) {
      // mark a random ~20% to receive a synchronous depolarizing kick that
      // pushes them toward (or over) threshold over the next ~2 ms.
      pulseAffected.fill(0);
      for (let i = 0; i < N; i++) {
        if (Math.random() < 0.20) pulseAffected[i] = 1;
      }
      pulseQueued = 4; // sim ms over which to deliver the kick
      lastClickAtMs = nowMs;
    }
  }

  // ----- simulation: advance STEPS_PER_FRAME * DT_MS ms of sim time -----
  const J = coupling;
  const decay = Math.exp(-DT_MS / TAU_MS);
  const drive = (1 - decay); // amount of (R*I_i) absorbed per step at equilibrium
  for (let step = 0; step < STEPS_PER_FRAME; step++) {
    // network input this step: sum of (J/N) over neurons that spiked.
    // We compute it after we know which neurons cross threshold, so first
    // integrate V, then detect spikes, then deliver the synaptic kick to
    // everyone simultaneously (delta-synapse, no transmission delay — fine
    // at this temporal resolution).
    let synKick = 0;
    // 1. exponential decay toward steady state V_inf = R*I_i, with R folded
    //    into I_i. Closed-form Euler-exact for the leak.
    for (let i = 0; i < N; i++) {
      if (refrac[i] > 0) {
        refrac[i] -= DT_MS;
        if (refrac[i] < 0) refrac[i] = 0;
        V[i] = V_RESET;
        continue;
      }
      V[i] = V[i] * decay + I[i] * drive;
    }
    // 2. external kick from pending click pulse
    if (pulseQueued > 0) {
      const kick = 0.35 * (DT_MS / 4); // total ~0.35 over the 4 ms window
      for (let i = 0; i < N; i++) if (pulseAffected[i] && refrac[i] === 0) V[i] += kick;
      pulseQueued -= DT_MS;
    }
    // 3. detect spikes
    let spikedThisStep = 0;
    for (let i = 0; i < N; i++) {
      if (refrac[i] === 0 && V[i] >= V_TH) {
        V[i] = V_RESET;
        refrac[i] = REFRAC_MS;
        spikedThisStep++;
        // Push onto the spike ring; if full, drop the oldest (will only
        // happen during a sustained high-J firestorm, which already gets
        // visually capped by the window cutoff below).
        spikeT[spikeHead] = nowMs;
        spikeN[spikeHead] = i;
        spikeHead = (spikeHead + 1) % SPIKE_CAP;
        if (spikeCount < SPIKE_CAP) spikeCount++;
        else spikeTail = (spikeTail + 1) % SPIKE_CAP;
      }
    }
    // 4. apply synaptic kick uniformly: each spike contributes J/N to every
    //    non-refractory neuron. (All-to-all coupling, mean-field style.)
    if (spikedThisStep > 0 && J > 0) {
      synKick = (J / N) * spikedThisStep;
      for (let i = 0; i < N; i++) if (refrac[i] === 0) V[i] += synKick;
    }
    nowMs += DT_MS;
  }

  // drop spikes older than the visible window — O(k) on the tail of the ring,
  // not O(n) like Array.shift().
  const cutoff = nowMs - WINDOW_MS;
  while (spikeCount > 0 && spikeT[spikeTail] < cutoff) {
    spikeTail = (spikeTail + 1) % SPIKE_CAP;
    spikeCount--;
  }

  // ----- rate histogram (binned spike count, smoothed) -----
  // We rebuild it from the spike ring each frame — cheap because N*window is bounded.
  rateHist.fill(0);
  const nBins = rateHist.length;
  for (let k = 0; k < spikeCount; k++) {
    const idx = (spikeTail + k) % SPIKE_CAP;
    const age = nowMs - spikeT[idx]; // 0..WINDOW_MS, 0 = newest
    const bin = nBins - 1 - Math.floor(age / RATE_BIN_MS);
    if (bin >= 0 && bin < nBins) rateHist[bin]++;
  }
  // simple 3-tap smoothing for display, into a reused buffer.
  const smoothed = rateSmoothed;
  for (let i = 0; i < nBins; i++) {
    const a = i > 0 ? rateHist[i - 1] : rateHist[i];
    const b = rateHist[i];
    const c = i < nBins - 1 ? rateHist[i + 1] : rateHist[i];
    smoothed[i] = (a + 2 * b + c) / 4;
  }

  // ----- layout -----
  const pad = 12;
  const titleH = 22;
  const hintH = 16;
  const rateH = Math.max(70, Math.min(110, H * 0.22));
  const gap = 10;
  const rasterTop = pad + titleH;
  const rasterBot = H - pad - rateH - gap - hintH;
  const rasterLeft = pad + 36; // axis labels
  const rasterRight = W - pad - 80; // legend column on the right
  const rasterW = rasterRight - rasterLeft;
  const rasterH = rasterBot - rasterTop;

  // ----- background -----
  ctx.fillStyle = "#08080c";
  ctx.fillRect(0, 0, W, H);

  // raster background
  ctx.fillStyle = "#0e0e16";
  ctx.fillRect(rasterLeft, rasterTop, rasterW, rasterH);
  ctx.strokeStyle = "#222";
  ctx.lineWidth = 1;
  ctx.strokeRect(rasterLeft + 0.5, rasterTop + 0.5, rasterW - 1, rasterH - 1);

  // raster dots
  const dotW = Math.max(1, Math.min(2, rasterW / (WINDOW_MS / DT_MS / 4)));
  const rowH = rasterH / N;
  for (let k = 0; k < spikeCount; k++) {
    const idx = (spikeTail + k) % SPIKE_CAP;
    const age = nowMs - spikeT[idx];
    if (age < 0 || age > WINDOW_MS) continue;
    const x = rasterLeft + (1 - age / WINDOW_MS) * rasterW;
    const y = rasterTop + spikeN[idx] * rowH;
    // newer spikes are brighter; recent click-affected pulse glows orange
    const recent = age < 30;
    if (recent) {
      ctx.fillStyle = "#ffd070";
    } else {
      ctx.fillStyle = "#a8d8ff";
    }
    ctx.fillRect(x, y, dotW, Math.max(1, rowH - 0.5));
  }

  // y-axis labels (neuron index)
  ctx.fillStyle = "#666";
  ctx.font = "10px monospace";
  ctx.textAlign = "right";
  ctx.fillText("n=1", rasterLeft - 4, rasterTop + 8);
  ctx.fillText(String(N), rasterLeft - 4, rasterBot - 2);
  ctx.save();
  ctx.translate(pad + 2, (rasterTop + rasterBot) / 2);
  ctx.rotate(-Math.PI / 2);
  ctx.textAlign = "center";
  ctx.fillStyle = "#888";
  ctx.fillText("neuron", 0, 0);
  ctx.restore();
  ctx.textAlign = "left";

  // x-axis time tick (now marker)
  ctx.strokeStyle = "#3a3a4a";
  ctx.beginPath();
  ctx.moveTo(rasterRight - 0.5, rasterTop);
  ctx.lineTo(rasterRight - 0.5, rasterBot);
  ctx.stroke();
  ctx.fillStyle = "#666";
  ctx.font = "10px monospace";
  ctx.fillText("now", rasterRight - 22, rasterBot + 11);
  ctx.fillText(`-${(WINDOW_MS / 1000).toFixed(1)}s`, rasterLeft + 2, rasterBot + 11);

  // ----- rate panel -----
  const rateTop = rasterBot + gap;
  const rateBot = rateTop + rateH;
  ctx.fillStyle = "#0e0e16";
  ctx.fillRect(rasterLeft, rateTop, rasterW, rateH);
  ctx.strokeStyle = "#222";
  ctx.strokeRect(rasterLeft + 0.5, rateTop + 0.5, rasterW - 1, rateH - 1);

  // scale: peak of smoothed (counts per bin) → top of panel, with floor
  let peak = 4;
  for (let i = 0; i < nBins; i++) if (smoothed[i] > peak) peak = smoothed[i];
  peak = Math.max(peak, 8);

  ctx.strokeStyle = "#7cf";
  ctx.lineWidth = 1.5;
  ctx.beginPath();
  for (let i = 0; i < nBins; i++) {
    const x = rasterLeft + (i / (nBins - 1)) * rasterW;
    const y = rateBot - (smoothed[i] / peak) * (rateH - 6) - 3;
    if (i === 0) ctx.moveTo(x, y); else ctx.lineTo(x, y);
  }
  ctx.stroke();
  // fill under the curve
  ctx.lineTo(rasterLeft + rasterW, rateBot);
  ctx.lineTo(rasterLeft, rateBot);
  ctx.closePath();
  ctx.fillStyle = "rgba(120,200,255,0.10)";
  ctx.fill();

  ctx.fillStyle = "#888";
  ctx.font = "10px monospace";
  ctx.fillText("population rate", rasterLeft + 6, rateTop + 12);

  // ----- right-side legend / regime indicator -----
  const legX = rasterRight + 12;
  ctx.fillStyle = "#ccc";
  ctx.font = "11px monospace";
  ctx.fillText("J", legX, rasterTop + 12);
  ctx.fillStyle = "#fc8";
  ctx.font = "14px monospace";
  ctx.fillText(coupling.toFixed(2), legX, rasterTop + 30);

  // regime label
  let regime;
  let regimeColor;
  if (coupling < 0.4) { regime = "async"; regimeColor = "#9af"; }
  else if (coupling < 1.4) { regime = "sparse bursts"; regimeColor = "#fc8"; }
  else if (coupling < 2.3) { regime = "synchrony"; regimeColor = "#f86"; }
  else { regime = "locked"; regimeColor = "#f48"; }
  ctx.fillStyle = regimeColor;
  ctx.font = "11px monospace";
  ctx.fillText(regime, legX, rasterTop + 50);

  // vertical J slider
  const sliderX = legX + 4;
  const sliderTop = rasterTop + 64;
  const sliderBot = rasterBot - 8;
  const sliderW = 6;
  ctx.fillStyle = "#1a1a26";
  ctx.fillRect(sliderX, sliderTop, sliderW, sliderBot - sliderTop);
  // marker
  const knobT = 1 - coupling / 3.5;
  const knobY = sliderTop + knobT * (sliderBot - sliderTop);
  ctx.fillStyle = "#fc8";
  ctx.fillRect(sliderX - 2, knobY - 2, sliderW + 4, 4);
  ctx.fillStyle = "#555";
  ctx.font = "9px monospace";
  ctx.fillText("3.5", sliderX + sliderW + 4, sliderTop + 4);
  ctx.fillText("0", sliderX + sliderW + 4, sliderBot);

  // mean rate (Hz, from total spikes in last second). Walk the ring backwards
  // from the newest spike until we leave the 1 s window.
  let recent1s = 0;
  for (let k = spikeCount - 1; k >= 0; k--) {
    const idx = (spikeTail + k) % SPIKE_CAP;
    if (nowMs - spikeT[idx] > 1000) break;
    recent1s++;
  }
  const meanRateHz = recent1s / N;
  ctx.fillStyle = "#aaa";
  ctx.font = "10px monospace";
  ctx.fillText(`${meanRateHz.toFixed(0)} Hz`, legX, rateBot - 4);

  // pulse flash indicator
  const sincePulse = nowMs - lastClickAtMs;
  if (sincePulse < 500) {
    ctx.fillStyle = `rgba(255,200,100,${0.6 * (1 - sincePulse / 500)})`;
    ctx.font = "10px monospace";
    ctx.fillText("kick!", legX, rasterTop + 78);
  }

  // ----- title + hint -----
  ctx.fillStyle = "#bbb";
  ctx.font = "12px monospace";
  ctx.textAlign = "left";
  ctx.fillText("LIF population — 100 neurons, all-to-all coupling J", pad, pad + 14);
  ctx.fillStyle = "#666";
  ctx.font = "10px monospace";
  ctx.fillText("drag Y: coupling J · click: kick 20% of neurons", pad, H - pad - 2);
}