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.
You are skimming: the title, the first figure, and the short version. Switch to Read in the header for the full page, or Deep to open every deep dive.
In 1952 Hodgkin and Huxley wrote down four differential equations for the squid axon 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.
- 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: C dV/dt = I − g_Na m³h (V − E_Na) − g_K n⁴ (V − E_K) − g_L (V − E_L).
Each gate: dx/dt = α_x(V)(1 − x) − β_x(V) x, with α and β the voltage-dependent opening and closing rates that Hodgkin and Huxley fit to their data. The squid parameters: g_Na = 120, g_K = 36, g_L = 0.3 mS/cm²; E_Na = 50, E_K = −77, E_L = −54.4 mV; C = 1 µF/cm². 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 VForward 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: τ dV/dt = −(V − E_L) + R I. 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 LIFLeaky integrate-and-fire neuronThe simplest useful neuron model: a leaky capacitor that fires and resets when it reaches a threshold. Glossary entry 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. dv/dt = 0.04v² + 5v + 140 − u + I; du/dt = a(bv − u); on spike, v ← c, u ← 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
Brian2Brian2A Python simulator for spiking neural networks in which you write the equations directly. Glossary entry 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.