In lesson 1 you looked at neurons firing on their own time. This lesson asks the question that most of systems neuroscience is built on: does this neuron respond to something in the world? You will align hundreds of trials to the moment a stimulus appeared, average them into the field’s most important plot, measure how the response scales with the stimulus, test whether it is real, and then run the same test on every neuron in seven brain areas to see which parts of the brain are listening.
flowchart LR A[Spike counts: neurons x trials x bins] --> B[Pick a neuron, trials and time windows] B --> C[Raster: every trial] B --> D[PSTH: average over trials] D --> E[Tuning curve: rate vs contrast] B --> F[Permutation test: is it real?] F --> G[Every neuron, every area]
The data
Same experiment as lesson 1, from Steinmetz and colleagues (2019), but a different session and a different shape. In the task, a mouse sat in front of two screens. On each trial a striped pattern appeared on the left screen, the right screen, both, or neither, at one of four contrast levels, and the mouse turned a wheel to indicate which side was brighter. The probes were in the left hemisphere, which processes the right visual field, so we will use the contrast of the right stimulus.
To keep the download small, I extracted one session into a 2.4 MB file: steinmetz_session11.zip. Save it next to your code. Do not unzip it: NumPy reads it directly. The full dataset, with all thirteen sessions, is available from Neuromatch Academy (88 MB).
The key difference from lesson 1 is that the spikes have already been cut into trials and counted in 10-millisecond bins. Instead of a list of spike times per neuron, we have a three-dimensional array: neurons by trials by time bins. Every trial is a 2.5-second window and the stimulus appears exactly 0.5 seconds in. This is the standard shape of trial-based neural data, and getting comfortable with it is half the lesson.
Step 1: load and orient yourself
import numpy as np
import matplotlib.pyplot as plt
dat = np.load("steinmetz_session11.zip", allow_pickle=True)
spks = dat["spks"] # neurons x trials x time bins: spike counts in 10 ms bins
area = dat["brain_area"] # one label per neuron
contrast = dat["contrast_right"] # stimulus contrast on the right screen, one value per trial
dt = float(dat["bin_size"]) # 0.01 seconds
t = np.arange(spks.shape[2]) * dt - 0.5 # time of each bin relative to stimulus onset
print(spks.shape, "neurons x trials x bins")
print("areas:", {a: int(np.sum(area == a)) for a in np.unique(area)})
print("contrast levels:", np.unique(contrast))
# (698, 340, 250) neurons x trials x bins
# areas: {'ACA': 16, 'CA1': 50, 'DG': 65, 'LGd': 11, 'LH': 18, 'MD': 126, 'MOs': 6, 'PL': 56, 'SUB': 105, 'VISam': 79, 'VISp': 66, 'root': 100}
# contrast levels: [0. 0.25 0.5 1. ]698 neurons, 340 trials, 250 bins. The area codes are from the Allen Brain Atlas: VISp is primary visual cortex, VISam a higher visual area, LGd the visual thalamus, CA1 and DG are hippocampus, MOs is secondary motor cortex, ACA and PL are prefrontal. The line that builds t is worth reading twice: it turns bin numbers into seconds and shifts them so that zero is the moment the stimulus appeared.
Step 2: find a neuron worth looking at
With 698 neurons we need a way to choose. We will define a response as the firing rate in a window shortly after the stimulus minus the rate just before it, on high-contrast trials, and pick the visual cortex neuron with the biggest one. Boolean arrays do the selecting: pre and post pick out time bins, high picks out trials, visp picks out neurons.
visp = np.where(area == "VISp")[0] # indices of primary visual cortex neurons
pre = (t >= -0.3) & (t < 0) # 300 ms before the stimulus
post = (t >= 0.05) & (t < 0.25) # 50 to 250 ms after it
high = contrast == 1.0 # full-contrast trials
resp = spks[:, high][:, :, post].mean(axis=(1, 2)) - spks[:, high][:, :, pre].mean(axis=(1, 2))
best = visp[np.argmax(resp[visp])]
print(f"most responsive VISp neuron: {best}")
print("one trial, spike counts per 10 ms bin:", spks[best, 3, 45:60])
# most responsive VISp neuron: 141
# one trial, spike counts per 10 ms bin: [0 0 0 0 0 0 0 0 0 0 1 0 1 1 1]Look at that last line. Bins 45 to 49 are the 50 milliseconds before the stimulus: silence. Bin 50 is the stimulus onset. Five bins later, the neuron starts firing. On a single trial that is suggestive. The rest of the lesson is about turning suggestive into certain.
Step 3: every trial at once
A raster plot with trials as rows, sorted so that the zero-contrast trials come first and the full-contrast trials last. np.nonzero finds every bin with at least one spike and returns its row and column.
order = np.argsort(contrast, kind="stable") # trial indices, sorted by contrast
fig, ax = plt.subplots(figsize=(10, 5))
rows, cols = np.nonzero(spks[best][order]) # (trial, bin) of every spike, in sorted order
ax.scatter(t[cols], rows, s=2, marker="|")
ax.axvline(0, color="orange") # stimulus onset
ax.set(xlabel="time from stimulus onset (s)", ylabel="trial (sorted by contrast)")
ax.invert_yaxis()
plt.show()
You can read the result off the picture before computing anything. Blank-screen trials at the top: nothing happens at time zero. Full-contrast trials at the bottom: a wall of spikes about 50 milliseconds after onset, on nearly every trial. In between, the response gets stronger as contrast increases. This is what a visual neuron looks like.
Step 4: the PSTH
The peri-stimulus time histogram is the raster averaged over trials: for each time bin, the mean spike count across trials, divided by the bin width to give spikes per second. It is the single most common plot in the field. We add a five-bin running average to smooth it, using np.convolve.
def psth(counts, dt, smooth=5):
"""Average spike count per bin across trials, in spikes/s, lightly smoothed."""
rate = counts.mean(axis=0) / dt
kernel = np.ones(smooth) / smooth
return np.convolve(rate, kernel, mode="same")
fig, ax = plt.subplots(figsize=(10, 3.8))
for lvl in [0.0, 0.25, 0.5, 1.0]:
ax.plot(t, psth(spks[best, contrast == lvl], dt), label=f"contrast {lvl}")
ax.axvline(0, color="orange")
ax.set(xlabel="time from stimulus onset (s)", ylabel="spikes / s", xlim=(-0.3, 1.0))
ax.legend()
plt.show()
Three things to notice. The response begins about 40 milliseconds after the stimulus, which is roughly how long it takes light hitting the retina to reach primary visual cortex in a mouse. The peak grows with contrast. And at zero contrast the trace is flat: this neuron does nothing when there is nothing on the right screen, even though the mouse is still doing the task.
Step 5: how much does it care? A tuning curve
To turn the picture into numbers, compute the rate in the response window on each trial, then the mean and the standard error for each contrast level. The standard error tells you how much to trust each mean.
levels = np.unique(contrast)
means, sems = [], []
for lvl in levels:
r = spks[best, contrast == lvl][:, post].sum(axis=1) / (post.sum() * dt) # rate per trial in the window
means.append(r.mean())
sems.append(r.std(ddof=1) / np.sqrt(len(r)))
print(dict(zip(levels, np.round(means, 1))))
# {0.0: 7.6, 0.25: 35.3, 0.5: 50.5, 1.0: 66.3}
fig, ax = plt.subplots(figsize=(6, 3.8))
ax.errorbar(levels, means, yerr=sems, marker="o", capsize=3)
ax.set(xlabel="stimulus contrast", ylabel="spikes / s, 50 to 250 ms after onset")
plt.show()
Step 6: is it real? A permutation test
A neuron with a big response on the plot could, in principle, be a fluke. The honest way to check is to ask how often chance alone would produce a difference this large. We do that by shuffling: pool the before and after rates, shuffle them so the labels are meaningless, recompute the difference, and repeat thousands of times. The fraction of shuffles that beat the real difference is the p-value. No formulas, no assumptions about the distribution, and you can read exactly what it is doing.
def permutation_p(before, after, n=5000, seed=0):
"""How often does shuffling the before/after labels give a difference as large as the real one?"""
rng = np.random.default_rng(seed)
observed = after.mean() - before.mean()
pooled = np.concatenate([before, after])
count = 0
for _ in range(n):
rng.shuffle(pooled)
shuffled_diff = pooled[len(before):].mean() - pooled[:len(before)].mean()
if shuffled_diff >= observed:
count += 1
return count / n
trials = spks[best, high]
before = trials[:, pre].sum(axis=1) / (pre.sum() * dt)
after = trials[:, post].sum(axis=1) / (post.sum() * dt)
print(f"before {before.mean():.1f} Hz, after {after.mean():.1f} Hz, p = {permutation_p(before, after):.4f}")
# before 5.1 Hz, after 66.3 Hz, p = 0.0000Not one of five thousand shuffles came close. For this neuron the answer to the lesson’s question is a definite yes. Now we have a test we can run on anything.
Step 7: which parts of the brain are listening?
Wrap the test in a function that takes a list of neurons and returns the fraction that respond, then run it on every area with a decent number of neurons. This takes a minute or two, because it is running a thousand shuffles for each of a few hundred neurons. That is fine. Real analysis takes time.
def fraction_responsive(idx, alpha=0.01):
hits = 0
for i in idx:
b = spks[i, high][:, pre].sum(axis=1)
a = spks[i, high][:, post].sum(axis=1)
if permutation_p(b, a, n=1000) < alpha:
hits += 1
return hits / len(idx)
for a in ["VISp", "VISam", "LGd", "CA1", "DG", "MOs", "ACA"]:
idx = np.where(area == a)[0]
print(f"{a:5s} {len(idx):3d} neurons, {fraction_responsive(idx)*100:4.0f}% respond to the stimulus")
# VISp 66 neurons, 47% respond to the stimulus
# VISam 79 neurons, 16% respond to the stimulus
# LGd 11 neurons, 9% respond to the stimulus
# CA1 50 neurons, 0% respond to the stimulus
# DG 65 neurons, 2% respond to the stimulus
# MOs 6 neurons, 17% respond to the stimulus
# ACA 16 neurons, 0% respond to the stimulus
Half of primary visual cortex responds within a quarter of a second of a stimulus in its visual field. The hippocampus, which sits a few millimetres away and is busy with memory and space, does not respond at all. Secondary motor cortex shows a weak, later bump, and if you look at the left panel you will see it: that is the mouse starting to move, not the stimulus itself. The visual thalamus number is low, but there are only eleven neurons there, and eleven is not enough to say much. Small samples are a fact of life in this field, and the first thing to check whenever a number surprises you.
You have just done, on real data, the analysis that a large fraction of systems neuroscience papers are built on: align, average, quantify, test, and compare across areas.
What you just did
- Worked with trial-structured data as a three-dimensional array and selected from it with Boolean masks.
- Made a trial-sorted raster and a peri-stimulus time histogram, the two standard plots of stimulus-evoked activity.
- Measured a tuning curve with error bars.
- Wrote a permutation test from scratch and understood every line of it.
- Applied it to hundreds of neurons and mapped which brain areas respond to a visual stimulus.
Exercises
- Repeat steps 3 and 4 using
contrast_leftinstead. The probes are in the left hemisphere. What do you expect, and what do you get? - Neuron 141 has a second, broader hump around 250 ms. Sort the raster by
dat["response"](which way the mouse turned the wheel: -1, 0 or 1) instead of by contrast. What is the second hump about? - Measure the response latency: the first time bin after onset where the high-contrast PSTH exceeds the pre-stimulus mean by three standard deviations. Do it for every VISp neuron and plot the distribution.
- Harder: the permutation test treats every trial as independent. Trials early and late in a session can differ, because the mouse gets tired. Split the trials into first half and second half and check whether the tuning curve changes.
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. The single-session extract hosted here contains the spike counts, brain areas, contrasts and responses from session 11 (mouse Lederberg, 5 December 2017) and nothing else.
Keep reading
The shape of a population: your first PCA
Foundations, lesson 4. Sixty-six neurons is a sixty-six-dimensional space. Write principal component analysis from scratch, find out how many dimensions visual cortex…
What did the mouse see? Your first neural decoder
Foundations, lesson 3. Turn 66 visual cortex neurons into a table, write logistic regression from scratch, and decode whether a stimulus was…
Your first neuron in Python
Foundations, lesson 1. Load a real recording of 734 neurons from a mouse brain, plot it, measure it, and build a working…