How a neural network learns: a playground

The figure below is a complete neural network: two inputs, a handful of hidden units you choose, one output, trained by gradient descent on whatever you point it at, running in your browser. It is the smallest network that shows the things that matter about large ones, and the point of this page is to let you see them happen rather than read about them. Each experiment is a couple of clicks. If the words neural network, hidden unit or training are new, read the beginner’s guide first; if you have done lesson 3, you have already built the version of this with no hidden units.

Loading the interactive figure…
The playground. Pick a dataset, choose how many hidden units, press train. Filled dots are training examples, hollow dots are the held-out test set, and the shading is the network’s current answer for every point in the plane.

What you are looking at

Each dot is an example with two measurements (its position) and a label (its colour). The network's job is to learn a rule that gives the right colour from the position. The big square shows the rule it currently has: for every point on the plane, blue or orange, stronger where it is more sure. The small diagram on the right is the network itself, with each weight drawn as a line whose thickness is its size and whose colour is its sign. The curves below it are the loss, a measure of how wrong the network is, on the training examples (solid) and on the examples it never sees during training (dashed).

Training is the same loop as lesson 3: show the network the training examples, measure how wrong it is, work out for every weight which direction would make it less wrong, and move each weight a little in that direction. The learning rate is how little. One press of "one step" does that once; "train" does it a few times per frame until you stop it.

Five experiments

  1. A straight line is enough. Choose two blobs and set hidden units to none. Press train. This is logistic regression, exactly lesson 3's decoder, and it can only draw a straight boundary. For two blobs a straight line is all you need, and both accuracies go to about 95 percent within a hundred epochs. Now watch the weights in the diagram: two lines, one from each input, and their thickness is the whole model.
  2. Rings need a bend. Choose rings, still with no hidden units, and train. It cannot do it: there is no straight line that separates an inner disc from an outer ring, and the accuracy stays near 50 percent however long you wait. Set hidden units to 1 and reset: a little better, since one unit can bend the boundary once. Set it to 3: the boundary closes into a loop and the accuracy climbs past 95 percent. Each hidden unit draws one straight line of its own; the output combines them, and three lines can enclose a region.
  3. The problem that stalled the field. Choose four corners: blue where both inputs have the same sign, orange where they differ. With no hidden units this is impossible, and in 1969 Minsky and Papert used exactly this example to argue that single-layer networks were a dead end, which helped stop research for a decade. Two hidden units solve it. Watch how: each learns one diagonal, and the output learns to say "one but not both".
  4. Overfitting, live. Choose spiral, set 8 hidden units and a learning rate of about 1.5, and let it train for a while. The training loss keeps falling. The dashed test loss falls, flattens, and then, often, begins to rise, while the shading grows strange tongues that fit single training dots. The network is memorising the examples rather than learning the rule, and the only thing that reveals it is the data you held out. This is why lesson 3 insisted on a test set, and why every serious model is judged on data it has not seen.
  5. Real neurons. Choose real: two neurons. The dots are the 340 trials from lesson 3, positioned by how many spikes neurons 141 and 190 fired, coloured by whether a stimulus was on the screen. With no hidden units the network reaches about 97 percent, which is lesson 3's result for these two cells. Now add hidden units. Nothing improves, because the two classes overlap where they overlap and no boundary, however bendy, can separate what the neurons themselves did not. More model does not fix less information. Lesson 3 got to 99.7 percent not with a cleverer network but with 64 more neurons.

Things to notice along the way

  • The learning rate is a dial between slow and unstable. Set it very low and training crawls. Set it to 3 on the spiral and the loss jumps around instead of falling: each step overshoots. Real training uses schedules that start high and shrink, for this reason.
  • Reset changes the answer. "Reset weights" starts from different random values. On the spiral, train twice from different starts and you will get two different boundaries with similar accuracy. There is no single right network, only a family of acceptable ones, which is one reason the weights of a trained model are hard to interpret.
  • Width buys shape, not understanding. Hidden units are straight lines; more of them means a boundary with more bends. That is the entire source of a network's flexibility, in two dimensions and, with more layers, in billions.
  • Nothing here knows what a ring is. The network has no concept of circles, corners or spirals. It has weights and a gradient. Everything that looks like understanding in the shading is the gradient having found a boundary that lowers the loss on these dots. Keep that in mind when a much larger version of the same machine produces a paragraph.

How it relates to the rest of the site

The network with no hidden units is the logistic regression of lesson 3, and the GLM primer explains why that is one member of a family. The hidden layer is new, and it is what turns a linear decoder into a universal one: a theorem from 1989 says that one hidden layer with enough units can approximate any reasonable boundary, which the rings and spiral experiments demonstrate in miniature. The train/test split is lesson 3's central rule, and the overfitting experiment is the reason for it. The real-data experiment is a reminder that in neuroscience the data, not the model, are usually the limit, which is the recurring lesson of lesson 4 and of the whole Foundations track.

The network is written in about a hundred lines of plain JavaScript inside this site's plugin: a forward pass, a backward pass computed by hand, and a loop. There is nothing in it that lesson 3 did not already contain except the hidden layer and the chain rule. The real dataset is neurons 141 and 190 from the Steinmetz recordings (Nature 2019, CC-BY, via Neuromatch).

How a brain works, how a language model works, and where they meet

This page is for someone who has heard that artificial intelligence was inspired by the brain, that language models are “neural networks”, and that nobody quite agrees on what either machine is doing. It explains, with no equations and no code, what a neuron does, what a brain does with billions of them, what a language model does, and where the two really overlap. The pictures move. Each one is a cartoon of the mechanism, drawn to be true in spirit rather than in detail, and the text says where the cartoon ends.

One neuron, one unit

Start with the smallest parts. A neuron is a cell. It receives signals from other neurons on its branches (dendrites), it has an electrical voltage across its skin, and when that voltage is pushed high enough it fires: a spike, a thousandth of a second long, that travels down its axon to the thousands of neurons it connects to. Between spikes the voltage leaks back to rest, so inputs that arrive close together count for more than inputs spread out in time. A neuron is a device for deciding, moment by moment, whether enough has happened recently.

An artificial unit, the building block of a language model, was invented in the 1940s and 50s as a cartoon of exactly that. It receives numbers. Each number is multiplied by a weight, which stands in for how strong a connection is. The products are added up and the sum is squashed into a fixed range, and that is the output. There is no voltage, no leak, no spike and no time: the unit is a formula, evaluated once.

Loading the interactive figure…
Watch them side by side. Left: spikes arrive, the voltage fills and leaks, and the neuron fires when it reaches threshold. Right: numbers arrive, meet their weights, are summed and squashed, and a number comes out.

Both are summarised by one idea: add up your inputs, weighted by how much you trust each source, and respond when the total is large. That idea is the whole inheritance from neuroscience to artificial intelligence, and lesson 1 builds the neuron version in thirty lines of Python so you can see the leak and the threshold do their work.

A brain: billions of them, in time

A brain is about 86 billion neurons wired into one another, each receiving from thousands and sending to thousands. Three things about how they work together are easy to miss from the single-neuron picture.

  • It is sparse. At any instant only a small fraction of neurons are firing; the typical neuron fires a few times a second, many fire less than once. Most of the brain is nearly silent most of the time, and that silence is informative: which cells are active, out of the many that could be, carries the message.
  • It runs in time. A spike is an event with a precise moment. Information is carried by when cells fire relative to each other, not just whether they fire, and rhythms (the alpha and theta waves of lesson 6) organise the timing across whole regions.
  • It is organised by what it does. Different regions handle vision, hearing, movement, memory and planning, each with its own wiring, connected in loops rather than in a single direction. There is no input end and output end; the system is always in the middle of doing something.

A language model: the next word, a trillion times

A large language model is a program with one job: given the words so far, produce a probability for every word that could come next. Type "The mouse saw the", and it reports that "cheese" is likely, "cat" fairly likely, "screen" possible, and "democracy" not. Pick a word from that list, add it to the text, ask again, and you have a machine that writes. Everything a chatbot does, answering, summarising, coding, is this loop.

Inside, the words are converted to lists of numbers, and those numbers pass through dozens of layers of the artificial units from the first figure, hundreds of billions of weights in all. One ingredient makes modern models different from earlier ones: attention. At each layer, each word position works out how much to draw on every earlier position. Reading "the" at the end of "The mouse saw the", the model attends strongly to "mouse" and "saw", because those are what make "cheese" and "cat" likely. Attention is how a model keeps track of what a sentence is about, and it is learned, not programmed.

Loading the interactive figure…
Left: a cartoon of a patch of cortex, with spikes rippling through a sparse network a few cells at a time. Right: a cartoon of a language model reading a sentence, attending to the earlier words, and producing probabilities for the next one. The sentences and numbers are illustrative.

How does it learn the weights? By reading. The model is shown trillions of words of text, one at a time, asked to predict the next, and told how wrong it was. A procedure called backpropagation works out, for every one of the billions of weights, which direction would have made the prediction less wrong, and nudges each by a tiny amount. Repeat for months on thousands of computers. No one writes the rules of grammar or the facts in; they are whatever settles into the weights because they make the next word more predictable.

Where they differ

A brainA large language model
Basic unitA neuron: a cell that fires brief electrical spikesAn artificial unit: a weighted sum of numbers, squashed
How manyAbout 86 billion neurons, each with thousands of connectionsHundreds of billions to trillions of weights, arranged in layers
SignalSpikes: all-or-nothing events, timed to the millisecondNumbers, passed through the layers once per word
ActivitySparse: most neurons are nearly silent most of the timeDense: every unit computes on every pass
TimeBuilt in: everything is a process unfolding at its own paceNot built in: a sequence is handled as a list of positions
LearningSynapses strengthen or weaken from local activity; nobody has found a global error signalEvery weight is nudged by a global error signal, backpropagation, over trillions of words
PowerAbout 20 wattsMegawatts to train; a data centre to run at scale
What it predictsWhat the senses will report next, and the outcome of actions, continuouslyThe next word, given the words so far

Two rows deserve emphasis. The brain learns without anything resembling backpropagation: a synapse changes based on what happens at that synapse, and no one has found the global error signal that would be needed to tell every connection how it contributed to a mistake. How brains learn as well as they do with only local rules is one of the central open questions of the field, and one of the places where neuroscience may still have something to give back to AI. And the brain runs on 20 watts, about a dim light bulb, because most of it is silent most of the time. A model that evaluates every unit on every pass cannot do that, and the gap in energy is not a detail; it is a sign that the two machines are organised on different principles.

Where they meet

Both are prediction machines. A brain is constantly predicting what its senses will report next and what its actions will bring about, and much of what it does can be read as correcting those predictions; a language model predicts the next word. Both store what they know in the strengths of connections rather than in a list of facts. Both are made of simple units that are individually dumb and collectively not. And both have turned out to be far more capable than their units would suggest, in ways their builders did not plan and do not fully understand.

Those similarities are why the two fields keep borrowing from each other: convolutional networks came from visual cortex, reinforcement learning from the dopamine system, and in the other direction, language models are now used as models of how the brain processes language, because their internal activity predicts brain recordings surprisingly well. The essay on the difference between AI and machine learning and the opening argument of this site go further into that exchange.

Words you will meet

  • Neuron: a brain cell that sends electrical spikes. Synapse: the junction where one neuron influences another; its strength is what learning changes.
  • Spike (action potential): the millisecond-long electrical pulse a neuron fires. Firing rate: spikes per second, the most common summary of what a neuron is doing.
  • EEG: the brain's electrical activity recorded from the scalp, which is the sum of millions of neurons at once.
  • Unit (artificial neuron): a weighted sum followed by a squashing function. Weight: a number standing for a connection strength; a model is its weights.
  • Token: the unit of text a language model reads, usually a word or part of one. Attention: the mechanism by which each position in the text draws on the others.
  • Backpropagation: the procedure that works out how to adjust every weight after a mistake. Training: doing that trillions of times.

Where to go from here

If this page made sense, you are ready for the learning path. Lesson 1 builds the neuron from the first figure in Python and loads a real recording of 734 of them. Lesson 3 builds the artificial unit from the same figure and uses it to read a mouse's mind. Lesson 5 shows what spikes look like coming off a real electrode, and lesson 7 turns brain rhythms into a control signal. None of it needs more than Python and an hour a week. And if you would rather see a network learn before you build one, the neural network playground trains one in your browser while you watch.

The animations on this page are cartoons of the mechanisms, built to be watched rather than measured. The neuron is a leaky integrate-and-fire model; the network is a random feed-forward sketch; the language-model sentences and probabilities are invented to show attention and next-word prediction. Real numbers for all of these are in the lessons.

Detecting a seizure: from scalp EEG to a Raspberry Pi monitor

This is the first article in the precision-medicine track, and it is a project rather than a lesson: take the tools from lesson 6 and lesson 7, point them at a clinical problem, and carry the result all the way to a device. The problem is seizure detection from scalp EEG. The data are real: seven seizures from one child in the CHB-MIT database. The detector is five features and a logistic regression, and it catches all seven with a five-second delay. The second half of the article turns it into a service that runs on a Raspberry Pi with an open-source EEG board, with the diagrams to build it. Throughout, it is specific about what this can and cannot detect, because that distinction is the whole of the difference between a project and a medical device.

Which seizures, exactly

“Seizure” covers very different events, and a detector built from scalp EEG can only find the ones that change the scalp EEG. The International League Against Epilepsy classifies seizures by where they start (focal or generalised), by awareness, and by whether there is movement. For a detector the more useful split is electrographic: does the seizure produce a rhythm at the scalp that is large, sustained and different from the background? The seizures in this article do. All seven start in the right hemisphere, under the right temporal chain of electrodes, as a rhythm that grows over seconds and spreads; they last 27 to 101 seconds; and they are what the CHB-MIT documentation describes as intractable (drug-resistant) seizures in a child being evaluated for surgery. In the ILAE language they are focal-onset seizures with a clear scalp EEG correlate, several of which spread to both hemispheres. Here is the honest scope of what the methods on this site can and cannot do with them:

EventWhat the scalp EEG showsThis detector
Focal seizure with a scalp signature (the seven here: right-hemisphere onset, 27 to 101 s)Rhythmic activity that starts in one region, grows, spreads and slows; large and sustainedYes. This is what it was built and tested on
Focal to bilateral tonic-clonic (a “grand mal” that starts focally)Same start, then very large, then muscle artefact swamps everythingYes, by the onset; the later part is mostly muscle, which also trips the features
Absence seizure (typical, childhood)3-per-second spike-and-wave, generalised, abrupt, 5 to 20 sProbably: large rhythmic delta/theta is exactly what the weights reward, but it is untested here; see the exercises
Focal aware seizure with little or no scalp changeOften nothing visible; the source is small or deepNo. Nothing to detect at the scalp
Myoclonic jerkA single brief burst, under a secondNo. Shorter than the window and the persistence rule
Non-epileptic event (fainting, a psychogenic seizure)No ictal rhythm; often movement artefactCannot tell; may false-alarm on the artefact
Which type a seizure isNeeds the clinical picture and usually videoNo. Detection is not diagnosis

Two more boundaries. Everything here is detection: raising an alarm while the seizure is happening. Prediction, raising it minutes before, is a different problem that decades of work have not solved reliably, and anyone selling it to you should be asked for their false-alarm rate. And the detector is patient-specific: trained on this child’s seizures and tested on her other seizures. A detector that works across patients it has never seen is a harder problem, and the state of the art there is still a human reading the trace.

The data

The CHB-MIT Scalp EEG Database was collected at Boston Children’s Hospital: 22 children with drug-resistant epilepsy, monitored for days while their medication was reduced to provoke seizures ahead of possible surgery, 23 channels at 256 Hz, with every seizure marked by a clinician. It is the dataset most seizure-detection papers of the last fifteen years were tested on. Patient chb01 is an eleven-year-old girl with seven seizures over the 40 hours of her recording. The channels are bipolar: each is the difference between two neighbouring electrodes, which cancels what they share, the same idea as lesson 7’s Laplacian. The file holds two chains of four, running front to back on the left (FP1-F7 to P7-O1) and on the right (FP2-F8 to P8-O2).

Step 1: look at one

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import butter, filtfilt

dat = np.load("chbmit_chb01.zip")
eeg = dat["eeg"].astype(float)                 # 8 channels x samples, microvolts
fs = float(dat["fs"])                          # 128 samples per second
names = list(dat["channels"])                  # bipolar pairs: each channel is the difference between two electrodes
sz_start, sz_end = dat["seizure_start"], dat["seizure_end"]      # sample indices of the seven seizures
seg_start = dat["segment_start"]               # where each recording segment begins in the array
print(eeg.shape, "=", eeg.shape[1] / fs / 60, "minutes;", len(sz_start), "seizures lasting", ((sz_end - sz_start) / fs).astype(int), "seconds")
print("channels:", names)
# (8, 494336) = 64.36666666666666 minutes; 7 seizures lasting [ 40  27  40  51  90  93 101] seconds
# channels: ['FP1-F7', 'F7-T7', 'T7-P7', 'P7-O1', 'FP2-F8', 'F8-T8', 'T8-P8', 'P8-O2']

k = 0                                          # the first seizure
t0 = sz_start[k] / fs
window = slice(int((t0 - 30) * fs), int((t0 + 60) * fs))
t = np.arange(eeg.shape[1]) / fs - t0
fig, ax = plt.subplots(figsize=(11, 6))
for i in range(8):
    ax.plot(t[window], eeg[i, window] - i * 400, lw=0.5)
ax.axvline(0, color="orange"); ax.axvline((sz_end[k] - sz_start[k]) / fs, color="orange", ls="--")
ax.set(yticks=-np.arange(8) * 400, yticklabels=names, xlabel="time from seizure onset (s)")
plt.show()
Seizure 1 of 7, on all eight channels. Thirty seconds of ordinary background, then at time zero a rhythm begins on the right-hand chain (blue), grows over ten seconds, spreads to the left chain, and stops forty seconds later, leaving a slower, flatter trace. The spectrogram of F8-T8 below shows the same thing as lesson 6 would: a burst of power that starts across many frequencies and slows as the seizure evolves. The horizontal line at 16 Hz is an equipment artefact present throughout.
Seizure 1 of 7, on all eight channels. Thirty seconds of ordinary background, then at time zero a rhythm begins on the right-hand chain (blue), grows over ten seconds, spreads to the left chain, and stops forty seconds later, leaving a slower, flatter trace. The spectrogram of F8-T8 below shows the same thing as lesson 6 would: a burst of power that starts across many frequencies and slows as the seizure evolves. The horizontal line at 16 Hz is an equipment artefact present throughout.

Everything a detector needs is visible by eye. The seizure is larger than the background, it is rhythmic, it evolves (fast at first, slowing toward the end, which is the signature that distinguishes a seizure from an artefact of fixed frequency), and it lasts. The features below are each a way of putting a number on one of those.

Step 2: five numbers per second

Cut the recording into two-second windows, one per second, as in lesson 7. For each window compute the power in four bands, delta through beta, the way lesson 6 did, and one more feature: line length, the average absolute difference between consecutive samples, which is large when the trace is both big and fast and is the single most used feature in the seizure-detection literature because it is cheap and it works. Take each feature on every channel and keep the loudest channel, so that a seizure confined to one side is not diluted by the quiet side. Five numbers per window. Label a window as seizure if its centre falls inside one of the clinician’s marked intervals.

def features(window_data, fs):
    """One row of numbers describing a 2-second window of multichannel EEG."""
    out = []
    for lo, hi in [(1, 4), (4, 8), (8, 13), (13, 30)]:                        # delta, theta, alpha, beta
        b, a = butter(3, [lo / (fs / 2), hi / (fs / 2)], "band")
        out.append(np.log(filtfilt(b, a, window_data, axis=1).var(axis=1) + 1e-6).max())   # band power, loudest channel
    out.append(np.log(np.abs(np.diff(window_data, axis=1)).mean(axis=1)).max())           # line length: how wiggly the trace is
    return np.array(out)

win, step = int(2 * fs), int(1 * fs)
starts = np.arange(0, eeg.shape[1] - win, step)                    # one window per second
X = np.array([features(eeg[:, s:s + win], fs) for s in starts])   # windows x 5 features
label = np.zeros(len(starts), dtype=int)
for s, e in zip(sz_start, sz_end):
    label[(starts + win / 2 >= s) & (starts + win / 2 <= e)] = 1  # a window is "seizure" if its centre is inside one
print(X.shape, "windows x features;", label.sum(), "seizure windows,", (label == 0).sum(), "normal windows")
# (3860, 5) windows x features; 449 seizure windows, 3411 normal windows
What the detector sees through seizure 1, one window per second. Top: delta and theta power on the loudest channel, which rise ten to fifty fold at onset. Middle: line length, which triples. Bottom: the classifier's score, with the threshold dashed and the alarm (red) seven seconds after onset. The two line-length spikes after the seizure are movement artefacts, and the score ignores them because their band powers do not match.
What the detector sees through seizure 1, one window per second. Top: delta and theta power on the loudest channel, which rise ten to fifty fold at onset. Middle: line length, which triples. Bottom: the classifier’s score, with the threshold dashed and the alarm (red) seven seconds after onset. The two line-length spikes after the seizure are movement artefacts, and the score ignores them because their band powers do not match.

Step 3: the classifier, tested honestly

The classifier is lesson 3‘s logistic regression on the five features, with a small penalty on the weights. The important part is how it is tested. The 449 seizure windows and 3,411 normal windows are not independent: windows from the same seizure are near-copies of each other, and a random split would put copies of every test seizure in the training set and report an accuracy that means nothing. So the split is by seizure: for each of the seven, train on the other six segments and score this one. Every score below was produced by a model that had never seen the seizure it was scoring. The seizure-free quarter hour is scored by a model trained on all seven.

def train_decoder(X, y, lr=0.1, steps=3000):
    w, b = np.zeros(X.shape[1]), 0.0
    for _ in range(steps):
        p = 1 / (1 + np.exp(-(X @ w + b)))
        w -= lr * (X.T @ (p - y) / len(y) + 0.01 * w)
        b -= lr * (p - y).mean()
    return w, b

def score(X, w, b):
    return 1 / (1 + np.exp(-(X @ w + b)))

segment_of = np.searchsorted(seg_start, starts, side="right") - 1    # which recording segment each window came from
scores = np.zeros(len(starts))
for seg in range(len(sz_start)):                                      # leave one seizure's segment out, test on it
    train = segment_of != seg
    mean, sd = X[train].mean(axis=0), X[train].std(axis=0)
    w, b = train_decoder((X[train] - mean) / sd, label[train])
    scores[segment_of == seg] = score((X[segment_of == seg] - mean) / sd, w, b)
mean, sd = X[segment_of < 7].mean(axis=0), X[segment_of < 7].std(axis=0)   # the seizure-free hour: scored by a model trained on all seven
w, b = train_decoder((X[segment_of < 7] - mean) / sd, label[segment_of < 7])
scores[segment_of == 7] = score((X[segment_of == 7] - mean) / sd, w, b)
print("weights (delta, theta, alpha, beta, line length):", w.round(2))
# weights (delta, theta, alpha, beta, line length): [1.21 1.38 0.49 0.4  0.42]

The weights are readable. Delta and theta power carry the most, line length and the faster bands less. That matches the figure: this patient’s seizures are slow rhythms, 2 to 8 Hz, that get large.

Step 4: from a score to an alarm

A score above 0.5 on a single window is not an alarm. Backgrounds have one-second bumps all the time, and a device that buzzed on each of them would be switched off within a day. Real detectors add two pieces of logic that have nothing to do with machine learning. Persistence: the score has to stay above threshold for several windows in a row. Refractory period: after an alarm, stay quiet for a while, so that a ninety-second seizure produces one alarm and not thirty. Both trade latency and sensitivity against false alarms, and the numbers that matter clinically are exactly those three: how many seizures were caught, how long after onset, and how many false alarms per hour.

The decision logic on a made-up score trace. A single window over the threshold is ignored. The third consecutive window fires the alarm. For the next thirty seconds the detector stays silent.
The decision logic on a made-up score trace. A single window over the threshold is ignored. The third consecutive window fires the alarm. For the next thirty seconds the detector stays silent.
Loading the interactive figure…
Tune the detector yourself. The sliders are the three settings from step 4; the numbers and alarms update as you move them, and clicking the timeline shows the EEG at that moment.
def alarms(scores, threshold=0.5, persistence=3, refractory=30):
    """Fire when the score has been above threshold for `persistence` windows in a row; then stay quiet for `refractory` windows."""
    fired, streak, quiet_until = [], 0, -1
    for i, s in enumerate(scores):
        streak = streak + 1 if s > threshold else 0
        if streak >= persistence and i > quiet_until:
            fired.append(i); quiet_until = i + refractory
    return np.array(fired)

def evaluate(fired):
    detected, latency = 0, []
    for s, e in zip(sz_start, sz_end):
        hits = [f for f in fired if s <= starts[f] + win / 2 <= e + 10 * fs]      # an alarm during the seizure (or within 10 s after)
        if hits:
            detected += 1; latency.append((starts[hits[0]] + win / 2 - s) / fs)
    false = [f for f in fired if label[max(0, f - 5):f + 6].sum() == 0]           # alarms with no seizure within 5 s
    hours_normal = (label == 0).sum() / 3600
    return detected, np.mean(latency) if latency else np.nan, len(false) / hours_normal

for persistence in [1, 3, 5]:
    d, lat, fa = evaluate(alarms(scores, persistence=persistence))
    print(f"persistence {persistence} windows: {d} of 7 seizures detected, mean latency {lat:4.1f} s, {fa:4.1f} false alarms per hour")
# persistence 1 windows: 7 of 7 seizures detected, mean latency  6.9 s, 23.2 false alarms per hour
# persistence 3 windows: 7 of 7 seizures detected, mean latency  5.1 s,  3.2 false alarms per hour
# persistence 5 windows: 7 of 7 seizures detected, mean latency  9.9 s,  1.1 false alarms per hour
All 64 minutes of the file: the classifier score, the clinician's seizure intervals in orange, alarms as triangles, false alarms in red. The seven seizures are each caught within seconds. The three false alarms all fall in the last quarter hour, the recording from a different hour of the day that no training segment resembled.
All 64 minutes of the file: the classifier score, the clinician's seizure intervals in orange, alarms as triangles, false alarms in red. The seven seizures are each caught within seconds. The three false alarms all fall in the last quarter hour, the recording from a different hour of the day that no training segment resembled.

All seven seizures caught at every setting. With persistence 3, the alarm comes 5 seconds after onset on average and there are 3.2 false alarms per hour; with persistence 5, latency doubles and false alarms drop to about one an hour. That is the trade-off every detector lives on, and where to sit on it is a question for the person wearing it, not the engineer. The odd result, that persistence 1 has longer latency than persistence 3, is the refractory period at work: a single-window false alarm just before a seizure silences the detector for thirty seconds and delays the real alarm.

Look at where the false alarms are. All three fall in the fifteen minutes from a different hour, a recording whose background the model had never seen. Within the segments around seizures, the background was learned and the false-alarm rate is zero; one hour later, the same patient, awake or asleep or moving differently, and the rate is twelve an hour. This is the central practical problem of seizure detection and of clinical machine learning generally: the data you test on must come from a different time than the data you trained on, or the number you report is fiction. It is lesson 3's train-test rule again, with time as the thing to hold out.

Step 5: save the model

Everything the device needs is five weights, a bias, the standardisation constants, and the three decision parameters. Write them to a file.

import json
model = {"fs": fs, "weights": w.tolist(), "bias": float(b), "mean": mean.tolist(), "sd": sd.tolist(),
         "threshold": 0.5, "persistence": 3, "refractory": 30}
json.dump(model, open("model.json", "w"), indent=1)
print("saved model.json:", {k: (np.round(v, 2).tolist() if isinstance(v, list) else v) for k, v in model.items()})
# saved model.json

Building the device

A detector that runs on a laptop against a file is an analysis. A detector that runs on a small computer against electrodes on a head, unattended, for days, is a system, and most of what makes it one is not the classifier. This section describes a build from open parts and gives the code. The pipeline is the one you have just run; what changes is that the data arrive one sample at a time and never stop.

The detector as a loop. Everything up to the classifier is lessons 6 and 7, the classifier is lesson 3, and the decision logic is step 4. On a device this loop runs once a second, forever.
The detector as a loop. Everything up to the classifier is lessons 6 and 7, the classifier is lesson 3, and the decision logic is step 4. On a device this loop runs once a second, forever.

Hardware

  • Electrodes. Eight gold-cup electrodes with conductive paste, placed by the 10-20 system at the positions the bipolar chains need (FP1, F7, T7, P7, O1 and the right-side equivalents), or a dry-electrode headband for a less fussy and noisier setup. Paste electrodes on a child for days is how the hospital did it, and it is a skill.
  • Amplifier. The signal is microvolts; the board that turns it into numbers is the one part that is hard to improvise. Two open-hardware options use the same Texas Instruments ADS1299 chip that research amplifiers use: the OpenBCI Cyton, an 8-channel board that streams at 250 Hz to a USB dongle and is supported by the BrainFlow library, and the PiEEG, an 8-channel shield that plugs directly onto a Raspberry Pi and talks over SPI. Both are open-source designs with published schematics.
  • Computer. A Raspberry Pi 4 or 5. The detector loop needs a few milliseconds a second; the Pi spends the rest idle. Power it from a battery bank if the wearer moves.
  • Outputs. A buzzer on a GPIO pin for the room, and a push notification for a phone elsewhere; the free ntfy service does the latter with one HTTP request and no account.
What talks to what. The amplifier board turns microvolts into numbers and hands them to the Pi over USB or SPI; one process fills a ring buffer, the detector reads the last two seconds from it once a second, and alarms go to a buzzer, a phone and a log. A systemd unit keeps it running.
What talks to what. The amplifier board turns microvolts into numbers and hands them to the Pi over USB or SPI; one process fills a ring buffer, the detector reads the last two seconds from it once a second, and alarms go to a buzzer, a phone and a log. A systemd unit keeps it running.

The service

One file. It has two sources that look identical to the rest of the program: a FileSource that replays a recording at real speed (or faster) and a BrainFlowSource that reads a Cyton. The Detector holds the last two seconds in a ring buffer, computes the same five features as step 2 on them, applies the saved weights, and runs the persistence and refractory logic. The main loop steps the detector once per second of data, not of wall-clock time, which is the kind of detail that only matters when you run it for real: a laggy USB read must not turn into alarms with the wrong timestamp.

"""monitor.py: a seizure detector as a long-running service.

Run on a laptop with a recording:   python monitor.py --source file --file chbmit_chb01.zip
Run on a Raspberry Pi with a board:  python monitor.py --source brainflow --port /dev/ttyUSB0
"""
import argparse, json, time, logging, urllib.request
import numpy as np
from collections import deque
from scipy.signal import butter, filtfilt

log = logging.getLogger("monitor")


# ---------- sources: anything with .read() that returns the newest samples as channels x n ----------
class FileSource:
    """Replays a recording at real speed, so the rest of the program cannot tell it from a live board."""
    def __init__(self, path, fs):
        d = np.load(path); self.eeg = d["eeg"].astype(float); self.fs = fs; self.pos = 0; self.t = time.time()
    def read(self):
        due = int((time.time() - self.t) * self.fs) - self.pos              # how many samples have "arrived" since last call
        chunk = self.eeg[:, self.pos:self.pos + due]; self.pos += chunk.shape[1]
        return chunk


class BrainFlowSource:
    """An OpenBCI Cyton (or any BrainFlow board) over USB. pip install brainflow."""
    def __init__(self, port, channels=8):
        from brainflow.board_shim import BoardShim, BrainFlowInputParams, BoardIds
        params = BrainFlowInputParams(); params.serial_port = port
        self.board = BoardShim(BoardIds.CYTON_BOARD, params)
        self.rows = BoardShim.get_eeg_channels(BoardIds.CYTON_BOARD)[:channels]
        self.fs = BoardShim.get_sampling_rate(BoardIds.CYTON_BOARD)         # 250
        self.board.prepare_session(); self.board.start_stream()
    def read(self):
        return self.board.get_board_data()[self.rows]                         # everything since the last call, in uV


# ---------- the detector: the same five features and the weights trained in the article ----------
class Detector:
    def __init__(self, model, fs, threshold=0.5, persistence=3, refractory=30):
        self.w, self.b = np.array(model["weights"]), model["bias"]
        self.mean, self.sd = np.array(model["mean"]), np.array(model["sd"])
        self.fs = fs; self.buffer = deque(maxlen=int(2 * fs))                 # the last two seconds, per sample
        self.threshold, self.persistence, self.refractory = threshold, persistence, refractory
        self.streak, self.quiet_until = 0, 0

    def features(self, window):
        out = []
        for lo, hi in [(1, 4), (4, 8), (8, 13), (13, 30)]:
            b, a = butter(3, [lo / (self.fs / 2), hi / (self.fs / 2)], "band")
            out.append(np.log(filtfilt(b, a, window, axis=1).var(axis=1) + 1e-6).max())
        out.append(np.log(np.abs(np.diff(window, axis=1)).mean(axis=1)).max())
        return (np.array(out) - self.mean) / self.sd

    def step(self, now):
        """Called once a second. Returns the score, and True if an alarm should fire."""
        window = np.array(self.buffer).T                                      # channels x samples
        if window.shape[1] < self.buffer.maxlen:
            return None, False
        s = 1 / (1 + np.exp(-(self.features(window) @ self.w + self.b)))
        self.streak = self.streak + 1 if s > self.threshold else 0
        fire = self.streak >= self.persistence and now > self.quiet_until
        if fire:
            self.quiet_until = now + self.refractory
        return s, fire


def notify(topic, text):
    """Push to a phone with ntfy.sh (free, no account): subscribe to the topic in the ntfy app."""
    try:
        urllib.request.urlopen(urllib.request.Request(f"https://ntfy.sh/{topic}", data=text.encode(), headers={"Priority": "high"}), timeout=5)
    except Exception as e:
        log.warning("notification failed: %s", e)


def main():
    ap = argparse.ArgumentParser()
    ap.add_argument("--source", choices=["file", "brainflow"], default="file")
    ap.add_argument("--file", default="chbmit_chb01.zip"); ap.add_argument("--port", default="/dev/ttyUSB0")
    ap.add_argument("--model", default="model.json"); ap.add_argument("--topic", default=None)
    ap.add_argument("--speed", type=float, default=1.0, help="replay faster than real time (file source only)")
    args = ap.parse_args()
    logging.basicConfig(level=logging.INFO, format="%(asctime)s %(levelname)s %(message)s")

    model = json.load(open(args.model))
    source = FileSource(args.file, model["fs"] * args.speed) if args.source == "file" else BrainFlowSource(args.port)
    det = Detector(model, model["fs"], model["threshold"], model["persistence"], model["refractory"])
    log.info("running; source=%s", args.source)
    fs = model["fs"]
    n_samples = 0                                                             # count data, not wall-clock time
    while True:
        chunk = source.read()
        for col in chunk.T:
            det.buffer.append(col)
            n_samples += 1
            if n_samples % int(fs) == 0:                                      # another second of data has arrived
                now = n_samples / fs
                s, fire = det.step(now)
                if fire:
                    log.warning("ALARM at %d s (score %.2f)", now, s)
                    if args.topic: notify(args.topic, f"Possible seizure at {time.strftime('%H:%M:%S')}")
        if args.source == "file" and source.pos >= source.eeg.shape[1]:
            log.info("end of recording"); break
        time.sleep(0.05)


if __name__ == "__main__":
    main()

Try it against the article's file, at fifty times real speed:

pip install numpy scipy
python monitor.py --source file --file chbmit_chb01.zip --speed 50
# 15:06:49 WARNING ALARM at 247 s (score 0.96)
# 15:06:51 WARNING ALARM at 644 s (score 0.97)
# 15:06:53 WARNING ALARM at 1032 s (score 0.77)
# ... one alarm per seizure within seconds of onset, three false alarms in the last quarter hour

The replay produces the same alarms as step 4, which is the test that the streaming code is equivalent to the offline analysis. That equivalence is worth checking every time you change either, because the two versions drift apart in practice: the offline version sees the whole recording and filters it in one pass; the online one sees two seconds and must get the same answer from them.

Packaging it

On the Pi, the service needs to start when the power comes on, restart if it crashes, and leave a log. That is what systemd is for, and the whole configuration is this file, installed with sudo cp seizure-monitor.service /etc/systemd/system/ && sudo systemctl enable --now seizure-monitor. The log is journalctl -u seizure-monitor -f.

[Unit]
Description=Seizure monitor (research build, not a medical device)
After=network-online.target

[Service]
User=pi
WorkingDirectory=/home/pi/seizure-monitor
ExecStart=/home/pi/seizure-monitor/venv/bin/python monitor.py --source brainflow --port /dev/ttyUSB0 --topic my-secret-topic-name
Restart=always
RestartSec=5

[Install]
WantedBy=multi-user.target

The model file it loads is the one step 5 wrote, so retraining on new seizures is a matter of replacing one small file:

{
 "fs": 128.0,
 "weights": [
  1.2117927709917629,
  1.3755688993410002,
  0.4899623463951435,
  0.40485118358151967,
  0.42444341307141376
 ],
 "bias": -4.152963406295196,
 "mean": [
  7.391718334533521,
  5.870881518587343,
  4.750163443873513,
  3.9429359195743587,
  2.5285384028505717
 ],
 "sd": [
  1.2656932215114118,
  0.9868093864383984,
  0.9293513440481552,
  1.1640161259724342,
  0.5585923924815547
 ],
 "threshold": 0.5,
 "persistence": 3,
 "refractory": 30
}

Three things the production version of this would add, in order of importance. An impedance check at startup, because an electrode that has come loose produces huge slow artefacts that look exactly like what the delta weight rewards; the ADS1299 can measure electrode impedance and both boards expose it. Artefact rejection before the classifier: the two line-length spikes after seizure 1 in the features figure are movements, and a blink detector from the ICA primer or a simple amplitude ceiling would catch most of them. Storage of every alarm's EEG, so that a human can later decide which alarms were real, which is how the threshold and the model get better over time, and how you would ever know the device's true false-alarm rate on its wearer.

What you built, and what it is worth

  • A seizure detector from scalp EEG: five features, a logistic regression, and two rules, that catches all seven of one patient's focal seizures within about five seconds at three false alarms an hour.
  • An honest evaluation: tested by leaving each seizure out, with the metrics clinicians use, and with the false alarms traced to the recording the model had not seen the background of.
  • A clear statement of scope: focal seizures with a scalp signature, yes; silent focal seizures and brief jerks, no; seizure type, never.
  • A streaming service that gives the same answers as the analysis, packaged to run unattended on a Raspberry Pi with an open-hardware amplifier.

The gap between this and a product is not the classifier; published clinical detectors use features much like these. It is the thousand hours of recordings across many patients needed to know the false-alarm rate, the artefact handling, the electrodes that stay on for a week, the battery, the failure modes when a wire breaks at 3 a.m., and the regulatory evidence that the whole thing does more good than harm. That is years of work by a team, and it is why those products cost what they do. What you have is the understanding of how they work, which is the part that cannot be bought.

Exercises

  1. Move the threshold from 0.5 to 0.3 and to 0.7 and redo step 4. Plot sensitivity and false alarms per hour against the threshold: this curve is how detectors are compared in papers.
  2. Download another CHB-MIT patient (chb03 and chb05 are good) and run steps 1 to 4 on them with the same code. Then test chb01's model on them without retraining. How well does one child's detector transfer to another?
  3. Add a sixth feature: the ratio of power on the right chain to the left chain. For this patient every seizure starts on the right. Does the detector get faster? What would happen on a patient whose seizures start on the left?
  4. Harder: the false alarms live in the seizure-free quarter hour. Download the full chb01_01.edf (42 MB) and chb01_02.edf from PhysioNet, add two hours of seizure-free background to the training set, and measure the false-alarm rate again on a third hour. This is what fixing a detector actually looks like.

Data: CHB-MIT Scalp EEG Database (Shoeb 2009; Goldberger et al. 2000) via PhysioNet, patient chb01, Open Data Commons Attribution licence. Line length as a seizure feature: Esteller et al. 2001. The analysis runs in about fifteen seconds; the monitor replay at 50x takes a minute and a half.

Your first neuron in Python

By the end of this lesson you will have loaded a real recording from a mouse brain, looked at how neurons actually fire, measured a few things about them, and built a working model of a neuron from scratch. All of it in Python, all of it in about thirty lines of code you write yourself. No prior programming experience is assumed.

The data

We will use a recording from Steinmetz and colleagues (2019), who inserted Neuropixels probes into the brains of mice performing a visual decision task. One session, 734 neurons, 45 minutes. Neuromatch Academy packaged it as a single file. Download it here: steinmetz_session.npz (42 MB). Save it in the folder where you will run your code.

Each neuron in the file is represented by a list of the moments, in seconds, at which it fired. That is all a spike train is: a list of times. Everything else in this lesson is built on that idea.

Step 1: load it

import numpy as np
import matplotlib.pyplot as plt

data = np.load("steinmetz_session.npz", allow_pickle=True)
spike_times = data["spike_times"]
print(len(spike_times), "neurons")
# 734 neurons

np.load opens the file. spike_times is now a collection of 734 arrays, one per neuron. len counts them. If you see 734 neurons, everything is working.

Step 2: look at one neuron

Let’s pick a neuron that fires reasonably often, so there is something to see. We compute every neuron’s firing rate (spikes divided by recording length) and choose the one closest to six spikes per second.

duration = max(st.max() for st in spike_times)      # length of the recording, in seconds
print(f"recording length: {duration/60:.1f} minutes")

rates = np.array([len(st) / duration for st in spike_times])   # spikes per second, one per neuron
idx = int(np.argmin(np.abs(rates - 6)))                          # the neuron closest to 6 Hz
neuron = spike_times[idx]

print(f"neuron {idx}: {len(neuron)} spikes, mean rate {rates[idx]:.2f} Hz")
print(neuron[:5])
# recording length: 45.0 minutes
# neuron 303: 16177 spikes, mean rate 5.99 Hz
# [0.011      0.03373333 0.07536667 0.1061     0.17953333]

Read the last line. Neuron 303 fired at 11 milliseconds, then 34, then 75, then 106, then 180. Those are real electrical events in a real brain, and you are looking at them in a Python list. Notice the gaps are uneven. Hold that thought.

Step 3: a raster plot

The standard picture of spiking activity is a raster: one row per neuron, one tick per spike. Here are forty neurons during the first minute.

fig, ax = plt.subplots(figsize=(10, 5))
for i in range(40):
    st = spike_times[i]
    st = st[st < 60]                              # keep only the first 60 seconds
    ax.vlines(st, i + 0.5, i + 1.5, lw=0.7)       # one small vertical tick per spike
ax.set(xlabel="time (s)", ylabel="neuron", xlim=(0, 60))
plt.show()
Forty neurons, one minute. Some rows are dense, some are almost empty, and none of them tick like a clock.
Forty neurons, one minute. Some rows are dense, some are almost empty, and none of them tick like a clock.

Two things stand out. Neurons differ enormously from each other: row 14 is busy, row 5 fires a handful of times in a minute. And no neuron is regular. The spacing between spikes looks random. That irregularity is not noise in the recording. It is how neurons behave, and we will come back to it.

Step 4: firing rate over time

A list of spike times is exact but hard to reason about. Neuroscientists usually convert it into a rate: how many spikes per second, in each second. np.histogram does the counting.

bin_size = 1.0                                       # seconds
edges = np.arange(0, 600 + bin_size, bin_size)       # the first ten minutes, in 1-second bins
counts, _ = np.histogram(neuron, bins=edges)
rate = counts / bin_size                             # spikes per second in each bin

fig, ax = plt.subplots(figsize=(10, 3.4))
ax.plot(edges[:-1], rate, lw=1)
ax.set(xlabel="time (s)", ylabel="spikes / s")
plt.show()
Neuron 303, ten minutes, one-second bins. The rate wanders between silence and twenty spikes a second as the mouse sees stimuli, decides, and moves.
Neuron 303, ten minutes, one-second bins. The rate wanders between silence and twenty spikes a second as the mouse sees stimuli, decides, and moves.

This is the first genuinely useful analysis in neuroscience: the rate over time is what you would align to a stimulus to ask whether a neuron responds to it. Every decoding algorithm and every brain-computer interface starts from something like this line.

Step 5: the gaps between spikes

Now to that irregularity. The gap between one spike and the next is called the inter-spike interval, or ISI. np.diff computes all of them at once. Dividing their standard deviation by their mean gives the coefficient of variation, or CV: a single number for how irregular the neuron is. A perfect clock has CV 0. A completely random process has CV 1.

isi = np.diff(neuron)                          # gap between consecutive spikes, in seconds
cv = isi.std() / isi.mean()
print(f"mean ISI {isi.mean()*1000:.0f} ms, CV = {cv:.2f}")
# mean ISI 167 ms, CV = 2.60

A CV above 2. This neuron is more irregular than random, because it fires in bursts: clusters of very short intervals separated by long pauses. Cortical neurons commonly have a CV near or above 1. Keep this number in mind. In a moment we will build a neuron and ask whether it can match it.

Step 6: the whole population

We computed a rate for every neuron in step 2. Let’s look at all 734 of them.

print(f"median rate {np.median(rates):.2f} Hz, max {rates.max():.1f} Hz")
print(f"{np.mean(rates < 1)*100:.0f}% of neurons fire below 1 Hz")
# median rate 1.90 Hz, max 46.4 Hz
# 33% of neurons fire below 1 Hz

fig, ax = plt.subplots(figsize=(10, 3.8))
ax.hist(rates, bins=np.logspace(-2, 2, 40))
ax.set_xscale("log")
ax.set(xlabel="mean firing rate (Hz, log scale)", ylabel="number of neurons")
plt.show()
Firing rates of all 734 neurons on a logarithmic axis. Most are nearly silent. A few do most of the talking.
Firing rates of all 734 neurons on a logarithmic axis. Most are nearly silent. A few do most of the talking.

The distribution is roughly a bell curve on a log axis, which means it is heavily skewed on a normal one: a small number of neurons fire ten or twenty times more than the typical one. This lognormal pattern shows up across brain areas and species. Nobody fully understands why, and it is one of the sharpest differences between real neural networks and artificial ones, where every unit is active on every pass. You have just reproduced a real finding with six lines of code.

Step 7: build a neuron

Time to make one. The simplest model that behaves like a neuron is the leaky integrate-and-fire neuron, and it has three ideas in it:

  • Integrate. Input current pushes the membrane voltage up.
  • Leak. Left alone, the voltage drifts back toward a resting level, like a bucket with a hole.
  • Fire. If the voltage reaches a threshold, the neuron spikes and the voltage resets.
flowchart TD

A[Start at rest, -70 mV] --> B[Add the input current]

B --> C[Leak back toward rest]

C --> D{Voltage above -50 mV?}

D -- no --> B

D -- yes --> E[Record a spike, reset to -65 mV]

E --> B
The leaky integrate-and-fire loop. Every time step does exactly this, and that is the whole model.

That is the whole model. Here it is as a function. The numbers are typical for a cortical neuron: resting at -70 mV, threshold at -50 mV, and a time constant of 20 ms, which sets how fast the leak works.

def lif(current, dt=1e-3, tau=0.02, R=1e8, v_rest=-0.070, v_thresh=-0.050, v_reset=-0.065):
    """Leaky integrate-and-fire neuron.
    current: input current in amps, one value per time step of dt seconds.
    Returns the membrane voltage over time and the spike times."""
    v = np.full(len(current), v_rest)
    spikes = []
    for t in range(1, len(current)):
        dv = (-(v[t-1] - v_rest) + R * current[t-1]) / tau   # leak toward rest, push from input
        v[t] = v[t-1] + dv * dt
        if v[t] >= v_thresh:                                   # threshold crossed
            spikes.append(t * dt)                              # record the spike
            v[t] = v_reset                                     # and reset
    return v, np.array(spikes)

Feed it a constant current and see what it does.

dt = 1e-3
t = np.arange(0, 2.0, dt)                        # two seconds in 1 ms steps
steady = np.full(len(t), 0.25e-9)                # 0.25 nanoamps, constant
v_steady, s_steady = lif(steady)
print(len(s_steady), "spikes, CV =", round(np.diff(s_steady).std() / np.diff(s_steady).mean(), 2))
# 71 spikes, CV = 0.0

Seventy-one spikes in two seconds, perfectly evenly spaced. CV of zero. It is a metronome. Real neuron 303 had a CV of 2.6. Our model is missing something important.

Step 8: make it real

What is missing is that a real neuron does not receive a steady current. It receives thousands of tiny, randomly timed inputs from other neurons. Their sum fluctuates wildly. Let’s give our model a fluctuating input instead: a mean that on its own is not enough to reach threshold, plus a lot of noise. Now only the random upward swings make it fire.

rng = np.random.default_rng(0)
noisy = 0.14e-9 + 0.55e-9 * rng.standard_normal(len(t))     # mean below threshold, big fluctuations
v_noisy, s_noisy = lif(noisy)
print(len(s_noisy), "spikes, CV =", round(np.diff(s_noisy).std() / np.diff(s_noisy).mean(), 2))
# 35 spikes, CV = 0.98

fig, axes = plt.subplots(2, 1, figsize=(10, 5.6), sharex=True)
for ax, v, s, title in [(axes[0], v_steady, s_steady, "constant input"),
                        (axes[1], v_noisy, s_noisy, "noisy input")]:
    ax.plot(t, v * 1000, lw=0.9)
    ax.vlines(s, -50, -20, lw=0.9)                # draw each spike as a tall line
    ax.set(ylabel="membrane potential (mV)", title=title)
axes[1].set_xlabel("time (s)")
plt.show()
Same neuron, two inputs. Top: constant current, perfectly regular spikes. Bottom: noisy current, and the spike timing looks like the raster in step 3.
Same neuron, two inputs. Top: constant current, perfectly regular spikes. Bottom: noisy current, and the spike timing looks like the raster in step 3.
Loading the interactive figure…
The same neuron, live. Drag the input current up and down, add noise, and watch the firing change; the two buttons jump to the settings of steps 7 and 8.

The CV jumped from 0 to about 1 with one change to the input. That is a real result from the 1990s, argued over in the literature for years: cortical neurons are irregular not because they are sloppy, but because they operate in a regime where noise, not the average input, decides when they fire. You just found it yourself in a dozen lines.

Run the noisy model for a full minute and compare its interval distribution with neuron 303:

t_long = np.arange(0, 60.0, dt)
_, s_long = lif(0.14e-9 + 0.55e-9 * rng.standard_normal(len(t_long)))
isi_model = np.diff(s_long)

fig, ax = plt.subplots(figsize=(10, 3.8))
ax.hist(isi[isi < 0.5] * 1000, bins=50, density=True, alpha=0.6, label="real neuron 303")
ax.hist(isi_model[isi_model < 0.5] * 1000, bins=50, density=True, alpha=0.6, label="noisy LIF model")
ax.set(xlabel="inter-spike interval (ms)", ylabel="density")
ax.legend()
plt.show()
Interval distributions. The model has the right overall shape: many short gaps, a long tail. It is missing the bursts that push the real neuron's CV above 2, which is a topic for a later lesson.
Interval distributions. The model has the right overall shape: many short gaps, a long tail. It is missing the bursts that push the real neuron's CV above 2, which is a topic for a later lesson.

What you just did

  • Loaded a real neural recording and understood what a spike train is.
  • Made a raster plot, the most common figure in systems neuroscience.
  • Converted spike times into a firing rate, the starting point of every decoding method.
  • Measured irregularity with the coefficient of variation.
  • Reproduced a real finding: neural firing rates are lognormally distributed.
  • Built a leaky integrate-and-fire neuron from scratch and discovered why real neurons are irregular.

If you had never written Python before today, you have now used arrays, loops, functions, conditionals and plotting, all in service of a real scientific question. That is the way this site teaches everything.

Exercises

  1. Change bin_size in step 4 to 0.1 and to 10. What changes, and which is more useful?
  2. Find the neuron with the highest CV in the whole population. Plot its raster. What does its firing look like?
  3. In step 8, slowly increase the mean current from 0.14 nA toward 0.25 nA while keeping the noise. What happens to the CV, and why?
  4. Harder: the model's CV tops out near 1, but neuron 303 reaches 2.6 by firing in bursts. What would you have to add to the model to make it burst?

Data: Steinmetz, Zatka-Haas, Carandini & Harris, "Distributed coding of choice, action and engagement across the mouse brain", Nature 2019, CC-BY 4.0, via Neuromatch Academy. Complete code for this lesson is in the blocks above and runs top to bottom as a single script.