Every lesson so far started from spike times, or spike counts, that someone else had already extracted. This is where they come from. A Neuropixels probe does not record neurons; it records voltage, thirty thousand times a second on each of 384 electrode sites, and somewhere in that voltage are the spikes. Turning one into the other is called spike detection and spike sorting, and in this lesson you will do both from scratch on ten seconds of a real recording: look at the raw trace, filter it, find the spikes, cut out their waveforms, use the PCA you wrote in lesson 4 plus a ten-line k-means to split them into neurons, and then check your answer against what a professional pipeline found on the same data. This is the first lesson of the Signals track.
flowchart LR A[Raw voltage: samples x sites] --> B[High-pass filter] B --> C[Threshold: find the troughs] C --> D[Cut out waveforms] D --> E[PCA + k-means: sort into units] E --> F[Footprints and refractory check] F --> G[Compare with the real sorter]
The data
Ten seconds of raw voltage from eight neighbouring sites of a Neuropixels probe in a mouse, recorded by the International Brain Laboratory (session KS023, 10 December 2019, from the Carandini and Harris lab at UCL) and shared under CC-BY as part of their brain-wide map. The IBL publishes the complete recording, 26 GB for this probe alone, so I cut out a slice: ibl_ks023_raw.zip. Save it next to your code and do not unzip it. Alongside the voltage, the file carries the spike times that the IBL’s own spike sorter found in the same ten seconds, which we will keep sealed until the last step.
The eight sites sit two to a row, with rows 20 micrometres apart along the probe. That spacing matters: a neuron’s spike is an electrical event in the tissue, and a site a few tens of micrometres away sees the same event, smaller. We will use that.
Step 1: look at raw voltage
The probe stores each sample as a 16-bit integer, and a number in the file header says how many microvolts one integer step is worth. Converting is one multiplication. Then plot one second on all eight sites, stacked so they do not overlap.
import numpy as np
import matplotlib.pyplot as plt
dat = np.load("ibl_ks023_raw.npz")
raw = dat["raw"] # samples x channels, int16, exactly as the probe wrote it
fs = float(dat["fs"]) # 30000 samples per second
uv = float(dat["uv_per_bit"]) # 2.34 microvolts per integer step
chan = dat["channel_ids"] # which of the probe's 384 electrode sites these are
t = np.arange(raw.shape[0]) / fs # time of each sample, in seconds
v = raw.astype(float) * uv # the same data in microvolts
print(raw.shape, raw.dtype, "=", raw.shape[0] / fs, "seconds on", raw.shape[1], "channels")
print("sites:", chan)
print("first samples on site 333:", raw[:6, 5], "=", v[:6, 5].round(1), "uV")
# (300000, 8) int16 = 10.0 seconds on 8 channels
# sites: [328 329 330 331 332 333 334 335]
# first samples on site 333: [54 52 49 51 49 50] = [126.6 121.9 114.8 119.5 114.8 117.2] uV
fig, ax = plt.subplots(figsize=(10, 6))
for i in range(v.shape[1]):
ax.plot(t, v[:, i] + i * 400, lw=0.5) # stack the channels 400 uV apart
ax.set(xlim=(0, 1), xlabel="time (s)", yticks=np.arange(8) * 400, yticklabels=chan, ylabel="electrode site")
plt.show()
Three things are visible before any analysis. There is a slow wobble shared by every site, which is the local field potential, the summed activity of thousands of neurons, and the subject of the next lesson. There is a fuzz of fast noise. And there are sharp downward blips, a few hundred microvolts deep and a millisecond wide, that appear on some sites and not others. Those are the spikes. A spike is a negative deflection at an extracellular electrode because the sodium rushing into the neuron is leaving the outside, where the electrode is.
Step 2: filter
The wobble and the spikes live at different speeds. The wobble changes over tens of milliseconds, which is to say below about 100 Hz; a spike is over in one millisecond, which is energy above 300 Hz. A high-pass filter keeps the fast part and discards the slow, and this one is a standard Butterworth design from SciPy. Two details: the cutoff is given as a fraction of half the sampling rate, which is how the function wants it, and filtfilt runs the filter forwards and then backwards, so the result lines up in time with the original. A filter run only forwards shifts everything slightly later, and in spike work a tenth of a millisecond matters.
from scipy.signal import butter, filtfilt
b, a = butter(3, 300 / (fs / 2), "high") # a high-pass filter: keep everything above 300 Hz
clean = filtfilt(b, a, v, axis=0) # filtfilt runs it forwards then backwards, so nothing shifts in time
fig, axes = plt.subplots(2, 1, figsize=(10, 5), sharex=True)
for ax, x, name in zip(axes, [v, clean], ["raw", "high-pass filtered"]):
ax.plot(t, x[:, 5], lw=0.6)
ax.set(ylabel="uV", title=name, xlim=(0.40, 0.50))
axes[-1].set_xlabel("time (s)")
plt.show()
Filtering is a large topic and the next lesson treats it properly, including what the filter threw away. For now the picture says enough: the baseline is flat, the spikes are untouched.
Step 3: find the spikes
A spike is a sample that goes much deeper than the noise. The question is how to measure the noise when the trace is full of spikes: the standard deviation would count them and come out too large. The trick, due to Quiroga and colleagues in 2004, is to take the median of the absolute voltage and divide by 0.6745, which gives the standard deviation of the noise alone, because the median hardly notices the rare large values. Threshold at five times that, find every downward crossing, locate the deepest point within a millisecond of it, and refuse a second detection within a millisecond of the first.
def detect_spikes(x, fs, n_sigma=5, refractory=0.001):
"""Times (in samples) of the troughs of negative deflections more than n_sigma noise levels deep."""
sigma = np.median(np.abs(x)) / 0.6745 # the noise level, estimated in a way spikes cannot inflate
threshold = -n_sigma * sigma
crossings = np.where((x[1:] < threshold) & (x[:-1] >= threshold))[0] + 1
peaks, last = [], -np.inf
for c in crossings:
p = c + np.argmin(x[c:c + int(0.001 * fs)]) # the deepest point within 1 ms of the crossing
if p - last > refractory * fs: # ignore a second crossing of the same spike
peaks.append(p)
last = p
return np.array(peaks), threshold
for i, site in enumerate(chan):
peaks, threshold = detect_spikes(clean[:, i], fs)
print(f"site {site}: noise {-threshold / 5:4.1f} uV, threshold {threshold:6.1f} uV, {len(peaks):3d} spikes = {len(peaks) / 10:4.1f} per second")
# site 328: noise 14.6 uV, threshold -73.0 uV, 363 spikes = 36.3 per second
# site 329: noise 15.4 uV, threshold -76.8 uV, 313 spikes = 31.3 per second
# site 330: noise 15.4 uV, threshold -76.9 uV, 481 spikes = 48.1 per second
# site 331: noise 15.3 uV, threshold -76.7 uV, 286 spikes = 28.6 per second
# site 332: noise 14.7 uV, threshold -73.5 uV, 409 spikes = 40.9 per second
# site 333: noise 14.8 uV, threshold -73.9 uV, 500 spikes = 50.0 per second
# site 334: noise 14.5 uV, threshold -72.4 uV, 401 spikes = 40.1 per second
# site 335: noise 14.9 uV, threshold -74.5 uV, 379 spikes = 37.9 per secondThirty to fifty detections per second on every site. In lesson 1 the typical neuron fired twice a second, so these are not single neurons. A site hears every neuron within something like fifty to a hundred micrometres, and in a dense region that is dozens. Detection finds spikes; it does not tell you whose they are. That is the job of sorting, and the rest of the lesson.
Step 4: cut out the waveforms
For each detected trough on site 333, cut out a window from one millisecond before to two after, on all eight sites at once. The result is a three-dimensional array, spikes by samples by sites, and as in lesson 2 the shape is half the lesson.
x = clean[:, 5] # site 333
peaks, threshold = detect_spikes(x, fs)
before, after = int(0.001 * fs), int(0.002 * fs) # 1 ms before the trough, 2 ms after
peaks = peaks[(peaks >= before) & (peaks < len(x) - after)]
waves = np.stack([clean[p - before:p + after] for p in peaks]) # spikes x samples x channels
tw = (np.arange(-before, after) / fs) * 1000 # waveform time axis, in ms
print(waves.shape, "spikes x samples x channels")
# (500, 90, 8) spikes x samples x channels
fig, axes = plt.subplots(1, 2, figsize=(11, 4.5))
axes[0].plot(tw, waves[:150, :, 5].T, color="gray", lw=0.4, alpha=0.5) # 150 spikes, site 333 only
axes[0].plot(tw, waves[:, :, 5].mean(axis=0), lw=2)
axes[0].axhline(threshold, color="orange", ls="--")
axes[0].set(xlabel="time from trough (ms)", ylabel="uV", title="site 333: 150 spikes and the average")
big = np.argmin(waves[:, before, 5]) # the single largest spike
for i in range(8):
axes[1].plot(tw, waves[big, :, i] + i * 150, lw=1) # the same spike on all eight sites
axes[1].set(xlabel="time from trough (ms)", yticks=np.arange(8) * 150, yticklabels=chan, title="one spike, eight sites")
plt.show()
The left panel is the reason sorting exists. These were all caught by the same threshold on the same site, but they are clearly not all the same thing: some troughs reach 100 microvolts, others 200, and the recoveries differ. The right panel is the clue that will separate them. A spike appears on several sites at once with an amplitude that falls off with distance, so each neuron has a spatial footprint on the probe set by where its body sits. Two neurons with similar waveforms on site 333 can still have different footprints.
Step 5: sort
Each spike is now 90 samples times 8 sites, 720 numbers. Finding groups in 720 dimensions is the problem lesson 4 solved: run PCA, keep a few components, and look. Then split the cloud with k-means, the simplest clustering algorithm there is: guess k centres, assign every point to its nearest centre, move each centre to the middle of its points, repeat until nothing moves. Both functions are below; pca is the one from lesson 4 without the sign-fixing loop.
def pca(data):
mean = data.mean(axis=0)
centered = data - mean
cov = centered.T @ centered / (len(data) - 1)
variance, components = np.linalg.eigh(cov)
order = np.argsort(variance)[::-1]
return components[:, order], variance[order], mean
def kmeans(points, k, seed=0):
"""Split points into k groups by repeatedly assigning each to the nearest centre and moving the centres."""
rng = np.random.default_rng(seed)
centres = points[rng.choice(len(points), k, replace=False)] # start from k random points
for _ in range(100):
distance = ((points[:, None, :] - centres[None, :, :]) ** 2).sum(axis=2) # every point to every centre
group = distance.argmin(axis=1) # nearest centre wins
new = np.array([points[group == j].mean(axis=0) for j in range(k)])
if np.allclose(new, centres): # nothing moved: done
break
centres = new
return group
features = waves.reshape(len(waves), -1) # each spike: 90 samples x 8 sites = 720 numbers
components, variance, mean = pca(features)
scores = (features - mean) @ components[:, :3] # keep three
print("variance in the first three components:", (variance[:3] / variance.sum()).round(2))
k = 5
group = kmeans(scores, k)
for j in range(k):
print(f"cluster {j}: {np.sum(group == j):3d} spikes, trough on site 333 {waves[group == j, before, 5].mean():6.0f} uV")
# variance in the first three components: [0.23 0.09 0.08]
# cluster 0: 145 spikes, trough on site 333 -89 uV
# cluster 1: 83 spikes, trough on site 333 -89 uV
# cluster 2: 109 spikes, trough on site 333 -123 uV
# cluster 3: 65 spikes, trough on site 333 -177 uV
# cluster 4: 98 spikes, trough on site 333 -157 uV
fig, ax = plt.subplots(figsize=(7, 6))
for j in range(k):
ax.scatter(scores[group == j, 0], scores[group == j, 1], s=10, alpha=0.7, label=f"cluster {j}")
ax.set(xlabel="component 1", ylabel="component 2")
ax.legend()
plt.show()
Two islands and a crowded mainland. The islands, clusters 3 and 4, have the deepest troughs and are the easy neurons. The mainland is three clusters that k-means split in the third component, which this plot cannot show. Whether those splits are real is exactly what the next two steps check. I chose five clusters after trying four and seeing two of the big units merged; choosing k is a judgement, and the exercises come back to it.
Step 6: footprints and the refractory test
Two checks that need no answer key. First, the footprint: average each cluster’s waveform on every site. If a cluster is one neuron, its average should peak on one site and fall off smoothly around it, and different clusters should peak in different places or with different shapes. Second, the refractory period: a neuron cannot fire twice within about two milliseconds, so if a cluster’s spike train has many intervals shorter than that, it contains more than one neuron.
fig, axes = plt.subplots(1, k, figsize=(13, 4.5), sharey=True)
for j, ax in enumerate(axes):
footprint = waves[group == j].mean(axis=0) # samples x channels: this cluster's average on every site
for i in range(8):
ax.plot(tw, footprint[:, i] + i * 100, lw=1.2)
ax.set(title=f"cluster {j}, n = {np.sum(group == j)}", xlabel="ms")
axes[0].set(yticks=np.arange(8) * 100, yticklabels=chan)
plt.show()
for j in range(k):
times = np.sort(peaks[group == j]) / fs
isi = np.diff(times)
print(f"cluster {j}: {np.mean(isi < 0.002):.1%} of intervals shorter than 2 ms")
# cluster 0: 0.0% of intervals shorter than 2 ms
# cluster 1: 2.4% of intervals shorter than 2 ms
# cluster 2: 0.9% of intervals shorter than 2 ms
# cluster 3: 0.0% of intervals shorter than 2 ms
# cluster 4: 0.0% of intervals shorter than 2 ms
Clusters 2, 3 and 4 are neurons. Each has a clean footprint with a single peak, on three different sites, and none of them has a meaningful number of refractory violations. Cluster 4 is the clearest case of why the footprint matters: its trough is on site 329, two rows below the site we detected it on, and on site 333 alone its waveform looks much like cluster 2’s. Cluster 0 passes the refractory test and peaks on 333, but with a shallower trough than cluster 2 and a broader shape; it is probably mostly one neuron. Cluster 1 has 2.4 percent of its intervals under two milliseconds, and its footprint looks like a blend of its neighbours; it is where the leftovers went.
Step 7: open the answer key
The file also holds what the IBL’s pipeline found: a version of Kilosort, run on all 384 sites of the full recording, followed by curation. For each spike we detected, look for a sorter spike within half a millisecond, and tally which sorter unit each of our clusters mostly corresponds to.
sorter_t = dat["sorter_samples"] # the professional pipeline's spike times, in samples
sorter_unit = dat["sorter_clusters"] # and which unit it assigned each one to
def nearest_sorter_unit(p, tolerance=15):
"""The unit the sorter assigned to the spike nearest to sample p, if there is one within half a millisecond."""
i = np.argmin(np.abs(sorter_t - p))
return sorter_unit[i] if abs(sorter_t[i] - p) <= tolerance else -1
sorter_label = np.array([nearest_sorter_unit(p) for p in peaks])
print("detected spikes the sorter also found:", f"{np.mean(sorter_label >= 0):.0%}")
for j in range(k):
units, counts = np.unique(sorter_label[group == j], return_counts=True)
top = np.argsort(counts)[::-1][:3]
print(f"cluster {j}: " + ", ".join(f"unit {u} x{c}" if u >= 0 else f"unmatched x{c}" for u, c in zip(units[top], counts[top])))
ids, site_of = dat["sorter_cluster_ids"], dat["sorter_cluster_channel"]
for u in [1302, 1325, 1310, 1327]:
print(f"unit {u}: peak site {site_of[ids == u][0]}, {np.sum(sorter_unit == u)} spikes in these 10 seconds")
# detected spikes the sorter also found: 93%
# cluster 0: unit 1327 x84, unmatched x25, unit 1325 x10
# cluster 1: unit 1333 x27, unit 1327 x18, unmatched x11
# cluster 2: unit 1325 x88, unit 1327 x7, unit 1336 x4
# cluster 3: unit 1310 x57, unit 1322 x2, unit 1325 x2
# cluster 4: unit 1302 x95, unit 1625 x1, unit 1333 x1
# unit 1302: peak site 329, 97 spikes in these 10 seconds
# unit 1325: peak site 333, 105 spikes in these 10 seconds
# unit 1310: peak site 331, 61 spikes in these 10 seconds
# unit 1327: peak site 333, 131 spikes in these 10 secondsOf the spikes we detected, 93 percent are spikes the sorter also found. Cluster 4 is unit 1302: 95 of its 98 spikes, and the sorter agrees it peaks on site 329. Cluster 3 is unit 1310 on site 331, cluster 2 is unit 1325 on site 333. Three neurons found correctly, from one site, with a threshold, PCA and k-means. Cluster 0 is mostly unit 1327, the fourth neuron on site 333, but 25 of its spikes have no match in the sorter’s output for the units near these sites. Cluster 1 is a mixture of unit 1333 from site 335 with leftovers of 1327, as its refractory violations suggested. And unit 1327 itself has 131 spikes in these ten seconds, of which we found 84 on site 333 and labelled another 18 as cluster 1: the smaller the spike, the more of them fall below the threshold or into the wrong cluster.
That is a fair summary of the difference between what you just wrote and a real spike sorter. Kilosort uses every site, not one, so a neuron is detected wherever it is loudest; it learns a template for each neuron and subtracts it, so two spikes that overlap in time can both be recovered; and it tracks the slow drift of the brain past the probe over an hour, which turns one neuron into two clusters in a naive sort. Those three things are most of what the sorting literature since 2015 has been about. The core idea, though, is the one you implemented: a spike is a point in waveform space, and neurons are clusters.
What you just did
- Read raw extracellular voltage straight from a Neuropixels probe and converted it to microvolts.
- High-pass filtered it with SciPy, zero-phase, and saw the spikes stand clear of the slow activity.
- Estimated the noise level robustly and detected spikes by threshold, with a refractory rule.
- Cut out multi-site waveforms and saw that a spike has a spatial footprint on the probe.
- Sorted spikes into neurons with your own PCA and your own k-means, and checked the result with footprints and refractory violations.
- Compared your sort with a professional pipeline and found three neurons exactly, and learned where the simple method fails and why.
Exercises
- Change the threshold to 4 and then 6 noise levels. How does the spike count change on site 333, and what happens to the islands in the component plot?
- Run k-means with k = 4. Which two sorter units end up in the same cluster, and what do their footprints have in common? Then try k = 6 and see what splits.
- A common step in real pipelines is to subtract the median across all sites at each sample, which removes noise shared by the whole probe. Try it with these eight sites:
clean = clean - np.median(clean, axis=1, keepdims=True). The noise drops, but so does something else. Look at figure 3’s right panel and explain why a median over eight neighbouring sites is a bad idea, and why a median over 384 is fine. - Harder: detect spikes on every site rather than just 333, then merge detections that fall within half a millisecond of each other on neighbouring sites, keeping the one with the deepest trough. How many distinct spikes are there in ten seconds, and how many neurons can you now sort out?
Data: International Brain Laboratory, brain-wide map of neural activity during complex behaviour, CC-BY 4.0; session KS023 2019-12-10 (cortexlab), probe 00, sites 328 to 335, seconds 1500 to 1510 of the recording, with the spike times from the IBL’s pykilosort pipeline for units peaking nearby. The code in this lesson runs top to bottom in a few seconds.
