Churchland et al. (2010) on dandiset 000139 (MC_Maze_Medium)¶
Paper. Churchland MM, Cunningham JP, Kaufman MT, Ryu SI, Shenoy KV (2010). Cortical preparatory activity: representation of movement or first cog in a dynamical machine? Neuron 68, 387-400. doi:10.1016/j.neuron.2010.09.015.
Dataset. DANDI:000139, "MC_Maze_Medium", version 0.220113.0408, streamed from the archive. This is the only published version of the dandiset (the other is the draft), so it is pinned; the paper cites no version. The dandiset holds one session of monkey Jenkins (2009-09-29) on the maze task, recorded with two Utah arrays, cut for the Neural Latents Benchmark into 250 training and validation trials (with behaviour and trial labels) and 100 test trials (700 ms segments, no condition labels, held-in units only).
Claims. claims/000139/10.1016-j.neuron.2010.09.015.yaml, written from the NIH author manuscript
(PMC2991102) with its supplement, and checked by a second reader. Methods the claims file leaves
open are taken from the paper (cited where used) or, where the paper is silent, run under every
reading the claims file names. This version answers the first and second verification rounds: the outputs they
asked for are added in the sections they concern, each marked as added in response to the verifier.
| Claim | Tested here |
|---|---|
| j-array-recordings | tested |
| fig4d-mean-correlation-j-array | tested |
| fig4abc-mean-correlation-other-datasets | blocked: data_missing (speed-task and single-electrode data not in the dandiset) |
| fig4-broad-unimodal | tested |
| fig4e-correlation-time-course | tested |
| suppfig1-population-psth-weak-movement-tuning | tested |
| fig6d-perimovement-best-j-array | tested |
| fig6-perimovement-best-each-dimensionality | tested |
| fig6-dimensionality-peak | tested |
| fig6-maze-task-spaces | tested |
| fig6abc-perimovement-best-other-datasets | blocked: data_missing |
| fig6e-summary-rank | blocked: data_missing |
| fig6-preparatory-dimensionality | tested |
| fig7a-shuffled-space-poor | tested |
| fig7b-window-start-slide | tested |
| ap-location-no-effect-j-array | tested (AP read as the array) |
| ap-location-no-effect-other-datasets | blocked: data_missing |
| fig7c-anterior-posterior-performance | tested (AP read as the array) |
| prep-activity-more-anterior | run; blocked after the analysis (metadata_insufficient: which array is PMd) |
| fig7d-pd-consistency | tested |
| some-neurons-better-elsewhere | tested |
| high-dimensional-task-spaces-underperform | tested |
| fig8-linear-dynamics-simulation | tested (simulation) |
import itertools
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy import stats
import diptest
import maze as M
pd.set_option("display.width", 200)
pd.set_option("display.max_columns", 40)
pd.set_option("display.max_rows", 200)
D = M.load(M.TRAIN_PATH)
TE = M.load(M.TEST_PATH)
tr = D["trials"]
units = D["units"]
assets = M.list_assets()
Data check¶
print("Assets at", M.VERSION)
for a in assets:
print(f" {a['path']} {a['size']} bytes")
print("Subject:", D["nwb"].subject.subject_id, D["nwb"].subject.species, "| session start", D["nwb"].session_start_time)
print("Train file: trials", len(tr), "| splits", tr.split.value_counts().to_dict(), "| success", tr.success.value_counts().to_dict())
print("Test file: trials", len(TE["trials"]), "columns", list(TE["trials"].columns), "| units", len(TE["units"]))
cond = (tr.trial_type * 10 + tr.trial_version).values
conds = np.unique(cond)
print("Conditions (trial_type x trial_version):", len(conds), "| maze_id values:", tr.maze_id.nunique(),
"| trial_version counts:", tr.trial_version.value_counts().sort_index().to_dict())
print("Each condition has one maze_id:", tr.groupby(["trial_type", "trial_version"]).maze_id.nunique().max() == 1)
print("Trials per condition: min", pd.Series(cond).value_counts().min(), "max", pd.Series(cond).value_counts().max())
print("num_targets by version:", tr.groupby("trial_version").num_targets.unique().to_dict(),
"| barriers by version:", tr.groupby("trial_version").num_barriers.apply(lambda x: sorted(x.unique())).to_dict())
print("Delay (ms): range", tr.delay.min(), "to", tr.delay.max(), "| trials with delay > 400 ms:", (tr.delay > 400).sum())
print("Delay column equals go - target times (ms): max |diff|",
np.abs((tr.go_cue_time - tr.target_on_time) * 1000 - tr.delay).max().round(6))
print("rt column equals move - go times (ms): max |diff|",
np.abs((tr.move_onset_time - tr.go_cue_time) * 1000 - tr.rt).max().round(6), "| rt range", tr.rt.min(), tr.rt.max(),
"| median", tr.rt.median())
print("Trial duration (s):", (tr.stop_time - tr.start_time).min().round(3), "to", (tr.stop_time - tr.start_time).max().round(3),
"| movement onset to trial end at least", (tr.stop_time - tr.move_onset_time).min().round(3), "s",
"| trial start to movement onset at least", (tr.move_onset_time - tr.start_time).min().round(3), "s")
print()
print("Units:", len(units), "| heldout", units.heldout.sum(), "| held-in", (~units.heldout).sum())
print("Spike times sorted in every unit:", all(np.all(np.diff(s) >= 0) for s in D["spikes"]))
print("Spikes per unit: min", min(len(s) for s in D["spikes"]), "max", max(len(s) for s in D["spikes"]))
print("Test-file units are the held-in train units:", set(TE["units"].id) == set(units.id[~units.heldout]))
el = D["electrodes"]
print("Electrode table: rows", len(el), "| groups", el.group_name.value_counts().to_dict(),
"| x, y, z all NaN:", bool(el[["x", "y", "z"]].isna().all().all()))
print("Electrode group referenced by units:", units.electrode_group_in_file.value_counts().to_dict())
print("Unit id is array*1000 + channel*10 + sort, with channel = electrode row + 1:",
bool(np.all(units.id_channel == units.electrode_row + 1)))
print("Units by id prefix (array):", units.id_array.value_counts().sort_index().to_dict(),
"| electrode rows referenced, prefix 1:", units.electrode_row[units.id_array == 1].min(), "to",
units.electrode_row[units.id_array == 1].max(), "; prefix 2:", units.electrode_row[units.id_array == 2].min(),
"to", units.electrode_row[units.id_array == 2].max())
# duplicate recordings: units with coincident spikes
st = D["spikes"]
coinc = []
for a, b in itertools.combinations(range(len(st)), 2):
x, y = st[a], st[b]
j = np.clip(np.searchsorted(y, x), 1, len(y) - 1)
coinc.append((np.minimum(np.abs(y[j] - x), np.abs(y[j - 1] - x)) < 2e-4).mean())
coinc = np.array(coinc)
print(f"Unit pairs: largest fraction of spikes within 0.2 ms of the other unit's {coinc.max():.3f} (99th percentile {np.percentile(coinc, 99):.4f})")
# test segments against train trials
hin = np.where(~units.heldout.values)[0]
test_ids = list(TE["units"].id)
order = [int(np.where(units.id.values == i)[0][0]) for i in test_ids]
def counts(sts, t0, t1):
return np.array([np.sum((s >= t0) & (s < t1)) for s in sts])
A = np.array([counts(TE["spikes"], r.move_onset_time - 0.25, r.move_onset_time + 0.45) for r in TE["trials"].itertuples()])
B = np.array([counts([st[i] for i in order], r.move_onset_time - 0.25, r.move_onset_time + 0.45) for r in tr.itertuples()])
L1 = np.abs(A[:, None, :] - B[None, :, :]).sum(2)
print("Test segments vs train trials (same 700 ms around onset): smallest summed |count difference|", L1.min(),
"| median", np.median(L1), "-> no test segment duplicates a train trial")
print()
tb = D["t_beh"]
print("Behaviour samples:", len(tb), "| sample interval within trials (ms):", np.round(np.median(np.diff(tb)) * 1000, 3),
"| hand x range (mm)", D["hand_pos"][:, 0].min().round(1), "to", D["hand_pos"][:, 0].max().round(1),
"| y", D["hand_pos"][:, 1].min().round(1), "to", D["hand_pos"][:, 1].max().round(1))
sp = np.hypot(D["hand_vel"][:, 0], D["hand_vel"][:, 1])
print("Hand speed (mm/s): median", np.median(sp).round(1), "| 99.9th percentile", np.percentile(sp, 99.9).round(1))
vnum = np.gradient(D["hand_pos"][:, 0], tb)
inside = np.abs(np.diff(tb, prepend=tb[0])) < 0.0015
print("hand_vel x against the derivative of hand_pos x: r =", np.corrcoef(vnum[inside], D["hand_vel"][inside, 0])[0, 1].round(3),
"| slope", np.polyfit(vnum[inside], D["hand_vel"][inside, 0], 1)[0].round(3), "(1 means hand_vel is in mm/s)")
Assets at 0.220113.0408
sub-Jenkins/sub-Jenkins_ses-medium_desc-test_ecephys.nwb 695928 bytes
sub-Jenkins/sub-Jenkins_ses-medium_desc-train_behavior+ecephys.nwb 76604764 bytes
Subject: Jenkins Macaca mulatta | session start 2009-09-29 00:00:00-07:00
Train file: trials 250 | splits {'train': 188, 'val': 62} | success {True: 250}
Test file: trials 100 columns ['start_time', 'stop_time', 'move_onset_time', 'split'] | units 114
Conditions (trial_type x trial_version): 27 | maze_id values: 9 | trial_version counts: {0: 81, 1: 86, 2: 83}
Each condition has one maze_id: True
Trials per condition: min 8 max 10
num_targets by version: {0: array([1]), 1: array([1]), 2: array([3])} | barriers by version: {0: [np.int64(0)], 1: [np.int64(8), np.int64(9)], 2: [np.int64(8), np.int64(9)]}
Delay (ms): range -3 to 998 | trials with delay > 400 ms: 235
Delay column equals go - target times (ms): max |diff| 0.0
rt column equals move - go times (ms): max |diff| 0.0 | rt range 247 989 | median 355.0
Trial duration (s): 2.066 to 4.066 | movement onset to trial end at least 0.959 s | trial start to movement onset at least 0.988 s
Units: 152 | heldout 38 | held-in 114
Spike times sorted in every unit: True
Spikes per unit: min 89 max 9891
Test-file units are the held-in train units: True
Electrode table: rows 192 | groups {'electrode_group_PMd': 96, 'electrode_group_M1': 96} | x, y, z all NaN: True
Electrode group referenced by units: {'electrode_group_PMd': 152}
Unit id is array*1000 + channel*10 + sort, with channel = electrode row + 1: True
Units by id prefix (array): {1: 75, 2: 77} | electrode rows referenced, prefix 1: 0 to 89 ; prefix 2: 14 to 94
Unit pairs: largest fraction of spikes within 0.2 ms of the other unit's 0.029 (99th percentile 0.0086)
Test segments vs train trials (same 700 ms around onset): smallest summed |count difference| 135 | median 240.0 -> no test segment duplicates a train trial Behaviour samples: 757450 | sample interval within trials (ms): 1.0 | hand x range (mm) -134.5 to 158.5 | y -172.6 to 99.5 Hand speed (mm/s): median 5.6 | 99.9th percentile 1008.1 hand_vel x against the derivative of hand_pos x: r = 1.0 | slope 1.0 (1 means hand_vel is in mm/s)
Clocks¶
Spikes, behaviour and the trial table share one clock if a task variable can be decoded from spikes after the event it follows and not before it. The maze condition is decoded (27 classes, nearest class mean, leave one trial out) from spike counts of all units in windows around target onset, with true labels and with labels permuted. Hand speed is shown around the dataset's movement onset.
def decode(X, y):
ok = 0
for i in range(len(y)):
m = np.ones(len(y), bool); m[i] = False
cs = np.unique(y[m])
cen = np.array([X[m & (y == c)].mean(0) for c in cs])
ok += cs[np.argmin(((cen - X[i]) ** 2).sum(1))] == y[i]
return ok / len(y)
rng_clk = np.random.default_rng(1)
rows = []
for name, ev in [("target", tr.target_on_time.values), ("go", tr.go_cue_time.values)]:
for w in [(-300, 0), (50, 350)]:
X = np.sqrt(np.array([counts(st, e + w[0] / 1000, e + w[1] / 1000) for e in ev]))
acc = decode(X, cond)
perm = np.mean([decode(X, rng_clk.permutation(cond)) for _ in range(5)])
rows.append(dict(event=name, window_ms=w, accuracy=round(acc, 3), permuted_labels=round(perm, 3)))
print(pd.DataFrame(rows).to_string(index=False), "| chance", round(1 / len(conds), 3))
vel_all = M.kin_traces(tb, D["hand_vel"], tr.move_onset_time.values, (-300, 300))
spd = np.hypot(vel_all[..., 0], vel_all[..., 1])
print("Mean hand speed around movement onset (mm/s): -300..-100 ms", spd[:, :200].mean().round(1),
"| -50..0", spd[:, 250:300].mean().round(1), "| 0..50", spd[:, 300:350].mean().round(1),
"| 100..300", spd[:, 400:].mean().round(1))
on_b = M.onset_from_speed(tb, D["hand_vel"], tr.go_cue_time.values)
dd = (on_b - tr.move_onset_time.values) * 1000
print("Onset from hand speed (10% of peak) minus the dataset's move_onset_time (ms): median", np.median(dd),
"| 5th to 95th percentile", np.percentile(dd, 5), "to", np.percentile(dd, 95))
event window_ms accuracy permuted_labels
target (-300, 0) 0.024 0.029
target (50, 350) 0.272 0.034
go (-300, 0) 0.336 0.026
go (50, 350) 0.452 0.035 | chance 0.037
Mean hand speed around movement onset (mm/s): -300..-100 ms 3.3 | -50..0 17.4 | 0..50 205.7 | 100..300 549.1
Onset from hand speed (10% of peak) minus the dataset's move_onset_time (ms): median 12.000000000000455 | 5th to 95th percentile 4.9999999999954525 to 27.000000000030866
Trial selection and condition means¶
Every analysis below starts from condition means under a reading of the trial selection: the movement onset (a, the dataset's; b, from hand speed) and the proto-trial rule (Supplemental Experimental Procedures, Task design: "the subset of trials (typically ~90% of all trials) where the horizontal and vertical velocity profiles correlated with those of a reference 'proto-trial' with r > 0.9"), proto (a) the condition mean or (b) a representative trial (the medoid), and the profiles compared (a) separately or (b) concatenated. Only trials with delays above 400 ms are kept (Methods, Task design). The main reading is onset (a), proto (a) the condition mean, and the profiles concatenated (b): comparing the profiles separately keeps 70 to 82% of trials, against the paper's "typically ~90%", while concatenating them keeps 93 to 96% (table below). The condition mean over the medoid is an assumption (a template is usually a mean); both are run.
RD = {}
def condition_data(onset="a", proto="mean", combine="separate"):
key = (onset, proto, combine)
if key in RD:
return RD[key]
mo = tr.move_onset_time.values if onset == "a" else on_b
mv = M.bin_spikes(st, mo, M.MOVE_WIN)
go = M.bin_spikes(st, tr.go_cue_time.values, M.GO_WIN)
pos = M.kin_traces(tb, D["hand_pos"], mo, M.MOVE_WIN)
vel = M.kin_traces(tb, D["hand_vel"], mo, M.MOVE_WIN)
r0 = M.MOVE_WIN[0]
keep_proto = M.proto_keep(vel[:, M.TMAT[0] - r0: M.TMAT[1] - r0], cond, proto, combine)
sel = (tr.delay.values > 400) & keep_proto
c = cond[sel]
out = dict(sel=sel, keep_proto=keep_proto,
go_rate=M.smooth_rate(M.cond_mean(go[sel].astype(np.float32), c, conds)),
mv_rate=M.smooth_rate(M.cond_mean(mv[sel].astype(np.float32), c, conds)),
pos=M.cond_mean(pos[sel], c, conds), vel=M.cond_mean(vel[sel], c, conds),
n_per_cond=pd.Series(c).value_counts().reindex(conds).values)
out["prep"] = M.window_mean(out["go_rate"], M.PREP, M.GO_WIN[0])
out["peri"] = M.window_mean(out["mv_rate"], M.PERI, M.MOVE_WIN[0])
RD[key] = out
return out
SEL_READINGS = [(o, p, c) for o in "ab" for p in ["mean", "medoid"] for c in ["separate", "concat"]]
rows = []
for r in SEL_READINGS:
d = condition_data(*r)
rows.append(dict(onset=r[0], proto=r[1], combine=r[2], proto_pass=int(d["keep_proto"].sum()),
proto_pass_pct=round(100 * d["keep_proto"].mean(), 1), kept=int(d["sel"].sum()),
min_per_cond=int(d["n_per_cond"].min()), max_per_cond=int(d["n_per_cond"].max())))
sel_table = pd.DataFrame(rows)
print(sel_table.to_string(index=False))
MAIN = ("a", "mean", "concat")
d0 = condition_data(*MAIN)
print("Preparatory rates (spikes/s): range", d0["prep"].min().round(2), "to", d0["prep"].max().round(2),
"| peri-movement", d0["peri"].min().round(2), "to", d0["peri"].max().round(2))
onset proto combine proto_pass proto_pass_pct kept min_per_cond max_per_cond
a mean separate 198 79.2 185 2 10
a mean concat 239 95.6 225 6 10
a medoid separate 175 70.0 164 2 9
a medoid concat 233 93.2 218 5 10
b mean separate 205 82.0 192 5 10
b mean concat 241 96.4 226 6 10
b medoid separate 175 70.0 163 2 9
b medoid concat 234 93.6 220 6 10
Preparatory rates (spikes/s): range 0.0 to 34.32 | peri-movement 0.0 to 37.84
Neurons¶
Preparatory tuning above 5 spikes/s (array threshold, Methods, Neural recordings and datasets), and for the correlations also peri-movement tuning at least 1.5 times as strong (Methods, Correlations). Tuning strength: reading (a) range across conditions, (b) standard deviation.
def neuron_sets(d, how):
tp = M.tuning_strength(d["prep"], how)
tq = M.tuning_strength(d["peri"], how)
prep_ok = tp > 5
both = prep_ok & (tq >= 1.5 * tp)
return prep_ok, both
rows = []
for r in SEL_READINGS:
for how in ["range", "sd"]:
p, b = neuron_sets(condition_data(*r), how)
rows.append(dict(onset=r[0], proto=r[1], combine=r[2], tuning=how, prep_pass=int(p.sum()), both_pass=int(b.sum())))
neuron_table = pd.DataFrame(rows)
print(neuron_table.to_string(index=False))
tp0 = M.tuning_strength(d0["prep"], "range")
print("Main reading, preparatory tuning range (spikes/s): median", np.median(tp0).round(2), "| max", tp0.max().round(2))
# split-half reliability of preparatory tuning: odd against even trials of each condition
sel0 = d0["sel"]
go0 = M.bin_spikes(st, tr.go_cue_time.values, M.PREP).sum(2) / ((M.PREP[1] - M.PREP[0]) / 1000)
halves = []
for h in (0, 1):
m = np.zeros(len(tr), bool)
for c in conds:
ii = np.where(sel0 & (cond == c))[0]
m[ii[h::2]] = True
halves.append(np.stack([go0[m & (cond == c)].mean(0) for c in conds]))
pdn = neuron_sets(d0, "range")[0]
rel = np.array([np.corrcoef(halves[0][:, n], halves[1][:, n])[0, 1] for n in np.where(pdn)[0]])
print("Split-half reliability of the preparatory rate vector (odd vs even trials), neurons over 5 spikes/s: median r",
np.round(np.median(rel), 3), "| Spearman-Brown corrected", np.round(2 * np.median(rel) / (1 + np.median(rel)), 3),
"| fraction with r > 0.5:", np.round((rel > 0.5).mean(), 3), "| n =", len(rel))
print("Main reading, preparatory tuning SD (spikes/s): median", np.median(M.tuning_strength(d0["prep"], "sd")).round(2),
"| max", M.tuning_strength(d0["prep"], "sd").max().round(2))
onset proto combine tuning prep_pass both_pass
a mean separate range 85 16
a mean separate sd 5 0
a mean concat range 76 17
a mean concat sd 5 0
a medoid separate range 86 15
a medoid separate sd 8 0
a medoid concat range 77 15
a medoid concat sd 6 0
b mean separate range 81 20
b mean separate sd 5 0
b mean concat range 74 18
b mean concat sd 6 0
b medoid separate range 89 15
b medoid separate sd 7 0
b medoid concat range 77 17
b medoid concat sd 6 0
Main reading, preparatory tuning range (spikes/s): median 4.93 | max 34.32
Split-half reliability of the preparatory rate vector (odd vs even trials), neurons over 5 spikes/s: median r 0.374 | Spearman-Brown corrected 0.544 | fraction with r > 0.5: 0.382 | n = 76 Main reading, preparatory tuning SD (spikes/s): median 1.19 | max 9.8
Noise reference for the 5 spikes/s threshold (added in response to the verifier, round 1). The paper accepts the lower array threshold "given the higher signal-to-noise ratio provided by the higher trial count" (Methods, Neural recordings and datasets); here there are 6 to 10 trials per condition. The same range measure on spike counts in the preparatory window, with condition labels shuffled among the kept trials (200 draws), shows how many units pass from noise alone, and a per-unit permutation p-value marks the units whose tuning is distinguishable from noise.
rng_thr = np.random.default_rng(5)
c_sel0 = cond[sel0]
cnt_prep = go0[sel0]
def range_counts(cs):
return M.tuning_strength(np.stack([cnt_prep[cs == c].mean(0) for c in conds]), "range")
rng_obs = range_counts(c_sel0)
rng_null = np.array([range_counts(rng_thr.permutation(c_sel0)) for _ in range(200)])
p_tune = ((rng_null >= rng_obs).sum(0) + 1) / (len(rng_null) + 1)
n_null = (rng_null > 5).sum(1)
print("Units with a preparatory range above 5 spikes/s (window counts, main selection): observed", int((rng_obs > 5).sum()),
"of", len(units), "| labels shuffled: median", int(np.median(n_null)), "range", n_null.min(), "to", n_null.max())
tuned_perm = pdn & (p_tune < 0.05)
print("PD neurons (range above 5 spikes/s):", int(pdn.sum()), "| of them with permutation p < 0.05:", int(tuned_perm.sum()),
"| all units with p < 0.05:", int((p_tune < 0.05).sum()), "(expected by chance about", round(0.05 * len(units), 1), ")")
print("Split-half reliability, PD neurons with permutation p < 0.05: median r",
np.round(np.median(rel[p_tune[np.where(pdn)[0]] < 0.05]), 3), "| the others:", np.round(np.median(rel[p_tune[np.where(pdn)[0]] >= 0.05]), 3))
Units with a preparatory range above 5 spikes/s (window counts, main selection): observed 75 of 152 | labels shuffled: median 59 range 51 to 67 PD neurons (range above 5 spikes/s): 76 | of them with permutation p < 0.05: 43 | all units with p < 0.05: 56 (expected by chance about 7.6 ) Split-half reliability, PD neurons with permutation p < 0.05: median r 0.611 | the others: 0.141
PD code on synthetic data¶
A check of maze.pd_performance_pca: 27 conditions, a T-matrix of rank 3 plus small noise, and
preparatory rates linear in its first three PCs plus noise. Performance should be near 1 from
k = 3 under both readings of S.
rng_syn = np.random.default_rng(3)
Z = rng_syn.standard_normal((27, 3))
Tsyn = Z @ rng_syn.standard_normal((3, 500)) + 0.01 * rng_syn.standard_normal((27, 500))
xs = 10 + Z @ np.array([2.0, -1.0, 0.5]) + 0.05 * rng_syn.standard_normal(27)
for R in "ab":
print("S reading", R, np.round(M.pd_performance_pca(xs, Tsyn @ Tsyn.T, [1, 2, 3, 6], R == "b"), 3))
print("Fixed space with the true Z:", round(M.pd_performance_fixed(xs, Z), 3))
S reading a [0.997 0.998 1. 1. ] S reading b [0.99 0.997 1. 1. ] Fixed space with the true Z: 1.0
j-array-recordings¶
The J-array dataset holds 146 single- and multi-unit isolations recorded with two 96-electrode arrays during the maze task, over 108 conditions, with a mean of 2166 trials per neuron.
n_units = len(units)
print("Units:", n_units, "(expected 146 +/- 14.6) -> within tolerance:", abs(n_units - 146) <= 14.6)
print("Units per array (id prefix): 1 ->", (units.id_array == 1).sum(), "| 2 ->", (units.id_array == 2).sum(),
"| electrodes per array in the table:", el.group_name.value_counts().to_dict())
print("Conditions:", len(conds), "(expected 108) | mazes", tr.maze_id.nunique(), "x versions", tr.trial_version.nunique())
print("maze_id values:", sorted(tr.maze_id.unique()))
print("Trials per neuron: train+val file", len(tr), "; with the 100 test segments", len(tr) + len(TE["trials"]),
"(paper: 2166 in the Results, 2155 in the Methods; reported only)")
print("Session description in the file:", D["nwb"].session_description)
gap = tr.start_time.values[1:] - tr.stop_time.values[:-1]
print("Trial timeline: first trial starts at", tr.start_time.min(), "s; next trial start minus trial end ranges from",
gap.min().round(3), "to", gap.max().round(3), "s, leaving no room for the failed or other trials a session holds, so the timeline was rebuilt",
"and cannot show which other trials or conditions the session held.")
Units: 152 (expected 146 +/- 14.6) -> within tolerance: True
Units per array (id prefix): 1 -> 75 | 2 -> 77 | electrodes per array in the table: {'electrode_group_PMd': 96, 'electrode_group_M1': 96}
Conditions: 27 (expected 108) | mazes 9 x versions 3
maze_id values: [np.int64(12), np.int64(13), np.int64(14), np.int64(20), np.int64(21), np.int64(22), np.int64(83), np.int64(84), np.int64(85)]
Trials per neuron: train+val file 250 ; with the 100 test segments 350 (paper: 2166 in the Results, 2155 in the Methods; reported only)
Session description in the file: Data from monkey Jenkins performing center-out delayed reaching task. This file contains trials from the full session on 2009-09-29 that can be used for training models for the Neural Latents Benchmark.
Trial timeline: first trial starts at 0.0 s; next trial start minus trial end ranges from 0.014 to 0.109 s, leaving no room for the failed or other trials a session holds, so the timeline was rebuilt and cannot show which other trials or conditions the session held.
Result: partially reproduced, resting on the 27 conditions being among the paper's 108. The dandiset holds 152 units (within 146 +/- 14.6) from two arrays, 75 and 77 by unit id, but 27 conditions (9 mazes in 3 versions), not 108, and 250 labelled trials per unit against the paper's 2166. The test file adds 100 unlabelled segments for the 114 held-in units only; the 38 held-out units have the 250 trials alone. The criterion's partial outcome applies when the units match and the conditions are a subset of the 108; that these 27 are among them is assumed, as nothing in the files shows it. The file describes its trials as taken from the session of 2009-09-29, and its rebuilt trial timeline cannot show whether that session also held other mazes. The paper's J-array dataset interleaved all 108 conditions, so if this session ran only these 27 it is a different recording and the match of the unit count says nothing about the paper's dataset; nothing in the files settles which.
fig4d-mean-correlation-j-array¶
For the J-array dataset, preparatory and peri-movement tuning are positively but only weakly correlated on average across neurons (mean r = 0.13, p = 0.03; mean Spearman r = 0.11).
Every reading: trial selection (8), tuning measure (2), test of the mean (3).
def corr_set(d, mask):
idx = np.where(mask)[0]
r = np.array([np.corrcoef(d["prep"][:, i], d["peri"][:, i])[0, 1] for i in idx])
rs = np.array([stats.spearmanr(d["prep"][:, i], d["peri"][:, i])[0] for i in idx])
return idx, r, rs
def mean_tests(r):
if len(r) < 3:
return dict(p_t=np.nan, p_wilcoxon=np.nan, p_fisher_t=np.nan)
return dict(p_t=stats.ttest_1samp(r, 0).pvalue, p_wilcoxon=stats.wilcoxon(r).pvalue,
p_fisher_t=stats.ttest_1samp(np.arctanh(np.clip(r, -0.999999, 0.999999)), 0).pvalue)
rows = []
for r_ in SEL_READINGS:
for how in ["range", "sd"]:
d = condition_data(*r_)
_, both = neuron_sets(d, how)
idx, r, rs = corr_set(d, both)
t = mean_tests(r)
row = dict(onset=r_[0], proto=r_[1], combine=r_[2], tuning=how, n=len(r),
mean_r=np.round(r.mean(), 3) if len(r) else np.nan,
mean_spearman=np.round(rs.mean(), 3) if len(rs) else np.nan,
**{k: np.round(v, 4) for k, v in t.items()})
pmin = np.nanmin(list(t.values())) if len(r) >= 3 else np.nan
row["within_tol_both"] = bool(len(r) and abs(row["mean_r"] - 0.13) <= 0.013 and abs(row["mean_spearman"] - 0.11) <= 0.011)
row["any_test_p<0.05"] = bool(pmin < 0.05) if len(r) >= 3 else False
row["in_0.054_0.231"] = bool(len(r) and 0.054 <= row["mean_r"] <= 0.231)
rows.append(row)
fig4d = pd.DataFrame(rows)
print(fig4d.to_string(index=False))
print()
print("Readings with both means in tolerance and p < 0.05:", int((fig4d.within_tol_both & fig4d["any_test_p<0.05"]).sum()))
print("Readings with mean r in 0.054..0.231 and p < 0.05:", int((fig4d["in_0.054_0.231"] & fig4d["any_test_p<0.05"]).sum()),
"of", int((fig4d.n >= 3).sum()), "readings with at least 3 neurons")
main_idx, main_r, main_rs = corr_set(d0, neuron_sets(d0, "range")[1])
print("Main reading", MAIN, "range: n =", len(main_r), "| mean r", main_r.mean().round(3),
"| SD", main_r.std(ddof=1).round(3), "| mean Spearman", main_rs.mean().round(3))
print("Per-array breakdown (main reading):", {a: (int((units.id_array.values[main_idx] == a).sum()),
np.round(main_r[units.id_array.values[main_idx] == a].mean(), 3)) for a in (1, 2)})
# threshold dependence (main selection, range measure)
rows = []
tpm = M.tuning_strength(d0["prep"], "range"); tqm = M.tuning_strength(d0["peri"], "range")
for thr in [3, 4, 5, 6, 8, 10]:
for ratio in [1.0, 1.5, 2.0]:
m = (tpm > thr) & (tqm >= ratio * tpm)
_, r, rs = corr_set(d0, m)
rows.append(dict(prep_thr=thr, peri_ratio=ratio, n=len(r), mean_r=np.round(r.mean(), 3) if len(r) else np.nan,
mean_spearman=np.round(rs.mean(), 3) if len(r) else np.nan,
p_t=np.round(stats.ttest_1samp(r, 0).pvalue, 4) if len(r) > 2 else np.nan))
print("Threshold dependence (main reading, range):")
print(pd.DataFrame(rows).to_string(index=False))
onset proto combine tuning n mean_r mean_spearman p_t p_wilcoxon p_fisher_t within_tol_both any_test_p<0.05 in_0.054_0.231
a mean separate range 16 0.163 0.161 0.0515 0.0833 0.0523 False False True
a mean separate sd 0 NaN NaN NaN NaN NaN False False False
a mean concat range 17 0.172 0.139 0.0422 0.0714 0.0443 False True True
a mean concat sd 0 NaN NaN NaN NaN NaN False False False
a medoid separate range 15 0.160 0.140 0.0489 0.0302 0.0501 False True True
a medoid separate sd 0 NaN NaN NaN NaN NaN False False False
a medoid concat range 15 0.170 0.119 0.0564 0.0637 0.0594 False False True
a medoid concat sd 0 NaN NaN NaN NaN NaN False False False
b mean separate range 20 0.163 0.157 0.0203 0.0296 0.0215 False True True
b mean separate sd 0 NaN NaN NaN NaN NaN False False False
b mean concat range 18 0.207 0.187 0.0130 0.0182 0.0147 False True True
b mean concat sd 0 NaN NaN NaN NaN NaN False False False
b medoid separate range 15 0.141 0.119 0.0826 0.0833 0.0798 True False True
b medoid separate sd 0 NaN NaN NaN NaN NaN False False False
b medoid concat range 17 0.183 0.159 0.0273 0.0348 0.0305 False True True
b medoid concat sd 0 NaN NaN NaN NaN NaN False False False
Readings with both means in tolerance and p < 0.05: 0
Readings with mean r in 0.054..0.231 and p < 0.05: 5 of 8 readings with at least 3 neurons
Main reading ('a', 'mean', 'concat') range: n = 17 | mean r 0.172 | SD 0.321 | mean Spearman 0.139
Per-array breakdown (main reading): {1: (9, np.float64(0.221)), 2: (8, np.float64(0.117))}
Threshold dependence (main reading, range):
prep_thr peri_ratio n mean_r mean_spearman p_t
3 1.0 54 0.143 0.118 0.0008
3 1.5 30 0.114 0.090 0.0405
3 2.0 16 -0.022 -0.020 0.7565
4 1.0 44 0.171 0.152 0.0007
4 1.5 23 0.140 0.126 0.0457
4 2.0 12 -0.010 0.015 0.9175
5 1.0 34 0.173 0.147 0.0029
5 1.5 17 0.172 0.139 0.0422
5 2.0 7 -0.006 -0.009 0.9621
6 1.0 25 0.177 0.162 0.0153
6 1.5 11 0.165 0.150 0.1639
6 2.0 4 -0.084 -0.114 0.7226
8 1.0 18 0.203 0.193 0.0236
8 1.5 6 0.210 0.192 0.2293
8 2.0 2 -0.133 -0.142 NaN
10 1.0 7 0.111 0.089 0.5579
10 1.5 3 0.135 0.089 0.7182
10 2.0 2 -0.133 -0.142 NaN
Result: partially reproduced, resting on the reading of the trial selection and of the tuning measure. Under the main reading 17 neurons pass both thresholds; the mean Pearson r is 0.172 and the mean Spearman r 0.139, both outside tolerance (0.13 +/- 0.013, 0.11 +/- 0.011), and the mean is above zero by the t-test (p = 0.0422) and the Fisher-z t-test (p = 0.0443), not by the Wilcoxon test (p = 0.0714). No reading puts both means in tolerance with p < 0.05. The mean r lies in 0.054 to 0.231 in all 8 selection readings, with p < 0.05 by some test in 5 of them, so the outcome depends on the selection reading. The SD measure of tuning leaves no neuron passing both thresholds. With these few trials per condition the counts move with the thresholds (threshold table).
fig4-broad-unimodal¶
Individual neurons span a broad range of preparatory to peri-movement correlations centred near zero, with little sign of separate populations of concordant (r near 1) and opposing (r near -1) tuning, except perhaps weakly for monkey B (Figure 4B).
rows = []
for r_ in SEL_READINGS:
d = condition_data(*r_)
_, r, _ = corr_set(d, neuron_sets(d, "range")[1])
if len(r) < 4:
continue
dip, p_dip = diptest.diptest(r)
rows.append(dict(onset=r_[0], proto=r_[1], combine=r_[2], n=len(r), min_r=r.min().round(3), max_r=r.max().round(3),
mean_r=r.mean().round(3), sd_r=r.std(ddof=1).round(3), dip=round(dip, 4), p_dip=round(p_dip, 3),
n_abs_gt_0_8=int((np.abs(r) > 0.8).sum()), n_abs_lt_0_5=int((np.abs(r) < 0.5).sum())))
f4u = pd.DataFrame(rows)
print(f4u.to_string(index=False))
print("SD-based tuning leaves too few neurons in every selection reading (see the neuron table).")
fig, ax = plt.subplots(figsize=(4, 3))
ax.hist(main_r, bins=np.linspace(-1, 1, 21), color="0.5")
ax.axvline(main_r.mean(), color="r")
ax.set_xlabel("preparatory vs peri-movement r"); ax.set_ylabel("neurons"); ax.set_title("main reading")
plt.show()
print("Histogram range: all", len(main_r), "values inside -1..1")
onset proto combine n min_r max_r mean_r sd_r dip p_dip n_abs_gt_0_8 n_abs_lt_0_5
a mean separate 16 -0.408 0.660 0.163 0.309 0.0563 0.972 0 14
a mean concat 17 -0.453 0.673 0.172 0.321 0.0665 0.788 0 15
a medoid separate 15 -0.397 0.674 0.160 0.288 0.0925 0.312 0 14
a medoid concat 15 -0.461 0.684 0.170 0.317 0.0597 0.960 0 13
b mean separate 20 -0.397 0.668 0.163 0.288 0.0888 0.195 0 18
b mean concat 18 -0.473 0.677 0.207 0.317 0.0617 0.855 0 16
b medoid separate 15 -0.359 0.537 0.141 0.291 0.0807 0.566 0 14
b medoid concat 17 -0.498 0.672 0.183 0.310 0.0457 0.993 0 15
SD-based tuning leaves too few neurons in every selection reading (see the neuron table).
Histogram range: all 17 values inside -1..1
Result: reproduced. In every selection reading the correlations run from negative (-0.36 to -0.50) to positive (0.54 to 0.68) values, the mean is 0.141 to 0.207, no neuron has |r| > 0.8, and the dip test does not reject unimodality (p = 0.195 to 0.993). Under the main reading n = 17, mean 0.172, SD 0.321, dip p = 0.788. The sample is small.
fig4e-correlation-time-course¶
The correlation between preparatory tuning and peri-movement activity at a single time point begins to decline about 200 ms before movement onset and is near zero by movement onset.
def corr_trace(prep, rate):
"""prep (c, n), rate (c, n, T) -> per-neuron correlation at each time (n, T)."""
a = prep - prep.mean(0)
b = rate - rate.mean(0)
num = (a[:, :, None] * b).sum(0)
den = np.sqrt((a ** 2).sum(0))[:, None] * np.sqrt((b ** 2).sum(0))
return num / den
rt_med = int(np.median(tr.rt[d0["sel"]]))
def fig4e_traces(d, mask):
idx = np.where(mask)[0]
r0 = M.MOVE_WIN[0]
mv_t = np.arange(-400, 101)
mv_c = corr_trace(d["prep"][:, idx], d["mv_rate"][:, idx, mv_t[0] - r0: mv_t[-1] + 1 - r0])
g0 = M.GO_WIN[0]
go_t = np.arange(-200, 101)
go_c = corr_trace(d["prep"][:, idx], d["go_rate"][:, idx, go_t[0] - g0: go_t[-1] + 1 - g0])
# reading (a): go-locked part mapped to movement time with the median RT, then movement-locked
go_tm = go_t - rt_med
keep_mv = mv_t > go_tm[-1]
sp_t = np.concatenate([go_tm, mv_t[keep_mv]])
sp_c = np.concatenate([go_c, mv_c[:, keep_mv]], axis=1)
return dict(b=(mv_t, mv_c), a=(sp_t, sp_c))
def judge_4e(t, c):
m = np.nanmean(c, 0)
def at(x):
return m[np.argmin(np.abs(t - x))]
v400 = at(-400)
pre = (t <= 0) & (t >= -400)
within = pre & (np.abs(m - v400) <= 0.1 * abs(v400))
decline = t[within].max() if within.any() else np.nan
on = c[:, np.argmin(np.abs(t))]
on = on[np.isfinite(on)]
t_ = mean_tests(on)
return dict(at_m400=round(v400, 3), at_m200=round(at(-200), 3), at_0=round(at(0), 3),
decline_ms=decline, p_t=round(t_["p_t"], 4), p_wil=round(t_["p_wilcoxon"], 4),
p_fisher=round(t_["p_fisher_t"], 4), n=len(on))
rows = []
for r_ in SEL_READINGS:
d = condition_data(*r_)
prep_ok, both = neuron_sets(d, "range")
for nset, mask in [("a both thresholds", both), ("b preparatory only", prep_ok)]:
tr_ = fig4e_traces(d, mask)
for rng_name in ["a", "b"]:
t, c = tr_[rng_name]
rows.append(dict(onset=r_[0], proto=r_[1], combine=r_[2], neurons=nset, time_range=rng_name, **judge_4e(t, c)))
f4e = pd.DataFrame(rows)
f4e["drop_200_to_0"] = f4e.at_m200 > f4e.at_0
f4e["near_zero_a"] = ~((f4e.at_0 > 0) & (f4e[["p_t", "p_wil", "p_fisher"]].min(1) < 0.05))
f4e["near_zero_b"] = f4e.at_0.abs() < 0.1
f4e["decline_in_-300_-100"] = f4e.decline_ms.between(-300, -100)
print(f4e.to_string(index=False))
ok4e = f4e.drop_200_to_0 & (f4e.near_zero_a | f4e.near_zero_b) & f4e["decline_in_-300_-100"]
print("Readings meeting the criterion:", int(ok4e.sum()), "of", len(f4e))
print("Median RT used for the splice (ms):", rt_med)
print("Go-locked part of the splice: go cue -200 to +100 ms, the preparatory window", M.PREP,
"| mapped to", -200 - rt_med, "to", 100 - rt_med, "ms from movement onset")
print("Readings meeting the criterion, of 8 selection readings, by neurons and time range:")
print(f4e.assign(meets=ok4e).groupby(["neurons", "time_range"]).meets.sum().to_string())
tr_main = fig4e_traces(d0, neuron_sets(d0, "range")[1])
fig, ax = plt.subplots(1, 2, figsize=(9, 3))
for a_, (nm, (t, c)) in zip(ax, tr_main.items()):
a_.plot(t, np.nanmean(c, 0), "k"); a_.axhline(0, color="0.7"); a_.axvline(0, color="0.7")
a_.set_title(f"time range ({nm})"); a_.set_xlabel("ms from movement onset"); a_.set_ylabel("mean r")
plt.show()
print("Plotted mean r range:", round(min(np.nanmean(c, 0).min() for t, c in tr_main.values()), 3), "to",
round(max(np.nanmean(c, 0).max() for t, c in tr_main.values()), 3))
onset proto combine neurons time_range at_m400 at_m200 at_0 decline_ms p_t p_wil p_fisher n drop_200_to_0 near_zero_a near_zero_b decline_in_-300_-100
a mean separate a both thresholds a 0.634 0.449 0.050 -266 0.4944 0.4954 0.4865 16 True True True True
a mean separate a both thresholds b 0.663 0.449 0.050 -312 0.4944 0.4954 0.4865 16 True True True False
a mean separate b preparatory only a 0.675 0.393 0.067 -278 0.0368 0.0808 0.0298 85 True False True True
a mean separate b preparatory only b 0.659 0.393 0.067 -299 0.0368 0.0808 0.0298 85 True False True True
a mean concat a both thresholds a 0.703 0.492 0.118 -263 0.1267 0.1202 0.1294 17 True True False True
a mean concat a both thresholds b 0.677 0.492 0.118 -307 0.1267 0.1202 0.1294 17 True True False False
a mean concat b preparatory only a 0.720 0.452 0.107 -278 0.0036 0.0071 0.0033 76 True False False True
a mean concat b preparatory only b 0.690 0.452 0.107 -279 0.0036 0.0071 0.0033 76 True False False True
a medoid separate a both thresholds a 0.674 0.348 0.017 -290 0.8135 0.8040 0.8096 15 True True True True
a medoid separate a both thresholds b 0.661 0.348 0.017 -300 0.8135 0.8040 0.8096 15 True True True True
a medoid separate b preparatory only a 0.683 0.352 0.054 -288 0.0856 0.1598 0.0668 86 True True True True
a medoid separate b preparatory only b 0.651 0.352 0.054 -298 0.0856 0.1598 0.0668 86 True True True True
a medoid concat a both thresholds a 0.701 0.518 0.151 -257 0.0702 0.0637 0.0751 15 True True False True
a medoid concat a both thresholds b 0.678 0.518 0.151 -300 0.0702 0.0637 0.0751 15 True True False True
a medoid concat b preparatory only a 0.714 0.443 0.108 -276 0.0027 0.0057 0.0026 77 True False False True
a medoid concat b preparatory only b 0.682 0.443 0.108 -281 0.0027 0.0057 0.0026 77 True False False True
b mean separate a both thresholds a 0.657 0.416 0.081 -269 0.2140 0.2305 0.2122 20 True True True True
b mean separate a both thresholds b 0.617 0.416 0.081 -298 0.2140 0.2305 0.2122 20 True True True True
b mean separate b preparatory only a 0.682 0.369 0.073 -275 0.0246 0.0582 0.0209 80 True False True True
b mean separate b preparatory only b 0.638 0.369 0.073 -287 0.0246 0.0582 0.0209 80 True False True True
b mean concat a both thresholds a 0.710 0.475 0.129 -277 0.0886 0.0814 0.0933 18 True True False True
b mean concat a both thresholds b 0.676 0.475 0.129 -325 0.0886 0.0814 0.0933 18 True True False False
b mean concat b preparatory only a 0.719 0.437 0.098 -276 0.0076 0.0137 0.0069 73 True False True True
b mean concat b preparatory only b 0.687 0.437 0.098 -290 0.0076 0.0137 0.0069 73 True False True True
b medoid separate a both thresholds a 0.609 0.362 0.114 -275 0.1296 0.1514 0.1254 15 True True False True
b medoid separate a both thresholds b 0.573 0.362 0.114 -302 0.1296 0.1514 0.1254 15 True True False False
b medoid separate b preparatory only a 0.660 0.346 0.070 -276 0.0187 0.0530 0.0151 88 True False True True
b medoid separate b preparatory only b 0.620 0.346 0.070 -283 0.0187 0.0530 0.0151 88 True False True True
b medoid concat a both thresholds a 0.701 0.494 0.152 -264 0.0395 0.0569 0.0441 17 True False False True
b medoid concat a both thresholds b 0.658 0.494 0.152 -315 0.0395 0.0569 0.0441 17 True False False False
b medoid concat b preparatory only a 0.715 0.416 0.098 -277 0.0059 0.0119 0.0055 76 True False True True
b medoid concat b preparatory only b 0.674 0.416 0.098 -285 0.0059 0.0119 0.0055 76 True False True True
Readings meeting the criterion: 22 of 32
Median RT used for the splice (ms): 354
Go-locked part of the splice: go cue -200 to +100 ms, the preparatory window (-200, 100) | mapped to -554 to -254 ms from movement onset
Readings meeting the criterion, of 8 selection readings, by neurons and time range:
neurons time_range
a both thresholds a 7
b 3
b preparatory only a 6
b 6
/tmp/ipykernel_2535/2232206344.py:7: RuntimeWarning: invalid value encountered in divide return num / den /tmp/ipykernel_2535/2232206344.py:51: Pandas4Warning: Starting with pandas version 4.0 all arguments of min will be keyword-only. f4e["near_zero_a"] = ~((f4e.at_0 > 0) & (f4e[["p_t", "p_wil", "p_fisher"]].min(1) < 0.05))
Plotted mean r range: 0.103 to 0.705
Result: reproduced, resting on the reading of the neurons and the time range; in the main selection it rests on the splice. With the neurons of Figure 4D and the spliced trace (readings a, a: the legend places 4E with 4A to 4D, and the Supplemental Experimental Procedures name the spliced trace for 4E), the main reading gives a mean r of 0.703 at 400 ms before onset, 0.492 at 200 ms before and 0.118 at onset, not above zero (t-test p = 0.1267), and the decline begins at -263 ms. These two readings meet the criterion in 7 of 8 selection readings.
The splice decides this result. Its go-locked part, go cue -200 to +100 ms, is the preparatory window itself, mapped to movement time with the median RT (354 ms, an assumption: the paper does not say where it splices) as -554 to -254 ms. The reference value at 400 ms before onset is therefore taken inside the window the preparatory rate is averaged over, and the decline at -263 ms is the end of that window shifted by the median RT, not a time read from the data. In the main selection no other combination meets the criterion: with the movement-locked trace alone the decline begins at -307 ms, just outside -300 to -100, and with all preparatory neurons the onset value (0.107) is above zero (p = 0.0036) and above 0.1. Without the splice the criterion is met by the movement-locked trace with all preparatory neurons in 6 of 8 selection readings (decline at -283 to -299 ms, onset value 0.054 to 0.098) and with the neurons of Figure 4D in 3 of 8 (decline at -298 to -300 ms, at the edge of the window). So the outcome holds under the criterion through those readings, with the decline at the early edge of -300 to -100 ms. 22 of 32 combinations meet the criterion.
suppfig1-population-psth-weak-movement-tuning¶
In population averages sorted by each neuron's preparatory preference, the preparatory differences carry weakly into the peri-movement epoch and show little relationship with preparatory tuning by movement onset.
def supp1(d, mask):
idx = np.where(mask)[0]
nc = len(conds)
ranks = [nc - 1, 2 * (nc - 1) // 3, (nc - 1) // 3, 0] # strongest, two intermediate, weakest
go_tr, mv_tr = [], []
for n in idx:
o = np.argsort(d["prep"][:, n])
mx = max(d["go_rate"][:, n].max(), d["mv_rate"][:, n].max())
go_tr.append(d["go_rate"][o[ranks], n] / mx)
mv_tr.append(d["mv_rate"][o[ranks], n] / mx)
go_tr = np.mean(go_tr, 0); mv_tr = np.mean(mv_tr, 0)
g0, m0 = M.GO_WIN[0], M.MOVE_WIN[0]
prep_diff = (go_tr[0, M.PREP[0] - g0: M.PREP[1] - g0] - go_tr[3, M.PREP[0] - g0: M.PREP[1] - g0]).mean()
onset_diff = mv_tr[0, -m0] - mv_tr[3, -m0]
return go_tr, mv_tr, prep_diff, onset_diff
rows = []
for r_ in SEL_READINGS:
d = condition_data(*r_)
for nset, mask in zip(["both thresholds", "preparatory only (variant)"], neuron_sets(d, "range")[::-1]):
_, _, pdf, odf = supp1(d, mask)
rows.append(dict(onset=r_[0], proto=r_[1], combine=r_[2], neurons=nset, n=int(mask.sum()), prep_diff=round(pdf, 3),
onset_diff=round(odf, 3), ratio=round(odf / pdf, 3)))
s1 = pd.DataFrame(rows)
s1["a_smaller"] = s1.onset_diff < s1.prep_diff
s1["b_less_than_half"] = s1.ratio < 0.5
print(s1.to_string(index=False))
go_tr, mv_tr, _, _ = supp1(d0, neuron_sets(d0, "range")[1])
fig, ax = plt.subplots(1, 2, figsize=(9, 3), sharey=True)
cols = ["r", "salmon", "yellowgreen", "g"]
for i in range(4):
ax[0].plot(np.arange(*M.GO_WIN), go_tr[i], color=cols[i]); ax[1].plot(np.arange(*M.MOVE_WIN), mv_tr[i], color=cols[i])
ax[0].set_xlabel("ms from go cue"); ax[1].set_xlabel("ms from movement onset"); ax[0].set_ylabel("normalized rate")
plt.show()
print("Plotted normalized rate range:", round(min(go_tr.min(), mv_tr.min()), 3), "to", round(max(go_tr.max(), mv_tr.max()), 3))
# Cross-validated variant (verifier, round 1): conditions ranked on the odd trials of each condition and the traces
# measured on the even trials, and the reverse, averaged; removes the regression to the mean of ranking and measuring
# on the same trials. Main selection only.
go_b = M.bin_spikes(st, tr.go_cue_time.values, M.GO_WIN)
mv_b = M.bin_spikes(st, tr.move_onset_time.values, M.MOVE_WIN)
def half_means(h):
m = np.zeros(len(tr), bool)
for c in conds:
ii = np.where(d0["sel"] & (cond == c))[0]
m[ii[h::2]] = True
go_r = M.smooth_rate(M.cond_mean(go_b[m].astype(np.float32), cond[m], conds))
mv_r = M.smooth_rate(M.cond_mean(mv_b[m].astype(np.float32), cond[m], conds))
return dict(go_rate=go_r, mv_rate=mv_r, prep=M.window_mean(go_r, M.PREP, M.GO_WIN[0]))
HALF = [half_means(0), half_means(1)]
def supp1_cv(rank_d, meas_d, idx):
nc = len(conds)
ranks = [nc - 1, 2 * (nc - 1) // 3, (nc - 1) // 3, 0]
go_tr, mv_tr = [], []
for n in idx:
o = np.argsort(rank_d["prep"][:, n])
mx = max(meas_d["go_rate"][:, n].max(), meas_d["mv_rate"][:, n].max())
go_tr.append(meas_d["go_rate"][o[ranks], n] / mx); mv_tr.append(meas_d["mv_rate"][o[ranks], n] / mx)
go_tr = np.mean(go_tr, 0); mv_tr = np.mean(mv_tr, 0)
g0, m0 = M.GO_WIN[0], M.MOVE_WIN[0]
return (go_tr[0, M.PREP[0] - g0: M.PREP[1] - g0] - go_tr[3, M.PREP[0] - g0: M.PREP[1] - g0]).mean(), mv_tr[0, -m0] - mv_tr[3, -m0]
rows = []
for nset, mask in zip(["both thresholds", "preparatory only (variant)"], neuron_sets(d0, "range")[::-1]):
idx = np.where(mask)[0]
for how, pairs in [("cross-validated", [(0, 1), (1, 0)]), ("same half (biased)", [(0, 0), (1, 1)])]:
res = np.array([supp1_cv(HALF[a], HALF[b], idx) for a, b in pairs])
pdf, odf = res.mean(0)
rows.append(dict(neurons=nset, n=len(idx), ranking=how, prep_diff=round(pdf, 3), onset_diff=round(odf, 3),
ratio=round(odf / pdf, 3), a_smaller=odf < pdf, b_less_than_half=odf / pdf < 0.5))
print("Half-trial variants, main selection:")
print(pd.DataFrame(rows).to_string(index=False))
onset proto combine neurons n prep_diff onset_diff ratio a_smaller b_less_than_half
a mean separate both thresholds 16 0.098 0.027 0.272 True True
a mean separate preparatory only (variant) 85 0.130 0.012 0.094 True True
a mean concat both thresholds 17 0.108 0.101 0.942 True False
a mean concat preparatory only (variant) 76 0.140 0.030 0.213 True True
a medoid separate both thresholds 15 0.089 0.032 0.357 True True
a medoid separate preparatory only (variant) 86 0.122 0.018 0.151 True True
a medoid concat both thresholds 15 0.111 0.100 0.902 True False
a medoid concat preparatory only (variant) 77 0.136 0.023 0.170 True True
b mean separate both thresholds 20 0.092 0.017 0.184 True True
b mean separate preparatory only (variant) 81 0.129 0.006 0.047 True True
b mean concat both thresholds 18 0.097 0.076 0.776 True False
b mean concat preparatory only (variant) 74 0.140 0.020 0.141 True True
b medoid separate both thresholds 15 0.083 0.029 0.350 True True
b medoid separate preparatory only (variant) 89 0.125 0.020 0.162 True True
b medoid concat both thresholds 17 0.097 0.067 0.692 True False
b medoid concat preparatory only (variant) 77 0.140 0.015 0.109 True True
Plotted normalized rate range: 0.01 to 0.183
Half-trial variants, main selection:
neurons n ranking prep_diff onset_diff ratio a_smaller b_less_than_half
both thresholds 17 cross-validated 0.040 0.045 1.115 False False
both thresholds 17 same half (biased) 0.112 0.068 0.607 True False
preparatory only (variant) 76 cross-validated 0.055 0.030 0.545 True False
preparatory only (variant) 76 same half (biased) 0.129 0.033 0.253 True True
Result: reproduced under the claims file's reading (a), not under (b). Under the main reading (17 neurons) the strongest-minus-weakest difference of normalized rate is 0.108 in the preparatory window and 0.101 at movement onset, a ratio of 0.942: smaller, but not less than half. Under the separate-profile selections the ratio is 0.184 to 0.357, and with all preparatory neurons 0.047 to 0.213. That the difference shows "little relationship" by onset holds only in those variants. Ranking and measuring on the same trials favours the claim: with conditions ranked on half the trials and measured on the other half, the preparatory difference of the 17 neurons falls to 0.040 and the onset difference (0.045) is no longer smaller (ratio 1.115); with all preparatory neurons it is 0.055 against 0.030 (ratio 0.545). So under the legend's method the claim holds, but on these data the result for the claim's neurons does not survive cross-validation.
PD performance in each space (Figure 6)¶
The PD model, leave-one-condition-out performance, the peri-movement, kinematic and task spaces
and the upper limit, as in the claims file's opening comment and maze.py. Neurons fitted: all
with preparatory tuning above 5 spikes/s (range), regardless of peri-movement activity (Methods,
Correlations, last sentence). The peri-movement and upper-limit T-matrices use every unit but the
one fitted. Readings run: S fixed (a) or recomputed in each fold (b); kinematic window end (a) 800
ms or (b) the end of the slowest condition's reach; kinematic signals (a) repeated to one block
per unit or (b) once; and, for the task spaces, segments equal in time or path length and peak
velocity read two ways (each listed as its own space).
KS = list(range(1, 15))
def fig6_all(d, how="range"):
prep_ok, _ = neuron_sets(d, how)
fit = np.where(prep_ok)[0]
X = d["prep"]
Gu = M.grams_per_unit(d["mv_rate"], M.MOVE_WIN[0])
Gt = Gu.sum(0)
Gp = M.prep_grams(X)
kin = M.condition_kinematics(d["pos"], d["vel"], M.MOVE_WIN[0])
v = tr[d["sel"]].groupby(cond[d["sel"]]).trial_version.first().reindex(conds).values
spaces = M.task_spaces(kin, (v > 0).astype(float), (v == 2).astype(float))
res = dict(fit=fit, kin=kin, spaces=spaces, curves={}, task={})
for rec in [False, True]:
R = "b" if rec else "a"
res["curves"][("peri", R)] = M.perf_curves_pca(X, Gt, Gu, fit, KS, rec)
res["curves"][("upper", R)] = M.perf_curves_pca(X, Gp.sum(0), Gp, fit, KS, rec)
res["curves"][("upper incl. neuron", R)] = M.perf_curves_pca(X, Gp.sum(0), None, fit, KS, rec)
for kw in ["a", "b"]:
end = M.TMAT[1] if kw == "a" else int(kin.reach_end_ms.max())
for rep in ["a", "b"]:
Gk = M.kinematic_gram(d["pos"], d["vel"], M.MOVE_WIN[0], len(units) - 1, end_ms=end, replicate=(rep == "a"))
res["curves"][(f"kinematic win {kw} rep {rep}", R)] = M.perf_curves_pca(X, Gk, None, fit, KS, rec)
for nm, S in spaces.items():
res["task"][nm] = np.array([M.pd_performance_fixed(X[:, n], S) for n in fit])
return res
F6 = fig6_all(d0)
print("Neurons fitted:", len(F6["fit"]), "| by array:", pd.Series(units.id_array.values[F6["fit"]]).value_counts().sort_index().to_dict())
print("Reach end (ms after onset) over conditions:", F6["kin"].reach_end_ms.min(), "to", F6["kin"].reach_end_ms.max())
curve_tab = pd.DataFrame({f"{k[0]} | S {k[1]}": v.mean(0) for k, v in F6["curves"].items()}, index=KS).T.round(3)
print("Mean performance over neurons, by k:")
print(curve_tab.to_string())
task_tab = pd.Series({k: v.mean() for k, v in F6["task"].items()}).round(3)
dims = {k: v.shape[1] for k, v in F6["spaces"].items()}
print()
print(pd.DataFrame({"dim": dims, "mean_perf": task_tab}).to_string())
Neurons fitted: 76 | by array: {1: 35, 2: 41}
Reach end (ms after onset) over conditions: 251 to 619
Mean performance over neurons, by k:
1 2 3 4 5 6 7 8 9 10 11 12 13 14
peri | S a -0.006 0.045 0.052 0.091 0.100 0.117 0.102 0.061 0.026 -0.008 -0.184 -0.396 -0.574 -0.980
upper | S a -0.014 0.054 0.026 0.056 0.068 0.061 0.050 -0.028 -0.099 -0.252 -0.281 -0.394 -0.549 -0.838
upper incl. neuron | S a 0.009 0.096 0.102 0.148 0.165 0.176 0.178 0.122 0.103 0.048 0.029 -0.016 -0.064 -0.097
kinematic win a rep a | S a 0.008 0.016 0.013 0.016 -0.005 -0.040 -0.070 -0.128 -0.176 -0.216 -0.328 -0.480 -0.696 -2.005
kinematic win a rep b | S a 0.008 0.016 0.013 0.016 -0.005 -0.040 -0.071 -0.128 -0.176 -0.216 -0.329 -0.481 -0.694 -2.011
kinematic win b rep a | S a 0.011 0.022 0.018 0.017 -0.004 -0.036 -0.070 -0.128 -0.175 -0.214 -0.216 -0.430 -0.506 -2.317
kinematic win b rep b | S a 0.011 0.022 0.018 0.017 -0.004 -0.036 -0.070 -0.128 -0.175 -0.214 -0.217 -0.431 -0.506 -2.319
peri | S b -0.004 0.056 0.065 0.108 0.129 0.145 0.195 0.203 0.205 0.205 0.221 0.225 0.226 0.226
upper | S b -0.015 0.058 0.042 0.080 0.103 0.115 0.144 0.184 0.172 0.165 0.157 0.151 0.138 0.098
upper incl. neuron | S b 0.004 0.092 0.091 0.151 0.181 0.205 0.236 0.277 0.275 0.282 0.289 0.297 0.300 0.299
kinematic win a rep a | S b 0.007 0.017 0.015 0.019 -0.002 -0.020 -0.029 -0.046 -0.082 -0.142 -0.105 -0.146 -0.188 -0.196
kinematic win a rep b | S b 0.007 0.017 0.015 0.019 -0.002 -0.020 -0.029 -0.046 -0.082 -0.142 -0.105 -0.146 -0.188 -0.196
kinematic win b rep a | S b 0.010 0.024 0.020 0.020 -0.001 -0.013 -0.028 -0.044 -0.085 -0.165 -0.093 -0.148 -0.181 -0.187
kinematic win b rep b | S b 0.009 0.024 0.020 0.020 -0.001 -0.013 -0.028 -0.044 -0.085 -0.165 -0.093 -0.148 -0.181 -0.187
dim mean_perf
endpoints + max speed 3 -0.045
initial velocity 2 0.030
endpoints + orth. excursion 3 0.001
endpoints + viapoint 4 -0.026
endpoints + max velocity (pk) 4 0.007
endpoints + max velocity (ex) 4 -0.005
endpoints + halfpoint 4 -0.012
segment dirs (t) 6 0.008
kitchen sink (t) 11 -0.244
segment dirs (l) 6 -0.027
kitchen sink (l) 11 -0.221
fig6d-perimovement-best-j-array¶
For the J-array dataset, the highest cross-validated PD performance over all tested spaces and dimensionalities is reached when the PD is expressed in the peri-movement space.
def fig6d_judge(F):
rows = []
best_task = max(v.mean() for v in F["task"].values())
best_task_name = max(F["task"], key=lambda k: F["task"][k].mean())
for R in "ab":
peri = F["curves"][("peri", R)].mean(0)
for kw in "ab":
for rep in "ab":
kin = F["curves"][(f"kinematic win {kw} rep {rep}", R)].mean(0)
rows.append(dict(S=R, kin_window=kw, kin_rep=rep, peri_best=peri.max().round(3), peri_best_k=KS[peri.argmax()],
kin_best=kin.max().round(3), kin_best_k=KS[kin.argmax()], task_best=round(best_task, 3),
task_best_space=best_task_name, upper_best=F["curves"][("upper", R)].mean(0).max().round(3),
upper_incl_best=F["curves"][("upper incl. neuron", R)].mean(0).max().round(3),
peri_highest=bool(peri.max() > max(kin.max(), best_task))))
return pd.DataFrame(rows)
f6d = fig6d_judge(F6)
print(f6d.to_string(index=False))
fig, ax = plt.subplots(1, 2, figsize=(10, 3.5), sharey=True)
for a_, R in zip(ax, "ab"):
for nm, colr in [("peri", "k"), ("upper", "0.6"), ("kinematic win a rep a", "b")]:
a_.plot(KS, F6["curves"][(nm, R)].mean(0), color=colr, label=nm)
for nm, v in F6["task"].items():
a_.plot(dims[nm], v.mean(), "g.", ms=8)
a_.axhline(0, color="0.8"); a_.set_title(f"S reading {R}"); a_.set_xlabel("dimensionality")
ax[0].set_ylabel("mean performance"); ax[0].legend(fontsize=7)
ax[0].set_ylim(-0.5, 1)
plt.show()
allv = np.concatenate([np.concatenate([v.mean(0) for v in F6["curves"].values()]), task_tab.values])
print("Points below the axis floor (-0.5):", int((allv < -0.5).sum()), "of", len(allv), "| min", allv.min().round(3))
S kin_window kin_rep peri_best peri_best_k kin_best kin_best_k task_best task_best_space upper_best upper_incl_best peri_highest a a a 0.117 6 0.016 2 0.03 initial velocity 0.068 0.178 True a a b 0.117 6 0.016 2 0.03 initial velocity 0.068 0.178 True a b a 0.117 6 0.022 2 0.03 initial velocity 0.068 0.178 True a b b 0.117 6 0.022 2 0.03 initial velocity 0.068 0.178 True b a a 0.226 13 0.019 4 0.03 initial velocity 0.184 0.300 True b a b 0.226 13 0.019 4 0.03 initial velocity 0.184 0.300 True b b a 0.226 13 0.024 2 0.03 initial velocity 0.184 0.300 True b b b 0.226 13 0.024 2 0.03 initial velocity 0.184 0.300 True
Points below the axis floor (-0.5): 12 of 207 | min -2.319
Robustness over the trial-selection readings (main tuning measure). The SD measure leaves too few neurons to fit.
F6_all = {MAIN: F6}
rows = []
for r_ in SEL_READINGS:
if r_ not in F6_all:
F6_all[r_] = fig6_all(condition_data(*r_))
t = fig6d_judge(F6_all[r_])
rows.append(dict(onset=r_[0], proto=r_[1], combine=r_[2], neurons=len(F6_all[r_]["fit"]),
peri_best_min=t.peri_best.min(), kin_best_max=t.kin_best.max(), task_best=t.task_best.iloc[0],
peri_highest_all=bool(t.peri_highest.all()), peri_highest_any=bool(t.peri_highest.any())))
f6d_sel = pd.DataFrame(rows)
print(f6d_sel.to_string(index=False))
by_arr = {a: F6["curves"][("peri", "a")][units.id_array.values[F6["fit"]] == a].mean(0).max().round(3) for a in (1, 2)}
print("Best peri-movement performance (S a) among neurons of array 1 and 2:", by_arr)
onset proto combine neurons peri_best_min kin_best_max task_best peri_highest_all peri_highest_any
a mean separate 85 0.000 -0.019 -0.015 True True
a mean concat 76 0.117 0.024 0.030 True True
a medoid separate 86 -0.042 -0.028 -0.028 False True
a medoid concat 77 0.076 0.019 0.027 True True
b mean separate 81 0.038 -0.003 -0.003 True True
b mean concat 74 0.120 0.033 0.040 True True
b medoid separate 89 -0.024 -0.025 -0.024 True True
b medoid concat 77 0.091 0.021 0.028 True True
Best peri-movement performance (S a) among neurons of array 1 and 2: {1: np.float64(0.291), 2: np.float64(-0.006)}
Chance level of each space (added in response to the verifier, round 1). The Results say "Points above zero indicate performance considerably better than chance". Each fitted neuron's preparatory rates are shuffled over conditions (10 draws) and the performance of every space is recomputed. Then the comparison is repeated on the PD neurons whose preparatory tuning is distinguishable from noise (permutation p < 0.05, above).
rng6 = np.random.default_rng(60)
X0, fit0 = d0["prep"], F6["fit"]
Gu0 = M.grams_per_unit(d0["mv_rate"], M.MOVE_WIN[0]); Gt0 = Gu0.sum(0)
Gk0 = M.kinematic_gram(d0["pos"], d0["vel"], M.MOVE_WIN[0], len(units) - 1)
NULL = {}
for draw in range(10):
Xs = X0.copy()
for n in fit0:
Xs[:, n] = X0[rng6.permutation(len(conds)), n]
for R in "ab":
NULL.setdefault(("peri", R), []).append(M.perf_curves_pca(Xs, Gt0, Gu0, fit0, KS, R == "b").mean(0))
NULL.setdefault(("kinematic win a rep a", R), []).append(M.perf_curves_pca(Xs, Gk0, None, fit0, KS, R == "b").mean(0))
for nm, S in F6["spaces"].items():
NULL.setdefault(nm, []).append(np.mean([M.pd_performance_fixed(Xs[:, n], S) for n in fit0]))
NULL = {k: np.array(v) for k, v in NULL.items()}
print("Chance: mean performance with conditions shuffled (10 draws), by k:")
print(pd.DataFrame({f"{k[0]} | S {k[1]}": v.mean(0) for k, v in NULL.items() if isinstance(k, tuple)}, index=KS).T.round(3).to_string())
print("Largest SD over draws of a curve point:", round(max(v.std(0).max() for k, v in NULL.items() if isinstance(k, tuple)), 3))
null_task = pd.Series({k: v.mean() for k, v in NULL.items() if not isinstance(k, tuple)})
print(pd.DataFrame({"dim": dims, "observed": task_tab, "chance": null_task.round(3)}).to_string())
for key in [("peri", "a"), ("peri", "b"), ("kinematic win a rep a", "a"), ("kinematic win a rep a", "b")]:
obs = F6["curves"][key].mean(0); ch = NULL[key].mean(0)
print(key, "| observed above chance at every k:", bool((obs > ch).all()), "| observed - chance, k = 1..14:", np.round(obs - ch, 3))
print("Task spaces above chance:", int((task_tab > null_task).sum()), "of", len(task_tab))
# variant: PD neurons whose tuning passes the permutation test
fitp = np.where(tuned_perm)[0]
rows = []
for R in "ab":
peri = M.perf_curves_pca(X0, Gt0, Gu0, fitp, KS, R == "b").mean(0)
kin = M.perf_curves_pca(X0, Gk0, None, fitp, KS, R == "b").mean(0)
task = {nm: np.mean([M.pd_performance_fixed(X0[:, n], S) for n in fitp]) for nm, S in F6["spaces"].items()}
bt = max(task, key=task.get)
rows.append(dict(S=R, neurons=len(fitp), peri_best=peri.max().round(3), peri_best_k=KS[peri.argmax()], peri_k14=peri[-1].round(3),
kin_best=kin.max().round(3), kin_best_k=KS[kin.argmax()], task_best=round(task[bt], 3), task_best_space=bt,
seg_dirs_t=round(task["segment dirs (t)"], 3), init_vel=round(task["initial velocity"], 3),
endpoints_max_speed=round(task["endpoints + max speed"], 3), peri_highest=bool(peri.max() > max(kin.max(), task[bt]))))
print("Variant, PD neurons with permutation p < 0.05 (kinematic window a, replication a):")
print(pd.DataFrame(rows).to_string(index=False))
Chance: mean performance with conditions shuffled (10 draws), by k:
1 2 3 4 5 6 7 8 9 10 11 12 13 14
peri | S a -0.120 -0.171 -0.222 -0.284 -0.363 -0.440 -0.519 -0.614 -0.717 -0.859 -1.134 -1.548 -1.953 -2.732
kinematic win a rep a | S a -0.119 -0.170 -0.227 -0.285 -0.353 -0.429 -0.513 -0.628 -0.754 -0.893 -1.097 -1.507 -1.730 -3.685
peri | S b -0.102 -0.132 -0.141 -0.164 -0.181 -0.214 -0.250 -0.260 -0.263 -0.277 -0.293 -0.297 -0.302 -0.308
kinematic win a rep a | S b -0.112 -0.168 -0.222 -0.278 -0.346 -0.399 -0.439 -0.489 -0.559 -0.682 -0.805 -0.930 -0.993 -1.068
Largest SD over draws of a curve point: 0.403
dim observed chance
endpoints + max speed 3 -0.045 -0.222
initial velocity 2 0.030 -0.172
endpoints + orth. excursion 3 0.001 -0.217
endpoints + viapoint 4 -0.026 -0.286
endpoints + max velocity (pk) 4 0.007 -0.291
endpoints + max velocity (ex) 4 -0.005 -0.298
endpoints + halfpoint 4 -0.012 -0.289
segment dirs (t) 6 0.008 -0.429
kitchen sink (t) 11 -0.244 -1.128
segment dirs (l) 6 -0.027 -0.441
kitchen sink (l) 11 -0.221 -0.978
('peri', 'a') | observed above chance at every k: True | observed - chance, k = 1..14: [0.115 0.217 0.274 0.375 0.463 0.557 0.621 0.675 0.743 0.852 0.949 1.151
1.379 1.751]
('peri', 'b') | observed above chance at every k: True | observed - chance, k = 1..14: [0.098 0.188 0.206 0.272 0.31 0.359 0.445 0.463 0.467 0.483 0.513 0.522
0.528 0.534]
('kinematic win a rep a', 'a') | observed above chance at every k: True | observed - chance, k = 1..14: [0.127 0.186 0.24 0.301 0.348 0.389 0.443 0.499 0.578 0.677 0.769 1.027
1.034 1.68 ]
('kinematic win a rep a', 'b') | observed above chance at every k: True | observed - chance, k = 1..14: [0.119 0.185 0.236 0.297 0.344 0.38 0.41 0.444 0.476 0.54 0.7 0.784
0.805 0.872]
Task spaces above chance: 11 of 11
Variant, PD neurons with permutation p < 0.05 (kinematic window a, replication a): S neurons peri_best peri_best_k peri_k14 kin_best kin_best_k task_best task_best_space seg_dirs_t init_vel endpoints_max_speed peri_highest a 43 0.336 7 -0.164 0.116 3 0.143 segment dirs (t) 0.143 0.102 0.02 True b 43 0.415 14 0.415 0.117 3 0.143 segment dirs (t) 0.143 0.102 0.02 True
Result: reproduced, resting on the tuning measure. Under the main reading the best peri-movement performance is 0.117 (k = 6) with S fixed and 0.226 (k = 13) with S recomputed, above the best kinematic space (0.016 to 0.024) and the best task space (initial velocity, 0.030) under all 8 combinations of S, kinematic window and replication. Over the selection readings the peri-movement space is highest under every combination in 7 of 8 and under some in the eighth. All performances are low: the split-half reliability of the preparatory tuning is a median r of 0.374. The paper's upper limit, whose T-matrix includes the fitted neuron, is 0.178 (S fixed) and 0.300 (recomputed); the variant without the fitted neuron (0.068, 0.184) is below the peri-movement best and so is no upper limit here. Neurons of array 1 carry most of the performance (best 0.291, array 2 -0.006).
Chance was computed for the peri-movement space, the kinematic space with window a and replication a (both readings of S) and the task spaces, and each is above its chance level at every dimensionality: with conditions shuffled, the peri-movement space gives -0.120 to -2.732 (S fixed) and -0.102 to -0.308 (recomputed), the task spaces -0.172 to -1.128, so the paper's "above zero is better than chance" holds for these spaces. Chance was not computed for kinematic window b, replication b or the upper limits. With labels shuffled a median of 59 of 152 units pass the 5 spikes/s range, against 75 observed, so noise alone matches most of the passing count; 43 of the 76 PD neurons have tuning distinguishable from noise (permutation p < 0.05). On those 43 the result is stronger: peri-movement best 0.336 (S fixed) and 0.415 (recomputed), against 0.116 to 0.117 for the kinematic space and 0.143 for the best task space (segment directions, time segments).
fig6-perimovement-best-each-dimensionality¶
Even taken one dimensionality at a time, PD performance is almost always higher in the peri-movement space than in other spaces of the same dimensionality.
TASK_SETS = {"pk, t": ["endpoints + max speed", "initial velocity", "endpoints + orth. excursion", "endpoints + viapoint",
"endpoints + max velocity (pk)", "endpoints + halfpoint", "segment dirs (t)", "kitchen sink (t)"],
"ex, l": ["endpoints + max speed", "initial velocity", "endpoints + orth. excursion", "endpoints + viapoint",
"endpoints + max velocity (ex)", "endpoints + halfpoint", "segment dirs (l)", "kitchen sink (l)"]}
rows = []
for R in "ab":
peri = F6["curves"][("peri", R)].mean(0)
for kw in "ab":
for rep in "ab":
kin = F6["curves"][(f"kinematic win {kw} rep {rep}", R)].mean(0)
for ts, names in TASK_SETS.items():
wins = [peri[i] > kin[i] for i in range(len(KS))]
wins += [peri[dims[nm] - 1] > F6["task"][nm].mean() for nm in names]
rows.append(dict(S=R, kin_window=kw, kin_rep=rep, task_variant=ts, comparisons=len(wins),
peri_wins=int(sum(wins)), fraction=round(np.mean(wins), 3)))
f6e = pd.DataFrame(rows)
print(f6e.to_string(index=False))
print("Fraction range over readings:", f6e.fraction.min(), "to", f6e.fraction.max())
S kin_window kin_rep task_variant comparisons peri_wins fraction a a a pk, t 22 21 0.955 a a a ex, l 22 21 0.955 a a b pk, t 22 21 0.955 a a b ex, l 22 21 0.955 a b a pk, t 22 20 0.909 a b a ex, l 22 20 0.909 a b b pk, t 22 20 0.909 a b b ex, l 22 20 0.909 b a a pk, t 22 21 0.955 b a a ex, l 22 21 0.955 b a b pk, t 22 21 0.955 b a b ex, l 22 21 0.955 b b a pk, t 22 21 0.955 b b a ex, l 22 21 0.955 b b b pk, t 22 21 0.955 b b b ex, l 22 21 0.955 Fraction range over readings: 0.909 to 0.955
Result: reproduced. The peri-movement space is higher in 20 to 21 of 22 comparisons (0.909 to 0.955) under every reading of S, the kinematic window and replication, and the task-space variants, above the 80% the criterion asks.
fig6-dimensionality-peak¶
Performance of the PCA-based spaces typically improves with dimensionality, peaks usually at a middle range, and declines at dimensionalities above 10, a sign of over-fitting.
rows = []
for (nm, R), v in F6["curves"].items():
if nm.startswith("upper incl"):
continue
m = v.mean(0)
kb = KS[m.argmax()]
rows.append(dict(space=nm, S=R, k1=m[0].round(3), best=m.max().round(3), best_k=kb, k14=m[-1].round(3),
above_k1=bool(m.max() > m[0]), best_k_2_10=bool(2 <= kb <= 10), k14_below_best=bool(m[-1] < m.max())))
f6p = pd.DataFrame(rows)
f6p["all_three"] = f6p.above_k1 & f6p.best_k_2_10 & f6p.k14_below_best
print(f6p.to_string(index=False))
for r_ in SEL_READINGS:
m = F6_all[r_]["curves"][("peri", "a")].mean(0); mb = F6_all[r_]["curves"][("peri", "b")].mean(0)
print(r_, "peri best k (S a):", KS[m.argmax()], "k14 below best:", m[-1] < m.max(), "| (S b):", KS[mb.argmax()], mb[-1] < mb.max())
for R in "ab":
ch = NULL[("peri", R)].mean(0)
print(f"S {R}: chance curve of the peri-movement space from {ch[0]:.3f} (k = 1) to {ch[-1]:.3f} (k = 14);",
f"fall from k = 6 to k = 14 {ch[5] - ch[-1]:.3f}")
space S k1 best best_k k14 above_k1 best_k_2_10 k14_below_best all_three
peri a -0.006 0.117 6 -0.980 True True True True
upper a -0.014 0.068 5 -0.838 True True True True
kinematic win a rep a a 0.008 0.016 2 -2.005 True True True True
kinematic win a rep b a 0.008 0.016 2 -2.011 True True True True
kinematic win b rep a a 0.011 0.022 2 -2.317 True True True True
kinematic win b rep b a 0.011 0.022 2 -2.319 True True True True
peri b -0.004 0.226 13 0.226 True False True False
upper b -0.015 0.184 8 0.098 True True True True
kinematic win a rep a b 0.007 0.019 4 -0.196 True True True True
kinematic win a rep b b 0.007 0.019 4 -0.196 True True True True
kinematic win b rep a b 0.010 0.024 2 -0.187 True True True True
kinematic win b rep b b 0.009 0.024 2 -0.187 True True True True
('a', 'mean', 'separate') peri best k (S a): 4 k14 below best: True | (S b): 11 True
('a', 'mean', 'concat') peri best k (S a): 6 k14 below best: True | (S b): 13 True
('a', 'medoid', 'separate') peri best k (S a): 2 k14 below best: True | (S b): 11 True
('a', 'medoid', 'concat') peri best k (S a): 6 k14 below best: True | (S b): 14 False
('b', 'mean', 'separate') peri best k (S a): 6 k14 below best: True | (S b): 14 False
('b', 'mean', 'concat') peri best k (S a): 6 k14 below best: True | (S b): 14 False
('b', 'medoid', 'separate') peri best k (S a): 2 k14 below best: True | (S b): 11 True
('b', 'medoid', 'concat') peri best k (S a): 6 k14 below best: True | (S b): 12 True
S a: chance curve of the peri-movement space from -0.120 (k = 1) to -2.732 (k = 14); fall from k = 6 to k = 14 2.292
S b: chance curve of the peri-movement space from -0.102 (k = 1) to -0.308 (k = 14); fall from k = 6 to k = 14 0.094
Result: reproduced, resting on the reading of S. With S fixed (adopted, below) the peri-movement curve rises from -0.006 at k = 1 to 0.117 at k = 6 and falls to -0.980 at k = 14; the kinematic curves peak at k = 2 and fall to about -2. With S recomputed in each fold the peri-movement curve peaks at k = 13 (0.226), outside 2 to 10, and is nearly flat above k = 10. S fixed is adopted because the Methods describe one matrix S per space and refit only the PD on the remaining conditions. The chance curves show that the two readings penalize extra dimensions differently: with conditions shuffled the peri-movement curve falls by 2.292 from k = 6 to k = 14 with S fixed and by 0.094 with S recomputed, so under the recomputed reading little decline is possible even for noise. With S fixed the best k is 2 to 6 in all 8 selection readings.
fig6-maze-task-spaces¶
Segment directions and initial velocity are among the task spaces that perform reasonably well, segment directions doing so for the maze task, while for the maze task the endpoints + maximum speed space does not perform well.
rows = []
for seg in ["segment dirs (t)", "segment dirs (l)"]:
for R in "ab":
sd_, iv, ems = F6["task"][seg].mean(), F6["task"]["initial velocity"].mean(), F6["task"]["endpoints + max speed"].mean()
p3 = F6["curves"][("peri", R)].mean(0)[2]
rows.append(dict(segments=seg, S=R, seg_dirs=round(sd_, 3), init_vel=round(iv, 3), endpoints_max_speed=round(ems, 3),
peri_3d=round(p3, 3), seg_above_0=sd_ > 0, iv_above_0=iv > 0, ems_below_peri3=ems < p3,
ems_below_seg_and_iv=ems < min(sd_, iv)))
f6m = pd.DataFrame(rows)
print(f6m.to_string(index=False))
for r_ in SEL_READINGS:
T = F6_all[r_]["task"]
print(r_, {k: round(T[k].mean(), 3) for k in ["segment dirs (t)", "segment dirs (l)", "initial velocity", "endpoints + max speed"]})
print("Chance (conditions shuffled):", {k: round(null_task[k], 3) for k in ["segment dirs (t)", "segment dirs (l)", "initial velocity", "endpoints + max speed"]})
print("Selection readings with segment dirs (t) > 0:", sum(F6_all[r_]["task"]["segment dirs (t)"].mean() > 0 for r_ in SEL_READINGS),
"| initial velocity > 0:", sum(F6_all[r_]["task"]["initial velocity"].mean() > 0 for r_ in SEL_READINGS), "of", len(SEL_READINGS))
segments S seg_dirs init_vel endpoints_max_speed peri_3d seg_above_0 iv_above_0 ems_below_peri3 ems_below_seg_and_iv
segment dirs (t) a 0.008 0.03 -0.045 0.052 True True True True
segment dirs (t) b 0.008 0.03 -0.045 0.065 True True True True
segment dirs (l) a -0.027 0.03 -0.045 0.052 False True True True
segment dirs (l) b -0.027 0.03 -0.045 0.065 False True True True
('a', 'mean', 'separate') {'segment dirs (t)': np.float64(-0.118), 'segment dirs (l)': np.float64(-0.089), 'initial velocity': np.float64(-0.015), 'endpoints + max speed': np.float64(-0.078)}
('a', 'mean', 'concat') {'segment dirs (t)': np.float64(0.008), 'segment dirs (l)': np.float64(-0.027), 'initial velocity': np.float64(0.03), 'endpoints + max speed': np.float64(-0.045)}
('a', 'medoid', 'separate') {'segment dirs (t)': np.float64(-0.156), 'segment dirs (l)': np.float64(-0.125), 'initial velocity': np.float64(-0.028), 'endpoints + max speed': np.float64(-0.106)}
('a', 'medoid', 'concat') {'segment dirs (t)': np.float64(-0.021), 'segment dirs (l)': np.float64(-0.052), 'initial velocity': np.float64(0.027), 'endpoints + max speed': np.float64(-0.053)}
('b', 'mean', 'separate') {'segment dirs (t)': np.float64(-0.109), 'segment dirs (l)': np.float64(-0.064), 'initial velocity': np.float64(-0.003), 'endpoints + max speed': np.float64(-0.065)}
('b', 'mean', 'concat') {'segment dirs (t)': np.float64(0.006), 'segment dirs (l)': np.float64(-0.021), 'initial velocity': np.float64(0.04), 'endpoints + max speed': np.float64(-0.04)}
('b', 'medoid', 'separate') {'segment dirs (t)': np.float64(-0.142), 'segment dirs (l)': np.float64(-0.106), 'initial velocity': np.float64(-0.024), 'endpoints + max speed': np.float64(-0.098)}
('b', 'medoid', 'concat') {'segment dirs (t)': np.float64(-0.02), 'segment dirs (l)': np.float64(-0.047), 'initial velocity': np.float64(0.028), 'endpoints + max speed': np.float64(-0.051)}
Chance (conditions shuffled): {'segment dirs (t)': np.float64(-0.429), 'segment dirs (l)': np.float64(-0.441), 'initial velocity': np.float64(-0.172), 'endpoints + max speed': np.float64(-0.222)}
Selection readings with segment dirs (t) > 0: 2 | initial velocity > 0: 4 of 8
Result: reproduced, resting on segments equal in time. The claims file names both segment divisions as readings, and the criterion holds under one: under the main reading segment directions with segments equal in time (0.008) and initial velocity (0.030) are above zero, and endpoints + maximum speed (-0.045) falls below the 3-D peri-movement space (0.052 and 0.065, reading a) and below both (reading b). With segments equal in path length segment directions are below zero (-0.027). All three are far above their chance levels (-0.429, -0.172 and -0.222). The result is fragile across the trial selection: segment directions (time) are above zero in 2 of 8 selection readings and initial velocity in 4. On the 43 PD neurons whose tuning passes the permutation test it is clearer: segment directions 0.143, initial velocity 0.102, endpoints + maximum speed 0.020.
high-dimensional-task-spaces-underperform¶
The high-dimensional task spaces, including the kitchen-sink space, always perform below the peri-movement space, so high dimensionality alone does not explain its performance.
rows = []
for R in "ab":
peri = F6["curves"][("peri", R)].mean(0)
for nm in ["kitchen sink (t)", "kitchen sink (l)", "segment dirs (t)", "segment dirs (l)"]:
v = F6["task"][nm].mean()
rows.append(dict(space=nm, dim=dims[nm], S=R, perf=round(v, 3), peri_same_k=peri[dims[nm] - 1].round(3),
peri_best=peri.max().round(3), below_same_k=v < peri[dims[nm] - 1], below_best=v < peri.max()))
f6h = pd.DataFrame(rows)
print(f6h.to_string(index=False))
for r_ in SEL_READINGS:
T, Cc = F6_all[r_]["task"], F6_all[r_]["curves"]
print(r_, "kitchen sink (t)", round(T["kitchen sink (t)"].mean(), 3), "| peri 11-D (S a)", Cc[("peri", "a")].mean(0)[10].round(3),
"| peri best (S a)", Cc[("peri", "a")].mean(0).max().round(3))
space dim S perf peri_same_k peri_best below_same_k below_best
kitchen sink (t) 11 a -0.244 -0.184 0.117 True True
kitchen sink (l) 11 a -0.221 -0.184 0.117 True True
segment dirs (t) 6 a 0.008 0.117 0.117 True True
segment dirs (l) 6 a -0.027 0.117 0.117 True True
kitchen sink (t) 11 b -0.244 0.221 0.226 True True
kitchen sink (l) 11 b -0.221 0.221 0.226 True True
segment dirs (t) 6 b 0.008 0.145 0.226 True True
segment dirs (l) 6 b -0.027 0.145 0.226 True True
('a', 'mean', 'separate') kitchen sink (t) -0.445 | peri 11-D (S a) -4.887 | peri best (S a) 0.0
('a', 'mean', 'concat') kitchen sink (t) -0.244 | peri 11-D (S a) -0.184 | peri best (S a) 0.117
('a', 'medoid', 'separate') kitchen sink (t) -0.482 | peri 11-D (S a) -9.046 | peri best (S a) -0.042
('a', 'medoid', 'concat') kitchen sink (t) -0.286 | peri 11-D (S a) -0.493 | peri best (S a) 0.076
('b', 'mean', 'separate') kitchen sink (t) -0.382 | peri 11-D (S a) -0.297 | peri best (S a) 0.038
('b', 'mean', 'concat') kitchen sink (t) -0.245 | peri 11-D (S a) -0.209 | peri best (S a) 0.12
('b', 'medoid', 'separate') kitchen sink (t) -0.498 | peri 11-D (S a) -4.114 | peri best (S a) -0.024
('b', 'medoid', 'concat') kitchen sink (t) -0.272 | peri 11-D (S a) -0.236 | peri best (S a) 0.091
Result: reproduced. Under the main reading the kitchen-sink spaces (-0.244, -0.221) and segment directions (0.008, -0.027) are below the peri-movement space at the same dimensionality and below its best (0.117 with S fixed, 0.226 recomputed), under both readings of the comparison and of S.
some-neurons-better-elsewhere¶
Although the peri-movement space performs best on average, some individual neurons are better captured in one of the other spaces.
rows = []
for R in "ab":
pc = F6["curves"][("peri", R)]
kbest = pc.mean(0).argmax()
peri_n = pc[:, kbest]
others = {}
kc = F6["curves"][("kinematic win a rep a", R)]
others["kinematic (win a rep a)"] = kc[:, kc.mean(0).argmax()]
for nm, v in F6["task"].items():
others[nm] = v
O = np.column_stack(list(others.values()))
better = O.max(1) > peri_n
rows.append(dict(S=R, peri_k=KS[kbest], neurons=len(peri_n), better_elsewhere=int(better.sum()),
fraction=round(better.mean(), 3)))
if R == "a":
print("Neurons better in each space (S a):", {k: int((v > peri_n).sum()) for k, v in others.items()})
print(pd.DataFrame(rows).to_string(index=False))
Neurons better in each space (S a): {'kinematic (win a rep a)': 28, 'endpoints + max speed': 22, 'initial velocity': 33, 'endpoints + orth. excursion': 26, 'endpoints + viapoint': 19, 'endpoints + max velocity (pk)': 25, 'endpoints + max velocity (ex)': 27, 'endpoints + halfpoint': 23, 'segment dirs (t)': 32, 'kitchen sink (t)': 17, 'segment dirs (l)': 24, 'kitchen sink (l)': 19}
S peri_k neurons better_elsewhere fraction
a 6 76 60 0.789
b 13 76 44 0.579
Result: reproduced. 60 of 76 neurons (0.789) with S fixed, and 44 of 76 (0.579) with S recomputed, perform better in at least one other space than in the peri-movement space at its best k.
fig6-preparatory-dimensionality¶
Preparatory activity needs 7 to 10 dimensions, depending on the dataset, to capture 90% of its variance.
rows = []
for r_ in SEL_READINGS:
d = condition_data(*r_)
prep_ok, _ = neuron_sets(d, "range")
for nset, mask in [("a PD neurons", prep_ok), ("b all units", np.ones(len(units), bool))]:
X = d["prep"][:, mask]
for norm in ["a none", "b unit range"]:
Xn = X if norm.startswith("a") else X / (X.max(0) - X.min(0) + 1e-12)
ev = M.explained(Xn @ Xn.T)
k90 = int(np.searchsorted(np.cumsum(ev), 0.9) + 1)
rows.append(dict(onset=r_[0], proto=r_[1], combine=r_[2], neurons=nset, n=int(mask.sum()), normalization=norm,
pcs_90=k90, in_6_3_to_11=6.3 <= k90 <= 11))
f6dim = pd.DataFrame(rows)
print(f6dim.to_string(index=False))
print("Bound c - 1 =", len(conds) - 1)
# noise reference: the same count when trial labels are shuffled within the main selection
rng_dim = np.random.default_rng(6)
sel = d0["sel"]
go_all = M.bin_spikes(st, tr.go_cue_time.values, (M.PREP[0], M.PREP[1]))
cnt = go_all[sel].sum(2).astype(float) / ((M.PREP[1] - M.PREP[0]) / 1000)
k_sh = []
for _ in range(20):
cs = rng_dim.permutation(cond[sel])
Xs = np.stack([cnt[cs == c].mean(0) for c in conds])[:, neuron_sets(d0, "range")[0]]
ev = M.explained(Xs @ Xs.T); k_sh.append(int(np.searchsorted(np.cumsum(ev), 0.9) + 1))
print("Reference: PCs for 90% with condition labels shuffled (main selection, PD neurons, 20 draws): median", np.median(k_sh),
"range", min(k_sh), "to", max(k_sh))
# Noise-corrected variant (verifier, round 1): PCs from the odd-trial condition means, their variance measured as the
# cross product of odd- and even-trial projections (and the reverse, averaged), so that noise independent between the
# halves drops out. Count of PCs for 90% of the summed cross-validated variance. Main selection, window counts.
def cv_dim(A, B):
Ac, Bc = A - A.mean(0), B - B.mean(0)
Vt = np.linalg.svd(Ac, full_matrices=False)[2]
return np.array([(Ac @ v) @ (Bc @ v) for v in Vt])
for nset, mask in [("a PD neurons", pdn), ("b all units", np.ones(len(units), bool))]:
A_, B_ = halves[0][:, mask], halves[1][:, mask]
sig = (cv_dim(A_, B_) + cv_dim(B_, A_)) / 2
cum = np.cumsum(sig) / sig.sum()
raw = [int(np.searchsorted(np.cumsum(M.explained(Z @ Z.T)), 0.9) + 1) for Z in (A_, B_)]
print(f"{nset}: cross-validated PCs for 90% = {int(np.searchsorted(cum, 0.9) + 1)}",
f"| first 10 shares {np.round(sig[:10] / sig.sum(), 3)}",
f"| plain count on each half {raw}")
onset proto combine neurons n normalization pcs_90 in_6_3_to_11
a mean separate a PD neurons 85 a none 13 False
a mean separate a PD neurons 85 b unit range 16 False
a mean separate b all units 152 a none 14 False
a mean separate b all units 152 b unit range 19 False
a mean concat a PD neurons 76 a none 12 False
a mean concat a PD neurons 76 b unit range 16 False
a mean concat b all units 152 a none 13 False
a mean concat b all units 152 b unit range 19 False
a medoid separate a PD neurons 86 a none 14 False
a medoid separate a PD neurons 86 b unit range 16 False
a medoid separate b all units 152 a none 15 False
a medoid separate b all units 152 b unit range 19 False
a medoid concat a PD neurons 77 a none 13 False
a medoid concat a PD neurons 77 b unit range 16 False
a medoid concat b all units 152 a none 14 False
a medoid concat b all units 152 b unit range 19 False
b mean separate a PD neurons 81 a none 13 False
b mean separate a PD neurons 81 b unit range 17 False
b mean separate b all units 152 a none 14 False
b mean separate b all units 152 b unit range 19 False
b mean concat a PD neurons 74 a none 12 False
b mean concat a PD neurons 74 b unit range 16 False
b mean concat b all units 152 a none 13 False
b mean concat b all units 152 b unit range 19 False
b medoid separate a PD neurons 89 a none 14 False
b medoid separate a PD neurons 89 b unit range 17 False
b medoid separate b all units 152 a none 15 False
b medoid separate b all units 152 b unit range 19 False
b medoid concat a PD neurons 77 a none 13 False
b medoid concat a PD neurons 77 b unit range 16 False
b medoid concat b all units 152 a none 14 False
b medoid concat b all units 152 b unit range 19 False
Bound c - 1 = 26
Reference: PCs for 90% with condition labels shuffled (main selection, PD neurons, 20 draws): median 17.0 range 17 to 17 a PD neurons: cross-validated PCs for 90% = 7 | first 10 shares [0.363 0.196 0.104 0.11 0.074 0.039 0.029 0.014 0.017 0.012] | plain count on each half [15, 15] b all units: cross-validated PCs for 90% = 7 | first 10 shares [0.358 0.195 0.101 0.111 0.078 0.04 0.032 0.014 0.015 0.007] | plain count on each half [16, 16]
Result: not reproduced. 90% of the variance of the preparatory matrix takes 12 to 19 PCs over every reading (12 for the main reading's PD neurons, unnormalized), outside 6.3 to 11. With condition labels shuffled the count is 17, so with 6 to 10 trials per condition the count reflects noise more than the structure of preparatory activity. A count that removes noise independent between halves of the trials (cross-validated PCA, odd against even trials) gives 7 PCs for both the PD neurons and all units, inside 6.3 to 11, while the plain count on each half gives 15 to 16. The criterion's count is the plain one, so the outcome stands, but the failure is likely due to the few trials per condition here.
fig7a-shuffled-space-poor¶
A peri-movement space built after randomly inverting each condition's peri-movement response on half the conditions performs poorly.
rng7a = np.random.default_rng(7)
X = d0["prep"]
r0 = M.MOVE_WIN[0]
ref = d0["mv_rate"][:, :, M.TMAT[0] - r0][:, :, None]
seg = d0["mv_rate"][:, :, M.TMAT[0] - r0: M.TMAT[1] - r0]
sh_curves = {"a": [], "b": []}
for draw in range(20):
flip = rng7a.random(seg.shape[:2]) < 0.5
segs = np.where(flip[:, :, None], 2 * ref - seg, seg)
Gu = M.grams_per_unit(segs, M.TMAT[0])
for R in "ab":
sh_curves[R].append(M.perf_curves_pca(X, Gu.sum(0), Gu, F6["fit"], KS, R == "b").mean(0))
rows = []
for R in "ab":
sh = np.mean(sh_curves[R], 0)
peri = F6["curves"][("peri", R)].mean(0)
kin = max(F6["curves"][(f"kinematic win {kw} rep {rp}", R)].mean(0).max() for kw in "ab" for rp in "ab")
task = task_tab.max()
rows.append(dict(S=R, shuffled_best=sh.max().round(3), shuffled_best_k=KS[sh.argmax()], peri_best=peri.max().round(3),
best_kin=round(kin, 3), best_task=round(task, 3), a_below_peri=sh.max() < peri.max(),
b_below_kin_and_task=sh.max() < min(kin, task),
draws_best_range=f"{np.min([c.max() for c in sh_curves[R]]):.3f} to {np.max([c.max() for c in sh_curves[R]]):.3f}"))
print(pd.DataFrame(rows).to_string(index=False))
print("Shuffled mean curve (S a):", np.round(np.mean(sh_curves["a"], 0), 3))
S shuffled_best shuffled_best_k peri_best best_kin best_task a_below_peri b_below_kin_and_task draws_best_range a -0.069 2 0.117 0.022 0.03 True True -0.099 to -0.022 b -0.010 6 0.226 0.024 0.03 True True -0.045 to 0.042 Shuffled mean curve (S a): [-0.08 -0.069 -0.084 -0.109 -0.156 -0.204 -0.269 -0.332 -0.412 -0.508 -0.618 -0.769 -0.928 -1.117]
Result: reproduced. The shuffled space's best performance, averaged over 20 draws, is -0.069 (S fixed) and -0.010 (recomputed), below the unshuffled peri-movement best (0.117, 0.226) and below the best kinematic and task spaces (readings a and b).
fig7b-window-start-slide¶
Sliding the start of the peri-movement window from 150 ms before movement onset to movement onset barely lowers PD performance, to 97% (monkey J) and 94% (monkeys A, B) of its original value.
Result: reproduced. With the window end fixed, performance at the best k falls from 0.117 to 0.107 (fraction 0.911) with S fixed and from 0.226 to 0.218 (0.962) with S recomputed, both within 0.97 +/- 0.097. With the whole window shifted the fractions are 0.853 and 0.937.
Anterior and posterior sites¶
The electrode table gives no positions (x, y, z are empty), so AP location can only be read as the array (reading (b) of the claims file): PMd anterior, M1 posterior. The table's location column puts PMd on electrode ids 1 to 96 and M1 on 97 to 192, but every unit references a PMd row; the array is read from the unit id (prefix 1 or 2). Prefix 1 is taken as PMd because the electrode table lists the PMd group first. This is an assumption that the files cannot check (cell below).
eg_df = D["nwb"].electrodes.to_dataframe()
print("Electrode table columns:", list(eg_df.columns), "| location by row:",
{loc: (int(eg_df.index[eg_df.location == loc].min()), int(eg_df.index[eg_df.location == loc].max())) for loc in eg_df.location.unique()},
"(electrode ids)")
for k, g in D["nwb"].electrode_groups.items():
print("Group", k, "| location:", g.location, "| description:", g.description, "| device:", g.device.description)
print("Units table columns:", list(D["nwb"].units.colnames), "| description:", D["nwb"].units.description)
print("Location of the electrodes the units reference:", units.electrode_group_in_file.value_counts().to_dict(),
"| prefix 1 rows", units.electrode_row[units.id_array == 1].min(), "to", units.electrode_row[units.id_array == 1].max(),
"; prefix 2 rows", units.electrode_row[units.id_array == 2].min(), "to", units.electrode_row[units.id_array == 2].max())
print("Nothing in the units or electrode tables ties id prefix 1 or 2 to PMd or M1.")
Electrode table columns: ['x', 'y', 'z', 'imp', 'location', 'filtering', 'group', 'group_name'] | location by row: {'PMd': (1, 96), 'M1': (97, 192)} (electrode ids)
Group electrode_group_M1 | location: Primary motor cortex | description: Electrodes in an implanted Utah array in M1 | device: 96-electrode Utah array in M1
Group electrode_group_PMd | location: Dorsal premotor cortex | description: Electrodes in implanted Utah Array in PMd | device: 96-electrode Utah array in PMd
Units table columns: ['heldout', 'spike_times', 'obs_intervals', 'electrodes'] | description: data on spiking units
Location of the electrodes the units reference: {'electrode_group_PMd': 152} | prefix 1 rows 0 to 89 ; prefix 2 rows 14 to 94
Nothing in the units or electrode tables ties id prefix 1 or 2 to PMd or M1.
ap-location-no-effect-j-array¶
For the J-array dataset the correlation between preparatory and peri-movement tuning does not depend on anterior-posterior location (r = 0.20, p = 0.27).
rows = []
for r_ in SEL_READINGS:
d = condition_data(*r_)
idx, r, _ = corr_set(d, neuron_sets(d, "range")[1])
ap = (units.id_array.values[idx] == 1).astype(float) # 1 = PMd (anterior)
if len(r) > 2 and 0 < ap.sum() < len(ap):
lr = stats.linregress(ap, r)
rows.append(dict(onset=r_[0], proto=r_[1], combine=r_[2], n=len(r), n_pmd=int(ap.sum()), r=round(lr.rvalue, 3),
p=round(lr.pvalue, 4), mean_r_pmd=round(r[ap == 1].mean(), 3), mean_r_m1=round(r[ap == 0].mean(), 3)))
fap = pd.DataFrame(rows)
fap["absence_holds_p>=0.05"] = fap.p >= 0.05
print(fap.to_string(index=False))
onset proto combine n n_pmd r p mean_r_pmd mean_r_m1 absence_holds_p>=0.05
a mean separate 16 8 -0.063 0.8172 0.145 0.182 True
a mean concat 17 9 0.165 0.5258 0.221 0.117 True
a medoid separate 15 8 -0.051 0.8562 0.147 0.176 True
a medoid concat 15 8 0.133 0.6367 0.209 0.127 True
b mean separate 20 12 0.082 0.7325 0.182 0.135 True
b mean concat 18 10 0.030 0.9073 0.215 0.197 True
b medoid separate 15 8 -0.126 0.6544 0.107 0.179 True
b medoid concat 17 10 0.114 0.6621 0.212 0.142 True
Result: reproduced, resting on the array as the AP location. Regressing each neuron's tuning correlation on its array gives r = 0.165, p = 0.5258 under the main reading (17 neurons, 9 on array 1), and p from 0.53 to 0.91 over the selection readings, so the paper's test finds no effect. The paper's J-array value is r = 0.20, p = 0.27.
prep-activity-more-anterior¶
Preparatory activity is more prevalent at anterior sites.
rows = []
for r_ in SEL_READINGS:
for how in ["range", "sd"]:
prep_ok, _ = neuron_sets(condition_data(*r_), how)
a1, a2 = units.id_array.values == 1, units.id_array.values == 2
tab = [[int(prep_ok[a1].sum()), int((~prep_ok[a1]).sum())], [int(prep_ok[a2].sum()), int((~prep_ok[a2]).sum())]]
rows.append(dict(onset=r_[0], proto=r_[1], combine=r_[2], tuning=how, pmd_pass=f"{tab[0][0]}/{a1.sum()}",
m1_pass=f"{tab[1][0]}/{a2.sum()}", frac_pmd=round(tab[0][0] / a1.sum(), 3),
frac_m1=round(tab[1][0] / a2.sum(), 3), fisher_p=round(stats.fisher_exact(tab)[1], 4)))
fpa = pd.DataFrame(rows)
fpa["a_higher_pmd"] = fpa.frac_pmd > fpa.frac_m1
fpa["b_fisher"] = fpa.a_higher_pmd & (fpa.fisher_p < 0.05)
print(fpa.to_string(index=False))
tp_ = M.tuning_strength(d0["prep"], "range")
print("Median preparatory tuning range (spikes/s), main reading: PMd", np.median(tp_[units.id_array == 1]).round(2),
"| M1", np.median(tp_[units.id_array == 2]).round(2), "| Mann-Whitney p",
round(stats.mannwhitneyu(tp_[units.id_array == 1], tp_[units.id_array == 2]).pvalue, 4))
onset proto combine tuning pmd_pass m1_pass frac_pmd frac_m1 fisher_p a_higher_pmd b_fisher
a mean separate range 37/75 48/77 0.493 0.623 0.1413 False False
a mean separate sd 4/75 1/77 0.053 0.013 0.2064 True False
a mean concat range 35/75 41/77 0.467 0.532 0.5165 False False
a mean concat sd 4/75 1/77 0.053 0.013 0.2064 True False
a medoid separate range 39/75 47/77 0.520 0.610 0.3263 False False
a medoid separate sd 5/75 3/77 0.067 0.039 0.4916 True False
a medoid concat range 35/75 42/77 0.467 0.545 0.4173 False False
a medoid concat sd 4/75 2/77 0.053 0.026 0.4391 True False
b mean separate range 36/75 45/77 0.480 0.584 0.2550 False False
b mean separate sd 4/75 1/77 0.053 0.013 0.2064 True False
b mean concat range 35/75 39/77 0.467 0.506 0.6308 False False
b mean concat sd 4/75 2/77 0.053 0.026 0.4391 True False
b medoid separate range 39/75 50/77 0.520 0.649 0.1382 False False
b medoid separate sd 4/75 3/77 0.053 0.039 0.7173 True False
b medoid concat range 35/75 42/77 0.467 0.545 0.4173 False False
b medoid concat sd 4/75 2/77 0.053 0.026 0.4391 True False
Median preparatory tuning range (spikes/s), main reading: PMd 4.69 | M1 5.12 | Mann-Whitney p 0.4091
Result: blocked (metadata_insufficient), after the analysis. The direction of the comparison depends on which array is PMd, and nothing in the files says (cell above). With the range measure 35 of 75 units on array 1 (0.467) and 41 of 77 on array 2 (0.532) pass the preparatory threshold under the main reading (Fisher p = 0.5165); array 1 has the lower fraction in every selection reading. With the SD measure array 1 is higher (4 of 75 against 1 of 77, p = 0.2064). If array 1 is PMd the range result goes against the claim; if it is M1 the result favours it, without reaching p < 0.05 in any reading. With labels shuffled 59 of 152 units pass the range threshold, so the fractions are largely noise.
fig7c-anterior-posterior-performance¶
The peri-movement space captures preparatory tuning at both anterior and posterior sites, with normalized performance only slightly different from that of randomly divided data, and slightly better at anterior sites.
rng7c = np.random.default_rng(70)
prep_ok0 = neuron_sets(d0, "range")[0]
Gu_all = M.grams_per_unit(d0["mv_rate"], M.MOVE_WIN[0])
def half_perf(members, R):
"""6-D peri-movement performance of the PD neurons in `members`, from a space built from
`members` only, and the upper limit from the same units at 6-D and at its best k. The upper
limit's T-matrix holds every neuron of the half, the fitted one included (Supplemental
Experimental Procedures, Estimating the upper limit: "each row containing the mean preparatory
activity for every neuron"); the variant without the fitted neuron is also returned."""
mem = np.where(members)[0]
fit = np.where(members & prep_ok0)[0]
Gu = Gu_all[mem]; Gt = Gu.sum(0)
Gp = M.prep_grams(X[:, mem])
loc = {n: i for i, n in enumerate(mem)}
peri = np.mean([M.pd_performance_pca(X[:, n], Gt - Gu[loc[n]], [6], R == "b")[0] for n in fit])
up_in = np.mean([M.pd_performance_pca(X[:, n], Gp.sum(0), KS, R == "b") for n in fit], 0)
up_ex = np.mean([M.pd_performance_pca(X[:, n], Gp.sum(0) - Gp[loc[n]], KS, R == "b") for n in fit], 0)
return peri, up_in[5], up_in.max(), up_ex[5], up_ex.max(), len(fit)
rows = []
for R in "ab":
res = {"all": half_perf(np.ones(len(units), bool), R),
"array 1": half_perf(units.id_array.values == 1, R),
"array 2": half_perf(units.id_array.values == 2, R)}
rand = []
for _ in range(20):
m = np.zeros(len(units), bool); m[rng7c.choice(len(units), len(units) // 2, replace=False)] = True
rand += [half_perf(m, R), half_perf(~m, R)]
rand = np.array(rand)
for nm, (p, ui6, uib, ue6, ueb, n) in res.items():
rows.append(dict(S=R, subset=nm, neurons=n, perf_6d=round(p, 3), upper_6d=round(ui6, 3), upper_best=round(uib, 3),
norm_a=round(p / ui6, 3), norm_b=round(p / uib, 3), upper_excl_6d=round(ue6, 3), norm_excl_6d=round(p / ue6, 3)))
rows.append(dict(S=R, subset="random halves (20 splits)", neurons=rand[:, 5].mean(), perf_6d=round(rand[:, 0].mean(), 3),
upper_6d=round(rand[:, 1].mean(), 3), upper_best=round(rand[:, 2].mean(), 3),
norm_a=round((rand[:, 0] / rand[:, 1]).mean(), 3), norm_b=round((rand[:, 0] / rand[:, 2]).mean(), 3),
upper_excl_6d=round(rand[:, 3].mean(), 3), norm_excl_6d=round((rand[:, 0] / rand[:, 3]).mean(), 3)))
print(f"S {R}, random halves: normalized (a) SD {(rand[:, 0] / rand[:, 1]).std():.3f}, range {(rand[:, 0] / rand[:, 1]).min():.3f} to "
f"{(rand[:, 0] / rand[:, 1]).max():.3f}; (b) SD {(rand[:, 0] / rand[:, 2]).std():.3f}")
f7c = pd.DataFrame(rows)
print(f7c.to_string(index=False))
for R in "ab":
for nrm in ["norm_a", "norm_b", "norm_excl_6d"]:
t = f7c[f7c.S == R].set_index("subset")[nrm]
a1, a2, rnd = t["array 1"], t["array 2"], t["random halves (20 splits)"]
print(f"S {R}, {nrm}: both > 0: {a1 > 0 and a2 > 0}; both within 10% of random ({rnd:.3f}): "
f"{abs(a1 - rnd) <= 0.1 * abs(rnd) and abs(a2 - rnd) <= 0.1 * abs(rnd)}; array 1 higher: {a1 > a2}")
S a, random halves: normalized (a) SD 0.153, range 0.067 to 0.869; (b) SD 0.118
S b, random halves: normalized (a) SD 0.103, range 0.352 to 0.736; (b) SD 0.066 S subset neurons perf_6d upper_6d upper_best norm_a norm_b upper_excl_6d norm_excl_6d a all 76.0 0.117 0.176 0.178 0.667 0.657 0.061 1.937 a array 1 35.0 0.271 0.356 0.575 0.761 0.472 0.128 2.111 a array 2 41.0 -0.060 0.150 0.161 -0.402 -0.374 -0.107 0.560 a random halves (20 splits) 38.0 0.095 0.211 0.277 0.455 0.350 -0.043 -0.711 b all 76.0 0.145 0.205 0.300 0.707 0.483 0.115 1.258 b array 1 35.0 0.270 0.352 0.616 0.766 0.438 0.170 1.581 b array 2 41.0 0.049 0.214 0.341 0.230 0.144 0.055 0.892 b random halves (20 splits) 38.0 0.140 0.250 0.417 0.564 0.336 0.084 -0.740 S a, norm_a: both > 0: False; both within 10% of random (0.455): False; array 1 higher: True S a, norm_b: both > 0: False; both within 10% of random (0.350): False; array 1 higher: True S a, norm_excl_6d: both > 0: True; both within 10% of random (-0.711): False; array 1 higher: True S b, norm_a: both > 0: True; both within 10% of random (0.564): False; array 1 higher: True S b, norm_b: both > 0: True; both within 10% of random (0.336): False; array 1 higher: True S b, norm_excl_6d: both > 0: True; both within 10% of random (-0.740): False; array 1 higher: True
Result: not reproduced. With the paper's upper limit (its T-matrix holds every neuron of the half, the fitted one included; Supplemental Experimental Procedures, Estimating the upper limit) neither array is within 10% of the random halves under any reading: at 6-D with S fixed array 1 gives 0.761, array 2 -0.402 and random halves 0.455 (SD 0.153 over splits); with S recomputed 0.766, 0.230 and 0.564; at the best k 0.472, -0.374 and 0.350 (S fixed) and 0.438, 0.144 and 0.336 (recomputed). Array 2 is below zero with S fixed. Array 1 is higher in every reading, which is "anterior higher" only if array 1 is PMd. The variant without the fitted neuron, which the round-1 notebook used, gives 6-D upper limits from -0.107 to 0.170, at or below zero in 2 of 8 (array 2 and random halves, S fixed), and ratios from -0.740 to 2.111. Unnormalized 6-D performance is 0.271 on array 1 and -0.060 on array 2 (S fixed; 0.270 and 0.049 recomputed): the peri-movement space captures tuning on array 1 but hardly on array 2. The within-10% test fails for both halves, so the outcome does not depend on which array is PMd.
fig7d-pd-consistency¶
The PD is most stable across random condition subsets when measured in the peri-movement space, for every comparison and all four datasets.
rng7d = np.random.default_rng(71)
Gt_all = Gu_all.sum(0)
kin_curve = F6["curves"][("kinematic win a rep a", "a")].mean(0)
k_kin = KS[kin_curve.argmax()]
Gk = M.kinematic_gram(d0["pos"], d0["vel"], M.MOVE_WIN[0], len(units) - 1)
S_kin = M.pca_scores(Gk, k_kin)
comps = [("endpoints + max speed", F6["spaces"]["endpoints + max speed"]),
("segment dirs (t)", F6["spaces"]["segment dirs (t)"]),
("segment dirs (l)", F6["spaces"]["segment dirs (l)"]),
(f"kinematic (best k = {k_kin})", S_kin)]
rows = []
for nm, S in comps:
k = S.shape[1]
other, peri = [], []
for n in F6["fit"]:
x = X[:, n]
Sp = M.pca_scores(Gt_all - Gu_all[n], k)
other.append(M.pd_consistency(x, S, rng7d))
peri.append(M.pd_consistency(x, Sp, rng7d))
other, peri = np.array(other), np.array(peri)
rows.append(dict(comparison=nm, dim=k, other_space=round(other.mean(), 3), peri_space=round(peri.mean(), 3),
peri_higher=peri.mean() > other.mean(), neurons_peri_higher=f"{int((peri > other).sum())}/{len(peri)}",
wilcoxon_p=round(stats.wilcoxon(peri, other).pvalue, 4)))
f7d = pd.DataFrame(rows)
print(f7d.to_string(index=False))
print("Kinematic best k from the Figure 6 curve (S a, window a, replication a):", k_kin)
comparison dim other_space peri_space peri_higher neurons_peri_higher wilcoxon_p
endpoints + max speed 3 0.597 0.698 True 56/76 0.0000
segment dirs (t) 6 0.605 0.613 True 39/76 0.8970
segment dirs (l) 6 0.542 0.615 True 45/76 0.0227
kinematic (best k = 2) 2 0.670 0.729 True 54/76 0.0004
Kinematic best k from the Figure 6 curve (S a, window a, replication a): 2
Result: reproduced. PD consistency is higher in the peri-movement space than in endpoints + maximum speed (0.698 against 0.597, Wilcoxon p < 0.0001), segment directions (0.613 against 0.605 for time segments, p = 0.897; 0.615 against 0.542 for path segments, p = 0.0227) and the kinematic space at its best k = 2 (0.729 against 0.670, p = 0.0004). The segment-directions comparison is close for time segments.
fig8-linear-dynamics-simulation¶
In a simulated two-unit linear dynamical system with 25 random initial states, preparatory activity correlates weakly with peri-movement activity at the single-unit level (r = 0.08), yet a PD in the PCA-based peri-movement space captures preparatory tuning perfectly.
x(t+1) = W x(t), W = [1, -2pi/360; 2pi/360, 1] (Figure 8 legend). Readings: initial states uniform in [-1, 1]^2 or standard Gaussian; the gray interval one full cycle (360 steps) or the first quarter cycle (90 steps) after release. 1000 draws each.
rng8 = np.random.default_rng(8)
eps = 2 * np.pi / 360
W = np.array([[1, -eps], [eps, 1]])
def simulate(x0, steps):
out = np.empty((len(x0), steps, 2)); x = x0.copy()
for t in range(steps):
x = x @ W.T; out[:, t] = x
return out
rows = []
for dist in ["uniform", "gaussian"]:
for steps in [360, 90]:
rs, perf = [], []
for _ in range(1000):
x0 = rng8.uniform(-1, 1, (25, 2)) if dist == "uniform" else rng8.standard_normal((25, 2))
tr8 = simulate(x0, steps)
rs.append(np.corrcoef(x0[:, 0], tr8[:, :, 0].mean(1))[0, 1])
T8 = tr8.reshape(25, -1)
S8 = M.pca_scores(T8 @ T8.T, 2)
A8 = np.column_stack([np.ones(25), S8])
coef, *_ = np.linalg.lstsq(A8, x0[:, 0], rcond=None)
perf.append(1 - ((A8 @ coef - x0[:, 0]) ** 2).mean() / x0[:, 0].var())
rs = np.array(rs)
lo, hi = np.percentile(rs, [2.5, 97.5])
rows.append(dict(initial=dist, interval_steps=steps, median_r=round(np.median(rs), 3), central95=f"{lo:.3f} to {hi:.3f}",
contains_0_08=lo <= 0.08 <= hi, min_perf=np.min(perf), max_perf=np.max(perf)))
f8 = pd.DataFrame(rows)
print(f8.to_string(index=False))
initial interval_steps median_r central95 contains_0_08 min_perf max_perf uniform 360 0.012 -0.397 to 0.387 True 1.0 1.0 uniform 90 0.703 0.496 to 0.847 False 1.0 1.0 gaussian 360 0.000 -0.400 to 0.404 True 1.0 1.0 gaussian 90 0.711 0.455 to 0.871 False 1.0 1.0
Result: reproduced under the full-cycle reading of the gray interval. Over one full cycle (360 steps) the central 95% of the correlation over 1000 draws is -0.397 to 0.387 (uniform initial states) and -0.400 to 0.404 (Gaussian), which contains 0.08; over a quarter cycle it is 0.496 to 0.847 and 0.455 to 0.871, which does not. The 2-D PCA peri-movement space fits unit 1's initial rates with performance 1.0 in every draw.
Summary¶
| Claim | Outcome |
|---|---|
| j-array-recordings | partially reproduced (152 units, but 27 conditions, assumed among the 108, and 250 trials) |
| fig4d-mean-correlation-j-array | partially reproduced, resting on the selection and tuning readings |
| fig4abc-mean-correlation-other-datasets | blocked: data missing |
| fig4-broad-unimodal | reproduced |
| fig4e-correlation-time-course | reproduced, resting on readings; in the main selection on the splice |
| suppfig1-population-psth-weak-movement-tuning | reproduced under reading (a) only; not under cross-validated ranking |
| fig6d-perimovement-best-j-array | reproduced |
| fig6-perimovement-best-each-dimensionality | reproduced |
| fig6-dimensionality-peak | reproduced, resting on S fixed across folds |
| fig6-maze-task-spaces | reproduced, resting on segments equal in time |
| fig6abc-perimovement-best-other-datasets | blocked: data missing |
| fig6e-summary-rank | blocked: data missing |
| fig6-preparatory-dimensionality | not reproduced |
| fig7a-shuffled-space-poor | reproduced |
| fig7b-window-start-slide | reproduced |
| ap-location-no-effect-j-array | reproduced, AP read as the array |
| ap-location-no-effect-other-datasets | blocked: data missing |
| fig7c-anterior-posterior-performance | not reproduced (paper's upper limit) |
| prep-activity-more-anterior | blocked: metadata insufficient (which array is PMd) |
| fig7d-pd-consistency | reproduced |
| some-neurons-better-elsewhere | reproduced |
| high-dimensional-task-spaces-underperform | reproduced |
| fig8-linear-dynamics-simulation | reproduced (full-cycle reading) |
A failure to reproduce here may be a failure of this attempt, not of the paper: see the reasoning and the differences from the paper's method.
How the analysis differs from the paper's. The paper's J-array dataset has 108 conditions and about 2166 trials per neuron; this dandiset has 27 conditions and 6 to 10 kept trials per condition, so condition means are noisy (split-half reliability of preparatory tuning, median r = 0.374) and every PD performance is low (best 0.117 with S fixed, 0.226 recomputed). The tuning strength is read as the range across conditions; the SD reading leaves almost no neuron above 5 spikes/s. The proto-trial is the condition-mean velocity profile with the two profiles concatenated (the reading closest to the paper's ~90% of trials kept); the other readings are shown. S is fixed across folds where the two readings differ. AP location is read as the array, and the array from the unit id. EMG, the speed-task and single-electrode datasets are not in the dandiset. With 6 to 10 trials per condition the 5 spikes/s threshold passes a median of 59 of 152 units on shuffled labels (75 observed); only 43 of the 76 PD neurons have tuning distinguishable from noise, and on them the Figure 6 results are stronger.
What was found about the dataset. The units' electrode references all point to PMd rows of the electrode table (whose location column puts M1 on ids 97 to 192), although the unit ids give 77 units on a second array whose references (rows 14 to 94) look like M1 channels without the offset of 96. The electrode table has no positions. The test file's 100 trials have no condition labels or go cue times. Nothing in the files links this session (2009-09-29) to the paper's J-array dataset or says which id prefix is PMd, and the 27 conditions are a quarter of the paper's 108.