Reproducing Sargolini et al. (2006) from dandiset 000582¶

Paper: Sargolini F, Fyhn M, Hafting T, McNaughton BL, Witter MP, Moser M-B, Moser EI (2006). Conjunctive representation of position, direction, and velocity in entorhinal cortex. Science 312:758–762. doi:10.1126/science.1125572

Data: DANDI:000582, version 0.251111.2151, streamed from the archive. Nothing is downloaded to disk.

What is tested¶

The paper is not open access and its full text was not read for this attempt. The three claims below are taken from the abstract, so they are at the level of the paper's conclusions, not of individual figure panels.

Claim Statement in the abstract
layer2-grid "Layer II was predominated by grid cells"
deep-hd-conjunctive "grid cells colocalized with head-direction cells and conjunctive grid × head-direction cells in the deeper layers"
speed "All cell types were modulated by running speed"

Because the methods were not available, the analysis uses common choices from the grid-cell literature (5 cm bins, Gaussian smoothing, grid score from the spatial autocorrelogram, mean vector length for head direction, chance levels from time-shifted spike trains). They are defined in sargolini.py next to this notebook.

In [1]:
import warnings

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from scipy import stats

import sargolini as sg

warnings.filterwarnings("ignore", category=RuntimeWarning)
plt.rcParams.update({"figure.dpi": 100, "axes.spines.top": False, "axes.spines.right": False})
LAYERS = ["II", "III", "V", "VI"]

Load the data¶

Each file is one recording session: the positions of one or two tracking LEDs on the rat's head and the spike times of the cells recorded in that session. Each cell carries the cortical layer it was recorded in.

In [2]:
sessions = sg.load_sessions()
cells = pd.DataFrame(
    [dict(subject=s.subject, layer=u.layer, two_leds=s.head_direction is not None,
          spikes=len(u.spike_times), minutes=(s.t[-1] - s.t[0]) / 60, box_m=s.box)
     for s in sessions for u in s.units]
)
print(f"{len(sessions)} sessions, {cells.subject.nunique()} rats, {len(cells)} cells")
summary = cells.assign(layer=cells.layer.replace("", "not given")).groupby("layer").agg(
    cells=("subject", "size"), rats=("subject", "nunique"),
    cells_with_two_leds=("two_leds", "sum"), median_spikes=("spikes", "median"))
summary

118 sessions, 15 rats, 620 cells
Out[2]:
cells rats cells_with_two_leds median_spikes
layer
II 71 5 43 1413.0
III 227 7 182 1053.0
V 118 3 79 922.5
VI 198 3 143 1341.0
not given 6 2 6 597.5
In [3]:
boxes = pd.DataFrame([dict(box_m=s.box, minutes=round((s.t[-1] - s.t[0]) / 60), cells=len(s.units))
                      for s in sessions])
boxes.groupby("box_m").agg(sessions=("cells", "size"), cells=("cells", "sum"),
                           median_minutes=("minutes", "median"))
Out[3]:
sessions cells median_minutes
box_m
1.0 95 464 10.0
1.5 3 30 20.0
2.0 19 120 20.0
2.5 1 6 20.0

Four things about the dataset matter for what follows:

  • The enclosure is not always 1 m. Every file's description says the rat ran in a 1 x 1 m enclosure, but in 23 sessions the tracked positions span 1.5 to 2.5 m. The distance between the two LEDs and the running speeds are the same as in the other sessions, so these are larger enclosures, not a scaling error. Rate maps here use the size read from the data.

  • Head direction needs two LEDs. It is computed from the line between them. Sessions with one LED can be used for position and speed but not for head direction.

  • Cells are listed per session. The files give no identifier that links a cell across sessions, so the same cell may be counted more than once. Counts below are of cell-sessions.

  • The dandiset is described as a sample. Its metadata says it "includes spike times for recorded grid cells" and that "the sample includes conjunctive cells and head direction cells from layers III and V". How the cells were selected from those in the paper is not stated.

Measures for every cell¶

For each cell: a firing-rate map and its grid score, the strength of head-direction tuning (mean vector length, where two LEDs are available), and the correlation of firing rate with running speed. Chance levels come from repeating the first two measures on spike trains shifted in time by a random amount, which keeps the cell's firing pattern and the rat's path but breaks the link between them.

One chance level is used for all cells: the 95th percentile of the shuffled values pooled over cells. About 5% of untuned cells pass it on average, though not exactly 5% for each cell.

In [4]:
N_SHUFFLES = 50
rng = np.random.default_rng(582)

rows, shuffled_grid, shuffled_hd, shuffled_speed = [], [], [], []
for s in sessions:
    speed = sg.running_speed(s)
    for u in s.units:
        samples = sg.spike_samples(s, u.spike_times)
        has_hd = s.head_direction is not None
        rows.append(dict(
            path=s.path, subject=s.subject, unit=u.name, layer=u.layer, two_leds=has_hd,
            grid_score=sg.grid_score(sg.autocorrelogram(sg.rate_map(s, samples))),
            hd_length=sg.head_direction_tuning(s, samples)[1] if has_hd else np.nan,
            speed_r=sg.speed_tuning(s, samples, speed)[1],
        ))
        for _ in range(N_SHUFFLES):
            shifted = sg.spike_samples(s, sg.shifted(s, u.spike_times, rng))
            shuffled_grid.append(sg.grid_score(sg.autocorrelogram(sg.rate_map(s, shifted))))
            shuffled_speed.append(sg.speed_tuning(s, shifted, speed)[1])
            if has_hd:
                shuffled_hd.append(sg.head_direction_tuning(s, shifted)[1])

df = pd.DataFrame(rows)
grid_chance = np.nanpercentile(shuffled_grid, 95)
hd_chance = np.nanpercentile(shuffled_hd, 95)
speed_chance = np.nanpercentile(shuffled_speed, 95)
df["is_grid"] = df.grid_score > grid_chance
df["is_hd"] = df.hd_length > hd_chance
print(f"Chance level (95th percentile of all shuffles pooled): grid score {grid_chance:.2f}, "
      f"head-direction vector length {hd_chance:.2f}, speed correlation {speed_chance:.2f}")

Chance level (95th percentile of all shuffles pooled): grid score 0.61, head-direction vector length 0.19, speed correlation 0.76

Examples¶

The cell with the highest grid score, the one with the strongest head-direction tuning, and the conjunctive cell with the highest grid score. Top: firing-rate map. Middle: spatial autocorrelogram. Bottom: firing rate by head direction.

In [5]:
by_path = {s.path: s for s in sessions}
both = df[df.is_grid & df.is_hd]
examples = {
    "Grid cell": df.loc[df[~df.is_hd.fillna(False)].grid_score.idxmax()],
    "Head-direction cell": df.loc[df[~df.is_grid].hd_length.idxmax()],
    "Conjunctive cell": both.loc[both.grid_score.idxmax()],
}
fig, axes = plt.subplots(3, 3, figsize=(8, 8))
angles = np.arange(-180 + sg.HD_BIN / 2, 180, sg.HD_BIN)
for column, (label, row) in enumerate(examples.items()):
    s = by_path[row.path]
    unit = next(u for u in s.units if u.name == row.unit)
    samples = sg.spike_samples(s, unit.spike_times)
    rate = sg.rate_map(s, samples)
    axes[0, column].imshow(rate.T, origin="lower", cmap="viridis")
    axes[0, column].set_title(f"{label}\nlayer {row.layer}, peak {np.nanmax(rate):.0f} Hz", fontsize=9)
    axes[1, column].imshow(sg.autocorrelogram(rate).T, origin="lower", cmap="RdBu_r", vmin=-1, vmax=1)
    axes[1, column].set_title(f"grid score {row.grid_score:.2f}", fontsize=9)
    for ax in axes[:2, column]:
        ax.set_xticks([]); ax.set_yticks([])
    if s.head_direction is not None:
        axes[2, column].plot(angles, sg.head_direction_tuning(s, samples)[0], color="k")
        axes[2, column].set_title(f"vector length {row.hd_length:.2f}", fontsize=9)
    else:
        axes[2, column].text(0.5, 0.5, "one LED:\nno head direction", ha="center", va="center",
                             transform=axes[2, column].transAxes)
    axes[2, column].set_xlabel("Head direction (deg)")
axes[2, 0].set_ylabel("Firing rate (Hz)")
plt.tight_layout()
No description has been provided for this image

Claim layer2-grid: "Layer II was predominated by grid cells"¶

All cells with a known layer and a grid score are used, since the grid score needs only position. (A few cells with too few spikes for an autocorrelogram have no score.)

In [6]:
known = df[df.layer.isin(LAYERS) & df.grid_score.notna()]

def fraction_table(frame, column):
    out = frame.groupby("layer")[column].agg(n="sum", cells="size")
    out = out.astype(int)
    ci = [stats.binomtest(int(n), int(c)).proportion_ci() for n, c in zip(out.n, out.cells)]
    out["fraction"] = (out.n / out.cells).round(2)
    out["95% CI"] = [f"{lo:.2f} to {hi:.2f}" for lo, hi in ci]
    return out.loc[LAYERS]

print(f"Grid cells: grid score above chance ({grid_chance:.2f})")
fraction_table(known, "is_grid")

Grid cells: grid score above chance (0.61)
Out[6]:
n cells fraction 95% CI
layer
II 22 70 0.31 0.21 to 0.44
III 76 225 0.34 0.28 to 0.40
V 20 114 0.18 0.11 to 0.26
VI 34 193 0.18 0.13 to 0.24
In [7]:
fig, axes = plt.subplots(1, 2, figsize=(9, 3.4))
bins = np.linspace(-1.3, 1.6, 30)
for layer in LAYERS:
    scores = known[known.layer == layer].grid_score.dropna()
    axes[0].hist(scores, bins=bins, histtype="step", density=True, label=f"Layer {layer} (n={len(scores)})")
axes[0].axvline(grid_chance, color="k", ls="--", lw=1)
axes[0].text(grid_chance, axes[0].get_ylim()[1], " chance", va="top", fontsize=8)
axes[0].set_xlabel("Grid score"); axes[0].set_ylabel("Density"); axes[0].legend(fontsize=8, frameon=False)

thresholds = np.linspace(-0.2, 1.0, 25)
for layer in LAYERS:
    scores = known[known.layer == layer].grid_score
    axes[1].plot(thresholds, [(scores > th).mean() for th in thresholds], label=f"Layer {layer}")
axes[1].axvline(grid_chance, color="k", ls="--", lw=1)
axes[1].axhline(0.5, color="0.7", lw=1)
axes[1].set_xlabel("Grid-score threshold"); axes[1].set_ylabel("Fraction of cells above threshold")
plt.tight_layout()

layer2 = known[known.layer == "II"].grid_score.dropna()
deeper = known[known.layer != "II"].grid_score.dropna()
test = stats.mannwhitneyu(layer2, deeper)
print(f"Median grid score: layer II {layer2.median():.2f}, deeper layers {deeper.median():.2f} "
      f"(Mann-Whitney p = {test.pvalue:.1e})")
for th in (0.0, 0.3):
    print(f"Fraction with grid score > {th}: "
          + ", ".join(f"{l} {(known[known.layer == l].grid_score > th).mean():.2f}" for l in LAYERS))

Median grid score: layer II 0.33, deeper layers 0.10 (Mann-Whitney p = 2.8e-03)
Fraction with grid score > 0.0: II 0.81, III 0.64, V 0.54, VI 0.52
Fraction with grid score > 0.3: II 0.54, III 0.49, V 0.27, VI 0.30
No description has been provided for this image

Cells from one rat are not independent, and the same cell may appear in several sessions, so the intervals and the p-value above are too narrow. The same comparison for each rat:

In [8]:
grid_per_rat = (known.assign(layers=np.where(known.layer == "II", "layer II", "layers III, V, VI"))
                .groupby(["layers", "subject"])
                .agg(cells=("grid_score", "size"), median_grid_score=("grid_score", "median"),
                     grid_cells=("is_grid", "sum")))
grid_per_rat["fraction"] = (grid_per_rat.grid_cells / grid_per_rat.cells).round(2)
grid_per_rat.round(2)
Out[8]:
cells median_grid_score grid_cells fraction
layers subject
layer II 10073 8 0.33 2 0.25
10697 9 0.53 2 0.22
10884 10 0.33 4 0.40
11265 26 0.75 14 0.54
11278 17 0.04 0 0.00
layers III, V, VI 10704 43 0.26 15 0.35
10884 37 0.28 8 0.22
10938 8 0.09 2 0.25
10962 11 -0.22 2 0.18
11016 12 0.03 3 0.25
11025 38 0.32 12 0.32
11084 52 -0.02 4 0.08
11138 44 -0.02 13 0.30
11207 189 0.03 30 0.16
11265 85 0.37 37 0.44
11340 10 0.07 4 0.40
11343 3 0.24 0 0.00

Result: partially reproduced. Layer II has a higher median grid score than the deeper layers. Against the chance level from shuffling, about a third of layer II cells are grid cells, the same as in layer III and about twice the fraction in layers V and VI.

At a lower threshold of 0.3, 54% of layer II cells pass, the most of any layer. That figure should be read with care: 0.3 is a common threshold for the standard grid score, but the score here is the best over several ring sizes and its chance level is about 0.6, so many cells pass 0.3 by chance; and 54% of 70 cells cannot be told apart from half, or from layer III's 49%.

The per-rat table shows how much rests on few animals. The layer II cells come from five rats, whose grid-cell fractions range from 0 of 17 to 14 of 26, so the layer II figure is not a stable estimate.

So grid cells are more common in layer II than in layers V and VI, and the next section shows they are the only tuned cell type found there above chance. Whether they are the majority of layer II cells depends on the criterion for a grid cell, and the paper's criterion was not available. Two more things limit the comparison: the dandiset is a sample whose selection is not described, and most sessions are 10 minutes long, which makes grid scores noisy.

Claim deep-hd-conjunctive: head-direction and conjunctive cells in the deeper layers¶

Only cells from sessions with two LEDs are used.

In [9]:
tracked = known[known.two_leds].copy()
tracked["kind"] = np.select(
    [tracked.is_grid & tracked.is_hd, tracked.is_grid, tracked.is_hd],
    ["conjunctive", "grid only", "head direction only"], "neither")
KINDS = ["grid only", "conjunctive", "head direction only", "neither"]
counts = pd.crosstab(tracked.layer, tracked.kind).reindex(index=LAYERS, columns=KINDS, fill_value=0)
counts["cells"] = counts.sum(axis=1)
counts
Out[9]:
kind grid only conjunctive head direction only neither cells
layer
II 12 2 3 26 43
III 27 39 65 49 180
V 1 16 54 6 77
VI 10 7 95 27 139
In [10]:
fractions = counts[KINDS].div(counts.cells, axis=0)
fig, axes = plt.subplots(1, 2, figsize=(9, 3.4))
bottom = np.zeros(len(LAYERS))
colors = {"grid only": "#3b6fb6", "conjunctive": "#8e5ea2", "head direction only": "#d1803a", "neither": "0.8"}
for kind in KINDS:
    axes[0].bar(LAYERS, fractions[kind], bottom=bottom, label=kind, color=colors[kind])
    bottom += fractions[kind].values
axes[0].set_xlabel("Layer"); axes[0].set_ylabel("Fraction of cells"); axes[0].legend(fontsize=8, frameon=False, loc="lower left")

axes[1].boxplot([tracked[tracked.layer == l].hd_length.dropna() for l in LAYERS], tick_labels=LAYERS, showfliers=False)
axes[1].axhline(hd_chance, color="k", ls="--", lw=1)
axes[1].set_xlabel("Layer"); axes[1].set_ylabel("Head-direction vector length")
plt.tight_layout()

hd_any = tracked.assign(hd=tracked.is_hd.astype(bool))
print("Fraction of cells with head-direction tuning above chance (about 5% expected by chance):")
print(fraction_table(hd_any, "hd").to_string())
ii = tracked[tracked.layer == "II"]; deep = tracked[tracked.layer != "II"]
table = [[int(ii.is_hd.sum()), int((~ii.is_hd).sum())], [int(deep.is_hd.sum()), int((~deep.is_hd).sum())]]
print(f"\nLayer II against deeper layers, Fisher exact p = {stats.fisher_exact(table)[1]:.1e}")
print(f"Conjunctive cells: layer II {int((ii.kind == 'conjunctive').sum())} of {len(ii)}, "
      f"deeper layers {int((deep.kind == 'conjunctive').sum())} of {len(deep)}")

Fraction of cells with head-direction tuning above chance (about 5% expected by chance):
         n  cells  fraction        95% CI
layer                                    
II       5     43      0.12  0.04 to 0.25
III    104    180      0.58  0.50 to 0.65
V       70     77      0.91  0.82 to 0.96
VI     102    139      0.73  0.65 to 0.81

Layer II against deeper layers, Fisher exact p = 9.5e-14
Conjunctive cells: layer II 2 of 43, deeper layers 62 of 396
No description has been provided for this image

The cells come from few rats, and the same cell may appear in several sessions, so the counts above are not independent and the intervals and p-value are too narrow. As a check that does not depend on that, the same fraction for each rat separately:

In [11]:
per_rat = (tracked.assign(layers=np.where(tracked.layer == "II", "layer II", "layers III, V, VI"))
           .groupby(["layers", "subject"]).is_hd.agg(head_direction="sum", cells="size"))
per_rat["fraction"] = (per_rat.head_direction / per_rat.cells).round(2)
per_rat
Out[11]:
head_direction cells fraction
layers subject
layer II 11265 3 26 0.12
11278 2 17 0.12
layers III, V, VI 11025 9 38 0.24
11084 22 27 0.81
11138 26 44 0.59
11207 150 189 0.79
11265 63 85 0.74
11340 4 10 0.40
11343 2 3 0.67

Result: reproduced. In layers III, V and VI most cells are head-direction tuned, and grid cells, head-direction cells and conjunctive cells are all found there together. In layer II, 12% of cells pass the head-direction threshold, which cannot be told apart from the roughly 5% expected by chance with this many cells, and 2 of 43 are conjunctive.

The difference holds within the one rat recorded in both layer II and the deeper layers (11265: 12% against 74%). Of the six other rats with deep-layer recordings, three have fractions of 59% to 81%, one has 24%, and two have ten cells or fewer.

Claim speed: "All cell types were modulated by running speed"¶

For each cell, the firing rate is computed in 5 cm/s bins of running speed, and correlated with speed. The correlation rests on at most ten bins, so it is noisy for a single cell; the chance level is the 95th percentile of the same correlation for time-shifted spike trains. The same two-LED cells and cell types as above are used; "neither" is what is left over, not a cell type of the paper.

In [12]:
fig, axes = plt.subplots(1, 2, figsize=(9, 3.4))
centres = (sg.SPEED_EDGES[:-1] + sg.SPEED_EDGES[1:]) / 2 * 100
speed_cache = {}
speed_rows = []
for kind in KINDS:
    curves = []
    for row in tracked[tracked.kind == kind].itertuples():
        s = by_path[row.path]
        if row.path not in speed_cache:
            speed_cache[row.path] = sg.running_speed(s)
        unit = next(u for u in s.units if u.name == row.unit)
        rate = sg.speed_tuning(s, sg.spike_samples(s, unit.spike_times), speed_cache[row.path])[0]
        if np.nanmean(rate) > 0:
            curves.append(rate / np.nanmean(rate))
    curves = np.array(curves)
    mean = np.nanmean(curves, axis=0)
    sem = np.nanstd(curves, axis=0) / np.sqrt(np.sum(~np.isnan(curves), axis=0))
    axes[0].errorbar(centres, mean, sem, label=kind, color=colors[kind] if kind != "neither" else "0.5")
    r = tracked[tracked.kind == kind].speed_r.dropna()
    above = int((r > speed_chance).sum())
    speed_rows.append(dict(
        kind=kind, cells=len(r), median_r=round(r.median(), 2),
        fraction_positive=round((r > 0).mean(), 2),
        fraction_above_chance=round(above / len(r), 2),
        p_above_chance=stats.binomtest(above, len(r), 0.05, alternative="greater").pvalue))
axes[0].set_xlabel("Running speed (cm/s)"); axes[0].set_ylabel("Firing rate / cell's mean rate")
axes[0].legend(fontsize=8, frameon=False)
axes[1].boxplot([tracked[tracked.kind == k].speed_r.dropna() for k in KINDS],
                tick_labels=[k.replace(" ", "\n") for k in KINDS], showfliers=False)
axes[1].axhline(0, color="k", lw=1)
axes[1].set_ylabel("Correlation of rate with speed")
plt.tight_layout()
print(f"Shuffled spike trains: median correlation {np.nanmedian(shuffled_speed):.2f}, "
      f"95th percentile {speed_chance:.2f}")
pd.DataFrame(speed_rows).set_index("kind")

Shuffled spike trains: median correlation -0.06, 95th percentile 0.76
Out[12]:
cells median_r fraction_positive fraction_above_chance p_above_chance
kind
grid only 50 0.58 0.80 0.30 1.287458e-08
conjunctive 64 0.58 0.78 0.20 1.443835e-05
head direction only 217 0.28 0.62 0.23 5.519612e-19
neither 108 0.54 0.76 0.33 2.345802e-20
No description has been provided for this image

Result: reproduced. In every cell type the firing rate rises with running speed. The median correlation is positive in each group, where time-shifted spike trains give a median near zero, and in each group four to seven times more cells pass the chance level than the 5% expected.

This is a statement about each type as a population, which is how the abstract's "all cell types" is read here. For a single cell the correlation is too noisy to say much, and most individual cells do not pass the chance level.

Summary¶

Claim Outcome Note
layer2-grid Partially reproduced Layer II has the highest grid scores and, at a grid-score threshold of 0.3, a small majority of grid cells. At the stricter shuffle-based threshold it has about a third, the same as layer III. The paper's criterion was not available.
deep-hd-conjunctive Reproduced Head-direction and conjunctive cells are found with grid cells in layers III, V and VI and are rare in layer II.
speed Reproduced Firing rate increases with running speed in every cell type, as a population.

What this attempt found about the dandiset¶

  • The enclosure size is wrong in the metadata of 23 sessions. Every file says 1 x 1 m, but positions span 1.5 to 2.5 m in those sessions. An analysis that trusts the description discards most of the data in them.
  • No identifier links a cell across sessions, so cells may be counted more than once.
  • The dandiset does not say how its cells were selected from those in the paper.
  • The second tracking LED is missing for 47 of the 118 sessions, so head direction cannot be computed for 167 of the 620 cell-sessions.
  • The paper is not open access, and the dandiset does not give the analysis parameters.