Installing requirements

The cell below installs every Python package needed to run this notebook, at fully pinned versions, using uv for fast resolution. In Colab the cell is collapsed by default — click the ▶ button to run it.

# install cell skipped during CI (deps preinstalled into system Python)

⚠️ Restart runtime after install

The install may upgrade packages already loaded in the kernel. Go to Runtime → Restart session, then Run all cells below (skip this install cell on re-run).

Figure 2: peri-head distance coding in the mouse brainstem

This notebook reproduces the Figure 2 panels of Xiao, Severson, et al. (2026), Peri-Head Distance Coding in the Mouse Brainstem, directly from Dandiset 001687, streaming the data rather than downloading it.

For each of six example PrV units it draws three panels:

Panel	Content
A. Raster	Spike times aligned to wall-pass onset, one row per trial, grouped and coloured by wall distance
B. PETH	Peri-event time histogram, one trace per wall distance, mean +/- SEM across trials
C. Tuning curve	The stored z-scored firing rate versus wall distance from processing/wall_tuning

The experiment. A head-fixed mouse sits while a motorised wall approaches and passes its head at a controlled distance. Second-order trigeminal (PrV) neurons respond to whiskers touching the wall, and report how far away that wall is. Two coding schemes appear in the data: proximity units that increase firing monotonically as the wall nears, and map units tuned to a preferred distance.

Citation.
Xiao, Severson, et al. (2026). Peri-head distance coding in the mouse brainstem, https://doi.org/10.1016/j.neuron.2026.05.027.
Dandiset 001687. https://doi.org/10.48324/dandi.001687/0.260805.1529.

1. What we stream, and why

The Dandiset holds two files per recording session:

• a processed file (no modality suffix) with sorted spike times, trials, tuning curves and PSTHs. These are small, a few MB each.
• a raw file (_ecephys+image) with the 30 kHz broadband, LFP, TTLs and a reference to the behaviour video. These are multi-GB.

This notebook only needs the processed files. It reads six of them, one per example unit, and never downloads a whole file: remfile fetches only the byte ranges that h5py actually asks for.

We pin the published version 0.260805.1529 rather than draft, so the figure is reproducible even if the Dandiset is later revised.

import warnings

import h5py
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import remfile
from dandi.dandiapi import DandiAPIClient
from matplotlib.colors import Normalize
from pynwb import NWBHDF5IO

warnings.filterwarnings("ignore", message=".*namespace.*already exists.*")

DANDISET_ID = "001687"
DANDISET_VERSION = "0.260805.1529"  # published snapshot; use "draft" for the latest

2. The six example units

Each row below identifies one Figure 2 unit by its DANDI asset path and its 0-based row in that file's units table. This mapping is also distributed as figure_unit_manifest.csv in the manuscript repository, reproduced inline here so the notebook stands alone.

tuning_type is the classification carried in the NWB units table itself.

FIGURE2_UNITS = pd.DataFrame(
    [
        # display, asset path,                                                      unit row, label,                          tuning
        ("Unit 1", "sub-Ephys2/sub-Ephys2_ses-20230208-sp-session-1.nwb",     4,  "20230208_sp_session_1_unit_5",  "proximity"),
        ("Unit 2", "sub-Ephys40/sub-Ephys40_ses-20230626-sp-session-1.nwb",   7,  "20230626_sp_session_1_unit_8",  "map"),
        ("Unit 3", "sub-Ephys88/sub-Ephys88_ses-20240905-sp-session-1.nwb",   3,  "20240905_sp_session_1_unit_4",  "map"),
        ("Unit 4", "sub-Ephys81/sub-Ephys81_ses-20240703-sp-session-1.nwb",   9,  "20240703_sp_session_1_unit_10", "map"),
        ("Unit 5", "sub-Ephys22/sub-Ephys22_ses-20221212-sp-session-2.nwb",  17,  "20221212_sp_session_2_unit_18", "map"),
        ("Unit 6", "sub-Ephys106/sub-Ephys106_ses-20240924-sp-session-1.nwb", 6,  "20240924_sp_session_1_unit_7",  "suppressed"),
    ],
    columns=["display_unit", "asset_path", "unit_index", "unit_label", "tuning_type"],
)
FIGURE2_UNITS

	display_unit	asset_path	unit_index	unit_label	tuning_type
0	Unit 1	sub-Ephys2/sub-Ephys2_ses-20230208-sp-session-...	4	20230208_sp_session_1_unit_5	proximity
1	Unit 2	sub-Ephys40/sub-Ephys40_ses-20230626-sp-sessio...	7	20230626_sp_session_1_unit_8	map
2	Unit 3	sub-Ephys88/sub-Ephys88_ses-20240905-sp-sessio...	3	20240905_sp_session_1_unit_4	map
3	Unit 4	sub-Ephys81/sub-Ephys81_ses-20240703-sp-sessio...	9	20240703_sp_session_1_unit_10	map
4	Unit 5	sub-Ephys22/sub-Ephys22_ses-20221212-sp-sessio...	17	20221212_sp_session_2_unit_18	map
5	Unit 6	sub-Ephys106/sub-Ephys106_ses-20240924-sp-sess...	6	20240924_sp_session_1_unit_7	suppressed

3. Stream each session

load_session opens one asset over HTTP and pulls out the three things the figure needs: the trials table, the units table, and the stored tuning curves.

def load_session(dandiset, asset_path):
    """Stream one NWB file and return (trials, units, tuning, distance_axis)."""
    asset = dandiset.get_asset_by_path(asset_path)
    url = asset.get_content_url(follow_redirects=1, strip_query=True)

    store = remfile.File(url)
    with h5py.File(store, "r") as h5, NWBHDF5IO(file=h5, load_namespaces=True) as io:
        nwb = io.read()
        trials = nwb.trials.to_dataframe()
        units = nwb.units.to_dataframe()

        tuning = None
        distances = None
        if "wall_tuning" in nwb.processing:
            module = nwb.processing["wall_tuning"]
            if "tuning_curves" in module.data_interfaces:
                tuning = module["tuning_curves"].to_dataframe()
            if "distance_axis" in module.data_interfaces:
                distances = module["distance_axis"].to_dataframe()["distance_mm"].to_numpy(float)

        return {
            "session_id": nwb.session_id or nwb.identifier,
            "trials": trials,
            "units": units,
            "tuning": tuning,
            "distances": distances,
        }


sessions = {}
with DandiAPIClient() as client:
    dandiset = client.get_dandiset(DANDISET_ID, DANDISET_VERSION)
    for asset_path in FIGURE2_UNITS["asset_path"].unique():
        sessions[asset_path] = load_session(dandiset, asset_path)
        info = sessions[asset_path]
        print(f"{info['session_id']:26s} {len(info['trials']):3d} trials  {len(info['units']):3d} units")

20230208_sp_session_1      200 trials   30 units

20230626_sp_session_1      200 trials   27 units

20240905_sp_session_1      200 trials   14 units

20240703_sp_session_1      200 trials   21 units

20221212_sp_session_2      255 trials   27 units

20240924_sp_session_1      200 trials   14 units

Each session has 200 wall passes. The wall distance for every pass is stored per trial in trials.wall_distance_mm, and the tuning analysis groups those into 2 mm bins spanning roughly 7 to 23 mm.

example = sessions[FIGURE2_UNITS.loc[0, "asset_path"]]
print("trials columns :", list(example["trials"].columns))
print("distance axis  :", np.round(example["distances"], 2))
example["trials"].head()

trials columns : ['start_time', 'stop_time', 'wall_distance_mm', 'trial_condition']
distance axis  : [ 6.74  8.74 10.74 12.74 14.74 16.74 18.74 20.74 22.74]

	start_time	stop_time	wall_distance_mm	trial_condition
id
0	40.755267	51.755267	10.742
1	55.474000	66.474000	22.742
2	69.976367	80.976367	24.742
3	84.794467	95.794467	14.742
4	99.580367	110.580367	18.742

4. Plotting helpers

Rasters and PETHs are computed here from the raw spike times rather than read from the stored PSTHs, so you can see exactly how they are built and change the window or binning. The tuning curves in panel C are the stored values.

Colour encodes wall distance throughout, near (red) to far (blue).

WINDOW = (0.0, 12.0)      # seconds relative to wall-pass onset
BIN_WIDTH = 0.2           # PETH bin, seconds
SMOOTH_BINS = 3           # boxcar smoothing width, in bins
TUNING_SCALE = 2.0        # matches the manuscript figure
PETH_TICK_STEPS = [10, 10, 20, 10, 2, 10]   # y-axis step per row


def distance_colors(distances):
    cmap = plt.get_cmap("turbo_r")
    norm = Normalize(vmin=float(np.nanmin(distances)), vmax=float(np.nanmax(distances)))
    return {float(d): cmap(norm(float(d))) for d in distances}


def smooth(values, width):
    if width <= 1:
        return values
    return np.convolve(values, np.ones(width) / width, mode="same")


def spikes_in_trial(spike_times, start_time, window):
    """Spike times relative to trial onset, using binary search on sorted times."""
    lo = np.searchsorted(spike_times, start_time + window[0], side="left")
    hi = np.searchsorted(spike_times, start_time + window[1], side="right")
    return spike_times[lo:hi] - start_time


def trials_at_distance(trials, distance):
    return trials[np.isclose(trials["wall_distance_mm"].to_numpy(float), distance)]

def plot_raster(ax, spike_times, trials, distances, colors, window):
    y = 0
    for distance in distances:
        for _, trial in trials_at_distance(trials, distance).iterrows():
            rel = spikes_in_trial(spike_times, float(trial["start_time"]), window)
            if rel.size:
                ax.vlines(rel, y + 0.08, y + 0.92,
                          color=colors[float(distance)], linewidth=0.25, alpha=0.85)
            y += 1
    ax.axvspan(3.0, 11.0, color="0.88", zorder=-10)   # wall in motion
    ax.set_xlim(window)
    ax.set_ylim(-1, max(y, 1))
    ax.set_yticks([])
    ax.tick_params(labelsize=8, length=2)
    for side in ("top", "right"):
        ax.spines[side].set_visible(False)


def plot_peth(ax, spike_times, trials, distances, colors, window, bin_width, smooth_bins):
    edges = np.arange(window[0], window[1] + bin_width, bin_width)
    centers = edges[:-1] + bin_width / 2
    for distance in distances:
        rows = trials_at_distance(trials, distance)
        if rows.empty:
            continue
        rates = np.asarray(
            [np.histogram(spikes_in_trial(spike_times, float(t["start_time"]), window),
                          bins=edges)[0] / bin_width for _, t in rows.iterrows()],
            dtype=float,
        )
        mean = smooth(np.nanmean(rates, axis=0), smooth_bins)
        sem = smooth(np.nanstd(rates, axis=0, ddof=1) / np.sqrt(max(len(rates), 1)), smooth_bins)
        color = colors[float(distance)]
        ax.plot(centers, mean, color=color, linewidth=1.1)
        ax.fill_between(centers, mean - sem, mean + sem, color=color, alpha=0.15, linewidth=0)
    ax.axvspan(3.0, 11.0, color="0.88", zorder=-10)
    ax.set_xlim(window)
    ax.set_ylabel("sp/s", fontsize=8)
    ax.tick_params(labelsize=8, length=2)
    for side in ("top", "right"):
        ax.spines[side].set_visible(False)


def plot_tuning(ax, session, unit_index, colors, scale):
    tuning, distances = session["tuning"], session["distances"]
    if tuning is None or distances is None:
        ax.text(0.5, 0.5, "no tuning", ha="center", va="center",
                transform=ax.transAxes, fontsize=8)
        ax.set_axis_off()
        return

    unit_label = str(session["units"].iloc[unit_index].get("unit_label", ""))
    match = tuning[tuning["unit_label"].astype(str) == unit_label] if "unit_label" in tuning else tuning.iloc[0:0]
    row = match.iloc[0] if not match.empty else tuning.iloc[unit_index]

    y = np.asarray(row["zfr_mean"], dtype=float) * scale
    yerr = np.asarray(row["zfr_sem"], dtype=float) * scale if "zfr_sem" in row else np.full_like(y, np.nan)
    x = distances[: len(y)]

    ax.plot(x, y, color="black", linewidth=0.8, zorder=1)
    for xi, yi, ei in zip(x, y, yerr):
        color = colors.get(float(xi), "black")
        ax.errorbar(xi, yi, yerr=None if np.isnan(ei) else ei, fmt="o",
                    color=color, ecolor=color, elinewidth=0.7, capsize=0,
                    markersize=3.5, zorder=2)
    ax.set_xlabel("Wall distance (mm)", fontsize=8)
    ax.set_ylabel("Firing rate (z-score)", fontsize=8)
    ticks = x[np.linspace(0, len(x) - 1, min(5, len(x)), dtype=int)]
    ax.set_xticks(ticks)
    ax.set_xticklabels([f"{v:.0f}" for v in ticks])
    ax.tick_params(labelsize=8, length=2)
    for side in ("top", "right"):
        ax.spines[side].set_visible(False)

5. Build the figure

One row per unit, ordered as in the manuscript. The shaded band marks the interval during which the wall is in motion past the animal's head.

def add_distance_legend(fig, distances, colors):
    ordered = np.array(sorted(distances), dtype=float)
    left, bottom, width, height = 0.10, 0.945, 0.20, 0.016
    ax = fig.add_axes([left, bottom, width, height])
    for i, distance in enumerate(ordered):
        ax.bar(i, 1, color=colors[float(distance)], width=1.0, align="edge")
    ax.set_xlim(0, len(ordered))
    ax.set_ylim(0, 1)
    ax.set_axis_off()
    fig.text(left, bottom + 0.020, "Wall distance (mm)", fontsize=8, ha="left")
    fig.text(left, bottom - 0.011, f"{ordered[0]:.0f}", fontsize=8, ha="left")
    fig.text(left + width, bottom - 0.011, f"{ordered[-1]:.0f}", fontsize=8, ha="right")


# Distances are shared across sessions; order far-to-near to match the manuscript.
all_distances = sorted({float(d) for s in sessions.values() for d in s["distances"]})
distances_desc = np.array(all_distances)[::-1]
colors = distance_colors(np.array(all_distances))

n_rows = len(FIGURE2_UNITS)
fig, axes = plt.subplots(
    n_rows, 3,
    figsize=(7.0, max(2.2, 1.45 * n_rows + 0.8)),
    sharex="col",
    gridspec_kw={"width_ratios": [1.0, 1.0, 0.9], "wspace": 0.35, "hspace": 0.18},
)

for row_idx, unit in FIGURE2_UNITS.iterrows():
    session = sessions[unit["asset_path"]]
    spike_times = np.asarray(session["units"].iloc[unit["unit_index"]]["spike_times"], dtype=float)
    row_distances = np.array(sorted(session["distances"]))[::-1]

    plot_raster(axes[row_idx, 0], spike_times, session["trials"], row_distances, colors, WINDOW)
    plot_peth(axes[row_idx, 1], spike_times, session["trials"], row_distances, colors,
              WINDOW, BIN_WIDTH, SMOOTH_BINS)

    step = PETH_TICK_STEPS[row_idx] if row_idx < len(PETH_TICK_STEPS) else 10
    upper = max(step, float(np.ceil(axes[row_idx, 1].get_ylim()[1] / step) * step))
    axes[row_idx, 1].set_ylim(0, upper)
    axes[row_idx, 1].set_yticks(np.arange(0, upper + step * 0.5, step))

    plot_tuning(axes[row_idx, 2], session, unit["unit_index"], colors, TUNING_SCALE)

    axes[row_idx, 0].set_ylabel(f"{unit['display_unit']}\n{unit['tuning_type']}",
                                fontsize=8, rotation=90, labelpad=12)
    axes[row_idx, 0].text(0.02, 0.94, unit["unit_label"],
                          transform=axes[row_idx, 0].transAxes,
                          fontsize=7, va="top", ha="left")

axes[0, 0].set_title("Rasters", fontsize=10)
axes[0, 1].set_title("PETHs", fontsize=10)
axes[0, 2].set_title("Tuning curves", fontsize=10)
axes[-1, 0].set_xlabel("Time (s)", fontsize=8)
axes[-1, 1].set_xlabel("Time (s)", fontsize=8)
fig.suptitle("Figure 2: peri-head distance coding in PrV", fontsize=11, x=0.99, y=0.988, ha="right")
add_distance_legend(fig, distances_desc, colors)
fig.subplots_adjust(top=0.915, left=0.10, right=0.98, bottom=0.06)
plt.show()

6. Reading the figure

• Unit 1 is a proximity unit: firing climbs steadily as the wall gets closer, so the whole tuning curve slopes in one direction.
• Units 2 to 5 are map units, each peaking at a preferred distance. Their PETH traces separate by colour, and the tuning curve has an interior maximum.
• Unit 6 is suppressed: the wall drives firing down rather than up.

Together these illustrate the paper's claim that PrV carries not one but two distance codes, a monotonic proximity signal and a distributed map.

7. Going further

• Other units. Every session's units table carries tuning_type, wall_responsive, pref_dist_mm and trial_include, so you can select units by classification instead of using the fixed list above.
• All sessions. client.get_dandiset(...).get_assets() lists all 191 assets. The processed files are the ones without a modality suffix.
• Raw data. The _ecephys+image assets hold the 30 kHz broadband, the 1 kHz LFP, the digital TTL record, and an ImageSeries pointing at the behaviour video, if you want to work from the unprocessed signals.