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 neuronsnp.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()
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()
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.60A 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()
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 --> BThat 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.0Seventy-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()
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()
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
- Change
bin_sizein step 4 to 0.1 and to 10. What changes, and which is more useful? - Find the neuron with the highest CV in the whole population. Plot its raster. What does its firing look like?
- 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?
- 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.
