ProjectTwo weekends$0Tier 0

Project A: Neurons from the equations up

Hodgkin–Huxley in Python from the four equations, then leaky integrate-and-fire, then Izhikevich, then a spiking network in Brian2. You heard a spike in Phase 1; now you make one.

AssumesWhat was that sound?Python for signals SpineComputational neuroscience

On this pageTwo weekends

Hodgkin–Huxley
Leaky integrate-and-fire
Izhikevich
A network in Brian2

In 1952 Hodgkin and Huxley wrote down four differential equations for the squid axonA quantitative description of membrane current and its application to conduction and excitation in nerveHodgkin, Huxley · J. Physiol. · 1952The four equations behind every spiking model in Phase 4, still readable in the original. and integrated them by hand on a mechanical calculator, and out came the action potential you heard from a cockroach leg. You can integrate them in twenty lines of Python and watch the same thing happen. Then you simplify them twice, to the models that networks are actually built from, and run a thousand neurons at once. The equations stop being symbols when you have watched them fire.

Hodgkin–Huxley

This interactive needs JavaScript. If you are reading a printout, the caption describes what it shows.

Figure 1. The model you are about to write, running live. Inject current and watch the membrane voltage and the three gating variables.

Try this

Press the brief pulse. One spike, then rest. Look at the gates: m rises fast, h falls, n rises slowly.
Raise the sustained current slowly. Below about 6 µA/cm² nothing; above it, a train of spikes appears at a nonzero rate straight away. That jump is a property of this model.
Raise it further. The rate rises, then the spikes shrink, then they stop. Too much current holds the sodium channels inactivated.

Hodgkin–Huxley

Four state variables: the membrane voltage V and three gating variables m, h, n between 0 and 1, which are the fraction of sodium activation gates open, sodium inactivation gates open, and potassium gates open.

The current balance, with $C$ the membrane capacitance, $I$ the injected current, and one leak term:

C\frac{dV}{dt} = I - g_{Na}m^3h,(V-E_{Na}) - g_K n^4,(V-E_K) - g_L,(V-E_L)

Each gate relaxes toward a voltage-dependent open fraction:

\frac{dx}{dt} = \alpha_x(V)(1-x) - \beta_x(V)x

with $\alpha$ and $\beta$ the voltage-dependent opening and closing rates that Hodgkin and Huxley fit to their data. Each gate has a steady value and a time constant worth knowing when checking an implementation:

x_\infty(V) = \frac{\alpha_x(V)}{\alpha_x(V)+\beta_x(V)}, \quad \tau_x(V) = \frac{1}{\alpha_x(V)+\beta_x(V)}

The squid parameters: $g_{Na} = 120$ , $g_K = 36$ , $g_L = 0.3\ \text{mS/cm}^2$ ; $E_{Na} = 50$ , $E_K = -77$ , $E_L = -54.4\ \text{mV}$ ; $C = 1\ \mu\text{F/cm}^2$ . The rate functions are in every textbook and in the site’s widget source.

import numpy as np
def hh(I_ext, T=100.0, dt=0.01):
    gNa, gK, gL, ENa, EK, EL, C = 120., 36., 0.3, 50., -77., -54.4, 1.
    an = lambda v: 0.01*(v+55)/(1-np.exp(-(v+55)/10)); bn = lambda v: 0.125*np.exp(-(v+65)/80)
    am = lambda v: 0.1*(v+40)/(1-np.exp(-(v+40)/10));  bm = lambda v: 4*np.exp(-(v+65)/18)
    ah = lambda v: 0.07*np.exp(-(v+65)/20);            bh = lambda v: 1/(1+np.exp(-(v+35)/10))
    n_steps = int(T/dt); V = np.empty(n_steps)
    v, m, h, n = -65., 0.05, 0.6, 0.32
    for i in range(n_steps):
        I = I_ext(i*dt)
        dv = (I - gNa*m**3*h*(v-ENa) - gK*n**4*(v-EK) - gL*(v-EL)) / C
        m += dt*(am(v)*(1-m) - bm(v)*m); h += dt*(ah(v)*(1-h) - bh(v)*h); n += dt*(an(v)*(1-n) - bn(v)*n)
        v += dt*dv; V[i] = v
    return V

Forward Euler with a 0.01 ms step is fine here. Plot V for a 2 ms pulse of 30 µA/cm²: one spike. Plot the gates alongside and you have reproduced Figure 1.

Then the classic experiments. Sweep the sustained current and count spikes per second: the frequency-current curve, with its jump at threshold. Deliver two pulses with a variable gap and find the shortest gap that produces a second spike: the refractory period. Block the potassium channels (set g_K to zero) and watch the spike fail to repolarize. Each is an afternoon and each is a figure for the build log.

Leaky integrate-and-fire

Throw away the gates. The membrane is a leaky capacitor:

\tau\frac{dV}{dt} = -(V-E_L) + RI

When $V$ crosses a threshold, record a spike and reset $V$ . Two lines of dynamics. It cannot produce a spike shape, but it captures when spikes happen, which for networks is what matters, and it runs ten thousand times faster. Reproduce the f-I curve and compare its shape to HH’s: the LIF curve rises smoothly from zero, the HH curve jumps.

Izhikevich

Two equations with four parameters that, by choosing the parameters, reproduce twenty different firing patterns real neurons show: regular spiking, bursting, chattering, fast spiking.

\frac{dv}{dt} = 0.04v^2 + 5v + 140 - u + I, \quad \frac{du}{dt} = a(bv-u)

On spike, $v \gets c$ and $u \gets u+d$ . With $a = 0.02$ , $b = 0.2$ , $c = -65$ , $d = 8$ it is a regular-spiking cortical neuron. Change $d$ to 2 and $c$ to −55 and it bursts. Reproduce the figure from the 2003 paper with all twenty patterns; it is the single most educational plot in computational neuroscience.

A network in Brian2

Brian2 lets you write the equations as strings and simulate thousands of neurons connected by synapses. The classic first network: 4000 LIF neurons, 80 percent excitatory and 20 percent inhibitory, randomly connected, driven by noise. Depending on the balance of excitation and inhibition it sits quietly, fires asynchronously, or falls into synchronized oscillations, and you can see population rhythms emerge from cells that individually have none. That emergence is where the oscillations explainer picks up.

from brian2 import *
eqs = 'dv/dt = (-(v - El) + I)/tau : volt (unless refractory)\nI : volt'
G = NeuronGroup(4000, eqs, threshold='v > -50*mV', reset='v = -65*mV', refractory=2*ms, method='euler')

Implement HH and reproduce a single spike with gates. Save the figure.
Sweep current for the f-I curve. Find the threshold current and the refractory period.
Implement LIF. Overlay its f-I curve on HH’s.
Implement Izhikevich and reproduce at least eight of the twenty firing patterns.
Install Brian2 and run the balanced random network. Vary the inhibitory weight and plot the population rate’s spectrum; find the setting where an oscillation appears.
Write the build log with all five figures.

What goes wrong

HH explodes to infinity. Time step too large, or a sign error in a rate function. Use 0.01 ms and check each α and β against the textbook at V = −65. No spikes at any current. Initial gate values wrong; start from the resting values given, or run 50 ms with no input first. Brian2 units errors. Every quantity has units in Brian2 and it refuses to add millivolts to milliseconds. This is a feature; read the error.

Stop here if you like. You have made an action potential from four equations and a rhythm from four thousand neurons. The explainers below are the biology behind the parameters.

ExhibitHodgkin–Huxley, LIF, and Izhikevich implemented from the equations with f-I curves and a firing-pattern figure, plus a Brian2 network showing an emergent oscillation.

Recall

What are the four state variables of the Hodgkin–Huxley model and what does each represent?

Membrane voltage V, and three gating variables: m (sodium activation), h (sodium inactivation), n (potassium activation), each the fraction of that gate type open.

Recall

What does the leaky integrate-and-fire model keep and throw away relative to Hodgkin–Huxley?

It keeps the membrane as a leaky capacitor and the timing of threshold crossings; it throws away the gating dynamics and therefore the spike shape and refractory mechanics.

Recall

Why is the Izhikevich model useful?

Two equations and four parameters reproduce most known cortical firing patterns at a tiny fraction of HH's cost, making large networks of realistic neurons feasible.