The Mathematical Beauty of a Neuron

A brain cell has a voltage, a small electrical difference between its inside and outside. Input from other cells pushes the voltage up. Left alone, it leaks back down to a resting level. If enough input arrives fast enough, the voltage crosses a tipping point and the cell fires: a quick electrical pulse that travels to other cells.

Each of those sentences is a rule about how fast the voltage is changing right now. Rules like that are called differential equations. They don't tell you where the voltage will be in a second. They only tell you which way it's heading at this instant. To find out where it ends up, you take a small step in that direction, look again, and repeat. Everything on this page is built on doing that well.

The state of a model neuron is a short list of numbers: its membrane voltage \(V\), maybe a slow adaptation current \(w\), and the currents arriving at its synapses. The model doesn't say what those numbers will be later. It says how fast each one is changing now, as a function of all of them:

That is an ordinary differential equation (ODE). For almost every model worth studying there is no formula for \(\mathbf{y}(t)\), so we integrate numerically: start from a known state, use \(\mathbf{f}\) to estimate where the state will be a short time \(h\) later, and repeat. The choice of how to make that estimate is the first decision in any neuron simulation, and it's the one people think about least. So that's where we start.

The simplest way to take a step: look at which way things are heading, and go that way for a bit. That's called Euler's method. It sounds fine, but there's a catch. The direction changes during the step, and Euler doesn't notice. On a curving path it always overshoots to the outside of the curve.

Below is a weight bouncing on a spring. Its real path is a perfect circle (the dashed line): it trades speed for stretch and back, forever. Euler (red) spirals outward, inventing energy from nothing. RK4 (blue) stays on the circle.

The simplest scheme is Euler's method: assume the slope stays constant over the step.

The error it makes on each step is \(O(h^2)\), so the error accumulated over a fixed time is \(O(h)\): halve the step, halve the error. That's slow. Worse, the errors aren't random. On a rotation, Euler always steps along the tangent, which points outside the circle, so every step adds a little energy. Below is the harmonic oscillator \(\ddot x = -x\), written as \(\dot x = v,\ \dot v = -x\). Its exact orbit is the unit circle.

So how does RK4 do so much better? Instead of looking once, it looks four times. It checks the direction where it is, then peeks halfway along, then peeks halfway again using what it just learned, then peeks at the end. It averages those four directions, giving the middle two double weight, and takes one step with the average. Four looks cost four times the work, but the answer can be thousands of times more accurate.

RK4 samples the slope four times inside each step and combines them:

Each probe uses the previous one to decide where to look. The weights 1, 2, 2, 1 are Simpson's rule: if \(\mathbf{f}\) depended only on \(t\), RK4 would be exactly Simpson integration of the slope. The figure below takes one step apart on \(y' = y\cos t\), whose exact solution \(y = e^{\sin t}\) lets us check every number.

In 1963, the weather scientist Edward Lorenz wrote down three simple rules for air rising and falling in a heated box. He found something strange. If he restarted his simulation from numbers rounded to three digits instead of six, the forecast came out completely different after a while. The tiniest difference grew and grew until it took over. This became known as the butterfly effect: the idea that a butterfly flapping its wings could change the weather weeks later.

Below, two paths start almost exactly at the same spot, differing by one part in 100 million. For a while they move together. Then they split up and loop around different wings of the butterfly. Yet both keep tracing the same butterfly shape. You can't predict where one path will be, but you can trust the shape.

The Lorenz system is three coupled equations, originally a drastically reduced model of convection:

At Lorenz's parameters \(\sigma = 10,\ \rho = 28,\ \beta = 8/3\), trajectories never settle and never repeat. They wind around two unstable fixed points, switching between the two wings in an order that looks random. Nearby trajectories separate exponentially, \(|\delta(t)| \approx |\delta_0|\, e^{\lambda t}\) with \(\lambda \approx 0.906\). That is what chaos means here: sensitive dependence on initial conditions, inside a bounded region of state space.

That makes Lorenz a hard test for a solver. Any numerical error acts like a perturbation and grows at the same exponential rate. Two trajectories that start \(10^{-8}\) apart below will lose each other after about \(\ln(10^{9})/0.906 \approx 23\) time units, no matter how good the integrator is.

The lesson for brain cells: when inputs are noisy, you can't predict the exact moment a cell will fire. What you can trust is the pattern, like how often it fires and how its connections change on average. That's what we'll watch.

The lesson to carry forward: with chaos, or with noise, individual trajectories can't be trusted for long, but their statistics can. A good solver keeps the shape of the attractor (where trajectories spend their time) even after the trajectory itself has diverged. Euler at a large step distorts the shape. RK4 at the same step preserves it. When we drive neurons with random input, the spike times will be different every run. The firing rates and the weight changes are what we'll read.

Now for the brain cell. Picture a bucket with a small hole in the bottom. Pour water in (that's the input) and the level rises. The hole lets it leak back out, so if you stop pouring, the level falls to empty. Draw a line near the top. When the water reaches the line, the cell fires: the bucket is dumped and it starts again.

That's the whole model, called the leaky integrate-and-fire neuron. "Integrate" means it adds up input. "Leaky" means it slowly forgets. "Fire" means it fires at the line. With a gentle trickle the leak wins and the level settles below the line: no firing at all. Pour harder and it fills and fires again and again, faster the harder you pour.

The leaky integrate-and-fire (LIF) neuron treats the membrane as a capacitor in parallel with a leak resistor. One ODE and one reset rule:

With \(C = 200\) pF and \(g_L = 10\) nS, the membrane time constant is \(\tau = C/g_L = 20\) ms. Under a constant current, \(V\) relaxes exponentially toward \(V_\infty = E_L + I/g_L\). If \(V_\infty\) stays below threshold, the neuron never fires. The smallest current that fires it is the rheobase, \(I_{\text{rh}} = g_L(\theta - E_L)\). Above it, the time from reset to threshold can be solved exactly, and with a refractory period \(t_{\text{ref}}\) after each spike, the firing rate is

That closed form is a gift: it lets us check the simulation. The dashed curve in the right panel is the formula. The dots are measured by running the RK4 neuron for a simulated second at each current.

The reset is a cliff, and RK4 can't see cliffs

RK4's accuracy comes from assuming the solution is smooth across the step. The reset breaks that assumption: \(V\) jumps from \(\theta\) to \(V_r\) instantly. If the code only checks for \(V \ge \theta\) at the end of each step and resets there, every spike is late by up to one step. Spike times carry the information downstream, and STDP in Section 6 depends on them at the millisecond level, so this matters. The fix is to find the crossing inside the step (here by linear interpolation), reset at that moment, and integrate the rest of the step from \(V_r\).

The bucket model is simple, maybe too simple. Its voltage creeps up to the line and then jumps back down: there's no actual spike, just a rule that says "fire now." Scientists have built better versions by changing one thing: how the voltage speeds up as it gets close to the tipping point. In a real cell, once the voltage gets high enough, it starts pushing itself up faster and faster. That runaway rise is the spike.

• LIF: the bucket. The voltage slows down as it approaches the line and is reset by a rule.
• QIF: near the tipping point, the voltage starts racing upward by itself.
• EIF: the same idea with an even sharper takeoff, which matches real cells well.
• AdEx: EIF plus tiredness. Each spike makes the next one a little harder.

The LIF's weakness shows in its equation: below threshold, \(dV/dt\) only ever pulls \(V\) back toward rest, harder the higher \(V\) goes. Nothing in the dynamics makes a spike. Real neurons have sodium channels that open as \(V\) rises, adding positive feedback, and that runaway upswing is the spike. The family of integrate-and-fire models differs mainly in one term: how \(dV/dt\) depends on \(V\).

In QIF and EIF the right-hand side turns upward and diverges, so \(V\) shoots off to infinity in finite time by itself. The simulation stops it at a peak and resets it. The exponential term in the EIF has a sharpness \(\Delta_T\). As \(\Delta_T \to 0\) it becomes a hard threshold and the EIF turns back into the LIF. Fits to cortical pyramidal cells give \(\Delta_T \approx 1\)–3 mV (Badel et al., 2008). AdEx adds a second variable \(w\), an adaptation current that jumps by \(b\) at every spike and decays with time constant \(\tau_w\).

Below, the faint dashed curve is the shape of a real spike, simulated with the full gate-by-gate model (Hodgkin–Huxley). Its parts are labeled: rest, tipping point, the fast rise, the top, the fall, the dip below rest, and a short time when it's hard to fire again. Pick a model to lay its spike on top and see which parts it gets.

The faint curve below is one action potential from the Hodgkin–Huxley model (squid axon parameters), integrated with the same RK4 code. Its phases are labeled. Each integrate-and-fire model is driven with a strong current step until its first spike, then the current is switched off so it can recover. All the spikes are aligned at the moment of the spike. The table separates what each model produces from its own dynamics from what is imposed by a reset rule or a cutoff.

The bottom panel is the key idea of this section: the shape of the upstroke is the shape of \(f(V)\). Where \(f\) crosses zero going downward, there is a stable rest state. Where it crosses going upward, there is an unstable point: push \(V\) past it and \(f\) keeps pushing. LIF's line never turns up, so its threshold has to be added as a rule. In QIF and EIF the threshold comes out of the equation, as the unstable fixed point. None of them repolarizes by itself. That needs a second, slower variable that pulls \(V\) back down, which is exactly what potassium channels do in Hodgkin–Huxley and what \(w\) does, crudely, in AdEx.

Real brain cells don't all fire the same way. Some tick like a metronome. Some start fast and tire. Some fire in bursts. The AdEx model can do all of these, just by changing a few numbers in its tiredness rule. Turn the tiredness off and every pattern turns back into the plain metronome.

One extra variable is enough to reproduce most of the firing patterns seen in cortical recordings. The presets below are parameter sets from Naud et al. (2008), each driven by the same kind of current step. The right panel is the \((V, w)\) phase plane, the same idea as the Lorenz picture but in two dimensions and with jumps. The dashed curves are nullclines, where \(\dot V = 0\) (green) or \(\dot w = 0\) (purple). Between spikes the trajectory follows the vector field. At each spike it leaves the right edge and is reinjected at \((V_r,\ w + b)\).

One cell on its own doesn't do much. Cells talk through connections called synapses. When cell A fires, it gives cell B a short kick of current, which raises B's voltage a little for a few milliseconds. One kick is usually not enough to make B fire. Several kicks close together add up.

Each connection also has a strength. A strong connection gives a big kick. A weak one barely registers. Learning, as we'll see next, is mostly about changing these strengths.

The simplest synapse model is a current that jumps when the presynaptic neuron spikes and decays exponentially:

\(I_{\text{syn}}\) feeds into neuron B's membrane equation as its input current. The state is now \([V_B, I_{\text{syn}}]\), integrated together with RK4. Filtering an exponential current through the membrane gives the familiar postsynaptic potential (PSP): a fast rise and a slower decay, here about 8 mV at peak for \(w = 1\). B's threshold is 15 mV above rest, so B fires only when several PSPs arrive within about one membrane time constant. This is temporal summation. The weight \(w\) is the number learning will change.

How does a connection know whether to get stronger? Through timing. If A fires just before B, A probably helped make B fire. That connection is useful, so it gets stronger. If A fires just after B, A had nothing to do with it, and the connection gets weaker. The closer together the two spikes, the bigger the change. This rule is called spike-timing-dependent plasticity, or STDP. It's often summed up as "fire together, wire together," but order matters.

How does a connection remember when its spikes happened? Each side keeps a fading memory of its last spike, like a footprint in sand. When B fires, the connection checks how fresh A's footprint is. When A fires, it checks B's.

Spike-timing-dependent plasticity, measured in hippocampal cultures by Bi & Poo (1998), changes a weight according to the relative timing \(\Delta t = t_{\text{post}} - t_{\text{pre}}\) of spike pairs:

A synapse doesn't store spike times. It keeps two decaying traces, which are ODEs just like everything else:

At each postsynaptic spike the weight gains \(A_+x\). At each presynaptic spike it loses \(A_-y\). Summed over a pair, that is exactly the window above. The figure runs the standard lab protocol: pair a pre and a post spike \(\Delta t\) apart, repeat 60 times at 1 Hz, and record the total change. The window isn't drawn in. It builds up from measurements.

The two rules differ in a way that matters for networks. In the additive rule a weight's change doesn't depend on the weight itself, so weights pushed in one direction keep going until they hit a bound. Under competition this gives a bimodal distribution: weights end up near 0 or near the maximum. The multiplicative rule adds a pull toward the middle that keeps weights unimodal and spread out (van Rossum, Bi & Turrigiano, 2000). The finale uses the additive rule because its winner-take-all behavior is easy to see.

Time to put everything together. One cell listens to 400 inputs. Every input fires about 20 times a second, at random. But there's a secret: half of the inputs are partly copying a hidden source. Now and then the source fires, and about one in five of those inputs fire at the same moment. Look at any single input and you can't tell which half it's in, since they all fire just as often. The only clue is that some of them tend to fire together.

The cell has no teacher and no score. It only has the timing rule from the last section. Press Play and watch: the connections from the inputs that fire together grow strong, and the rest fade away. The cell has found the hidden group by itself.

Now everything runs at once. One LIF neuron receives 400 Poisson inputs at 20 Hz each, through current-based exponential synapses, and every synapse runs the additive STDP rule. Inputs 1–200 share a hidden source: each source spike is copied by every one of them independently with probability \(c = 0.2\), and each makes up the rest of its 20 Hz independently. Inputs 201–400 are independent Poisson. Single-input rates are identical, so only correlations distinguish the groups. The window is slightly depression-dominated, \(A_-\tau_- = 1.05\,A_+\tau_+\). That is the regime Song, Miller & Abbott (2000) showed makes STDP competitive.

The full state, integrated together at \(h = 0.25\) ms, is \(V\), the synaptic current \(I\), the postsynaptic trace \(y\), 400 presynaptic traces \(x_i\), and 400 weights \(w_i\). The panel on the right lights up each equation as it acts. Notice the timescales: \(V\) and \(I\) move in milliseconds, the traces in tens of milliseconds, and the weights over tens of seconds.

Why does it work? When several inputs fire at the same moment, their kicks add up and often make the cell fire right after. So inputs from the hidden group often fire just before the cell, and the timing rule strengthens them. The inputs that fire on their own are as likely to fire just after the cell as just before. The rule is tuned to weaken slightly more than it strengthens, so on balance they fade. Strong connections make the cell fire even more in step with the group, so the effect feeds itself.

Why it works: a correlated input often fires together with others in its group, their PSPs sum, and the neuron tends to spike a few milliseconds later. So correlated inputs get more pre-before-post pairings than chance would give, and net LTP. An independent input's timing relative to the output is nearly random. It samples both lobes of the window equally, and because the LTD lobe is slightly larger, it drifts down. The weights that win make the output follow the source more closely, which strengthens the correlation they feed on: positive feedback, held in check by the hard bound \(w \le 1\). Depression dominance also keeps the output rate bounded. As weights grow, extra output spikes bring extra LTD.