In lesson 3 a decoder read the stimulus out of 66 neurons at once, and you never saw what it was looking at. Sixty-six neurons means sixty-six numbers per trial, a point in a sixty-six-dimensional space, and nobody can picture that. This lesson is about the standard way of squashing that space down to two or three axes so that you can. The method is principal component analysis, PCA, and you will write it in ten lines, apply it to two neurons where you can see exactly what it does, then to the whole population, then to the population’s activity over time, where the trial becomes a path through the space the neurons define. Along the way you will find out how many dimensions primary visual cortex actually uses, and fall into the single most common PCA trap so that you recognise it next time.

flowchart LR

A[Trials x neurons table] --> B[Two neurons: see the directions]

B --> C[Write PCA: centre, covariance, eigenvectors]

C --> D[All 66 neurons: the loud-neuron trap]

D --> E[Standardise: how many dimensions?]

E --> F[Every trial as a point]

F --> G[Every moment as a point: the trajectory]

Step 1: two neurons, one picture

Start where you can see everything. Take the trials-by-neurons table from lesson 3 and keep just two columns: neuron 141, the one we have followed since lesson 2, and neuron 190, which the decoder gave the second largest positive weight. Plot every trial as a point.

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
area = dat["brain_area"]
contrast = dat["contrast_right"]
dt = float(dat["bin_size"])
t = np.arange(spks.shape[2]) * dt - 0.5

post = (t >= 0.05) & (t < 0.25)
visp = np.where(area == "VISp")[0]
X = spks[visp][:, :, post].sum(axis=2).T  # trials x neurons, exactly as in lesson 3
y = (contrast > 0).astype(int)

a, b = list(visp).index(141), list(visp).index(190)   # columns of X for neurons 141 and 190
print(f"correlation between neurons 141 and 190: {np.corrcoef(X[:, a], X[:, b])[0, 1]:.2f}")
# correlation between neurons 141 and 190: 0.79

fig, ax = plt.subplots(figsize=(6, 5.5))
ax.scatter(X[y == 0, a], X[y == 0, b], s=14, alpha=0.6, label="blank screen")
ax.scatter(X[y == 1, a], X[y == 1, b], s=14, alpha=0.6, label="stimulus")
ax.set(xlabel="neuron 141 spikes", ylabel="neuron 190 spikes")
ax.legend()
plt.show()

The two neurons are strongly correlated, 0.79, so the cloud is an elongated ellipse. Now the key idea. This cloud has a long axis and a short axis, and the long axis is not “neuron 141” or “neuron 190”. It is a mixture, roughly 0.8 of neuron 141 plus 0.6 of neuron 190. If you had to describe every trial with one number instead of two, position along that long axis is the number that loses the least. That axis is the first principal component. The short axis, at right angles to it, is the second. PCA is nothing more than finding those axes, in any number of dimensions, and reporting how much of the cloud’s spread each one carries.

Step 2: write PCA

Ten lines, using one function from NumPy’s linear algebra module that you have not met before. First move the cloud so its centre is at zero. Then compute the covariance matrix: one row and one column per neuron, where each entry says whether two neurons go up and down together across trials (positive), in opposition (negative), or independently (zero). The diagonal holds each neuron’s own variance. The directions of the cloud are the eigenvectors of this matrix, and the variance along each direction is the matching eigenvalue. You do not need to know how np.linalg.eigh finds them, only that it does, and that it returns them smallest first, so we reverse the order.

def pca(data):
    """Principal component analysis of a samples x features array.
    Returns the components (one per column, sorted by variance), the variance along each, and the mean."""
    mean = data.mean(axis=0)
    centered = data - mean                                  # put the cloud's centre at zero
    cov = centered.T @ centered / (len(data) - 1)           # features x features covariance matrix
    variance, components = np.linalg.eigh(cov)              # directions of the cloud and the variance along each
    order = np.argsort(variance)[::-1]                      # eigh returns smallest first; we want largest first
    variance, components = variance[order], components[:, order]
    for j in range(components.shape[1]):                    # the sign of a component is arbitrary:
        if components[np.abs(components[:, j]).argmax(), j] < 0:
            components[:, j] *= -1                          # make its biggest entry positive, for consistent plots
    return components, variance, mean

two = X[:, [a, b]]
components, variance, mean = pca(two)
print("component 1:", components[:, 0].round(2), f"carries {variance[0] / variance.sum():.0%} of the variance")
print("component 2:", components[:, 1].round(2), f"carries {variance[1] / variance.sum():.0%} of the variance")
# component 1: [0.8 0.6] carries 90% of the variance
# component 2: [-0.6  0.8] carries 10% of the variance

The first component is 0.8 of neuron 141 and 0.6 of neuron 190, which is the orange arrow, and it carries 90 percent of the variance. That is what “reduce two dimensions to one” means: keep the position along the orange arrow, throw away the position along the teal one, and you have thrown away a tenth of the information. The for loop at the end of the function handles a small annoyance: an axis has no preferred direction, so an eigenvector and its negative are equally valid, and different computers can return either. Flipping each one so that its largest entry is positive makes the plots come out the same way every time.

Step 3: all 66 neurons, and a trap

The function does not care how many columns the table has. Run it on the full table.

components, variance, mean = pca(X)
fraction = variance / variance.sum()
print("variance carried by the first five components:", fraction[:5].round(3))
print("neurons with the biggest weight in component 1:", visp[np.argsort(-np.abs(components[:, 0]))[:4]])
print("standard deviation of those neurons' counts:", X.std(axis=0)[np.argsort(-np.abs(components[:, 0]))[:4]].round(1))
print("median standard deviation across all 66:", np.median(X.std(axis=0)).round(1))
# variance carried by the first five components: [0.724 0.07  0.031 0.026 0.022]
# neurons with the biggest weight in component 1: [141 190 184 120]
# standard deviation of those neurons' counts: [6.  4.8 3.5 3.6]
# median standard deviation across all 66: 0.4

72 percent of the variance in one component, out of 66. That sounds like a spectacular result, and it is the trap. Look at which neurons carry that component: 141, 190, 184 and 120, and look at their standard deviations, 6, 4.8, 3.5 and 3.6 spikes, against a median across the population of 0.4. PCA finds the directions of largest variance, and variance is measured in spikes squared, so a neuron that fires fifteen times per window has more of it than fifty quiet neurons put together. The first component is not “what the population does”. It is “what the four loudest neurons do”, and the other 62 barely got a vote.

Whether that is a problem depends on the question. If you believe loud neurons matter more, the raw analysis is right. Almost nobody believes that, and the field’s standard fix is to put every neuron on the same scale before asking about directions.

Step 4: standardise, then ask how many dimensions

Subtract each neuron’s mean and divide by its standard deviation, so that every column has mean 0 and spread 1. This is the same z-scoring you would apply before most machine learning, and it has one prerequisite: a neuron that never fired in the window has a standard deviation of zero and no direction at all, so drop it first. Three of the 66 go.

active = X.std(axis=0) > 0                              # drop neurons that never fired in the window
X, visp = X[:, active], visp[active]
Z = (X - X.mean(axis=0)) / X.std(axis=0)                 # every neuron: mean 0, spread 1
print(Z.shape[1], "neurons kept")

components, variance, mean = pca(Z)
fraction = variance / variance.sum()
cumulative = np.cumsum(fraction)
print("variance carried by the first five components:", fraction[:5].round(3))
print("components needed for half the variance:", np.searchsorted(cumulative, 0.5) + 1)
print("components needed for 80% of the variance:", np.searchsorted(cumulative, 0.8) + 1)
# 63 neurons kept
# variance carried by the first five components: [0.19  0.053 0.042 0.035 0.029]
# components needed for half the variance: 12
# components needed for 80% of the variance: 30

fig, ax = plt.subplots(figsize=(10, 3.8))
ax.bar(np.arange(1, 21), fraction[:20] * 100)
ax.plot(np.arange(1, 21), cumulative[:20] * 100, marker="o", color="orange")
ax.set(xlabel="component", ylabel="% of variance", xticks=np.arange(1, 21))
plt.show()

Now the first component carries 19 percent, and it takes twelve components to account for half the variance and thirty to account for 80 percent of it. This is the honest number, and it is worth sitting with. Primary visual cortex, in a 200 ms window, in a task with one stimulus that varies along one axis, is not a one-dimensional place. The population is doing many things at once, most of which have nothing to do with the contrast of the stimulus. There is a live argument in the field about what those dimensions are, and Steinmetz and colleagues’ answer, from this very dataset, is that a great many of them are the animal’s own movements.

Step 5: every trial as a point

Even so, the first two components are the best two-dimensional picture of the population that exists, so draw it. Projecting a trial onto a component is a dot product: multiply each neuron’s standardised count by the component’s weight for that neuron and add up. @ does it for all trials and both components at once.

scores = (Z - mean) @ components[:, :2]                 # every trial's position along components 1 and 2

fig, ax = plt.subplots(figsize=(7, 6))
for lvl in [0.0, 0.25, 0.5, 1.0]:
    pick = contrast == lvl
    ax.scatter(scores[pick, 0], scores[pick, 1], s=14, alpha=0.7, label=f"contrast {lvl}")
ax.set(xlabel="component 1", ylabel="component 2")
ax.legend()
plt.show()

for lvl in [0.0, 0.25, 0.5, 1.0]:
    print(f"contrast {lvl}: mean position along component 1 = {scores[contrast == lvl, 0].mean():+.1f}")
print(f"guess 'stimulus' whenever component 1 is positive: {np.mean((scores[:, 0] > 0) == y):.1%} correct, with no labels used")
# contrast 0.0: mean position along component 1 = -2.5
# contrast 0.25: mean position along component 1 = -0.5
# contrast 0.5: mean position along component 1 = +2.1
# contrast 1.0: mean position along component 1 = +4.6
# guess 'stimulus' whenever component 1 is positive: 84.7% correct, with no labels used

PCA was never told which trials had a stimulus. It was told nothing at all; it looked only at how the neurons co-vary. And yet the first thing it found is the stimulus: the four contrast levels line up along component 1 in order, and drawing a line at zero classifies 85 percent of trials correctly. Compare that with lesson 3, where the decoder was given the labels and reached 99.7 percent. The gap between 85 and 99.7 is the difference between the direction of most variance and the direction of most information about the stimulus. They are related here because the stimulus is the biggest thing happening to visual cortex. In other areas, or for subtler variables, they can be entirely different directions, which is why unsupervised methods like PCA and supervised ones like the decoder answer different questions.

Step 6: every moment as a point

So far each trial has been a single point, a 200 ms snapshot. The last step changes what a point is. Average the population’s firing rate over trials as in lesson 2, but for every neuron at once, so that each 10 ms bin of the trial is now a row of the table, a snapshot of what all 63 neurons were doing at that moment. Do it separately for stimulus and blank trials, stack the two, standardise, and run PCA. Now a component is a pattern of activity across the population, and the trial itself becomes a path through the space.

def population_rates(trials, smooth=7):
    """Trial-averaged, lightly smoothed firing rate of every active VISp neuron: neurons x time bins, in spikes/s."""
    rate = spks[visp][:, trials].mean(axis=1) / dt
    kernel = np.ones(smooth) / smooth
    return np.array([np.convolve(r, kernel, mode="same") for r in rate])

R_stim, R_blank = population_rates(y == 1), population_rates(y == 0)
both = np.concatenate([R_stim, R_blank], axis=1).T      # 500 time points x 63 neurons
scale = both.std(axis=0)                                # again, put every neuron on the same footing
components, variance, mean = pca(both / scale)
print("variance carried by the first three components:", (variance[:3] / variance.sum()).round(3))
# variance carried by the first three components: [0.525 0.075 0.046]

path_stim = (R_stim.T / scale - mean) @ components[:, :2]     # 250 time points x 2
path_blank = (R_blank.T / scale - mean) @ components[:, :2]

show = slice(30, 151)                                   # -0.2 to +1.0 s
marks = [50, 55, 60, 70, 80, 100, 150]                   # bins: 0, 50, 100, 200, 300, 500, 1000 ms after onset

fig, axes = plt.subplots(1, 2, figsize=(11, 4.6))
axes[0].plot(path_blank[show, 0], path_blank[show, 1], color="gray", label="blank screen")
axes[0].plot(path_stim[show, 0], path_stim[show, 1], label="stimulus")
axes[0].scatter(path_stim[marks, 0], path_stim[marks, 1], s=30, zorder=3)
for m in marks:
    axes[0].annotate(f"{round(t[m] * 1000)} ms", (path_stim[m, 0], path_stim[m, 1]), fontsize=8, xytext=(4, 4), textcoords="offset points")
axes[0].set(xlabel="component 1", ylabel="component 2")
axes[0].legend()

for lvl in [0.0, 0.25, 0.5, 1.0]:
    path = (population_rates(contrast == lvl).T / scale - mean) @ components[:, 0]
    axes[1].plot(t, path, label=f"contrast {lvl}")
axes[1].axvline(0, color="orange")
axes[1].set(xlabel="time from stimulus onset (s)", ylabel="component 1", xlim=(-0.3, 1.0))
axes[1].legend()
plt.show()

Read the left panel like a map. Before the stimulus both paths sit in the same small patch. Then the stimulus path leaves: by 50 ms it is on its way, at 100 ms it is as far from home as it will get, and by 200 ms it is swinging back, before a second, smaller excursion around 300 ms and a slow return that is still not complete at one second. That whole loop is the transient response and the second hump from lesson 2’s PSTH, but now for the population rather than one neuron, and now as a shape rather than a curve. The blank path goes nowhere along component 1 and drifts up component 2, which is the population doing something slow and unrelated to vision while the mouse waits and moves. This picture, a trial as a trajectory through a low-dimensional space, is how a great deal of modern systems neuroscience thinks about population activity.

The right panel asks how the path depends on the stimulus, and the answer is the cleanest result in the lesson: it is the same path at four different scales. Contrast does not send the population in a different direction. It sends it further along the same one, and the first component, with 52 percent of the time-course variance, is essentially a contrast axis.

What you just did

• Saw what a principal component is on two neurons: the long axis of the cloud of trials.
• Wrote PCA from scratch with a covariance matrix and np.linalg.eigh, and understood every line.
• Ran it on 66 neurons, fell into the loud-neuron trap, and climbed out by standardising.
• Measured the dimensionality of visual cortex activity: 12 components for half the variance, not one.
• Found the stimulus without labels, and understood why unsupervised and supervised methods find different directions.
• Turned a trial into a trajectory through population space and read the response off it as a shape.

Exercises

1. Single trials. Step 6 averaged over trials before running PCA. Project a handful of individual stimulus trials onto the same components (smooth them more heavily, say 15 bins) and plot them over the average path. How much do single trials wobble around it?
2. Repeat step 4 for VISam, MD and CA1. Which area needs the most components to reach half its variance, and how does that fit with what lessons 2 and 3 found there?
3. Compare component 1 from step 4 with the decoder weights from lesson 3 (train the decoder on all trials of the standardised table). Compute the correlation between the two 63-element vectors. Are they the same direction?
4. Harder: the blank-screen path in step 6 drifts along component 2. Sort the blank trials by dat["response"], average each group separately, and project them. Is component 2 about what the mouse does?

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 code in this lesson runs top to bottom as a single script in a few seconds.

In lesson 2 you asked whether a neuron responds to the stimulus. This lesson turns the question around: given only the neurons, can you tell what the stimulus was? You will turn a recording into a table, guess the stimulus from one neuron, discover why that guess cannot be trusted, then train a classifier on 66 neurons at once and find that it gets the answer right on 339 trials out of 340. Along the way you will build the three habits that separate machine learning from wishful thinking: hold out test data, cross-validate, and check what luck alone would score. At the end you will point the decoder at eight brain areas and at every moment of the trial, and watch the information appear.

This is the first step toward a brain-computer interface. A BCI is a loop (here is the map), and the stage in the middle of it, where neural activity is turned into a guess about what the user wants, is exactly what you are about to build.

flowchart LR

A[Spike counts: neurons x trials x bins] --> B[Table: one row per trial, one column per neuron]

B --> C[One neuron + a threshold]

C --> D[Train / test split]

D --> E[Population decoder: 66 neurons]

E --> F[Cross-validate]

F --> G[Shuffle the labels: what does luck score?]

G --> H[Every area, every moment]

The data

Same file as lesson 2: steinmetz_session11.zip, one session from Steinmetz and colleagues (2019), with 698 neurons, 340 trials and spike counts in 10-millisecond bins. On each trial a striped pattern appeared on the right screen at one of four contrasts, or the screen stayed blank. The question for the decoder is the simplest one possible: was there anything on the right screen, or not?

Step 1: turn the recording into a table

Every classifier ever built wants the same thing: a table with one row per example and one column per measurement, plus a list of the right answers. In machine learning the table is called X and the answers are called y. Here, an example is a trial, a measurement is one neuron’s spike count in the 50 to 250 ms window from lesson 2, and the answer is 1 if the right screen showed a stimulus and 0 if it was blank. Two lines of NumPy build the whole table.

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

post = (t >= 0.05) & (t < 0.25)           # the response window from lesson 2
visp = np.where(area == "VISp")[0]        # primary visual cortex neurons

X = spks[visp][:, :, post].sum(axis=2).T  # trials x neurons: each neuron's spike count in the window
y = (contrast > 0).astype(int)            # 1 if there was a stimulus on the right, 0 if the screen was blank

print(X.shape, "trials x neurons")
print("stimulus on", y.sum(), "of", len(y), "trials")
print("trial 0:", X[0, :12], "... label", y[0])
print("trial 3:", X[3, :12], "... label", y[3])
# (340, 66) trials x neurons
# stimulus on 173 of 340 trials
# trial 0: [0 0 0 0 0 0 0 0 0 0 0 0] ... label 0
# trial 3: [2 2 4 2 0 0 2 0 0 0 2 0] ... label 1

340 rows, 66 columns. Trial 0 was a blank screen and the first dozen visual cortex neurons were silent. Trial 3 had a stimulus and they were not. The .T at the end of the X line transposes the array so that trials are rows, which is the convention every machine learning tool expects. Everything from here on works on X and y; the raw recording is not needed again until step 7.

Step 2: one neuron, one threshold

Start with the neuron we know best. Neuron 141 fired 66 spikes per second to a full-contrast stimulus and 8 to a blank screen. So here is a decoder: count its spikes in the window, and if there are more than some number, say the stimulus was there. Before choosing the number, look at the two piles of trials.

one = spks[141][:, post].sum(axis=1)      # neuron 141's spike count in the window, on every trial

fig, ax = plt.subplots(figsize=(10, 3.8))
bins = np.arange(one.max() + 2) - 0.5     # one bar per whole number of spikes
ax.hist(one[y == 0], bins=bins, alpha=0.6, label="blank screen")
ax.hist(one[y == 1], bins=bins, alpha=0.6, label="stimulus")
ax.set(xlabel="spikes from neuron 141, 50 to 250 ms after onset", ylabel="number of trials")
ax.legend()
plt.show()

A decoder is a rule that turns a row of the table into a guess, and its accuracy is the fraction of trials it gets right. Try every threshold from 1 to 9.

def accuracy(guess, truth):
    """Fraction of trials where the guess matches the truth."""
    return np.mean(guess == truth)

for k in range(1, 10):
    guess = one > k                       # True where the neuron fired more than k spikes
    print(f"more than {k} spikes means stimulus: {accuracy(guess, y):.1%} correct")
# more than 1 spikes means stimulus: 81.8% correct
# more than 2 spikes means stimulus: 85.0% correct
# more than 3 spikes means stimulus: 87.9% correct
# more than 4 spikes means stimulus: 87.9% correct
# more than 5 spikes means stimulus: 88.2% correct
# more than 6 spikes means stimulus: 85.3% correct
# more than 7 spikes means stimulus: 82.6% correct
# more than 8 spikes means stimulus: 79.4% correct
# more than 9 spikes means stimulus: 76.5% correct

88 percent from a single neuron and a single number. That is genuinely good. It is also, as it stands, slightly dishonest, and the next step is about why.

Step 3: never grade yourself on the questions you studied

We picked the threshold by trying all of them on the 340 trials and keeping the best. Then we reported the accuracy on the same 340 trials. That is grading yourself on the exam you used to revise. With one number to choose it barely matters, but with 66 weights to choose, as in the next step, a decoder can memorise the quirks of the trials it was trained on and score brilliantly on them while knowing nothing that transfers to a new trial. The fix is a rule so important that it is the one thing to take away from this lesson if you take away nothing else: choose the rule on some trials, and measure it on different ones.

rng.permutation shuffles the trial numbers, and slicing gives us 240 training trials and 100 test trials that the decoder never sees until it is judged.

rng = np.random.default_rng(0)
shuffled = rng.permutation(len(y))                  # the trial numbers 0..339 in random order
train, test = shuffled[:240], shuffled[240:]        # 240 trials to learn from, 100 to be examined on

scores = [accuracy(one[train] > k, y[train]) for k in range(15)]
best_k = int(np.argmax(scores))                     # the threshold that did best on the training trials

print(f"rule learned from training trials: more than {best_k} spikes means stimulus")
print(f"training trials: {scores[best_k]:.1%} correct")
print(f"test trials:     {accuracy(one[test] > best_k, y[test]):.1%} correct")
# rule learned from training trials: more than 3 spikes means stimulus
# training trials: 88.8% correct
# test trials:     86.0% correct

The rule chosen on the training trials scores 88.8 percent on them and 86 percent on the held-out test trials. The drop is small here because the rule is simple. Watch for it in everything you do from now on: the gap between training and test accuracy is the size of the lie you would have told yourself.

Step 4: 66 neurons, and a model neuron to read them

Neuron 141 is one of 66 in primary visual cortex, and every one of them saw the stimulus. Combining them needs a rule with more than one number: give each neuron a weight, multiply its spike count by that weight, add everything up, and say “stimulus” if the sum is above zero. Neurons that fire more for the stimulus should get positive weights, neurons that fire less should get negative ones, and neurons that do not care should get weights near zero. The only question is how to find the 66 weights, and the answer is to learn them from the training trials, one small correction at a time.

Look at the shape of that rule before reading the code. Inputs arrive, each is multiplied by a weight, they are summed, and the sum is pushed through a threshold. That is a neuron. Not the leaky integrate-and-fire neuron of lesson 1 but the other kind, the one Frank Rosenblatt built out of motors and potentiometers in 1958 and called a perceptron, and which is still the basic unit of every deep network. We are going to decode 66 real neurons with one artificial one.

def train_decoder(X, y, lr=0.01, steps=1000):
    """Logistic regression, trained by gradient descent.
    X: trials x neurons spike counts. y: 0 or 1 per trial.
    Returns one weight per neuron and a bias."""
    w = np.zeros(X.shape[1])                        # start with every weight at zero
    b = 0.0
    for _ in range(steps):
        p = 1 / (1 + np.exp(-(X @ w + b)))          # weighted sum of the counts, squashed to a probability
        error = p - y                               # how wrong it was on each trial, and in which direction
        w -= lr * (X.T @ error) / len(y)            # nudge each weight against its share of the error
        b -= lr * error.mean()
    return w, b

def predict(X, w, b):
    """1 where the weighted sum says stimulus, 0 where it says blank."""
    return (X @ w + b > 0).astype(int)

w, b = train_decoder(X[train], y[train])
print(f"training trials: {accuracy(predict(X[train], w, b), y[train]):.1%} correct")
print(f"test trials:     {accuracy(predict(X[test], w, b), y[test]):.1%} correct")
# training trials: 99.6% correct
# test trials:     100.0% correct

Read train_decoder as a loop of three moves, repeated a thousand times. X @ w + b computes the weighted sum for every trial at once; @ is matrix multiplication, and it does in one symbol what would otherwise be a loop over trials inside a loop over neurons. The 1 / (1 + np.exp(-...)) wrapper squashes each sum to a number between 0 and 1, the decoder’s probability that the stimulus was there. error is the gap between that probability and the truth. And the update line moves each neuron’s weight in the direction that would have shrunk the error, by an amount proportional to how much that neuron fired. The learning rate lr keeps the steps small.

That update rule is worth a second look, because it is local: the change to a neuron’s weight depends only on that neuron’s own activity and the error. A synapse could implement it. This method has a name, logistic regression, and a longer history than the name suggests. In a library like scikit-learn it is one line. We wrote it out so that there is nothing hidden.

99.6 percent on the training trials and 100 percent on the 100 test trials. The population knows something no single neuron knows. To see how confident it is, look at the weighted sum itself, on the test trials only.

evidence = X[test] @ w + b                          # the decoder's weighted sum on each test trial

fig, ax = plt.subplots(figsize=(10, 3.8))
bins = np.linspace(evidence.min(), evidence.max(), 40)
ax.hist(evidence[y[test] == 0], bins=bins, alpha=0.6, label="blank screen")
ax.hist(evidence[y[test] == 1], bins=bins, alpha=0.6, label="stimulus")
ax.axvline(0, color="orange")                       # the decision boundary
ax.set(xlabel="decoder's weighted sum", ylabel="number of test trials")
ax.legend()
plt.show()

Step 5: cross-validation

One split of 240 and 100 gives one number, and 100 test trials is a small exam: one lucky trial moves the score by a full percentage point. Cross-validation fixes this by dividing the trials into five groups and letting each group be the test set in turn, training on the other four. Every trial gets tested exactly once, by a decoder that never saw it, and the five scores are averaged. This is the standard way to report a decoder’s accuracy, and the function below is the one we will use for the rest of the lesson.

def cross_validate(X, y, folds=5, seed=0):
    """Average test accuracy when every trial takes one turn in the test set."""
    rng = np.random.default_rng(seed)
    parts = np.array_split(rng.permutation(len(y)), folds)   # five random, equal groups of trials
    scores = []
    for part in parts:
        is_test = np.zeros(len(y), dtype=bool)
        is_test[part] = True                                 # this group is the test set this time
        w, b = train_decoder(X[~is_test], y[~is_test])       # train on everything else
        scores.append(accuracy(predict(X[is_test], w, b), y[is_test]))
    return np.mean(scores)

print(f"VISp population, cross-validated: {cross_validate(X, y):.1%} correct")
print(f"neuron 141 alone, cross-validated: {cross_validate(X[:, visp == 141], y):.1%} correct")
# VISp population, cross-validated: 99.7% correct
# neuron 141 alone, cross-validated: 87.9% correct

99.7 percent: 339 of 340 trials. The single-neuron rule from step 2, run through the same machinery, gets 87.9. The next two steps make sure that 99.7 means what we think it means.

Step 6: what would luck score?

Fifty percent is not the only kind of chance. A decoder can find structure in things that have nothing to do with the stimulus: a slow drift in firing rates over the session, say, or a neuron that fires more in the second half of the experiment when more of the stimulus trials happened to be. The permutation test from lesson 2 handles this here too. Shuffle the labels so that each trial keeps its spike counts but is assigned a random answer, run the whole cross-validation, and see what accuracy comes out. Do it 200 times. This is the slow part of the lesson; it takes a minute or two.

rng = np.random.default_rng(1)
chance = np.array([cross_validate(X, rng.permutation(y)) for _ in range(200)])

print(f"shuffled labels: mean {chance.mean():.1%}, best of 200 shuffles {chance.max():.1%}")
print(f"95% of shuffles score below {np.percentile(chance, 95):.1%}")
# shuffled labels: mean 49.6%, best of 200 shuffles 58.8%
# 95% of shuffles score below 54.7%

With the labels scrambled the decoder averages 49.6 percent, and its best result in 200 attempts is 58.8. Our real score is 99.7. Anything under about 55 percent could be luck; that line is going to matter in the next step, where the numbers are not so clear-cut.

Step 7: which areas know?

Wrap the whole pipeline in a function that takes a list of neurons, and run it on each brain area. The seven areas from lesson 2 are here, plus MD, the mediodorsal thalamus, which has the most neurons of any area in this recording.

def decode(neurons, window=post):
    """Cross-validated accuracy of a decoder that reads these neurons in this time window."""
    counts = spks[neurons][:, :, window].sum(axis=2).T
    return cross_validate(counts, y)

names = ["VISp", "VISam", "LGd", "MD", "CA1", "DG", "MOs", "ACA"]
by_area = []
for a in names:
    idx = np.where(area == a)[0]
    by_area.append(decode(idx))
    print(f"{a:5s} {len(idx):3d} neurons, {by_area[-1]:.0%} correct")

fig, ax = plt.subplots(figsize=(10, 3.8))
ax.bar(names, by_area)
ax.axhline(np.percentile(chance, 95), color="orange")       # anything below this line could be luck
ax.set(ylabel="decoding accuracy", ylim=(0, 1))
plt.show()
# VISp   66 neurons, 100% correct
# VISam  79 neurons, 71% correct
# LGd    11 neurons, 57% correct
# MD    126 neurons, 77% correct
# CA1    50 neurons, 52% correct
# DG     65 neurons, 63% correct
# MOs     6 neurons, 57% correct
# ACA    16 neurons, 61% correct

Primary visual cortex is essentially perfect. VISam, a higher visual area, gets 71 percent, and CA1, which in lesson 2 had not one neuron responding to the stimulus, is at chance. So far the map from lesson 2 is confirmed. Then there is MD at 77 percent, the second best area in the recording, and DG at 63, even though lesson 2 found that only two percent of its neurons respond to the stimulus. A decoder pools weak signals that a single-neuron test misses, and part of the explanation is that. But look at what else these trials differ in. On stimulus trials the mouse almost always turns the wheel; on blank trials it often does nothing. An area that carries no visual information at all, but knows what the mouse is about to do, will decode the stimulus above chance, because in this task the two are correlated. Restrict the analysis to trials where the mouse made the same response, and MD’s advantage over simply guessing the commoner label shrinks to a few points.

This is the most important caveat in decoding, and it applies to every result in the field, including the ones in brain-computer interfaces. A decoder tells you that the information is present in an area. It does not tell you why it is there, or what the area is doing with it, or whether it would still be there if the animal’s behaviour were different. Steinmetz and colleagues built their whole paper around this problem, and found that signals related to movement are present nearly everywhere in the mouse brain.

Step 8: when does the brain know?

So far the decoder has read one fixed window. Slide a 100 ms window across the trial instead, training a fresh decoder at each position, and you get the accuracy as a function of time. This is the neural equivalent of watching the information arrive.

starts = np.arange(10, 140, 5)                      # window start, in bins: every 50 ms
acc_time = []
for s in starts:
    window = np.zeros(len(t), dtype=bool)
    window[s:s + 10] = True                         # a 100 ms window starting at bin s
    acc_time.append(decode(visp, window))

fig, ax = plt.subplots(figsize=(10, 3.8))
ax.plot(t[starts] + 0.05, acc_time, marker="o")     # plot each window at its centre
ax.axvline(0, color="orange")                       # stimulus onset
ax.axhline(0.5, color="gray", ls=":")               # chance
ax.set(xlabel="centre of 100 ms window, time from stimulus onset (s)", ylabel="decoding accuracy", ylim=(0.4, 1))
plt.show()

Before the stimulus, chance. That is a sanity check as much as a result: a decoder that scores well before the stimulus appears has found a leak, and you should go looking for it. In the first 100 ms after onset the accuracy jumps to 97 percent, which is the transient response from lesson 2 doing its work. It peaks at 99 in the 50 to 150 ms window, dips around 200 ms, and then stays above 90 percent for most of the next half second, partly because the stimulus is still on the screen and partly, as step 7 warned, because the mouse is now moving.

What you just did

• Turned a neural recording into the trials-by-features table that every machine learning method starts from.
• Built a one-neuron decoder and learned why accuracy must be measured on held-out trials.
• Wrote logistic regression from scratch, trained it by gradient descent, and understood it as a model neuron with 66 learnable synapses.
• Cross-validated it, and established what luck alone would score by shuffling the labels.
• Decoded the stimulus from eight brain areas and from every moment of the trial, and met the central caveat of the field: information present is not the same as information used.

Every brain-computer interface, from the cursor-control systems in clinical trials to the speech decoders in the news, is this pipeline with more neurons, a fancier classifier, and a loop that feeds the guess back to the user. You have now built the part in the middle.

Exercises

1. Decode the left stimulus instead: y = (dat["contrast_left"] > 0).astype(int). The probes are in the left hemisphere. What does VISp score now, and what does that tell you?
2. Train the decoder on all 340 VISp trials and sort the neurons by the size of their weight. Is neuron 141 at the top? The neuron with the largest weight has a negative sign. Plot its PSTH for blank and stimulus trials, as in lesson 2, and work out what the decoder is using it for.
3. How many neurons do you need? Pick random subsets of 1, 2, 5, 10, 20 and 40 VISp neurons, cross-validate each, and plot accuracy against the number of neurons. Where does it saturate?
4. Harder, and the real first step of a BCI: decode what the mouse is about to do. Keep only the trials where it turned the wheel (dat["response"] != 0), label them by direction, and use a window from 250 to 750 ms. Try VISp, VISam and MD. Which area knows which way the mouse will turn, and how early can you decode it?

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 code in this lesson runs top to bottom as a single script; the shuffle test in step 6 is the only slow part.