Churchland et al. (2010) on dandiset 000140 (MC_Maze_Small)¶
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.
Claims. claims/000140/10.1016-j.neuron.2010.09.015.yaml, written from the NIH author manuscript
in PubMed Central (PMC2991102) and its supplement (NIHMS240235-supplement-01.pdf), and checked by a
second reader. The claims file decides what is tested; method details it lacks are taken from the
paper and cited where they are used, here and in maze_prep.py.
Dataset. DANDI:000140, version 0.220113.0408 (doi:10.48324/dandi.000140/0.220113.0408), streamed. The paper cites no version of the dandiset. 0.220113.0408 is the only published version; the draft (modified 2023-06-20) has the same two assets and size. A published version is fixed, so it is pinned.
The dandiset is the Neural Latents Benchmark's MC_Maze_Small: one session of monkey J (Jenkins, 2009-09-28) with two 96-electrode arrays, 100 train/val trials over 27 maze conditions, and 100 test segments without conditions. The paper's J-array dataset has 108 conditions and 2166 trials per neuron; the dandiset holds 27 of its conditions (or of a session like it) with 1 to 4 analysable trials each. The single-electrode datasets of monkeys A, B and J and all EMG are not in it.
| Claim | Tested here |
|---|---|
| fig4d-mean-correlation-j-array | tested |
| fig4abc-mean-correlation-other-datasets | blocked: data_missing (single-electrode A, B, J) |
| fig4-broad-near-zero | tested |
| fig4e-correlation-time-course | tested |
| suppfig1-population-averages-lose-tuning | tested |
| fig6d-perimovement-best-j-array | tested |
| fig6abc-perimovement-best-other-datasets | blocked: data_missing (single-electrode data, EMG) |
| fig6-perimovement-best-each-dimensionality | tested |
| fig6-maze-endpoint-speed-space-poor | tested |
| fig6-dimensionality-trend | tested |
| fig6e-summary-rank | blocked: data_missing (A, B, J and EMG); J-array variant shown |
| prep-dimensionality-7-to-10 | tested |
| fig7a-shuffled-space-poor | tested |
| fig7b-window-start-slide | tested (monkey J reading (a), J-array); A and B not in the dandiset |
| ap-location-no-effect-j-array | tested on array identity (assumption) |
| ap-location-no-effect-other-datasets | blocked: data_missing |
| fig7c-anterior-posterior-performance | run on array identity; blocked: other (the normalisation is undefined) |
| prep-activity-more-anterior | run on array identity; blocked: metadata_insufficient (which array is anterior) |
| fig7d-pd-consistency | tested |
| fig8-linear-dynamics-simulation | tested (simulation, no data) |
import os
os.environ.setdefault("OMP_NUM_THREADS", "1")
os.environ.setdefault("OPENBLAS_NUM_THREADS", "1")
import time
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy import stats
import maze_prep as mp
T_START = time.time()
pd.set_option("display.width", 200)
pd.set_option("display.max_columns", 40)
urls = mp.asset_urls()
print("dandiset", mp.DANDISET_ID, "version", mp.DANDISET_VERSION)
for p, u in urls.items():
print(" ", p, u)
data = mp.load_train(urls[mp.TRAIN_PATH])
test = mp.load_test(urls[mp.TEST_PATH])
trials, units = data["trials"], data["units"]
print(f"loaded in {time.time() - T_START:.1f} s")
dandiset 000140 version 0.220113.0408 sub-Jenkins/sub-Jenkins_ses-small_desc-test_ecephys.nwb https://dandiarchive.s3.amazonaws.com/blobs/118/a0f/118a0fea-39e2-4e85-b01d-a72df466f266 sub-Jenkins/sub-Jenkins_ses-small_desc-train_behavior+ecephys.nwb https://dandiarchive.s3.amazonaws.com/blobs/7f4/e60/7f4e60f3-2184-4e68-ace1-504a2e564939
loaded in 18.9 s
Data check¶
What was loaded, and each check of the description against the values.
meta = data["meta"]
print("subject", meta["subject"], "session start", meta["session_start"])
print("train file trials:", len(trials), "splits:", trials.split.value_counts().to_dict())
print("test file segments:", len(test["trials"]), "columns:", test["colnames"], "behaviour:", test["has_behavior"])
print("test segment length (s):", np.unique(np.round(test["trials"].stop_time - test["trials"].start_time, 3)),
"move onset re start (s):", np.unique(np.round(test["trials"].move_onset_time - test["trials"].start_time, 3)))
print("all train trials successful:", bool(trials.success.all()))
print("conditions (maze_id x trial_version):", trials.groupby(["maze_id", "trial_version"]).ngroups,
"| mazes:", trials.maze_id.nunique(), "| trial_type values:", trials.trial_type.nunique())
print("trial_version vs targets/barriers:")
print(pd.crosstab(trials.trial_version, [trials.num_targets, trials.num_barriers]))
print("delay (ms): min %d median %d max %d; trials with delay > 400 ms: %d"
% (trials.delay.min(), trials.delay.median(), trials.delay.max(), (trials.delay > mp.MIN_DELAY_MS).sum()))
kept_delay = trials[trials.delay > mp.MIN_DELAY_MS]
per_cond = kept_delay.groupby("cond").size()
print("trials per condition after the delay rule: min %d, median %.0f, max %d; conditions with 1 trial: %d"
% (per_cond.min(), per_cond.median(), per_cond.max(), (per_cond == 1).sum()))
print("delay column = go cue - target onset (max abs diff, ms):",
float(np.abs((trials.go_cue_time - trials.target_on_time) * 1000 - trials.delay).max()))
print("rt column = move onset - go cue (max abs diff, ms):",
float(np.abs((trials.move_onset_time - trials.go_cue_time) * 1000 - trials.rt).max()))
print("rt (ms): min %d median %d max %d" % (trials.rt.min(), trials.rt.median(), trials.rt.max()))
print("target onset - trial start (s): %.3f to %.3f; trial stop - move onset (s): %.3f to %.3f"
% ((trials.target_on_time - trials.start_time).min(), (trials.target_on_time - trials.start_time).max(),
(trials.stop_time - trials.move_onset_time).min(), (trials.stop_time - trials.move_onset_time).max()))
print("so go cue - 300 ms and move onset + 900 ms lie inside every trial:",
bool(((trials.go_cue_time - 0.3) > trials.start_time).all() and ((trials.move_onset_time + 0.9) < trials.stop_time).all()))
print()
print("units:", len(units), "| heldout:", int(units.heldout.sum()))
print("unit id prefix (array) counts:", units.array_prefix.value_counts().to_dict())
print("electrode group that units/electrodes points to:", units.electrode_group.value_counts().to_dict())
print("electrode index == channel - 1 for every unit:", bool((units.electrode == units.channel - 1).all()))
print("electrodes table groups:", meta["electrode_groups"], "| x, y, z all NaN:", meta["electrode_xyz_all_nan"])
oi = meta["obs_intervals"].reshape(len(units), len(trials), 2)
print("obs_intervals of every unit equal the trial intervals:",
bool(np.allclose(oi, np.c_[trials.start_time, trials.stop_time][None])))
n_spk = np.array([len(s) for s in data["spikes"]])
dur = (trials.stop_time - trials.start_time).sum()
print("spikes per unit: min %d median %d max %d; mean rate over trials %.2f to %.2f spikes/s"
% (n_spk.min(), np.median(n_spk), n_spk.max(), n_spk.min() / dur, n_spk.max() / dur))
print("spike times sorted in every unit:", all(np.all(np.diff(s) >= 0) for s in data["spikes"]))
inside = np.zeros(len(n_spk), int)
for u, s in enumerate(data["spikes"]):
inside[u] = sum(((s >= a) & (s <= b)).sum() for a, b in zip(trials.start_time, trials.stop_time))
print("spikes outside the trial intervals:", int((n_spk - inside).sum()))
beh = data["beh"]
print()
print("hand_pos:", meta["pos_description"], "| stored unit", meta["pos_unit"], "conversion", meta["pos_conversion"])
print("hand_vel:", meta["vel_description"], "| stored unit", meta["vel_unit"], "conversion", meta["vel_conversion"])
print("hand position range x %.1f to %.1f, y %.1f to %.1f (values, mm)" % (beh["pos"][:, 0].min(), beh["pos"][:, 0].max(),
beh["pos"][:, 1].min(), beh["pos"][:, 1].max()))
dt = np.diff(beh["t"])
print("behaviour sampling: %.0f%% of intervals are 1 ms; gaps (between trials) up to %.3f s" % (100 * np.mean(np.isclose(dt, 0.001)), dt.max()))
# velocity is the derivative of position (description); check the scale on the 1 ms samples
dpos = np.diff(beh["pos"], axis=0) / dt[:, None]
ok = np.isclose(dt, 0.001)
slope = np.polyfit(beh["vel"][1:][ok, 0], dpos[ok, 0], 1)[0]
print("d(position)/dt against hand_vel, slope %.3f (1 if both are mm and mm/s)" % slope)
sp = np.hypot(beh["vel"][:, 0], beh["vel"][:, 1])
print("hand speed 99.9th percentile %.0f mm/s" % np.percentile(sp, 99.9))
print("target positions (mm) x %d to %d, y %d to %d" % (trials.target_x.min(), trials.target_x.max(), trials.target_y.min(), trials.target_y.max()))
subject Jenkins session start 2009-09-28T00:00:00-07:00
train file trials: 100 splits: {'train': 75, 'val': 25}
test file segments: 100 columns: ['start_time', 'stop_time', 'move_onset_time', 'split'] behaviour: False
test segment length (s): [0.7] move onset re start (s): [0.25]
all train trials successful: True
conditions (maze_id x trial_version): 27 | mazes: 9 | trial_type values: 9
trial_version vs targets/barriers:
num_targets 1 3
num_barriers 0 9 9
trial_version
0 32 0 0
1 0 34 0
2 0 0 34
delay (ms): min 14 median 615 max 999; trials with delay > 400 ms: 78
trials per condition after the delay rule: min 1, median 3, max 4; conditions with 1 trial: 2
delay column = go cue - target onset (max abs diff, ms): 4.774847184307873e-11
rt column = move onset - go cue (max abs diff, ms): 5.1386450650170445e-11
rt (ms): min 241 median 321 max 758
target onset - trial start (s): 0.659 to 0.952; trial stop - move onset (s): 0.979 to 1.626
so go cue - 300 ms and move onset + 900 ms lie inside every trial: True
units: 142 | heldout: 35
unit id prefix (array) counts: {1: 72, 2: 70}
electrode group that units/electrodes points to: {'electrode_group_PMd': 142}
electrode index == channel - 1 for every unit: True
electrodes table groups: {'electrode_group_PMd': 96, 'electrode_group_M1': 96} | x, y, z all NaN: True
obs_intervals of every unit equal the trial intervals: True
spikes per unit: min 40 median 654 max 3660; mean rate over trials 0.14 to 12.72 spikes/s
spike times sorted in every unit: True
spikes outside the trial intervals: 0
hand_pos: Monkey hand position in mm. First column x, second column y | stored unit meters conversion 0.001
hand_vel: Monkey hand velocity in mm/s, calculated from monkey hand position. First column x, second column y | stored unit m/s conversion 0.001
hand position range x -135.0 to 153.3, y -159.8 to 98.0 (values, mm)
behaviour sampling: 100% of intervals are 1 ms; gaps (between trials) up to 0.110 s
d(position)/dt against hand_vel, slope 1.000 (1 if both are mm and mm/s)
hand speed 99.9th percentile 1057 mm/s
target positions (mm) x -118 to 133, y -83 to 82
The positions are stored as numbers in mm with unit meters and conversion 0.001, which is
consistent; the derivative of position matches hand_vel with slope about 1, so velocity is in mm/s.
Every unit's electrodes entry points to the PMd group, but the unit ids carry two array prefixes
(1 and 2) and the electrode index equals the channel number minus one for both, so units of the
second array point to PMd electrodes as well. The electrodes table has no positions. Both are findings.
Clock check¶
Spikes and the trial table: decode the side of the reach target (target x > 0) from population spike counts in 100 ms windows around target onset, leave-one-trial-out, nearest class mean on z-scored counts. Before target onset the target is unknown and decoding should be at chance; after it, above chance. The same is shown with the labels of the next trial (a shifted label sequence). Hand and trial table: hand speed around the stored movement onset.
rng_clock = np.random.default_rng(0)
def window_counts(t_ref, w0, w1):
X = np.zeros((len(trials), len(units)))
for u, s in enumerate(data["spikes"]):
lo = np.searchsorted(s, t_ref + w0)
hi = np.searchsorted(s, t_ref + w1)
X[:, u] = hi - lo
return X
def loo_centroid_accuracy(X, y):
X = (X - X.mean(0)) / (X.std(0) + 1e-9)
hits = 0
for i in range(len(y)):
tr_ = np.arange(len(y)) != i
m1 = X[tr_ & (y == 1)].mean(0); m0 = X[tr_ & (y == 0)].mean(0)
pred = int(np.sum((X[i] - m1) ** 2) < np.sum((X[i] - m0) ** 2))
hits += pred == y[i]
return hits / len(y)
y_side = (trials.target_x.values > 0).astype(int)
y_next = np.r_[y_side[1:], y_side[:1]]
rows = []
for w0 in np.arange(-0.4, 0.6, 0.1):
X = window_counts(trials.target_on_time.values, w0, w0 + 0.1)
rows.append({"window re target onset (ms)": f"{int(round(w0*1000))} to {int(round((w0+0.1)*1000))}",
"accuracy, true labels": loo_centroid_accuracy(X, y_side),
"accuracy, next trial's labels": loo_centroid_accuracy(X, y_next)})
clock = pd.DataFrame(rows)
null = [loo_centroid_accuracy(window_counts(trials.target_on_time.values, 0.2, 0.3), rng_clock.permutation(y_side)) for _ in range(50)]
print(clock.round(2).to_string(index=False))
print("chance: majority class %.2f; permuted labels at 200 to 300 ms, 95th percentile %.2f" % (max(y_side.mean(), 1 - y_side.mean()), np.percentile(null, 95)))
speed_at = lambda lag: np.array([np.hypot(*beh["vel"][np.searchsorted(beh["t"], t + lag)]) for t in trials.move_onset_time])
s_before, s_on, s_after = speed_at(-0.1), speed_at(0.0), speed_at(0.1)
print("hand speed (mm/s) median at move onset -100 ms %.1f, at onset %.1f, +100 ms %.1f; trials faster at +100 than -100: %d of %d"
% (np.median(s_before), np.median(s_on), np.median(s_after), (s_after > s_before).sum(), len(trials)))
window re target onset (ms) accuracy, true labels accuracy, next trial's labels
-400 to -300 0.54 0.52
-300 to -200 0.48 0.57
-200 to -100 0.49 0.58
-100 to 0 0.48 0.57
0 to 100 0.68 0.26
100 to 200 0.79 0.50
200 to 300 0.78 0.45
300 to 400 0.81 0.51
400 to 500 0.73 0.57
500 to 600 0.79 0.45
chance: majority class 0.53; permuted labels at 200 to 300 ms, 95th percentile 0.63
hand speed (mm/s) median at move onset -100 ms 2.6, at onset 47.0, +100 ms 744.2; trials faster at +100 than -100: 100 of 100
Duplicated recordings¶
The test file holds 0.7 s segments with movement onset at 250 ms. Are any of them copies of train trials? Each test segment's spike counts of the held-in units are compared with the counts of every train trial in the window 250 ms before to 450 ms after its movement onset. Train trials are also compared with each other.
held = units.index[~units.heldout].values
assert np.array_equal(units.id.values[held], test["ids"])
tt = test["trials"]
Xtest = np.zeros((len(tt), len(held)))
for k, u in enumerate(held):
s = test["spikes"][k]
Xtest[:, k] = np.searchsorted(s, tt.stop_time.values) - np.searchsorted(s, tt.start_time.values)
Xtrain = window_counts(trials.move_onset_time.values, -0.25, 0.45)[:, held]
same_test_train = sum(np.any(np.all(Xtrain == x, axis=1)) for x in Xtest)
dist = np.abs(Xtest[:, None, :] - Xtrain[None]).sum(-1)
Xfull = window_counts(trials.start_time.values, 0, 2.0)
same_train = sum(np.all(Xfull[i] == Xfull[j]) for i in range(len(trials)) for j in range(i + 1, len(trials)))
print("test segments identical to a train trial window:", same_test_train, "of", len(tt))
print("smallest summed count difference between a test segment and any train window: %d (median over segments %d)"
% (dist.min(), np.median(dist.min(1))))
print("pairs of train trials with identical spike counts in their first 2 s:", same_train)
print("total spikes, test file %d, train file %d" % (sum(len(s) for s in test["spikes"]), n_spk.sum()))
test segments identical to a train trial window: 0 of 100 smallest summed count difference between a test segment and any train window: 117 (median over segments 156) pairs of train trials with identical spike counts in their first 2 s: 0 total spikes, test file 32897, train file 131669
Shared computations¶
Condition means under each reading of the open details the claims file names: movement onset (O3a,
the dataset's move_onset_time; O3b, hand speed reaching 15% of its trial peak) and the proto-trial
filter (O4a, b, c). Only trials with delay > 400 ms are used. A condition is a maze_id and
trial_version pair (27). Spike trains are smoothed with a 20 ms Gaussian before averaging
(Supplemental Experimental Procedures). The 10 ms time step of the T matrices is O5.
t0 = time.time()
keep = np.flatnonzero(trials.delay.values > mp.MIN_DELAY_MS)
conds = trials.cond.values[keep]
C = len(data["conditions"])
GO_L = np.arange(mp.GO_SPAN[0], mp.GO_SPAN[1] + 1)
MV_L = np.arange(mp.MV_SPAN[0], mp.MV_SPAN[1] + 1)
onsets = {"a": trials.move_onset_time.values,
"b": np.array([mp.speed_onset(beh, trials.iloc[i]) for i in range(len(trials))])}
d_on = (onsets["b"] - onsets["a"])[keep] * 1000
print("O3b minus O3a onset (ms): min %.0f median %.0f max %.0f" % (d_on.min(), np.median(d_on), d_on.max()))
RATES, KINS = {}, {}
for o in "ab":
RATES[o] = mp.trial_rates(data, onsets[o], keep)
KINS[o] = [mp.trial_kinematics(beh, onsets[o][i]) for i in keep]
TSEL = (MV_L >= mp.T_WIN[0]) & (MV_L <= mp.T_WIN[1]) & ((MV_L - mp.T_WIN[0]) % mp.STEP_MS == 0)
CONFIGS = {}
for o in "ab":
vp = np.stack([k["vel"][(k["lags_ms"] >= mp.T_WIN[0]) & (k["lags_ms"] <= mp.T_WIN[1])] for k in KINS[o]])
for f in "abc":
kp = mp.proto_filter(vp, conds, f)
Rg, Rm = RATES[o]
cfg = {"kept": kp,
"CMg": mp.condition_means(Rg[kp], conds[kp], C),
"CMm": mp.condition_means(Rm[kp], conds[kp], C)}
cfg["P"] = mp.epoch_mean(cfg["CMg"], GO_L, mp.PREP_WIN)
cfg["M"] = mp.epoch_mean(cfg["CMm"], MV_L, mp.PERI_WIN)
kins = [KINS[o][i] for i in np.flatnonzero(kp)]
cfg["KIN"] = {k: mp.condition_means(np.stack([x[k] for x in kins]), conds[kp], C) for k in ("pos", "vel", "acc")}
cfg["KIN"]["lags_ms"] = kins[0]["lags_ms"]
for seg in ("length", "time"):
feats = [mp.task_features(x, seg) for x in kins]
cfg["TS_" + seg] = mp.task_space_matrices(feats, trials.iloc[keep[kp]], conds[kp], C)
cfg["reach_end"] = np.array([mp.task_features(x)["reach_end_ms"][0] for x in kins])
CONFIGS[(o, f)] = cfg
rows = []
for (o, f), cfg in CONFIGS.items():
n = np.bincount(conds[cfg["kept"]], minlength=C)
rows.append({"onset": o, "filter": f, "trials kept": int(cfg["kept"].sum()), "of": len(keep),
"trials per condition min": n.min(), "median": np.median(n), "max": n.max()})
print(pd.DataFrame(rows).to_string(index=False))
re_ = CONFIGS[("a", "a")]["reach_end"]
print("reach end (ms after onset, 10%% of peak speed): min %d median %d max %d; at the 800 ms limit: %d trials"
% (re_.min(), np.median(re_), re_.max(), (re_ >= 800).sum()))
print(f"shared computations {time.time() - t0:.1f} s")
MAIN = ("a", "a")
O3b minus O3a onset (ms): min 8 median 25 max 42
onset filter trials kept of trials per condition min median max
a a 66 78 1 2.0 4
a b 71 78 1 3.0 4
a c 78 78 1 3.0 4
b a 66 78 1 2.0 4
b b 73 78 1 3.0 4
b c 78 78 1 3.0 4
reach end (ms after onset, 10% of peak speed): min 236 median 382 max 800; at the 800 ms limit: 1 trials
shared computations 3.0 s
Reliability of condition means¶
The paper's J-array dataset averaged about 20 trials per condition (2166 trials per neuron over 108 conditions); here there are 1 to 4. As a reference for every claim below, the split-half reliability of each unit's preparatory and peri-movement tuning: over the conditions with at least two trials, the correlation between the condition means of one trial (the first) and of the others.
Rg, Rm = RATES["a"]
Pt = mp.epoch_mean(Rg, GO_L, mp.PREP_WIN) # trials x units
Mt = mp.epoch_mean(Rm, MV_L, mp.PERI_WIN)
multi = [c for c in range(C) if (conds == c).sum() >= 2]
first = np.array([np.flatnonzero(conds == c)[0] for c in multi])
rest_P = np.stack([Pt[np.flatnonzero(conds == c)[1:]].mean(0) for c in multi])
rest_M = np.stack([Mt[np.flatnonzero(conds == c)[1:]].mean(0) for c in multi])
rel_P = mp.rowwise_corr(Pt[first], rest_P)
rel_M = mp.rowwise_corr(Mt[first], rest_M)
pt_main = mp.tuning(CONFIGS[MAIN]["P"], "range")
print("conditions with >= 2 trials:", len(multi))
for name, r in (("preparatory", rel_P), ("peri-movement", rel_M)):
r = r[np.isfinite(r)]
print(f"{name} split-half r over units: median {np.median(r):.2f}, IQR {np.percentile(r,25):.2f} to {np.percentile(r,75):.2f}")
sel5 = pt_main > mp.PREP_TUNING_MIN
print("same, units with > 5 spikes/s preparatory range: preparatory median %.2f, peri-movement median %.2f"
% (np.nanmedian(rel_P[sel5]), np.nanmedian(rel_M[sel5])))
sel_sd = mp.tuning(CONFIGS[MAIN]["P"], "sd") > mp.PREP_TUNING_MIN
print("same, units with > 5 spikes/s preparatory SD (n = %d): preparatory median %.2f (values %s), peri-movement median %.2f"
% (sel_sd.sum(), np.nanmedian(rel_P[sel_sd]), np.round(np.sort(rel_P[sel_sd]), 2), np.nanmedian(rel_M[sel_sd])))
conditions with >= 2 trials: 25 preparatory split-half r over units: median 0.20, IQR -0.03 to 0.40 peri-movement split-half r over units: median 0.28, IQR 0.07 to 0.53 same, units with > 5 spikes/s preparatory range: preparatory median 0.26, peri-movement median 0.34 same, units with > 5 spikes/s preparatory SD (n = 8): preparatory median 0.57 (values [0.03 0.11 0.38 0.53 0.62 0.63 0.7 0.99]), peri-movement median 0.49
/home/runner/work/encore/encore/entries/000140-churchland-2010/maze_prep.py:270: RuntimeWarning: invalid value encountered in divide return (a * b).sum(0) / np.sqrt((a ** 2).sum(0) * (b ** 2).sum(0))
Readings adopted where readings differ (author, step 1.6)¶
Every claim below is run under every reading the claims file names. Where they give different outcomes, these are adopted as assumptions, and the outcome is marked as resting on them:
- O1, "spikes/s of preparatory tuning" (Methods, Neural recordings and datasets): the range of the
preparatory rate across conditions (O1a). The Methods speak of a rate difference ("Firing rate
modulation was typically lower for the array recordings ... We thus used a lower inclusion criterion:
5 spikes/s of preparatory tuning"), which the range measures directly. The SD reading (O1b) is shown everywhere.
- O3, movement onset: the dataset's
move_onset_time(O3a), computed by the data's authors ("calculated offline using a robust algorithm", its description). The paper defers to Churchland et al. (2006c). - O4, the proto-trial: a trial of the condition, the one best correlated with the condition mean (O4a), since the Supplemental Experimental Procedures call it a "reference 'proto-trial'".
- Units in the peri-movement T matrix ("all neurons except the neuron whose preparatory activity we were trying to capture", Methods): all recorded units (main) and, as a second reading, the units passing the inclusion rule; both are run in fig6d. Under O1a they give the same outcome.
fig4d-mean-correlation-j-array¶
For the J-array dataset, the correlation across conditions between a unit's preparatory and peri-movement firing rates is positive on average but small (mean Pearson r = 0.13, p = 0.03; mean Spearman r = 0.11).
Criterion: units with > 5 spikes/s of preparatory tuning and peri-movement tuning at least 1.5 times as strong, under O1 (range, SD) and O2 (i epoch means, ii largest tuning at any time); O3 and O4 as above; the mean against zero by (a) one-sample t-test, (b) t-test on Fisher z, (c) Wilcoxon. Expected mean Pearson 0.13 ± 0.013 and Spearman 0.11 ± 0.011.
def fig4_units(cfg, O1, O2):
pt = mp.tuning(cfg["P"], O1)
qt = mp.peri_tuning(cfg["CMm"], MV_L, cfg["M"], O1, O2)
return (pt > mp.PREP_TUNING_MIN) & (qt >= mp.PERI_RATIO_MIN * pt)
def fig4_corrs(cfg, sel):
r = mp.rowwise_corr(cfg["P"][:, sel], cfg["M"][:, sel])
rs = np.array([stats.spearmanr(cfg["P"][:, j], cfg["M"][:, j])[0] for j in np.flatnonzero(sel)])
return r, rs
rows = []
for (o, f), cfg in CONFIGS.items():
for O1 in ("range", "sd"):
for O2 in ("i", "ii"):
sel = fig4_units(cfg, O1, O2)
row = {"O3": o, "O4": f, "O1": O1, "O2": O2, "n units": int(sel.sum())}
if sel.sum() >= 3:
r, rs = fig4_corrs(cfg, sel)
row.update({"mean r": r.mean(), "mean Spearman": rs.mean(),
"p t": stats.ttest_1samp(r, 0).pvalue,
"p t(z)": stats.ttest_1samp(np.arctanh(np.clip(r, -0.999, 0.999)), 0).pvalue,
"p Wilcoxon": stats.wilcoxon(r).pvalue})
row["within tol"] = abs(r.mean() - 0.13) <= 0.013 and abs(rs.mean() - 0.11) <= 0.011
row["positive, p<0.05 (any test)"] = r.mean() > 0 and min(row["p t"], row["p t(z)"], row["p Wilcoxon"]) < 0.05
rows.append(row)
fig4_table = pd.DataFrame(rows)
print(fig4_table.round(3).to_string(index=False))
print("units selected under O1 range / O1 SD (preparatory > 5 spikes/s, main O3/O4):",
int((mp.tuning(CONFIGS[MAIN]["P"], "range") > 5).sum()), "/", int((mp.tuning(CONFIGS[MAIN]["P"], "sd") > 5).sum()), "of", len(units))
ok = fig4_table.dropna(subset=["mean r"])
print("readings with >= 3 units:", len(ok), "| within tolerance and p < 0.05:",
int((ok["within tol"] & ok["positive, p<0.05 (any test)"]).sum()),
"| positive at p < 0.05:", int(ok["positive, p<0.05 (any test)"].sum()))
print("mean r over readings: min %.3f max %.3f; mean Spearman min %.3f max %.3f"
% (ok["mean r"].min(), ok["mean r"].max(), ok["mean Spearman"].min(), ok["mean Spearman"].max()))
sel_main = fig4_units(CONFIGS[MAIN], "range", "i")
r_main, rs_main = fig4_corrs(CONFIGS[MAIN], sel_main)
print("main reading (O1 range, O2 i, O3a, O4a): n = %d, mean r = %.3f, mean Spearman = %.3f, conditions = %d"
% (sel_main.sum(), r_main.mean(), rs_main.mean(), C))
for arr in ("PMd", "M1"):
m = (units.array.values == arr)[sel_main]
print(f" {arr} (by unit id prefix): n = {m.sum()}, mean r = {r_main[m].mean():.3f}")
# the units of the widest selection (O2 ii) for reference
sel_ii = fig4_units(CONFIGS[MAIN], "range", "ii")
r_ii, rs_ii = fig4_corrs(CONFIGS[MAIN], sel_ii)
print("O2 ii (main otherwise): n = %d, mean r = %.3f (t-test p = %.3f), mean Spearman %.3f" % (sel_ii.sum(), r_ii.mean(), stats.ttest_1samp(r_ii, 0).pvalue, rs_ii.mean()))
O3 O4 O1 O2 n units mean r mean Spearman p t p t(z) p Wilcoxon within tol positive, p<0.05 (any test) a a range i 19 -0.079 -0.020 0.086 0.085 0.169 False False a a range ii 88 0.011 0.043 0.662 0.587 0.973 False False a a sd i 0 NaN NaN NaN NaN NaN NaN NaN a a sd ii 5 0.095 -0.053 0.607 0.582 0.625 False False a b range i 13 -0.148 -0.053 0.023 0.026 0.017 False False a b range ii 85 0.013 0.056 0.606 0.542 0.842 False False a b sd i 0 NaN NaN NaN NaN NaN NaN NaN a b sd ii 4 -0.024 -0.044 0.901 0.909 1.000 False False a c range i 15 -0.106 -0.032 0.078 0.077 0.169 False False a c range ii 85 0.035 0.061 0.199 0.189 0.265 False False a c sd i 0 NaN NaN NaN NaN NaN NaN NaN a c sd ii 4 -0.029 -0.020 0.891 0.900 1.000 False False b a range i 18 -0.041 -0.000 0.379 0.366 0.442 False False b a range ii 89 -0.003 0.025 0.897 0.993 0.554 False False b a sd i 0 NaN NaN NaN NaN NaN NaN NaN b a sd ii 5 0.060 -0.061 0.745 0.709 0.812 False False b b range i 15 -0.072 -0.013 0.162 0.162 0.389 False False b b range ii 86 0.006 0.041 0.815 0.740 0.838 False False b b sd i 0 NaN NaN NaN NaN NaN NaN NaN b b sd ii 4 -0.045 -0.053 0.803 0.806 1.000 False False b c range i 15 -0.063 -0.003 0.219 0.212 0.489 False False b c range ii 83 0.011 0.035 0.686 0.667 0.608 False False b c sd i 0 NaN NaN NaN NaN NaN NaN NaN b c sd ii 3 0.069 0.061 0.753 0.743 0.750 False False units selected under O1 range / O1 SD (preparatory > 5 spikes/s, main O3/O4): 93 / 8 of 142 readings with >= 3 units: 18 | within tolerance and p < 0.05: 0 | positive at p < 0.05: 0 mean r over readings: min -0.148 max 0.095; mean Spearman min -0.061 max 0.061 main reading (O1 range, O2 i, O3a, O4a): n = 19, mean r = -0.079, mean Spearman = -0.020, conditions = 27 PMd (by unit id prefix): n = 10, mean r = -0.060 M1 (by unit id prefix): n = 9, mean r = -0.099 O2 ii (main otherwise): n = 88, mean r = 0.011 (t-test p = 0.662), mean Spearman 0.043
Noise reference for fig4d, fig4-broad-near-zero, fig4e and suppfig1¶
With 1 to 4 trials per condition, the preparatory and peri-movement means of a condition share the noise of the same few trials, and the preparatory window ends 100 ms after the go cue, about 220 ms before movement onset (median reaction time 321 ms). Two references are computed on the main reading (O3a, O4a, O1a):
- condition-shuffled null: the condition labels of the kept trials are permuted (counts kept), and the whole pipeline (selection, correlations, time course, population averages) is run on the shuffled means, 20 times. There is no tuning in these data, only trial noise.
- disjoint trials: within each condition with at least two kept trials, alternate trials go to two sets; the preparatory means come from one set and the movement-locked means from the other, both ways. Single-trial conditions are dropped. Units are those of fig4d on all trials. Shared trial noise cannot correlate the two means here.
The time-course and population-average functions used by fig4e and suppfig1 are defined here.
LAGS_MV = np.arange(-400, 401, 10)
LAGS_GO = np.arange(-200, 201, 10)
def time_course(cfg, sel):
P = cfg["P"][:, sel]
mv = np.array([np.nanmean(mp.rowwise_corr(P, cfg["CMm"][:, sel, np.flatnonzero(MV_L == t)[0]])) for t in LAGS_MV])
go = np.array([np.nanmean(mp.rowwise_corr(P, cfg["CMg"][:, sel, np.flatnonzero(GO_L == t)[0]])) for t in LAGS_GO])
return go, mv
def judge_tc(mv):
at0 = mv[LAGS_MV == 0][0]; at200 = mv[LAGS_MV == -200][0]
drops = mv[:-5] - mv[5:] # drop over 50 ms starting at each lag
start = LAGS_MV[:-5][np.argmax(drops)]
return {"r at -200": at200, "r at 0": at0, "largest 50 ms drop starts": start,
"pass": at0 < 0.1 and at0 < 0.5 * at200 and -300 <= start <= 0,
"partial": at0 >= 0.1 and at0 < at200}
def pop_avg(cfg, sel, normalise=True):
P = cfg["P"][:, sel]
n_c = P.shape[0]
order = np.argsort(P, axis=0)
ranks = [order[-1], order[int(round((n_c - 1) * 2 / 3))], order[int(round((n_c - 1) / 3))], order[0]]
idx = np.flatnonzero(sel)
go_tr, mv_tr = [], []
for rk in ranks:
g = np.stack([cfg["CMg"][c, u] for c, u in zip(rk, idx)])
m = np.stack([cfg["CMm"][c, u] for c, u in zip(rk, idx)])
if normalise:
mx = np.maximum(cfg["CMg"][:, idx].max((0, 2)), cfg["CMm"][:, idx].max((0, 2)))
g, m = g / mx[:, None], m / mx[:, None]
go_tr.append(g.mean(0)); mv_tr.append(m.mean(0))
return np.array(go_tr), np.array(mv_tr)
def judge_pop(go_tr, mv_tr):
prep = (GO_L >= mp.PREP_WIN[0]) & (GO_L <= mp.PREP_WIN[1])
spread_prep = (go_tr[0] - go_tr[3])[prep].mean()
spread_0 = (mv_tr[0] - mv_tr[3])[MV_L == 0][0]
top150 = mv_tr[:, MV_L == -150][:, 0]
return {"spread prep": spread_prep, "spread at onset": spread_0, "ratio": spread_0 / spread_prep,
"strongest on top at -150": bool(np.argmax(top150) == 0),
"pass": spread_0 < 0.5 * spread_prep and np.argmax(top150) == 0}
def means_from(o, idx_prep, idx_peri, labels, cond_list):
"""Condition means over `cond_list`, go-locked from trials idx_prep and movement-locked from
idx_peri (positions in `keep`), with the given condition label of every kept trial."""
Rg, Rm = RATES[o]
cm = lambda R, idx: np.stack([R[idx[labels[idx] == c]].mean(0) for c in cond_list])
cfg = {"CMg": cm(Rg, idx_prep), "CMm": cm(Rm, idx_peri)}
cfg["P"] = mp.epoch_mean(cfg["CMg"], GO_L, mp.PREP_WIN)
cfg["M"] = mp.epoch_mean(cfg["CMm"], MV_L, mp.PERI_WIN)
return cfg
def noise_row(cfg, sel=None):
sel_b = mp.tuning(cfg["P"], "range") > mp.PREP_TUNING_MIN
sel = fig4_units(cfg, "range", "i") if sel is None else sel
r, _ = fig4_corrs(cfg, sel)
q1, q3 = np.nanpercentile(r, [25, 75])
ja, jb = judge_tc(time_course(cfg, sel)[1]), judge_tc(time_course(cfg, sel_b)[1])
pa = judge_pop(*pop_avg(cfg, sel, True))
return {"n > 5 spikes/s": int(sel_b.sum()), "n fig4": int(sel.sum()), "mean r": np.nanmean(r), "IQR r": q3 - q1,
"fig4e(a) r -200": ja["r at -200"], "r 0": ja["r at 0"], "drop": ja["largest 50 ms drop starts"], "pass (a)": ja["pass"],
"fig4e(b) r -200": jb["r at -200"], "r 0 (b)": jb["r at 0"], "pass (b)": jb["pass"],
"suppfig1 ratio": pa["ratio"], "top at -150": pa["strongest on top at -150"], "suppfig1 pass": pa["pass"]}
rng_null = np.random.default_rng(7)
kp_main = np.flatnonzero(CONFIGS[MAIN]["kept"])
rows = [{"set": "data (main reading)", **noise_row(CONFIGS[MAIN])}]
for s in range(20):
lab = conds.copy(); lab[kp_main] = rng_null.permutation(conds[kp_main])
rows.append({"set": f"shuffle {s + 1}", **noise_row(means_from("a", kp_main, kp_main, lab, range(C)))})
null_table = pd.DataFrame(rows)
print("condition-shuffled null, main reading")
print(null_table.round(3).to_string(index=False))
nt = null_table.iloc[1:]
print("over 20 shuffles: units passing > 5 spikes/s range %d to %d (data %d of %d); fig4e passes (a) %d, (b) %d; suppfig1 passes %d;"
% (nt["n > 5 spikes/s"].min(), nt["n > 5 spikes/s"].max(), null_table["n > 5 spikes/s"].iloc[0], len(units),
nt["pass (a)"].sum(), nt["pass (b)"].sum(), nt["suppfig1 pass"].sum()))
print(" mean r %.3f to %.3f; IQR of r %.3f to %.3f; fig4e(b) r at -200 %.2f to %.2f and at onset %.3f to %.3f"
% (nt["mean r"].min(), nt["mean r"].max(), nt["IQR r"].min(), nt["IQR r"].max(),
nt["fig4e(b) r -200"].min(), nt["fig4e(b) r -200"].max(), nt["r 0 (b)"].min(), nt["r 0 (b)"].max()))
rows = []
for key, cfg in CONFIGS.items():
kp = np.flatnonzero(cfg["kept"])
A, B, both = [], [], []
for c in range(C):
idx = kp[conds[kp] == c]
if len(idx) >= 2:
A += list(idx[0::2]); B += list(idx[1::2]); both.append(c)
A, B = np.array(A), np.array(B)
sel = fig4_units(cfg, "range", "i")[None]
for name, (X, Y) in (("prep A, move B", (A, B)), ("prep B, move A", (B, A))):
d = means_from(key[0], X, Y, conds, both)
rows.append({"O3": key[0], "O4": key[1], "split": name, "conditions": len(both), "trials": f"{len(X)}/{len(Y)}",
**noise_row(d, sel[0])})
disj_table = pd.DataFrame(rows).drop(columns=["n > 5 spikes/s"])
print()
print("disjoint trials for the preparatory and movement-locked means (units of fig4d, O1 range, O2 i)")
print(disj_table.round(3).to_string(index=False))
print("disjoint readings: fig4e passes (a) %d, (b) %d of %d; suppfig1 passes %d; fig4e(b) r at -200 %.2f to %.2f"
% (disj_table["pass (a)"].sum(), disj_table["pass (b)"].sum(), len(disj_table), disj_table["suppfig1 pass"].sum(),
disj_table["fig4e(b) r -200"].min(), disj_table["fig4e(b) r -200"].max()))
/home/runner/work/encore/encore/entries/000140-churchland-2010/maze_prep.py:270: RuntimeWarning: invalid value encountered in divide return (a * b).sum(0) / np.sqrt((a ** 2).sum(0) * (b ** 2).sum(0))
condition-shuffled null, main reading
set n > 5 spikes/s n fig4 mean r IQR r fig4e(a) r -200 r 0 drop pass (a) fig4e(b) r -200 r 0 (b) pass (b) suppfig1 ratio top at -150 suppfig1 pass
data (main reading) 93 19 -0.079 0.201 0.295 -0.078 -250 True 0.388 -0.001 True 0.186 True True
shuffle 1 92 10 0.026 0.239 0.188 -0.062 -250 True 0.334 -0.028 True 0.393 False False
shuffle 2 95 13 -0.002 0.101 0.319 -0.123 120 False 0.280 0.001 True -0.355 True True
shuffle 3 82 9 -0.081 0.169 0.214 -0.009 -130 True 0.314 -0.001 True 0.193 True True
shuffle 4 89 13 0.038 0.330 0.208 -0.003 -280 True 0.301 -0.023 True 0.323 False False
shuffle 5 92 13 0.013 0.069 0.308 -0.021 -190 True 0.328 0.015 True 0.264 True True
shuffle 6 87 15 -0.060 0.291 0.196 -0.008 -230 True 0.340 0.000 True -0.396 False False
shuffle 7 90 15 -0.024 0.226 0.273 0.053 -280 True 0.332 0.009 True 0.146 True True
shuffle 8 87 13 0.019 0.184 0.206 -0.062 -250 True 0.292 0.004 True 0.207 False False
shuffle 9 96 11 0.000 0.422 0.201 -0.102 -240 True 0.308 -0.015 True -0.753 True True
shuffle 10 90 13 0.012 0.285 0.317 0.012 -230 True 0.369 -0.003 True 0.014 True True
shuffle 11 86 19 -0.005 0.366 0.258 -0.001 -240 True 0.367 0.006 True 0.048 False False
shuffle 12 88 13 0.008 0.098 0.228 -0.015 -260 True 0.342 0.017 True 0.236 True True
shuffle 13 87 16 0.028 0.227 0.345 0.022 -210 True 0.332 0.034 True -0.041 False False
shuffle 14 89 16 0.074 0.158 0.268 0.040 -230 True 0.308 0.031 True -0.496 False False
shuffle 15 88 15 -0.038 0.101 0.284 0.017 -250 True 0.308 -0.023 True -0.314 True True
shuffle 16 95 13 -0.000 0.269 0.246 0.009 -260 True 0.275 0.034 True -0.204 False False
shuffle 17 93 15 0.025 0.150 0.298 0.009 -210 True 0.327 0.010 True -0.163 True True
shuffle 18 90 13 -0.066 0.260 0.167 -0.046 -260 True 0.333 0.023 True 0.288 False False
shuffle 19 92 18 -0.106 0.287 0.295 -0.047 -240 True 0.334 0.004 True -0.140 True True
shuffle 20 92 10 0.043 0.146 0.387 -0.004 170 False 0.369 0.035 True -0.112 True True
over 20 shuffles: units passing > 5 spikes/s range 82 to 96 (data 93 of 142); fig4e passes (a) 18, (b) 20; suppfig1 passes 11;
mean r -0.106 to 0.074; IQR of r 0.069 to 0.422; fig4e(b) r at -200 0.28 to 0.37 and at onset -0.028 to 0.035
/home/runner/work/encore/encore/entries/000140-churchland-2010/maze_prep.py:270: RuntimeWarning: invalid value encountered in divide return (a * b).sum(0) / np.sqrt((a ** 2).sum(0) * (b ** 2).sum(0))
disjoint trials for the preparatory and movement-locked means (units of fig4d, O1 range, O2 i) O3 O4 split conditions trials n fig4 mean r IQR r fig4e(a) r -200 r 0 drop pass (a) fig4e(b) r -200 r 0 (b) pass (b) suppfig1 ratio top at -150 suppfig1 pass a a prep A, move B 23 35/27 19 -0.090 0.232 0.134 -0.003 180 False 0.171 0.030 True 0.282 False False a a prep B, move A 23 27/35 19 -0.048 0.313 0.090 -0.077 -360 False 0.137 -0.021 False -0.176 False False a b prep A, move B 24 38/30 13 -0.100 0.171 0.100 -0.032 200 False 0.177 0.024 True 0.208 False False a b prep B, move A 24 30/38 13 -0.118 0.386 0.106 -0.049 -60 True 0.146 -0.008 True -0.434 False False a c prep A, move B 25 44/32 15 0.025 0.162 0.115 0.045 250 False 0.184 0.024 True 0.017 False False a c prep B, move A 25 32/44 15 -0.153 0.274 0.140 -0.104 -200 True 0.161 0.009 False -0.610 False False b a prep A, move B 23 34/28 18 -0.049 0.215 0.043 0.029 210 False 0.126 0.017 True 0.325 False False b a prep B, move A 23 28/34 18 -0.062 0.215 0.134 -0.095 90 False 0.122 -0.019 False -0.306 True True b b prep A, move B 24 40/30 15 0.036 0.311 0.022 0.052 210 False 0.128 0.013 True 0.327 True True b b prep B, move A 24 30/40 15 -0.128 0.216 0.145 -0.090 -200 True 0.164 -0.007 True -0.746 False False b c prep A, move B 25 44/32 15 0.045 0.227 0.029 0.062 220 False 0.146 0.015 True 0.371 False False b c prep B, move A 25 32/44 15 -0.121 0.261 0.094 -0.099 -220 True 0.159 0.003 True -0.732 False False disjoint readings: fig4e passes (a) 4, (b) 9 of 12; suppfig1 passes 2; fig4e(b) r at -200 0.12 to 0.18
Noise reference. With condition labels shuffled, 82 to 96 units pass the > 5 spikes/s range rule, against 93 with the true labels: under O1a the selection admits units whose range comes from trial noise. The paper justified the 5 spikes/s threshold by its trial count ("acceptable given the higher signal-to-noise ratio provided by the higher trial count", Methods), which this dandiset does not have. The conclusions of fig4d, fig4-broad-near-zero, fig4e and suppfig1 below refer to both references.
Result: not reproduced. Under none of the 18 readings with units is the mean correlation within 0.13 ± 0.013 or positive at p < 0.05: mean r ranges from -0.148 to 0.095 and mean Spearman r from -0.061 to 0.061 (main reading -0.079 and -0.020, 19 units, 27 conditions); the one reading at p < 0.05 has a negative mean. With condition labels shuffled the mean r is -0.106 to 0.074, so the main value is within what trial noise gives. The condition means here rest on 1 to 4 trials, and their split-half reliability (above) is low, so a weak true correlation such as 0.13 could not be resolved; the outcome is a property of these 100 trials, not a test of the paper's J-array data.
fig4abc-mean-correlation-other-datasets¶
Blocked, data_missing: the dandiset holds no single-electrode recordings of monkeys A, B or J (one
subject, Jenkins, two Utah arrays).
fig4-broad-near-zero¶
Correlations between preparatory and peri-movement tuning span a broad range centred near zero, with no separate populations of concordant (r near 1) and opposing (r near -1) tuning; most neurons have weakly correlated tuning.
Criterion: on the units and correlations of fig4d, median r in [-0.3, 0.3], IQR at least 0.3, more than half of units with |r| < 0.5, fewer than 15% above 0.8 and fewer than 15% below -0.8.
def broad(r):
q1, q3 = np.percentile(r, [25, 75])
out = {"n": len(r), "median": np.median(r), "IQR": q3 - q1, "frac |r|<0.5": np.mean(np.abs(r) < 0.5),
"frac r>0.8": np.mean(r > 0.8), "frac r<-0.8": np.mean(r < -0.8)}
out["pass"] = (-0.3 <= out["median"] <= 0.3 and out["IQR"] >= 0.3 and out["frac |r|<0.5"] > 0.5
and out["frac r>0.8"] < 0.15 and out["frac r<-0.8"] < 0.15)
return out
rows = []
for (o, f), cfg in CONFIGS.items():
for O1 in ("range", "sd"):
for O2 in ("i", "ii"):
sel = fig4_units(cfg, O1, O2)
if sel.sum() >= 3:
rows.append({"O3": o, "O4": f, "O1": O1, "O2": O2, **broad(fig4_corrs(cfg, sel)[0])})
broad_table = pd.DataFrame(rows)
print(broad_table.round(3).to_string(index=False))
print("readings passing:", int(broad_table["pass"].sum()), "of", len(broad_table))
print("main reading r range: %.2f to %.2f" % (r_main.min(), r_main.max()))
fig, ax = plt.subplots(1, 2, figsize=(8, 3))
for a, r, t in ((ax[0], r_main, "O2 i (main)"), (ax[1], r_ii, "O2 ii")):
a.hist(r, bins=np.arange(-1, 1.01, 0.1), color="0.4")
a.axvline(r.mean(), color="r")
a.set_xlabel("r, preparatory vs peri-movement"); a.set_title(f"{t}, n = {len(r)}")
ax[0].set_ylabel("units")
plt.tight_layout(); plt.show()
O3 O4 O1 O2 n median IQR frac |r|<0.5 frac r>0.8 frac r<-0.8 pass a a range i 19 -0.066 0.201 0.947 0.000 0.0 False a a range ii 88 -0.018 0.262 0.955 0.000 0.0 False a a sd ii 5 0.095 0.691 0.800 0.000 0.0 True a b range i 13 -0.098 0.246 0.923 0.000 0.0 False a b range ii 85 -0.001 0.285 0.953 0.000 0.0 False a b sd ii 4 -0.064 0.422 1.000 0.000 0.0 True a c range i 15 -0.035 0.226 0.933 0.000 0.0 False a c range ii 85 0.008 0.286 0.965 0.000 0.0 False a c sd ii 4 -0.069 0.434 1.000 0.000 0.0 True b a range i 18 -0.008 0.193 0.944 0.000 0.0 False b a range ii 89 -0.024 0.244 0.966 0.000 0.0 False b a sd ii 5 -0.033 0.646 0.800 0.000 0.0 True b b range i 15 -0.010 0.189 0.933 0.000 0.0 False b b range ii 86 0.002 0.272 0.977 0.012 0.0 False b b sd ii 4 -0.093 0.419 1.000 0.000 0.0 True b c range i 15 -0.010 0.173 1.000 0.000 0.0 False b c range ii 83 0.016 0.250 0.964 0.000 0.0 False b c sd ii 3 0.052 0.329 1.000 0.000 0.0 True readings passing: 6 of 18 main reading r range: -0.53 to 0.18
Result: not reproduced, resting on the assumption that tuning is the range across conditions (O1a). Under O1a (12 readings) the median correlation is near zero (-0.098 to 0.016), 92% to 100% of units have |r| < 0.5, at most 1.2% lie above 0.8 and none below -0.8, but the interquartile range (0.173 to 0.286) is below the 0.3 the criterion asks for a broad range in every O1a reading. The criterion has no partial rule, and the condition that fails is the one that tests the breadth the claim states. The conditions that hold are also met by noise: with condition labels shuffled the IQR is 0.069 to 0.422 and the mean r -0.106 to 0.074. Under O1b (SD; 3 to 5 units) all 6 readings pass.
fig4e-correlation-time-course¶
When peri-movement activity is taken at single time points, its correlation with preparatory tuning starts to fall about 200 ms before movement onset and is near zero by movement onset.
Criterion: the rate at single time points of the 20 ms smoothed condition means, swept along the go-locked and movement-locked means, correlated per unit with the preparatory vector, averaged over units ((a) the fig4 units, (b) all units with > 5 spikes/s preparatory tuning). Reproduced when the mean at movement onset is below 0.1 and below half its value at -200 ms, and the largest drop over a 50 ms step begins between -300 ms and onset.
rows, curves = [], {}
for (o, f), cfg in CONFIGS.items():
for O1 in ("range", "sd"):
for O2 in ("i", "ii"):
for units_reading in ("a", "b"):
sel = fig4_units(cfg, O1, O2) if units_reading == "a" else mp.tuning(cfg["P"], O1) > mp.PREP_TUNING_MIN
if sel.sum() < 3 or (units_reading == "b" and O2 == "ii"):
continue
go, mv = time_course(cfg, sel)
curves[(o, f, O1, O2, units_reading)] = (go, mv)
rows.append({"O3": o, "O4": f, "O1": O1, "O2": O2 if units_reading == "a" else "-",
"units": units_reading, "n": int(sel.sum()), "r at go-200": go[0], "r at go": go[LAGS_GO == 0][0],
**judge_tc(mv)})
tc_table = pd.DataFrame(rows)
print(tc_table.round(3).to_string(index=False))
print("readings passing:", int(tc_table["pass"].sum()), "of", len(tc_table))
go_m, mv_m = curves[("a", "a", "range", "i", "a")]
go_b, mv_b = curves[("a", "a", "range", "i", "b")]
print("main, fig4 units: r from %.2f to %.2f over -400..400 ms re onset" % (mv_m.min(), mv_m.max()))
print("main, all prep-tuned units: r at -400 %.2f, -200 %.2f, -100 %.2f, 0 %.2f, +200 %.2f"
% tuple(mv_b[np.isin(LAGS_MV, [-400, -200, -100, 0, 200])]))
fig, ax = plt.subplots(1, 2, figsize=(8, 3), sharey=True, gridspec_kw={"width_ratios": [1, 2]})
for go, mv, lab in ((go_m, mv_m, "fig4 units"), (go_b, mv_b, "all > 5 spikes/s")):
ax[0].plot(LAGS_GO, go, label=lab); ax[1].plot(LAGS_MV, mv, label=lab)
ax[0].set_xlabel("ms re go cue"); ax[1].set_xlabel("ms re movement onset"); ax[0].set_ylabel("mean r with preparatory rates")
for a in ax: a.axhline(0, color="k", lw=0.5)
ax[1].legend(); plt.tight_layout(); plt.show()
/home/runner/work/encore/encore/entries/000140-churchland-2010/maze_prep.py:270: RuntimeWarning: invalid value encountered in divide return (a * b).sum(0) / np.sqrt((a ** 2).sum(0) * (b ** 2).sum(0))
O3 O4 O1 O2 units n r at go-200 r at go r at -200 r at 0 largest 50 ms drop starts pass partial a a range i a 19 0.341 0.620 0.295 -0.078 -250 True False a a range - b 93 0.410 0.605 0.388 -0.001 -220 True False a a range ii a 88 0.404 0.595 0.376 -0.011 -220 True False a a sd - b 8 0.621 0.787 0.555 0.106 -130 False True a a sd ii a 5 0.632 0.801 0.567 0.091 90 False False a b range i a 13 0.350 0.643 0.298 -0.075 -240 True False a b range - b 91 0.436 0.611 0.399 0.019 -220 True False a b range ii a 85 0.423 0.596 0.381 0.008 -210 True False a b sd - b 8 0.637 0.803 0.593 0.126 -130 False True a b sd ii a 4 0.565 0.782 0.543 0.147 -220 False True a c range i a 15 0.334 0.604 0.330 -0.065 -230 True False a c range - b 91 0.452 0.614 0.414 0.022 -220 True False a c range ii a 85 0.442 0.599 0.395 0.012 -210 True False a c sd - b 7 0.691 0.834 0.661 0.172 -130 False True a c sd ii a 4 0.625 0.766 0.561 0.142 90 False True b a range i a 18 0.306 0.578 0.216 -0.018 -270 True False b a range - b 93 0.416 0.606 0.299 0.021 -250 True False b a range ii a 89 0.404 0.597 0.288 0.015 -230 True False b a sd - b 8 0.619 0.784 0.482 0.180 -160 False True b a sd ii a 5 0.629 0.797 0.492 0.210 -230 False True b b range i a 15 0.328 0.606 0.207 -0.026 -260 True False b b range - b 91 0.438 0.608 0.295 0.036 -250 True False b b range ii a 86 0.426 0.594 0.276 0.026 -250 True False b b sd - b 8 0.646 0.794 0.468 0.180 -150 False True b b sd ii a 4 0.585 0.765 0.349 0.177 -240 False True b c range i a 15 0.307 0.602 0.200 -0.043 -250 True False b c range - b 91 0.452 0.614 0.302 0.021 -250 True False b c range ii a 83 0.436 0.599 0.277 0.007 -250 True False b c sd - b 7 0.691 0.834 0.514 0.244 -150 False True b c sd ii a 3 0.663 0.770 0.444 0.240 60 False True readings passing: 18 of 30 main, fig4 units: r from -0.13 to 0.64 over -400..400 ms re onset main, all prep-tuned units: r at -400 0.62, -200 0.39, -100 0.12, 0 -0.00, +200 -0.02
Result: reproduced, resting on the assumption that tuning is the range across conditions (O1a); trial noise also meets the criterion, so the pass carries little evidence. Under every reading with O1a (18 of 18) the mean correlation at movement onset is below 0.1 and below half its value at 200 ms before onset, and the largest 50 ms drop starts between 270 and 210 ms before onset; in the main reading it falls from 0.295 at -200 ms to -0.078 at onset, and for all units with preparatory tuning from 0.39 to -0.00. With condition labels shuffled the criterion is met in 18 of 20 shuffles for the fig4 units and in 20 of 20 for all units with > 5 spikes/s (r at -200 ms 0.28 to 0.37, at onset -0.028 to 0.035): the preparatory and movement-locked means share the noise of the same 1 to 4 trials, and the preparatory window ends about 220 ms before onset. With the two means taken from disjoint trials, the correlation at -200 ms is only 0.12 to 0.18 for all units, and the criterion is met in 4 of 12 readings (fig4 units) and 9 of 12 (all units). Under O1b (SD, 3 to 8 units) the correlation at onset stays between 0.09 and 0.24, partially reproduced in 11 of 12 readings; in the twelfth (O3a, O4a, O2 ii, 5 units) r at onset is 0.091 but the largest drop starts at +90 ms, so it is neither a pass nor a partial.
suppfig1-population-averages-lose-tuning¶
Population-average PSTHs grouped by each neuron's preparatory preference show weak peri-movement differences that shrink with time, with little relationship to preparatory tuning by movement onset.
Criterion: per fig4 unit, the conditions with the strongest and weakest preparatory rate and those at the 1/3 and 2/3 ranks; each unit normalised by its maximum rate; averaged. Reproduced when the spread between the strongest and weakest traces at movement onset is less than half the spread in the preparatory epoch, with the strongest still on top at -150 ms. Unnormalised average as a variant.
rows = []
for (o, f), cfg in CONFIGS.items():
for O1 in ("range", "sd"):
for O2 in ("i", "ii"):
sel = fig4_units(cfg, O1, O2)
if sel.sum() < 3:
continue
for norm in (True, False):
rows.append({"O3": o, "O4": f, "O1": O1, "O2": O2, "n": int(sel.sum()), "normalised": norm,
**judge_pop(*pop_avg(cfg, sel, norm))})
pop_table = pd.DataFrame(rows)
print(pop_table.round(3).to_string(index=False))
print("normalised readings passing:", int(pop_table[pop_table.normalised]["pass"].sum()), "of", int(pop_table.normalised.sum()),
"| unnormalised:", int(pop_table[~pop_table.normalised]["pass"].sum()), "of", int((~pop_table.normalised).sum()))
g_tr, m_tr = pop_avg(CONFIGS[MAIN], sel_main, True)
fig, ax = plt.subplots(1, 2, figsize=(8, 3), sharey=True, gridspec_kw={"width_ratios": [1, 2]})
cols = ["#b2182b", "#ef8a62", "#67a9cf", "#1a9850"]
for k in range(4):
ax[0].plot(GO_L, g_tr[k], color=cols[k]); ax[1].plot(MV_L, m_tr[k], color=cols[k])
ax[0].set_xlabel("ms re go cue"); ax[1].set_xlabel("ms re movement onset"); ax[0].set_ylabel("normalised rate")
plt.tight_layout(); plt.show()
print("plotted normalised rates range %.3f to %.3f" % (min(g_tr.min(), m_tr.min()), max(g_tr.max(), m_tr.max())))
O3 O4 O1 O2 n normalised spread prep spread at onset ratio strongest on top at -150 pass a a range i 19 True 0.170000 0.032 0.186 True True a a range i 19 False 9.333000 2.125 0.228 True True a a range ii 88 True 0.266000 0.006 0.022 True True a a range ii 88 False 12.421000 0.520 0.042 True True a a sd ii 5 True 0.400000 -0.019 -0.048 True True a a sd ii 5 False 25.027000 -1.245 -0.050 True True a b range i 13 True 0.163000 0.019 0.116 True True a b range i 13 False 9.425000 1.129 0.120 True True a b range ii 85 True 0.270000 0.028 0.103 True True a b range ii 85 False 12.257000 1.327 0.108 True True a b sd ii 4 True 0.385000 0.016 0.042 True True a b sd ii 4 False 24.746000 0.746 0.030 True True a c range i 15 True 0.170000 -0.026 -0.155 True True a c range i 15 False 9.027000 -1.579 -0.175 True True a c range ii 85 True 0.270000 0.022 0.083 True True a c range ii 85 False 11.564000 1.199 0.104 True True a c sd ii 4 True 0.386000 0.071 0.183 True True a c sd ii 4 False 23.652000 3.812 0.161 True True b a range i 18 True 0.187000 0.050 0.270 True True b a range i 18 False 10.141000 2.083 0.205 True True b a range ii 89 True 0.269000 0.017 0.063 True True b a range ii 89 False 12.579000 1.177 0.094 True True b a sd ii 5 True 0.394000 0.078 0.199 True True b a sd ii 5 False 24.879000 5.177 0.208 True True b b range i 15 True 0.179000 0.029 0.164 True True b b range i 15 False 9.831000 -0.527 -0.054 True True b b range ii 86 True 0.268000 0.020 0.074 True True b b range ii 86 False 12.039000 1.217 0.101 True True b b sd ii 4 True 0.384000 0.108 0.281 True True b b sd ii 4 False 24.162001 6.411 0.265 True True b c range i 15 True 0.172000 0.023 0.134 True True b c range i 15 False 8.596000 1.016 0.118 True True b c range ii 83 True 0.265000 0.005 0.021 True True b c range ii 83 False 11.481000 0.777 0.068 True True b c sd ii 3 True 0.364000 0.107 0.295 True True b c sd ii 3 False 23.391001 6.358 0.272 True True normalised readings passing: 18 of 18 | unnormalised: 18 of 18
plotted normalised rates range 0.000 to 0.257
Result: reproduced; trial noise also meets the criterion in about half the shuffles, so the pass carries little evidence. In all 36 readings (normalised and not) the spread between the strongest and weakest preparatory conditions at movement onset is less than half the preparatory spread (ratio 0.186 in the main normalised reading, -0.175 to 0.295 over readings), with the strongest-preparatory trace on top at 150 ms before onset. With condition labels shuffled the criterion is met in 11 of 20 shuffles. With the preparatory and movement-locked means from disjoint trials it is met in 2 of 12 readings. This suggests that the strongest trace being on top at -150 ms, the part of the result that noise does not give every time, comes mostly from noise shared within the same trials, but it does not show it: each disjoint half uses about half the trials, and the added noise could also hide a real effect.
fig6d-perimovement-best-j-array¶
For the J-array dataset, the cross-validated performance of a preferred direction in capturing preparatory tuning is highest when the PD is expressed in the PCA-based peri-movement space, compared with the kinematic space and the task spaces.
Criterion: units with > 5 spikes/s preparatory tuning (O1); leave-one-condition-out regression on a space plus offset; performance 1 - MSE / variance, averaged over units. Peri-movement space: PC scores of the condition means from -150 to +800 ms of all other units, not normalised. Kinematic space: PC scores of x, y position, velocity and acceleration from -150 ms to the end of the reach (600 or 800 ms), each scaled to unit range; reading (a) repeats the six variables to about n columns and (b) uses them once. The eight maze task spaces of the Methods. k = 1 to 14 and 1 to 20. Reproduced when under each reading the best peri-movement performance over k exceeds the best kinematic one and every task space; partially when it exceeds all but one.
Readings (a) and (b) of the kinematic matrix give the same scores up to a common scale (each variable repeated the same number of times multiplies the Gram matrix by that number), and regression with an offset is unchanged by scaling, so one is computed and the identity is checked below. Segment directions are computed with segments of equal path length and, as a variant, of equal duration.
KMAX = 20
t0 = time.time()
def perf_matrix(P, spaces_per_unit, idx, kmax=KMAX):
return np.array([mp.loo_performance_curve(P[:, j], spaces_per_unit[j], kmax) for j in idx])
PERF = {}
for key, cfg in CONFIGS.items():
P = cfg["P"]
peri = mp.leave_one_out_scores(cfg["CMm"][:, :, TSEL], KMAX)
upper = mp.leave_one_out_scores(P[:, :, None], KMAX)
allidx = np.arange(len(units))
res = {"peri": perf_matrix(P, peri, allidx), "upper": perf_matrix(P, upper, allidx)}
for end in (600, 800):
ks = mp.pca_scores(mp.kinematic_matrix(cfg["KIN"], end), KMAX)
res[f"kin{end}"] = perf_matrix(P, [ks] * len(units), allidx)
for seg in ("length", "time"):
for name, S in cfg["TS_" + seg].items():
res[f"{seg}:{name}"] = np.array([mp.loo_performance(P[:, j], S) for j in allidx])
res["peri_scores"] = peri
PERF[key] = res
print(f"performance for all units and readings in {time.time() - t0:.1f} s")
# check: repeating the six kinematic variables equally leaves performance unchanged
K = mp.kinematic_matrix(CONFIGS[MAIN]["KIN"], 800)
Krep = np.tile(K, (1, 23))
p1 = mp.loo_performance_curve(CONFIGS[MAIN]["P"][:, 0], mp.pca_scores(K, 10), 10)
p2 = mp.loo_performance_curve(CONFIGS[MAIN]["P"][:, 0], mp.pca_scores(Krep, 10), 10)
print("kinematic reading (a) vs (b), max abs difference in performance (unit 0, k 1-10): %.2e" % np.nanmax(np.abs(p1 - p2)))
def fig6_summary(key, O1, kmax, end, seg):
cfg, res = CONFIGS[key], PERF[key]
sel = mp.tuning(cfg["P"], O1) > mp.PREP_TUNING_MIN
peri = np.nanmean(res["peri"][sel, :kmax], 0)
kin = np.nanmean(res[f"kin{end}"][sel, :kmax], 0)
tasks = {name: np.nanmean(res[f"{seg}:{name}"][sel]) for name in mp.TASK_SPACES}
return sel, peri, kin, tasks
rows = []
for key in CONFIGS:
for O1 in ("range", "sd"):
for kmax in (14, 20):
for end in (600, 800):
for seg in ("length", "time"):
sel, peri, kin, tasks = fig6_summary(key, O1, kmax, end, seg)
best_p = peri.max()
others = {"kinematic": kin.max(), **tasks}
beaten = [k for k, v in others.items() if v >= best_p]
rows.append({"O3": key[0], "O4": key[1], "O1": O1, "kmax": kmax, "kin end": end, "segments": seg,
"n": int(sel.sum()), "best peri": best_p, "best peri k": int(np.argmax(peri) + 1),
"best kin": kin.max(), "best task": max(tasks.values()),
"best task space": max(tasks, key=tasks.get), "spaces >= peri": len(beaten),
"outcome": "reproduced" if not beaten else ("partial" if len(beaten) == 1 else "no")})
fig6_table = pd.DataFrame(rows)
print(fig6_table.round(3).to_string(index=False))
print(fig6_table.outcome.value_counts().to_dict())
sel6, peri6, kin6, tasks6 = fig6_summary(MAIN, "range", 20, 800, "length")
upper6 = np.nanmean(PERF[MAIN]["upper"][sel6], 0)
print("main reading: n = %d units; mean performance by k:" % sel6.sum())
print(pd.DataFrame({"k": np.arange(1, KMAX + 1), "peri": peri6, "kinematic": kin6, "upper limit": upper6}).round(3).to_string(index=False))
print("task spaces:", {k: round(v, 3) for k, v in tasks6.items()})
med = {k: np.nanmedian(PERF[MAIN][k][sel6, :6], 0) for k in ("peri", "kin800", "upper")}
print("median over units, k = 1..6:", {k: np.round(v, 3) for k, v in med.items()})
print("units with peri performance > 0 at their best k <= 14: %d of %d" % ((np.nanmax(PERF[MAIN]["peri"][sel6, :14], 1) > 0).sum(), sel6.sum()))
sd_sel = mp.tuning(CONFIGS[MAIN]["P"], "sd") > 5
print("O1 SD units (n = %d): peri best %.3f, kin best %.3f, upper best %.3f"
% (sd_sel.sum(), np.nanmean(PERF[MAIN]["peri"][sd_sel, :14], 0).max(), np.nanmean(PERF[MAIN]["kin800"][sd_sel, :14], 0).max(),
np.nanmean(PERF[MAIN]["upper"][sd_sel, :14], 0).max()))
fig, ax = plt.subplots(figsize=(5, 3.5))
ks = np.arange(1, KMAX + 1)
ax.plot(ks, peri6, "k-o", ms=3, label="peri-movement"); ax.plot(ks, kin6, "-o", color="tab:purple", ms=3, label="kinematic")
ax.plot(ks, upper6, color="0.6", label="upper limit")
dims = {name: S.shape[1] for name, S in CONFIGS[MAIN]["TS_length"].items()}
ax.scatter([dims[n] for n in tasks6], list(tasks6.values()), color="tab:green", label="task spaces", zorder=3)
ax.set_ylim(-1.5, 0.2); ax.set_xlabel("dimensionality"); ax.set_ylabel("mean performance"); ax.legend(fontsize=7)
plt.tight_layout(); plt.show()
allv = np.r_[peri6, kin6, upper6, list(tasks6.values())]
print("points below the axis (< -1.5): %d of %d; plotted values range %.2f to %.2f" % ((allv < -1.5).sum(), len(allv), allv.min(), allv.max()))
/home/runner/work/encore/encore/entries/000140-churchland-2010/maze_prep.py:285: RuntimeWarning: invalid value encountered in scalar divide return 1.0 - np.mean(press ** 2) / np.var(x)
performance for all units and readings in 4.2 s
kinematic reading (a) vs (b), max abs difference in performance (unit 0, k 1-10): 1.19e-14
O3 O4 O1 kmax kin end segments n best peri best peri k best kin best task best task space spaces >= peri outcome
a a range 14 600 length 93 -0.044 2 -0.052 -0.048 initial vel. 0 reproduced
a a range 14 600 time 93 -0.044 2 -0.052 -0.048 initial vel. 0 reproduced
a a range 14 800 length 93 -0.044 2 -0.053 -0.048 initial vel. 0 reproduced
a a range 14 800 time 93 -0.044 2 -0.053 -0.048 initial vel. 0 reproduced
a a range 20 600 length 93 -0.044 2 -0.052 -0.048 initial vel. 0 reproduced
a a range 20 600 time 93 -0.044 2 -0.052 -0.048 initial vel. 0 reproduced
a a range 20 800 length 93 -0.044 2 -0.053 -0.048 initial vel. 0 reproduced
a a range 20 800 time 93 -0.044 2 -0.053 -0.048 initial vel. 0 reproduced
a a sd 14 600 length 8 0.209 4 0.138 0.133 endpoints + halfpoint 0 reproduced
a a sd 14 600 time 8 0.209 4 0.138 0.135 segment dirs. 0 reproduced
a a sd 14 800 length 8 0.209 4 0.136 0.133 endpoints + halfpoint 0 reproduced
a a sd 14 800 time 8 0.209 4 0.136 0.135 segment dirs. 0 reproduced
a a sd 20 600 length 8 0.209 4 0.138 0.133 endpoints + halfpoint 0 reproduced
a a sd 20 600 time 8 0.209 4 0.138 0.135 segment dirs. 0 reproduced
a a sd 20 800 length 8 0.209 4 0.136 0.133 endpoints + halfpoint 0 reproduced
a a sd 20 800 time 8 0.209 4 0.136 0.135 segment dirs. 0 reproduced
a b range 14 600 length 91 -0.045 2 -0.048 -0.041 initial vel. 1 partial
a b range 14 600 time 91 -0.045 2 -0.048 -0.041 initial vel. 1 partial
a b range 14 800 length 91 -0.045 2 -0.050 -0.041 initial vel. 1 partial
a b range 14 800 time 91 -0.045 2 -0.050 -0.041 initial vel. 1 partial
a b range 20 600 length 91 -0.045 2 -0.048 -0.041 initial vel. 1 partial
a b range 20 600 time 91 -0.045 2 -0.048 -0.041 initial vel. 1 partial
a b range 20 800 length 91 -0.045 2 -0.050 -0.041 initial vel. 1 partial
a b range 20 800 time 91 -0.045 2 -0.050 -0.041 initial vel. 1 partial
a b sd 14 600 length 8 0.202 4 0.135 0.117 endpoints + halfpoint 0 reproduced
a b sd 14 600 time 8 0.202 4 0.135 0.117 endpoints + halfpoint 0 reproduced
a b sd 14 800 length 8 0.202 4 0.133 0.117 endpoints + halfpoint 0 reproduced
a b sd 14 800 time 8 0.202 4 0.133 0.117 endpoints + halfpoint 0 reproduced
a b sd 20 600 length 8 0.202 4 0.135 0.117 endpoints + halfpoint 0 reproduced
a b sd 20 600 time 8 0.202 4 0.135 0.117 endpoints + halfpoint 0 reproduced
a b sd 20 800 length 8 0.202 4 0.133 0.117 endpoints + halfpoint 0 reproduced
a b sd 20 800 time 8 0.202 4 0.133 0.117 endpoints + halfpoint 0 reproduced
a c range 14 600 length 91 -0.023 2 -0.040 -0.036 initial vel. 0 reproduced
a c range 14 600 time 91 -0.023 2 -0.040 -0.036 initial vel. 0 reproduced
a c range 14 800 length 91 -0.023 2 -0.043 -0.036 initial vel. 0 reproduced
a c range 14 800 time 91 -0.023 2 -0.043 -0.036 initial vel. 0 reproduced
a c range 20 600 length 91 -0.023 2 -0.040 -0.036 initial vel. 0 reproduced
a c range 20 600 time 91 -0.023 2 -0.040 -0.036 initial vel. 0 reproduced
a c range 20 800 length 91 -0.023 2 -0.043 -0.036 initial vel. 0 reproduced
a c range 20 800 time 91 -0.023 2 -0.043 -0.036 initial vel. 0 reproduced
a c sd 14 600 length 7 0.247 2 0.186 0.228 kitchen sink 0 reproduced
a c sd 14 600 time 7 0.247 2 0.186 0.186 kitchen sink 0 reproduced
a c sd 14 800 length 7 0.247 2 0.182 0.228 kitchen sink 0 reproduced
a c sd 14 800 time 7 0.247 2 0.182 0.186 kitchen sink 0 reproduced
a c sd 20 600 length 7 0.247 2 0.186 0.228 kitchen sink 0 reproduced
a c sd 20 600 time 7 0.247 2 0.186 0.186 kitchen sink 0 reproduced
a c sd 20 800 length 7 0.247 2 0.182 0.228 kitchen sink 0 reproduced
a c sd 20 800 time 7 0.247 2 0.182 0.186 kitchen sink 0 reproduced
b a range 14 600 length 93 -0.046 2 -0.051 -0.042 initial vel. 1 partial
b a range 14 600 time 93 -0.046 2 -0.051 -0.042 initial vel. 1 partial
b a range 14 800 length 93 -0.046 2 -0.053 -0.042 initial vel. 1 partial
b a range 14 800 time 93 -0.046 2 -0.053 -0.042 initial vel. 1 partial
b a range 20 600 length 93 -0.046 2 -0.051 -0.042 initial vel. 1 partial
b a range 20 600 time 93 -0.046 2 -0.051 -0.042 initial vel. 1 partial
b a range 20 800 length 93 -0.046 2 -0.053 -0.042 initial vel. 1 partial
b a range 20 800 time 93 -0.046 2 -0.053 -0.042 initial vel. 1 partial
b a sd 14 600 length 8 0.140 4 0.135 0.128 endpoints + halfpoint 0 reproduced
b a sd 14 600 time 8 0.140 4 0.135 0.128 endpoints + halfpoint 0 reproduced
b a sd 14 800 length 8 0.140 4 0.133 0.128 endpoints + halfpoint 0 reproduced
b a sd 14 800 time 8 0.140 4 0.133 0.128 endpoints + halfpoint 0 reproduced
b a sd 20 600 length 8 0.140 4 0.135 0.128 endpoints + halfpoint 0 reproduced
b a sd 20 600 time 8 0.140 4 0.135 0.128 endpoints + halfpoint 0 reproduced
b a sd 20 800 length 8 0.140 4 0.133 0.128 endpoints + halfpoint 0 reproduced
b a sd 20 800 time 8 0.140 4 0.133 0.128 endpoints + halfpoint 0 reproduced
b b range 14 600 length 91 -0.043 2 -0.046 -0.033 initial vel. 1 partial
b b range 14 600 time 91 -0.043 2 -0.046 -0.033 initial vel. 1 partial
b b range 14 800 length 91 -0.043 2 -0.049 -0.033 initial vel. 1 partial
b b range 14 800 time 91 -0.043 2 -0.049 -0.033 initial vel. 1 partial
b b range 20 600 length 91 -0.043 2 -0.046 -0.033 initial vel. 1 partial
b b range 20 600 time 91 -0.043 2 -0.046 -0.033 initial vel. 1 partial
b b range 20 800 length 91 -0.043 2 -0.049 -0.033 initial vel. 1 partial
b b range 20 800 time 91 -0.043 2 -0.049 -0.033 initial vel. 1 partial
b b sd 14 600 length 8 0.129 2 0.134 0.121 initial vel. 1 partial
b b sd 14 600 time 8 0.129 2 0.134 0.121 initial vel. 1 partial
b b sd 14 800 length 8 0.129 2 0.132 0.121 initial vel. 1 partial
b b sd 14 800 time 8 0.129 2 0.132 0.121 initial vel. 1 partial
b b sd 20 600 length 8 0.129 2 0.134 0.121 initial vel. 1 partial
b b sd 20 600 time 8 0.129 2 0.134 0.121 initial vel. 1 partial
b b sd 20 800 length 8 0.129 2 0.132 0.121 initial vel. 1 partial
b b sd 20 800 time 8 0.129 2 0.132 0.121 initial vel. 1 partial
b c range 14 600 length 91 -0.026 2 -0.039 -0.029 initial vel. 0 reproduced
b c range 14 600 time 91 -0.026 2 -0.039 -0.029 initial vel. 0 reproduced
b c range 14 800 length 91 -0.026 2 -0.043 -0.029 initial vel. 0 reproduced
b c range 14 800 time 91 -0.026 2 -0.043 -0.029 initial vel. 0 reproduced
b c range 20 600 length 91 -0.026 2 -0.039 -0.029 initial vel. 0 reproduced
b c range 20 600 time 91 -0.026 2 -0.039 -0.029 initial vel. 0 reproduced
b c range 20 800 length 91 -0.026 2 -0.043 -0.029 initial vel. 0 reproduced
b c range 20 800 time 91 -0.026 2 -0.043 -0.029 initial vel. 0 reproduced
b c sd 14 600 length 7 0.249 2 0.187 0.217 kitchen sink 0 reproduced
b c sd 14 600 time 7 0.249 2 0.187 0.175 initial vel. 0 reproduced
b c sd 14 800 length 7 0.249 2 0.183 0.217 kitchen sink 0 reproduced
b c sd 14 800 time 7 0.249 2 0.183 0.175 initial vel. 0 reproduced
b c sd 20 600 length 7 0.249 2 0.187 0.217 kitchen sink 0 reproduced
b c sd 20 600 time 7 0.249 2 0.187 0.175 initial vel. 0 reproduced
b c sd 20 800 length 7 0.249 2 0.183 0.217 kitchen sink 0 reproduced
b c sd 20 800 time 7 0.249 2 0.183 0.175 initial vel. 0 reproduced
{'reproduced': 64, 'partial': 32}
main reading: n = 93 units; mean performance by k:
k peri kinematic upper limit
1 -0.055 -0.053 -0.045
2 -0.044 -0.058 -0.067
3 -0.104 -0.081 -0.086
4 -0.275 -0.091 -0.114
5 -0.981 -0.126 -0.159
6 -1.369 -0.191 -0.192
7 -2.285 -0.257 -0.343
8 -3.124 -0.324 -0.476
9 -3.487 -0.372 -0.528
10 -4.069 -0.428 -0.589
11 -4.399 -0.553 -0.670
12 -5.042 -0.778 -0.825
13 -5.716 -0.882 -1.111
14 -8.266 -1.222 -1.610
15 -10.835 -1.462 -2.015
16 -16.163 -2.050 -2.646
17 -20.155 -2.120 -3.323
18 -28.106 -2.503 -4.075
19 -34.954 -2.849 -5.626
20 -44.316 -3.805 -6.514
task spaces: {'initial vel.': np.float64(-0.048), 'endpoints + max. spd.': np.float64(-0.143), 'endpoints + orth. exc.': np.float64(-0.126), 'endpoints + viapoint': np.float64(-0.138), 'endpoints + max. vel.': np.float64(-0.185), 'endpoints + halfpoint': np.float64(-0.126), 'segment dirs.': np.float64(-0.176), 'kitchen sink': np.float64(-0.573)}
median over units, k = 1..6: {'peri': array([-0.104, -0.103, -0.141, -0.238, -0.73 , -1.091]), 'kin800': array([-0.088, -0.103, -0.138, -0.141, -0.177, -0.219]), 'upper': array([-0.077, -0.096, -0.097, -0.092, -0.116, -0.161])}
units with peri performance > 0 at their best k <= 14: 47 of 93
O1 SD units (n = 8): peri best 0.209, kin best 0.136, upper best 0.123
points below the axis (< -1.5): 26 of 68; plotted values range -44.32 to -0.04
Spread of the fig6d comparisons, a second reading of the T matrix, and a noise reference.
- Per reading, the peri-movement performance per unit (at the k best on average) is compared with the best other space per unit by a paired t-test, with SEs across units.
- The Methods say each row of the T matrix "contained data from all neurons except the neuron whose preparatory activity we were trying to capture", and that "Analysis of preparatory responses was restricted to neurons with at least 10 spikes/s of preparatory tuning" (5 for the arrays). "All neurons" may be all recorded units (main, 142) or the included units only; both are run.
- Condition-shuffled null (main reading, 10 shuffles): the neural condition labels are permuted, the kinematic and task spaces are left intact, and the same comparison is made.
rows = []
for key in CONFIGS:
res = PERF[key]
for O1 in ("range", "sd"):
sel = mp.tuning(CONFIGS[key]["P"], O1) > mp.PREP_TUNING_MIN
for kmax in (14, 20):
for end in (600, 800):
for seg in ("length", "time"):
peri = res["peri"][sel, :kmax]
pu = peri[:, np.argmax(np.nanmean(peri, 0))]
kin = res[f"kin{end}"][sel, :kmax]
others = {"kinematic": kin[:, np.argmax(np.nanmean(kin, 0))],
**{name: res[f"{seg}:{name}"][sel] for name in mp.TASK_SPACES}}
best_o = max(others, key=lambda k: np.nanmean(others[k]))
d = pu - others[best_o]
rows.append({"O3": key[0], "O4": key[1], "O1": O1, "kmax": kmax, "kin end": end, "segments": seg,
"peri": np.nanmean(pu), "SE peri": np.nanstd(pu) / np.sqrt(len(pu)), "best other": best_o,
"other": np.nanmean(others[best_o]), "peri - other": np.nanmean(d),
"SE diff": np.nanstd(d) / np.sqrt(len(d)), "p paired t": stats.ttest_rel(pu, others[best_o], nan_policy="omit").pvalue})
f6_tests = pd.DataFrame(rows)
print(f6_tests.drop_duplicates(subset=["O3", "O4", "O1", "best other", "peri", "other"]).round(3).to_string(index=False))
print("(rows identical in all values but kmax, kin end and segments are shown once)")
print("readings where the peri-movement space is better than the best other space at p < 0.05: %d of %d; p range %.2f to %.2f"
% (((f6_tests["p paired t"] < 0.05) & (f6_tests["peri - other"] > 0)).sum(), len(f6_tests), f6_tests["p paired t"].min(), f6_tests["p paired t"].max()))
rows = []
for key, cfg in CONFIGS.items():
for O1 in ("range", "sd"):
sel = mp.tuning(cfg["P"], O1) > mp.PREP_TUNING_MIN
sc = mp.leave_one_out_scores(cfg["CMm"][:, :, TSEL], KMAX, members=sel)
pm = np.nanmean(perf_matrix(cfg["P"], sc, np.flatnonzero(sel)), 0)
kin = np.nanmean(PERF[key]["kin800"][sel], 0).max()
tasks = {n: np.nanmean(PERF[key][f"length:{n}"][sel]) for n in mp.TASK_SPACES}
best_t = max(tasks, key=tasks.get)
n_beat = int(kin >= pm.max()) + sum(v >= pm.max() for v in tasks.values())
rows.append({"O3": key[0], "O4": key[1], "O1": O1, "n": int(sel.sum()),
"peri best, T of all units": np.nanmean(PERF[key]["peri"][sel], 0).max(),
"peri best, T of included units": pm.max(), "k": int(np.argmax(pm) + 1), "kin best": kin,
"best task": tasks[best_t], "spaces >= peri (included)": n_beat,
"outcome (included)": "reproduced" if n_beat == 0 else ("partial" if n_beat == 1 else "no")})
memb_table = pd.DataFrame(rows)
print()
print("peri-movement T matrix from all recorded units or from the included units only (kmax 20, kin end 800, segments by length)")
print(memb_table.round(3).to_string(index=False))
rng_f6 = np.random.default_rng(8)
ks_main = mp.pca_scores(mp.kinematic_matrix(CONFIGS[MAIN]["KIN"], 800), KMAX)
rows = []
for s in range(10):
lab = conds.copy(); lab[kp_main] = rng_f6.permutation(conds[kp_main])
d = means_from("a", kp_main, kp_main, lab, range(C))
sel = mp.tuning(d["P"], "range") > mp.PREP_TUNING_MIN
sc = mp.leave_one_out_scores(d["CMm"][:, :, TSEL], 14)
pm = np.nanmean(perf_matrix(d["P"], sc, np.flatnonzero(sel), 14), 0)
km = np.nanmean(perf_matrix(d["P"], [ks_main] * len(units), np.flatnonzero(sel), 14), 0)
iv = np.nanmean([mp.loo_performance(d["P"][:, j], CONFIGS[MAIN]["TS_length"]["initial vel."]) for j in np.flatnonzero(sel)])
rows.append({"shuffle": s + 1, "n": int(sel.sum()), "peri best": pm.max(), "k": int(np.argmax(pm) + 1), "kin best": km.max(), "initial vel.": iv})
f6_null = pd.DataFrame(rows)
print()
print("condition-shuffled null, main reading (data: peri %.3f, kinematic %.3f, initial vel. %.3f)" % (peri6.max(), kin6.max(), tasks6["initial vel."]))
print(f6_null.round(3).to_string(index=False))
print("null peri best %.3f to %.3f, kinematic best %.3f to %.3f" % (f6_null["peri best"].min(), f6_null["peri best"].max(),
f6_null["kin best"].min(), f6_null["kin best"].max()))
O3 O4 O1 kmax kin end segments peri SE peri best other other peri - other SE diff p paired t a a range 14 600 length -0.044 0.021 initial vel. -0.048 0.004 0.018 0.813 a a sd 14 600 length 0.209 0.097 kinematic 0.138 0.072 0.125 0.608 a a sd 14 800 length 0.209 0.097 kinematic 0.136 0.073 0.124 0.596 a b range 14 600 length -0.045 0.022 initial vel. -0.041 -0.004 0.018 0.815 a b sd 14 600 length 0.202 0.091 kinematic 0.135 0.068 0.147 0.680 a b sd 14 800 length 0.202 0.091 kinematic 0.133 0.070 0.146 0.669 a c range 14 600 length -0.023 0.022 initial vel. -0.036 0.013 0.018 0.492 a c sd 14 600 length 0.247 0.053 kitchen sink 0.228 0.018 0.079 0.836 a c sd 14 600 time 0.247 0.053 kitchen sink 0.186 0.060 0.110 0.628 b a range 14 600 length -0.046 0.021 initial vel. -0.042 -0.004 0.018 0.820 b a sd 14 600 length 0.140 0.132 kinematic 0.135 0.005 0.123 0.968 b a sd 14 800 length 0.140 0.132 kinematic 0.133 0.007 0.123 0.958 b b range 14 600 length -0.043 0.023 initial vel. -0.033 -0.010 0.018 0.572 b b sd 14 600 length 0.129 0.125 kinematic 0.134 -0.005 0.104 0.965 b b sd 14 800 length 0.129 0.125 kinematic 0.132 -0.003 0.103 0.982 b c range 14 600 length -0.026 0.022 initial vel. -0.029 0.003 0.018 0.885 b c sd 14 600 length 0.249 0.058 kitchen sink 0.217 0.032 0.078 0.720 b c sd 14 600 time 0.249 0.058 kinematic 0.187 0.062 0.081 0.504 b c sd 14 800 time 0.249 0.058 kinematic 0.183 0.066 0.079 0.470 (rows identical in all values but kmax, kin end and segments are shown once) readings where the peri-movement space is better than the best other space at p < 0.05: 0 of 96; p range 0.47 to 0.98
peri-movement T matrix from all recorded units or from the included units only (kmax 20, kin end 800, segments by length) O3 O4 O1 n peri best, T of all units peri best, T of included units k kin best best task spaces >= peri (included) outcome (included) a a range 93 -0.044 -0.044 2 -0.053 -0.048 0 reproduced a a sd 8 0.209 -0.007 2 0.136 0.133 7 no a b range 91 -0.045 -0.045 2 -0.050 -0.041 1 partial a b sd 8 0.202 0.024 2 0.133 0.117 7 no a c range 91 -0.023 -0.026 2 -0.043 -0.036 0 reproduced a c sd 7 0.247 0.085 1 0.182 0.228 6 no b a range 93 -0.046 -0.046 2 -0.053 -0.042 1 partial b a sd 8 0.140 -0.022 2 0.133 0.128 8 no b b range 91 -0.043 -0.043 2 -0.049 -0.033 1 partial b b sd 8 0.129 0.042 2 0.132 0.121 6 no b c range 91 -0.026 -0.029 2 -0.043 -0.029 1 partial b c sd 7 0.249 0.101 2 0.183 0.217 6 no
condition-shuffled null, main reading (data: peri -0.044, kinematic -0.053, initial vel. -0.048)
shuffle n peri best k kin best initial vel.
1 90 -0.214 1 -0.108 -0.158
2 89 -0.192 1 -0.108 -0.160
3 94 -0.140 1 -0.130 -0.167
4 92 -0.169 1 -0.121 -0.163
5 95 -0.126 1 -0.113 -0.153
6 91 -0.188 1 -0.133 -0.184
7 90 -0.208 1 -0.118 -0.167
8 91 -0.137 1 -0.129 -0.173
9 86 -0.092 1 -0.138 -0.178
10 91 -0.148 1 -0.123 -0.177
null peri best -0.214 to -0.092, kinematic best -0.138 to -0.108
Result: partially reproduced, resting on the assumption O1a. In the main reading the best peri-movement performance (-0.044 at k = 2) is above the best kinematic (-0.053) and every task space (best, initial velocity, -0.048). Over the 96 readings 64 meet the criterion and 32 are partial: under O1a with the condition mean as proto-trial (O4b), or with the speed-based onset (O3b) and a filter, and under O1b with O3b and O4b, the initial-velocity or kinematic space matches or beats it. The criterion asks that it hold under each reading, so the outcome is partial. The differences are small against their spread: in no reading is the peri-movement space better than the best other space at p < 0.05 (paired t p 0.47 to 0.98; main reading difference 0.004, SE 0.018). With the T matrix built from the included units only, O1a gives the same (reproduced or partial in each of 6 settings), while under O1b the peri-movement space falls to -0.022 to 0.101 and 6 to 8 spaces match or beat it. Under O1a every space performs below zero, worse than predicting the mean rate, but above the condition-shuffled null (peri-movement -0.214 to -0.092, kinematic -0.138 to -0.108), so the data carry some tuning these spaces capture.
fig6abc-perimovement-best-other-datasets¶
Blocked, data_missing: single-electrode datasets and EMG are not in the dandiset.
fig6-perimovement-best-each-dimensionality¶
Taken one dimensionality at a time, performance is almost always higher in the peri-movement space than in other spaces of the same dimensionality.
Criterion: for each k under the fig6d readings, the peri-movement space against the kinematic space at the same k and every task space of that dimensionality; at least 80% of comparisons won.
rows = []
for key in CONFIGS:
for O1 in ("range", "sd"):
for kmax in (14, 20):
for end in (600, 800):
for seg in ("length", "time"):
sel, peri, kin, tasks = fig6_summary(key, O1, kmax, end, seg)
wins = total = 0
for k in range(1, kmax + 1):
total += 1; wins += peri[k - 1] > kin[k - 1]
for name, v in tasks.items():
if dims[name] == k:
total += 1; wins += peri[k - 1] > v
rows.append({"O3": key[0], "O4": key[1], "O1": O1, "kmax": kmax, "kin end": end, "segments": seg,
"wins": wins, "comparisons": total, "fraction": wins / total})
each_table = pd.DataFrame(rows)
print(each_table.round(3).to_string(index=False))
print("readings with fraction >= 0.8:", int((each_table.fraction >= 0.8).sum()), "of", len(each_table),
"| fraction range %.2f to %.2f" % (each_table.fraction.min(), each_table.fraction.max()))
O3 O4 O1 kmax kin end segments wins comparisons fraction a a range 14 600 length 4 22 0.182 a a range 14 600 time 4 22 0.182 a a range 14 800 length 4 22 0.182 a a range 14 800 time 4 22 0.182 a a range 20 600 length 4 28 0.143 a a range 20 600 time 4 28 0.143 a a range 20 800 length 4 28 0.143 a a range 20 800 time 4 28 0.143 a a sd 14 600 length 7 22 0.318 a a sd 14 600 time 7 22 0.318 a a sd 14 800 length 8 22 0.364 a a sd 14 800 time 8 22 0.364 a a sd 20 600 length 7 28 0.250 a a sd 20 600 time 7 28 0.250 a a sd 20 800 length 8 28 0.286 a a sd 20 800 time 8 28 0.286 a b range 14 600 length 2 22 0.091 a b range 14 600 time 2 22 0.091 a b range 14 800 length 2 22 0.091 a b range 14 800 time 2 22 0.091 a b range 20 600 length 2 28 0.071 a b range 20 600 time 2 28 0.071 a b range 20 800 length 2 28 0.071 a b range 20 800 time 2 28 0.071 a b sd 14 600 length 7 22 0.318 a b sd 14 600 time 7 22 0.318 a b sd 14 800 length 8 22 0.364 a b sd 14 800 time 8 22 0.364 a b sd 20 600 length 7 28 0.250 a b sd 20 600 time 7 28 0.250 a b sd 20 800 length 8 28 0.286 a b sd 20 800 time 8 28 0.286 a c range 14 600 length 5 22 0.227 a c range 14 600 time 5 22 0.227 a c range 14 800 length 5 22 0.227 a c range 14 800 time 5 22 0.227 a c range 20 600 length 5 28 0.179 a c range 20 600 time 5 28 0.179 a c range 20 800 length 5 28 0.179 a c range 20 800 time 5 28 0.179 a c sd 14 600 length 6 22 0.273 a c sd 14 600 time 6 22 0.273 a c sd 14 800 length 6 22 0.273 a c sd 14 800 time 6 22 0.273 a c sd 20 600 length 6 28 0.214 a c sd 20 600 time 6 28 0.214 a c sd 20 800 length 6 28 0.214 a c sd 20 800 time 6 28 0.214 b a range 14 600 length 3 22 0.136 b a range 14 600 time 3 22 0.136 b a range 14 800 length 3 22 0.136 b a range 14 800 time 3 22 0.136 b a range 20 600 length 3 28 0.107 b a range 20 600 time 3 28 0.107 b a range 20 800 length 3 28 0.107 b a range 20 800 time 3 28 0.107 b a sd 14 600 length 7 22 0.318 b a sd 14 600 time 7 22 0.318 b a sd 14 800 length 8 22 0.364 b a sd 14 800 time 8 22 0.364 b a sd 20 600 length 7 28 0.250 b a sd 20 600 time 7 28 0.250 b a sd 20 800 length 8 28 0.286 b a sd 20 800 time 8 28 0.286 b b range 14 600 length 2 22 0.091 b b range 14 600 time 2 22 0.091 b b range 14 800 length 2 22 0.091 b b range 14 800 time 2 22 0.091 b b range 20 600 length 2 28 0.071 b b range 20 600 time 2 28 0.071 b b range 20 800 length 2 28 0.071 b b range 20 800 time 2 28 0.071 b b sd 14 600 length 4 22 0.182 b b sd 14 600 time 4 22 0.182 b b sd 14 800 length 4 22 0.182 b b sd 14 800 time 4 22 0.182 b b sd 20 600 length 4 28 0.143 b b sd 20 600 time 4 28 0.143 b b sd 20 800 length 4 28 0.143 b b sd 20 800 time 4 28 0.143 b c range 14 600 length 5 22 0.227 b c range 14 600 time 5 22 0.227 b c range 14 800 length 6 22 0.273 b c range 14 800 time 6 22 0.273 b c range 20 600 length 5 28 0.179 b c range 20 600 time 5 28 0.179 b c range 20 800 length 6 28 0.214 b c range 20 800 time 6 28 0.214 b c sd 14 600 length 6 22 0.273 b c sd 14 600 time 6 22 0.273 b c sd 14 800 length 6 22 0.273 b c sd 14 800 time 6 22 0.273 b c sd 20 600 length 7 28 0.250 b c sd 20 600 time 7 28 0.250 b c sd 20 800 length 6 28 0.214 b c sd 20 800 time 6 28 0.214 readings with fraction >= 0.8: 0 of 96 | fraction range 0.07 to 0.36
Result: not reproduced. Under no reading does the peri-movement space win 80% of the same-dimension comparisons; the fraction is 0.07 to 0.36 over 96 readings. Its performance falls steeply with k (from -0.104 at k = 3 to -1.369 at k = 6 in the main reading) while the kinematic space falls slowly.
fig6-maze-endpoint-speed-space-poor¶
In the maze task the endpoints + maximum speed space, which did well for the speed task, does not perform well.
Criterion: the 3-D endpoints + max. speed performance below the 3-D peri-movement performance by at least 0.1 and below half the best peri-movement performance over k.
rows = []
for key in CONFIGS:
for O1 in ("range", "sd"):
for kmax in (14, 20):
sel, peri, kin, tasks = fig6_summary(key, O1, kmax, 800, "length")
e = tasks["endpoints + max. spd."]
d = PERF[key]["peri"][sel, 2] - PERF[key]["length:endpoints + max. spd."][sel]
rows.append({"O3": key[0], "O4": key[1], "O1": O1, "kmax": kmax, "endpoints+max spd": e, "peri k=3": peri[2],
"diff": peri[2] - e, "SE diff": np.nanstd(d) / np.sqrt(len(d)), "p paired t": stats.ttest_rel(PERF[key]["peri"][sel, 2], PERF[key]["length:endpoints + max. spd."][sel], nan_policy="omit").pvalue,
"best peri": peri.max(), "below peri3 by >= 0.1": e <= peri[2] - 0.1, "below half best": e < 0.5 * peri.max(),
"pass": (e <= peri[2] - 0.1) and (e < 0.5 * peri.max())})
ems_table = pd.DataFrame(rows)
print(ems_table.round(3).to_string(index=False))
print("readings passing:", int(ems_table["pass"].sum()), "of", len(ems_table))
# O1 SD with the peri-movement T matrix of the included units only (see fig6d)
rows = []
for key, cfg in CONFIGS.items():
sel = mp.tuning(cfg["P"], "sd") > mp.PREP_TUNING_MIN
sc = mp.leave_one_out_scores(cfg["CMm"][:, :, TSEL], 14, members=sel)
pm = np.nanmean(perf_matrix(cfg["P"], sc, np.flatnonzero(sel), 14), 0)
e = np.nanmean(PERF[key]["length:endpoints + max. spd."][sel])
rows.append({"O3": key[0], "O4": key[1], "n": int(sel.sum()), "endpoints+max spd": e, "peri k=3 (included)": pm[2],
"best peri (included)": pm.max(), "pass": (e <= pm[2] - 0.1) and (e < 0.5 * pm.max())})
print("O1 SD, T matrix of the included units:")
print(pd.DataFrame(rows).round(3).to_string(index=False))
O3 O4 O1 kmax endpoints+max spd peri k=3 diff SE diff p paired t best peri below peri3 by >= 0.1 below half best pass a a range 14 -0.143 -0.104 0.039 0.025 0.115 -0.044 False True False a a range 20 -0.143 -0.104 0.039 0.025 0.115 -0.044 False True False a a sd 14 -0.029 0.110 0.140 0.122 0.321 0.209 True True True a a sd 20 -0.029 0.110 0.140 0.122 0.321 0.209 True True True a b range 14 -0.138 -0.127 0.011 0.029 0.707 -0.045 False True False a b range 20 -0.138 -0.127 0.011 0.029 0.707 -0.045 False True False a b sd 14 -0.029 0.101 0.130 0.129 0.379 0.202 True True True a b sd 20 -0.029 0.101 0.130 0.129 0.379 0.202 True True True a c range 14 -0.126 -0.087 0.039 0.023 0.098 -0.023 False True False a c range 20 -0.126 -0.087 0.039 0.023 0.098 -0.023 False True False a c sd 14 0.018 0.202 0.185 0.088 0.099 0.247 True True True a c sd 20 0.018 0.202 0.185 0.088 0.099 0.247 True True True b a range 14 -0.147 -0.093 0.054 0.024 0.028 -0.046 False True False b a range 20 -0.147 -0.093 0.054 0.024 0.028 -0.046 False True False b a sd 14 -0.040 0.113 0.153 0.122 0.279 0.140 True True True b a sd 20 -0.040 0.113 0.153 0.122 0.279 0.140 True True True b b range 14 -0.136 -0.122 0.013 0.030 0.657 -0.043 False True False b b range 20 -0.136 -0.122 0.013 0.030 0.657 -0.043 False True False b b sd 14 -0.032 0.114 0.147 0.132 0.333 0.129 True True True b b sd 20 -0.032 0.114 0.147 0.132 0.333 0.129 True True True b c range 14 -0.126 -0.061 0.065 0.020 0.002 -0.026 False True False b c range 20 -0.126 -0.061 0.065 0.020 0.002 -0.026 False True False b c sd 14 0.017 0.209 0.192 0.093 0.107 0.249 True True True b c sd 20 0.017 0.209 0.192 0.093 0.107 0.249 True True True readings passing: 12 of 24 O1 SD, T matrix of the included units: O3 O4 n endpoints+max spd peri k=3 (included) best peri (included) pass a a 8 -0.029 -0.088 -0.007 False a b 8 -0.029 -0.048 0.024 False a c 7 0.018 0.026 0.085 False b a 8 -0.040 -0.068 -0.022 False b b 8 -0.032 -0.053 0.042 False b c 7 0.017 0.045 0.101 False
Result: not reproduced, resting on the assumption O1a. Under O1a (12 readings) the endpoints + max. speed space (-0.126 to -0.147) is below the 3-D peri-movement space by 0.011 to 0.065, less than the 0.1 required (main: -0.143 against -0.104, SE of the difference 0.025). Under O1b (12 readings) it passes with the T matrix of all recorded units: -0.040 to 0.018 against 0.101 to 0.209. With the T matrix of the 7 or 8 included units, the 3-D peri-movement space gives -0.088 to 0.045 and no O1b reading passes.
fig6-dimensionality-trend¶
Performance of the PCA-based spaces typically improves with dimensionality and then declines above about 10 dimensions.
Criterion: for the peri-movement space over k = 1 to 20, performance at the best k above that at k = 2, best k between 3 and 14, and performance at k = 20 below that at k = 10; partially when only the rise holds.
rows = []
for key in CONFIGS:
for O1 in ("range", "sd"):
sel, peri, kin, tasks = fig6_summary(key, O1, 20, 800, "length")
kb = int(np.argmax(peri) + 1)
rise = peri.max() > peri[1] and 3 <= kb <= 14
fall = peri[19] < peri[9]
rows.append({"O3": key[0], "O4": key[1], "O1": O1, "best k": kb, "perf best": peri.max(), "perf k=2": peri[1],
"perf k=10": peri[9], "perf k=20": peri[19], "rise": rise, "decline": fall,
"kin best k": int(np.argmax(kin) + 1)})
trend_table = pd.DataFrame(rows)
print(trend_table.round(3).to_string(index=False))
print("readings with rise and decline:", int((trend_table.rise & trend_table.decline).sum()), "| rise only:",
int((trend_table.rise & ~trend_table.decline).sum()), "| no rise:", int((~trend_table.rise).sum()), "of", len(trend_table))
O3 O4 O1 best k perf best perf k=2 perf k=10 perf k=20 rise decline kin best k a a range 2 -0.044 -0.044 -4.069 -44.316 False True 1 a a sd 4 0.209 0.122 -1.397 -21.284 True True 3 a b range 2 -0.045 -0.045 -2.879 -35.254 False True 1 a b sd 4 0.202 0.119 -1.976 -26.724 True True 3 a c range 2 -0.023 -0.023 -2.421 -39.978 False True 1 a c sd 2 0.247 0.247 -0.836 -19.801 False True 3 b a range 2 -0.046 -0.046 -3.905 -61.254 False True 1 b a sd 4 0.140 0.123 -1.480 -37.652 True True 3 b b range 2 -0.043 -0.043 -2.783 -42.729 False True 1 b b sd 2 0.129 0.129 -2.032 -23.285 False True 3 b c range 2 -0.026 -0.026 -2.979 -45.220 False True 1 b c sd 2 0.249 0.249 -2.714 -18.377 False True 3 readings with rise and decline: 3 | rise only: 0 | no rise: 9 of 12
Result: not reproduced, resting on the assumption O1a. Under O1a the best k is 2 in all 6 readings, so performance does not rise with dimensionality; the decline above 10 holds (main: -4.069 at k = 10, -44.316 at k = 20). Under O1b the rise and the decline both hold in 3 of 6 readings (best k = 4).
fig6e-summary-rank¶
Blocked, data_missing: Figure 6E combines monkeys A, B and J with EMG spaces, none of which are in
the dandiset. As the criterion asks, the same ranking without EMG is shown on the J-array data: for
limits 3, 6 and 12, the best PCA-based performance with k up to the limit and each task space at or
under it, with SEs across units and paired tests against the best space.
res, sel = PERF[MAIN], mp.tuning(CONFIGS[MAIN]["P"], "range") > 5
for lim in (3, 6, 12):
vals = {"peri-movement": res["peri"][sel, :lim], "kinematic": res["kin800"][sel, :lim]}
per_unit = {}
for name, v in vals.items():
kb = np.argmax(np.nanmean(v, 0))
per_unit[name] = v[:, kb]
for name in mp.TASK_SPACES:
if dims[name] <= lim:
per_unit[name] = res[f"length:{name}"][sel]
means = {k: np.nanmean(v) for k, v in per_unit.items()}
best = max(means, key=means.get)
out = []
for k, v in per_unit.items():
p_t = stats.ttest_rel(per_unit[best], v, nan_policy="omit").pvalue if k != best else np.nan
p_w = stats.wilcoxon(per_unit[best], v).pvalue if k != best else np.nan
out.append({"space": k, "mean": means[k], "SE": np.nanstd(v) / np.sqrt(np.isfinite(v).sum()), "p paired t": p_t, "p Wilcoxon": p_w})
print(f"limit {lim}: best space {best}")
print(pd.DataFrame(out).round(4).to_string(index=False))
limit 3: best space peri-movement
space mean SE p paired t p Wilcoxon
peri-movement -0.0437 0.0211 NaN NaN
kinematic -0.0530 0.0135 0.6192 0.7030
initial vel. -0.0479 0.0179 0.8133 0.8466
endpoints + max. spd. -0.1435 0.0208 0.0000 0.0000
endpoints + orth. exc. -0.1260 0.0219 0.0003 0.0002
limit 6: best space peri-movement
space mean SE p paired t p Wilcoxon
peri-movement -0.0437 0.0211 NaN NaN
kinematic -0.0530 0.0135 0.6192 0.7030
initial vel. -0.0479 0.0179 0.8133 0.8466
endpoints + max. spd. -0.1435 0.0208 0.0000 0.0000
endpoints + orth. exc. -0.1260 0.0219 0.0003 0.0002
endpoints + viapoint -0.1379 0.0251 0.0000 0.0000
endpoints + max. vel. -0.1851 0.0278 0.0000 0.0000
endpoints + halfpoint -0.1256 0.0254 0.0003 0.0002
segment dirs. -0.1756 0.0333 0.0000 0.0000
limit 12: best space peri-movement
space mean SE p paired t p Wilcoxon
peri-movement -0.0437 0.0211 NaN NaN
kinematic -0.0530 0.0135 0.6192 0.7030
initial vel. -0.0479 0.0179 0.8133 0.8466
endpoints + max. spd. -0.1435 0.0208 0.0000 0.0000
endpoints + orth. exc. -0.1260 0.0219 0.0003 0.0002
endpoints + viapoint -0.1379 0.0251 0.0000 0.0000
endpoints + max. vel. -0.1851 0.0278 0.0000 0.0000
endpoints + halfpoint -0.1256 0.0254 0.0003 0.0002
segment dirs. -0.1756 0.0333 0.0000 0.0000
kitchen sink -0.5729 0.0653 0.0000 0.0000
prep-dimensionality-7-to-10¶
Preparatory activity is moderately high dimensional; depending on the dataset 7 to 10 dimensions capture 90% of its variance.
Criterion: PCA of the c x n matrix of condition-mean preparatory rates of the fig6d units; the number of components reaching 90% of the variance; expected 8.65 ± 2.35 (6.3 to 11). Unnormalised main, z-scored variant; bounded by c - 1 = 26.
rows = []
for key, cfg in CONFIGS.items():
for O1 in ("range", "sd"):
sel = mp.tuning(cfg["P"], O1) > mp.PREP_TUNING_MIN
Psel = cfg["P"][:, sel]
n90, ev = mp.n_components_for(Psel)
z = (Psel - Psel.mean(0)) / Psel.std(0)
n90z, _ = mp.n_components_for(z)
rows.append({"O3": key[0], "O4": key[1], "O1": O1, "n units": int(sel.sum()), "PCs for 90%": n90,
"z-scored": n90z, "PC1 share": ev[0], "within 6.3-11": 6.3 <= n90 <= 11})
dim_table = pd.DataFrame(rows)
print(dim_table.round(3).to_string(index=False))
# reference: the same count for a matrix of pure trial noise (conditions shuffled across trials)
rng_dim = np.random.default_rng(1)
cfg = CONFIGS[MAIN]; kp = cfg["kept"]
Pt_k = Pt[kp]; sel = mp.tuning(cfg["P"], "range") > 5
null_n = []
for _ in range(50):
perm = rng_dim.permutation(conds[kp])
null_n.append(mp.n_components_for(mp.condition_means(Pt_k, perm, C)[:, sel])[0])
print("same count with condition labels shuffled across trials (50 shuffles): median %d, range %d to %d"
% (np.median(null_n), min(null_n), max(null_n)))
O3 O4 O1 n units PCs for 90% z-scored PC1 share within 6.3-11 a a range 93 16 18 0.195 False a a sd 8 6 6 0.349 False a b range 91 16 18 0.209 False a b sd 8 5 6 0.363 False a c range 91 16 18 0.222 False a c sd 7 5 5 0.405 False b a range 93 17 18 0.194 False b a sd 8 6 6 0.347 False b b range 91 16 18 0.214 False b b sd 8 5 6 0.367 False b c range 91 16 18 0.222 False b c sd 7 5 5 0.405 False same count with condition labels shuffled across trials (50 shuffles): median 17, range 16 to 17
Result: not reproduced. 90% of the variance needs 16 or 17 components under O1a (91 to 93 units) and 5 or 6 under O1b (7 or 8 units), outside 6.3 to 11 in all 12 readings. With the condition labels shuffled across trials the count is 16 to 17, the same as with the true labels, so under O1a the count measures trial noise in 1 to 4 trials per condition, not the dimensionality of preparatory tuning.
fig7a-shuffled-space-poor¶
A control peri-movement space in which each condition's peri-movement response is inverted at random on half the conditions performs poorly.
Criterion: for each unit and condition with probability 0.5 the activity after -150 ms is reflected around the rate at -150 ms; the shuffled T matrix is the space; per condition for all units as a variant; 20 shuffles. Reproduced when the best shuffled performance over k is below half the best intact performance and below the best kinematic performance.
rng_shuf = np.random.default_rng(2)
t0 = time.time()
def shuffled_perf(cfg, sel, mode, rng, kmax=14):
B = cfg["CMm"][:, :, TSEL].copy() # c x n x t, first sample is -150 ms
flip = rng.random((C, len(units))) < 0.5 if mode == "unit" else np.repeat(rng.random((C, 1)) < 0.5, len(units), 1)
B[flip] = 2 * B[flip][:, :1] - B[flip]
sc = mp.leave_one_out_scores(B, kmax)
return np.nanmean(perf_matrix(cfg["P"], sc, np.flatnonzero(sel), kmax), 0)
rows = []
for key in CONFIGS:
cfg = CONFIGS[key]
for O1 in ("range", "sd"):
sel = mp.tuning(cfg["P"], O1) > mp.PREP_TUNING_MIN
intact = np.nanmean(PERF[key]["peri"][sel, :14], 0).max()
kin = np.nanmean(PERF[key]["kin800"][sel, :14], 0).max()
for mode in ("unit", "condition"):
best = [shuffled_perf(cfg, sel, mode, rng_shuf).max() for _ in range(20)]
rows.append({"O3": key[0], "O4": key[1], "O1": O1, "coin per": mode, "intact best": intact, "kin best": kin,
"shuffled best (mean of 20)": np.mean(best), "min": np.min(best), "max": np.max(best),
"pass": np.mean(best) < 0.5 * intact and np.mean(best) < kin})
shuf_table = pd.DataFrame(rows)
print(shuf_table.round(3).to_string(index=False))
print(f"{time.time() - t0:.1f} s")
O3 O4 O1 coin per intact best kin best shuffled best (mean of 20) min max pass a a range unit -0.044 -0.053 -0.097 -0.153 -0.072 True a a range condition -0.044 -0.053 -0.101 -0.134 -0.068 True a a sd unit 0.209 0.136 -0.056 -0.331 0.070 True a a sd condition 0.209 0.136 -0.065 -0.191 0.033 True a b range unit -0.045 -0.050 -0.090 -0.115 -0.054 True a b range condition -0.045 -0.050 -0.103 -0.131 -0.066 True a b sd unit 0.202 0.133 -0.020 -0.098 0.078 True a b sd condition 0.202 0.133 -0.037 -0.130 0.062 True a c range unit -0.023 -0.043 -0.086 -0.112 -0.061 True a c range condition -0.023 -0.043 -0.103 -0.124 -0.081 True a c sd unit 0.247 0.182 0.041 -0.022 0.112 True a c sd condition 0.247 0.182 -0.011 -0.130 0.139 True b a range unit -0.046 -0.053 -0.086 -0.111 -0.055 True b a range condition -0.046 -0.053 -0.107 -0.135 -0.071 True b a sd unit 0.140 0.133 -0.041 -0.116 0.044 True b a sd condition 0.140 0.133 -0.043 -0.104 0.101 True b b range unit -0.043 -0.049 -0.099 -0.169 -0.070 True b b range condition -0.043 -0.049 -0.101 -0.136 -0.065 True b b sd unit 0.129 0.132 -0.014 -0.255 0.090 True b b sd condition 0.129 0.132 -0.052 -0.122 0.040 True b c range unit -0.026 -0.043 -0.100 -0.174 -0.064 True b c range condition -0.026 -0.043 -0.102 -0.134 -0.070 True b c sd unit 0.249 0.183 0.007 -0.110 0.086 True b c sd condition 0.249 0.183 -0.047 -0.109 0.039 True 19.7 s
Result: reproduced, weakly. In all 24 readings (six O3/O4 settings, O1a and O1b, the coin per unit or per condition) the shuffled space's best performance (mean of 20 shuffles) is below half the intact best and below the best kinematic performance: main reading -0.097 against -0.044 intact and -0.053 kinematic. Under O1a the intact space is itself below zero, so "below half the intact best" is met trivially; the evidence is the O1b readings (7 or 8 units), where the shuffled mean is -0.065 to 0.041 against an intact 0.129 to 0.249.
fig7b-window-start-slide¶
Sliding the start of the peri-movement window later costs little; with the window starting at movement onset, performance is still 97% (monkey J) and 94% (monkeys A, B) of its original value.
Criterion: performance averaged across units at the best k for each start time, window end fixed; the J value under reading (a), J-array; expected ratio 0.97 ± 0.097. Monkeys A and B are not in the dandiset.
Result: not reproduced, resting on the SD reading of tuning (O1b), the only reading under which the ratio is defined. Under O1b the ratio of performance with the window starting at movement onset to that at -150 ms is 0.471 to 0.761 in the 6 readings, outside 0.97 ± 0.097. Under O1a (6 readings) performance at -150 ms is below zero (-0.046 to -0.023), so the ratio is undefined; from there it falls to -0.059 to -0.046 at onset. The monkeys A and B value cannot be computed from this dandiset.
Anterior and posterior sites¶
The electrodes table has no positions, and every unit's electrodes entry points to the PMd group.
The three claims below use array identity as the only anterior-posterior information: PMd anterior,
M1 posterior (Supplemental Figure 3 legend: "the first in caudal PMd ... the second in surface M1").
Which array is which is not in the file: the electrode index of both prefixes points into the PMd
group. Prefix 1 is read as PMd only because the PMd group comes first in the electrodes table. Each
claim's result is shown under both mappings: the regression r with the arrays swapped
(ap-location-no-effect-j-array), the raw performance of each half under each mapping (fig7c), and the
outcome under each mapping (prep-activity-more-anterior).
ap-location-no-effect-j-array¶
For J-array, the preparatory to peri-movement tuning correlation does not depend on anterior-posterior recording location (r = 0.20, p = 0.27).
Criterion: a claim of absence; the fig4d correlations regressed on AP location, reproduced when p >= 0.05. With array identity only, the two-group comparison (two-sided t-test) is the variant the criterion names; regression on a 0/1 location is the same test.
rows = []
for (o, f), cfg in CONFIGS.items():
for O1 in ("range", "sd"):
for O2 in ("i", "ii"):
sel = fig4_units(cfg, O1, O2)
if sel.sum() < 4:
continue
r, _ = fig4_corrs(cfg, sel)
ant = (units.array.values == "PMd")[sel].astype(float)
if ant.min() == ant.max():
continue
lr = stats.linregress(ant, r)
tt_ = stats.ttest_ind(r[ant == 1], r[ant == 0])
rows.append({"O3": o, "O4": f, "O1": O1, "O2": O2, "n PMd": int(ant.sum()), "n M1": int((1 - ant).sum()),
"mean r PMd": r[ant == 1].mean(), "mean r M1": r[ant == 0].mean(),
"regression r": lr.rvalue, "p": lr.pvalue, "t-test p": tt_.pvalue})
ap_table = pd.DataFrame(rows)
ap_table["regression r, arrays swapped"] = -ap_table["regression r"] # the 0/1 location is reversed; p is unchanged
print(ap_table.round(3).to_string(index=False))
print("readings with p >= 0.05:", int((ap_table.p >= 0.05).sum()), "of", len(ap_table), "under either mapping of the arrays")
O3 O4 O1 O2 n PMd n M1 mean r PMd mean r M1 regression r p t-test p regression r, arrays swapped a a range i 10 9 -0.060 -0.099 0.105 0.668 0.668 -0.105 a a range ii 43 45 0.016 0.006 0.020 0.854 0.854 -0.020 a a sd ii 4 1 0.095 0.095 0.000 1.000 1.000 -0.000 a b range i 7 6 -0.124 -0.176 0.131 0.671 0.671 -0.131 a b range ii 42 43 0.006 0.020 -0.030 0.787 0.787 0.030 a b sd ii 3 1 -0.060 0.085 -0.207 0.793 0.793 0.207 a c range i 8 7 -0.104 -0.109 0.013 0.965 0.965 -0.013 a c range ii 41 44 0.036 0.035 0.003 0.981 0.981 -0.003 a c sd ii 2 2 0.126 -0.184 0.462 0.538 0.538 -0.462 b a range i 11 7 -0.056 -0.017 -0.101 0.691 0.691 0.101 b a range ii 44 45 -0.003 -0.003 -0.001 0.995 0.995 0.001 b a sd ii 4 1 0.083 -0.033 0.134 0.830 0.830 -0.134 b b range i 9 6 -0.072 -0.073 0.004 0.989 0.989 -0.004 b b range ii 41 45 -0.005 0.016 -0.045 0.679 0.679 0.045 b b sd ii 3 1 -0.081 0.064 -0.221 0.779 0.779 0.221 b c range i 9 6 -0.079 -0.039 -0.107 0.705 0.705 0.107 b c range ii 41 42 0.011 0.010 0.001 0.990 0.990 -0.001 readings with p >= 0.05: 17 of 17 under either mapping of the arrays
Result: reproduced, resting on the assumption that array identity gives AP location. In all 17 readings with units in both arrays the tuning correlation does not differ between the two arrays (p 0.538 to 1.000; main reading r = 0.105, p = 0.668, 10 and 9 units), under either mapping of the arrays: swapping them changes only the sign of r. With 19 units and two locations the test has little power.
ap-location-no-effect-other-datasets¶
Blocked, data_missing: single-electrode penetrations of monkeys A, B and J are not in the dandiset.
fig7c-anterior-posterior-performance¶
The peri-movement space captures preparatory tuning at both anterior and posterior sites, with normalised performance slightly better at anterior sites.
Criterion: a 6-D peri-movement space; units split by AP (here by array); performance normalised by its upper limit (the PCA space of the preparatory matrix) and compared with random halves; the space built from the half's units or from all units. Reproduced when both halves reach 70% of the random-half normalised performance and the anterior half is higher.
rng_ap = np.random.default_rng(3)
cfg = CONFIGS[MAIN]; P = cfg["P"]
prep_sel = mp.tuning(P, "range") > mp.PREP_TUNING_MIN
def half_perf(members, space_from, prep_sel=prep_sel):
"""Mean 6-D peri performance and upper limit (best k) over prep-tuned units in `members`."""
mem_space = members if space_from == "half" else np.ones(len(units), bool)
sc = mp.leave_one_out_scores(cfg["CMm"][:, :, TSEL], 6, members=mem_space)
up = mp.leave_one_out_scores(P[:, :, None], 14, members=mem_space)
idx = np.flatnonzero(members & prep_sel)
perf6 = np.nanmean([mp.loo_performance(P[:, j], sc[j][:, :6]) for j in idx])
upper = np.nanmean(perf_matrix(P, up, idx, 14), 0).max()
return perf6, upper, len(idx)
rows = []
pmd = units.array.values == "PMd"
for space_from in ("half", "all"):
for lab, mem in (("all units", np.ones(len(units), bool)), ("PMd (anterior)", pmd), ("M1 (posterior)", ~pmd)):
p6, up, n = half_perf(mem, space_from)
rows.append({"space from": space_from, "subset": lab, "n tuned": n, "perf k=6": p6, "upper limit": up,
"normalised": p6 / up if up > 0 else np.nan})
rand = []
for _ in range(20):
mem = np.zeros(len(units), bool); mem[rng_ap.choice(len(units), pmd.sum(), replace=False)] = True
p6, up, n = half_perf(mem, space_from)
rand.append((p6, up))
rand = np.array(rand)
rows.append({"space from": space_from, "subset": "random halves (20)", "n tuned": np.nan, "perf k=6": rand[:, 0].mean(),
"upper limit": rand[:, 1].mean(), "normalised": np.nan})
sd_sel = mp.tuning(P, "sd") > mp.PREP_TUNING_MIN
for lab, mem in (("all units, O1 SD units", np.ones(len(units), bool)), ("PMd (anterior), O1 SD units", pmd), ("M1 (posterior), O1 SD units", ~pmd)):
p6, up, n = half_perf(mem, "all", sd_sel)
rows.append({"space from": "all", "subset": lab, "n tuned": n, "perf k=6": p6, "upper limit": up,
"normalised": p6 / up if up > 0 else np.nan})
rand = []
for _ in range(20):
mem = np.zeros(len(units), bool); mem[rng_ap.choice(len(units), pmd.sum(), replace=False)] = True
if (mem & sd_sel).sum():
p6, up, n = half_perf(mem, "all", sd_sel)
rand.append((p6, up, p6 / up if up > 0 else np.nan))
rand = np.array(rand)
rows.append({"space from": "all", "subset": f"random halves ({len(rand)}), O1 SD units", "n tuned": np.nan, "perf k=6": rand[:, 0].mean(),
"upper limit": rand[:, 1].mean(), "normalised": np.nanmean(rand[:, 2])})
ap7_table = pd.DataFrame(rows)
print(ap7_table.round(3).to_string(index=False))
print("upper limits <= 0 make the normalisation undefined:", int((ap7_table["upper limit"] <= 0).sum()), "of", len(ap7_table))
for a_lab, b_lab, name in (("PMd (anterior)", "M1 (posterior)", "prefix 1 = PMd"), ("M1 (posterior)", "PMd (anterior)", "arrays swapped (prefix 2 = PMd)")):
t_ = ap7_table[ap7_table["space from"] == "half"].set_index("subset")["perf k=6"]
print(f"{name}: raw 6-D performance, anterior {t_[a_lab]:.3f}, posterior {t_[b_lab]:.3f} (space from the half's units)")
space from subset n tuned perf k=6 upper limit normalised
half all units 93.0 -1.369 -0.045 NaN
half PMd (anterior) 44.0 -1.535 -0.034 NaN
half M1 (posterior) 49.0 -1.087 -0.097 NaN
half random halves (20) NaN -1.008 -0.059 NaN
all all units 93.0 -1.369 -0.045 NaN
all PMd (anterior) 44.0 -1.454 -0.019 NaN
all M1 (posterior) 49.0 -1.293 -0.069 NaN
all random halves (20) NaN -1.400 -0.046 NaN
all all units, O1 SD units 8.0 -1.144 0.123 -9.277
all PMd (anterior), O1 SD units 4.0 -1.751 0.145 -12.048
all M1 (posterior), O1 SD units 4.0 -0.537 0.184 -2.915
all random halves (20), O1 SD units NaN -0.952 0.156 -7.282
upper limits <= 0 make the normalisation undefined: 8 of 12
prefix 1 = PMd: raw 6-D performance, anterior -1.535, posterior -1.087 (space from the half's units)
arrays swapped (prefix 2 = PMd): raw 6-D performance, anterior -1.087, posterior -1.535 (space from the half's units)
Result: blocked (other): the normalised measure the criterion uses does not exist here. Under O1a the upper-limit performance is at or below zero for every subset (-0.019 to -0.097), so performance cannot be normalised by it. The raw 6-D performance is far below zero at both arrays: -1.535 (prefix 1) and -1.087 (prefix 2) with the space from the half's units, -1.454 and -1.293 from all units, against -1.008 and -1.400 for random halves. Under O1b (4 and 4 units) the upper limits are positive (0.145 and 0.184) but the 6-D performance is not (normalised -12.048 and -2.915, random halves -7.282), so a threshold at 70% of the random-half value has no meaning either. Under O1b the upper limit of all 8 units (0.123, fig6d cell) is also below the best peri-movement performance (0.209), so it is not a ceiling. Which half is anterior is also not known: with prefix 1 as PMd the anterior half is lower in raw performance (-1.535 against -1.087), with the arrays swapped it is higher.
prep-activity-more-anterior¶
Preparatory activity is more prevalent at anterior sites.
Criterion: the fraction of units passing > 5 spikes/s preparatory tuning (O1) higher anterior than posterior, two-sided Fisher exact test at p < 0.05; partially when higher without significance.
rows = []
for (o, f), cfg in CONFIGS.items():
for O1 in ("range", "sd"):
passed = mp.tuning(cfg["P"], O1) > mp.PREP_TUNING_MIN
a_ = passed[pmd].sum(); b_ = passed[~pmd].sum()
p = stats.fisher_exact([[a_, pmd.sum() - a_], [b_, (~pmd).sum() - b_]]).pvalue
fa, fb = a_ / pmd.sum(), b_ / (~pmd).sum()
judge = lambda ant, post: "reproduced" if ant > post and p < 0.05 else ("partial" if ant > post else "not reproduced")
rows.append({"O3": o, "O4": f, "O1": O1, "prefix 1 passing": f"{a_}/{pmd.sum()}", "prefix 2 passing": f"{b_}/{(~pmd).sum()}",
"frac prefix 1": fa, "frac prefix 2": fb, "Fisher p": p,
"outcome if prefix 1 = PMd": judge(fa, fb), "outcome if prefix 2 = PMd": judge(fb, fa)})
prev_table = pd.DataFrame(rows)
print(prev_table.round(4).to_string(index=False))
print(prev_table.groupby("O1")[["outcome if prefix 1 = PMd", "outcome if prefix 2 = PMd"]].agg(lambda x: x.value_counts().to_dict()))
# the same with a measure less sensitive to noise: preparatory range relative to the condition-shuffled range
cfg = CONFIGS[MAIN]; kp = cfg["kept"]
rng_prev = np.random.default_rng(4)
null_rng = np.stack([np.ptp(mp.condition_means(Pt[kp], rng_prev.permutation(conds[kp]), C), axis=0) for _ in range(200)])
z_tune = (np.ptp(cfg["P"], 0) - null_rng.mean(0)) / null_rng.std(0)
print("preparatory range above the 95th percentile of its condition-shuffled null: PMd %d/%d, M1 %d/%d (Fisher p = %.3g); units with no spikes in the epoch give no z"
% ((z_tune[pmd] > 1.645).sum(), pmd.sum(), (z_tune[~pmd] > 1.645).sum(), (~pmd).sum(),
stats.fisher_exact([[(z_tune[pmd] > 1.645).sum(), (z_tune[pmd] <= 1.645).sum()], [(z_tune[~pmd] > 1.645).sum(), (z_tune[~pmd] <= 1.645).sum()]]).pvalue))
print("median preparatory range (spikes/s): PMd %.2f, M1 %.2f; mean rate in the preparatory epoch: PMd %.2f, M1 %.2f"
% (np.median(np.ptp(cfg["P"], 0)[pmd]), np.median(np.ptp(cfg["P"], 0)[~pmd]), cfg["P"].mean(0)[pmd].mean(), cfg["P"].mean(0)[~pmd].mean()))
O3 O4 O1 prefix 1 passing prefix 2 passing frac prefix 1 frac prefix 2 Fisher p outcome if prefix 1 = PMd outcome if prefix 2 = PMd
a a range 44/72 49/70 0.6111 0.7000 0.2931 not reproduced partial
a a sd 4/72 4/70 0.0556 0.0571 1.0000 not reproduced partial
a b range 44/72 47/70 0.6111 0.6714 0.4877 not reproduced partial
a b sd 4/72 4/70 0.0556 0.0571 1.0000 not reproduced partial
a c range 43/72 48/70 0.5972 0.6857 0.2977 not reproduced partial
a c sd 3/72 4/70 0.0417 0.0571 0.7167 not reproduced partial
b a range 44/72 49/70 0.6111 0.7000 0.2931 not reproduced partial
b a sd 4/72 4/70 0.0556 0.0571 1.0000 not reproduced partial
b b range 43/72 48/70 0.5972 0.6857 0.2977 not reproduced partial
b b sd 4/72 4/70 0.0556 0.0571 1.0000 not reproduced partial
b c range 43/72 48/70 0.5972 0.6857 0.2977 not reproduced partial
b c sd 3/72 4/70 0.0417 0.0571 0.7167 not reproduced partial
outcome if prefix 1 = PMd outcome if prefix 2 = PMd
O1
range {'not reproduced': 6} {'partial': 6}
sd {'not reproduced': 6} {'partial': 6}
preparatory range above the 95th percentile of its condition-shuffled null: PMd 7/72, M1 6/70 (Fisher p = 1); units with no spikes in the epoch give no z
median preparatory range (spikes/s): PMd 6.58, M1 10.19; mean rate in the preparatory epoch: PMd 2.40, M1 3.73
/tmp/ipykernel_2266/939575653.py:19: RuntimeWarning: invalid value encountered in divide z_tune = (np.ptp(cfg["P"], 0) - null_rng.mean(0)) / null_rng.std(0)
Result: blocked (metadata_insufficient): the outcome depends on which array is anterior, which the file does not say. Under O1a 44 of 72 prefix-1 units (0.611) and 49 of 70 prefix-2 units (0.700) pass the > 5 spikes/s rule (Fisher p = 0.293). With prefix 1 as PMd the anterior fraction is lower, not reproduced in all 12 readings; with prefix 2 as PMd it is higher without significance, partially reproduced in all 12 (Fisher p 0.293 to 1.000). Under O1b 4 of 72 and 4 of 70 pass. Units whose range is above the 95th percentile of its condition-shuffled null are 7 of 72 and 6 of 70, and with condition labels shuffled 82 to 96 units pass the O1a rule (noise reference), so under O1a these counts mostly measure noise.
fig7d-pd-consistency¶
The PD measured on two random sets of 25 conditions is more stable (higher dot product) in the peri-movement space than in the task, kinematic and EMG spaces, for every comparison and every dataset.
Criterion: per unit, 25 conditions drawn with replacement, PD fitted, another 25 drawn, PD refitted; dot product averaged over units; 100 repeats. Maze comparison spaces: endpoints + max. spd. (3-D), segment dirs. (6-D), and the kinematic space at its best k; the peri-movement space at the same dimensionality. The task-space axes have mixed units (mm, mm/s), and the PD's direction depends on their scale, which the paper does not give: raw units and z-scored axes are both run.
rng_pd = np.random.default_rng(5)
rows = []
for key, cfg in CONFIGS.items():
P = cfg["P"]
ks_kin = mp.pca_scores(mp.kinematic_matrix(cfg["KIN"], 800), KMAX)
for O1 in ("range", "sd"):
sel = np.flatnonzero(mp.tuning(P, O1) > mp.PREP_TUNING_MIN)
kin_k = int(np.argmax(np.nanmean(PERF[key]["kin800"][sel, :14], 0)) + 1)
comps = {"endpoints + max. spd.": cfg["TS_length"]["endpoints + max. spd."],
"segment dirs.": cfg["TS_length"]["segment dirs."], "kinematic (best k)": ks_kin[:, :kin_k]}
for name, S in comps.items():
k = S.shape[1]
for zs in (False, True):
c_other = np.nanmean([mp.pd_consistency(P[:, j], S, rng_pd, standardize=zs) for j in sel])
c_peri = np.nanmean([mp.pd_consistency(P[:, j], PERF[key]["peri_scores"][j][:, :k], rng_pd, standardize=zs) for j in sel])
rows.append({"O3": key[0], "O4": key[1], "O1": O1, "space": name, "k": k, "z-scored axes": zs, "n units": len(sel),
"consistency other": c_other, "consistency peri": c_peri, "peri higher": c_peri > c_other})
pd_table = pd.DataFrame(rows)
print(pd_table.round(3).to_string(index=False))
print("comparisons with peri higher:", int(pd_table["peri higher"].sum()), "of", len(pd_table))
print(pd_table.groupby(["O1", "space"])["peri higher"].agg(["sum", "count"]))
O3 O4 O1 space k z-scored axes n units consistency other consistency peri peri higher
a a range endpoints + max. spd. 3 False 93 0.523 0.617 True
a a range endpoints + max. spd. 3 True 93 0.531 0.630 True
a a range segment dirs. 6 False 93 0.525 0.416 False
a a range segment dirs. 6 True 93 0.527 0.449 False
a a range kinematic (best k) 1 False 93 0.563 0.587 True
a a range kinematic (best k) 1 True 93 0.548 0.591 True
a a sd endpoints + max. spd. 3 False 8 0.576 0.780 True
a a sd endpoints + max. spd. 3 True 8 0.631 0.776 True
a a sd segment dirs. 6 False 8 0.615 0.545 False
a a sd segment dirs. 6 True 8 0.629 0.607 False
a a sd kinematic (best k) 3 False 8 0.672 0.765 True
a a sd kinematic (best k) 3 True 8 0.740 0.784 True
a b range endpoints + max. spd. 3 False 91 0.525 0.617 True
a b range endpoints + max. spd. 3 True 91 0.531 0.631 True
a b range segment dirs. 6 False 91 0.526 0.452 False
a b range segment dirs. 6 True 91 0.526 0.470 False
a b range kinematic (best k) 1 False 91 0.573 0.581 True
a b range kinematic (best k) 1 True 91 0.573 0.586 True
a b sd endpoints + max. spd. 3 False 8 0.601 0.823 True
a b sd endpoints + max. spd. 3 True 8 0.608 0.827 True
a b sd segment dirs. 6 False 8 0.641 0.622 False
a b sd segment dirs. 6 True 8 0.654 0.660 True
a b sd kinematic (best k) 3 False 8 0.666 0.814 True
a b sd kinematic (best k) 3 True 8 0.740 0.839 True
a c range endpoints + max. spd. 3 False 91 0.531 0.634 True
a c range endpoints + max. spd. 3 True 91 0.560 0.644 True
a c range segment dirs. 6 False 91 0.547 0.480 False
a c range segment dirs. 6 True 91 0.536 0.511 False
a c range kinematic (best k) 1 False 91 0.578 0.562 False
a c range kinematic (best k) 1 True 91 0.573 0.581 True
a c sd endpoints + max. spd. 3 False 7 0.646 0.843 True
a c sd endpoints + max. spd. 3 True 7 0.674 0.865 True
a c sd segment dirs. 6 False 7 0.697 0.622 False
a c sd segment dirs. 6 True 7 0.712 0.694 False
a c sd kinematic (best k) 3 False 7 0.726 0.838 True
a c sd kinematic (best k) 3 True 7 0.799 0.858 True
b a range endpoints + max. spd. 3 False 93 0.508 0.618 True
b a range endpoints + max. spd. 3 True 93 0.516 0.625 True
b a range segment dirs. 6 False 93 0.515 0.401 False
b a range segment dirs. 6 True 93 0.512 0.433 False
b a range kinematic (best k) 1 False 93 0.554 0.591 True
b a range kinematic (best k) 1 True 93 0.570 0.576 True
b a sd endpoints + max. spd. 3 False 8 0.602 0.813 True
b a sd endpoints + max. spd. 3 True 8 0.617 0.800 True
b a sd segment dirs. 6 False 8 0.645 0.512 False
b a sd segment dirs. 6 True 8 0.645 0.589 False
b a sd kinematic (best k) 3 False 8 0.653 0.792 True
b a sd kinematic (best k) 3 True 8 0.727 0.788 True
b b range endpoints + max. spd. 3 False 91 0.526 0.622 True
b b range endpoints + max. spd. 3 True 91 0.529 0.638 True
b b range segment dirs. 6 False 91 0.525 0.459 False
b b range segment dirs. 6 True 91 0.530 0.465 False
b b range kinematic (best k) 1 False 91 0.566 0.566 False
b b range kinematic (best k) 1 True 91 0.566 0.580 True
b b sd endpoints + max. spd. 3 False 8 0.593 0.814 True
b b sd endpoints + max. spd. 3 True 8 0.612 0.834 True
b b sd segment dirs. 6 False 8 0.659 0.626 False
b b sd segment dirs. 6 True 8 0.672 0.646 False
b b sd kinematic (best k) 3 False 8 0.650 0.817 True
b b sd kinematic (best k) 3 True 8 0.744 0.825 True
b c range endpoints + max. spd. 3 False 91 0.533 0.654 True
b c range endpoints + max. spd. 3 True 91 0.548 0.656 True
b c range segment dirs. 6 False 91 0.546 0.485 False
b c range segment dirs. 6 True 91 0.544 0.501 False
b c range kinematic (best k) 1 False 91 0.576 0.590 True
b c range kinematic (best k) 1 True 91 0.597 0.566 False
b c sd endpoints + max. spd. 3 False 7 0.644 0.837 True
b c sd endpoints + max. spd. 3 True 7 0.669 0.881 True
b c sd segment dirs. 6 False 7 0.701 0.624 False
b c sd segment dirs. 6 True 7 0.697 0.678 False
b c sd kinematic (best k) 3 False 7 0.714 0.854 True
b c sd kinematic (best k) 3 True 7 0.794 0.875 True
comparisons with peri higher: 46 of 72
sum count
O1 space
range endpoints + max. spd. 12 12
kinematic (best k) 9 12
segment dirs. 0 12
sd endpoints + max. spd. 12 12
kinematic (best k) 12 12
segment dirs. 1 12
Result: not reproduced. Over the six O3/O4 settings, both O1 readings and raw or z-scored axes (72 comparisons) the peri-movement space is more consistent in 46: than endpoints + max. speed in 24 of 24, than the kinematic space at its best k in 21 of 24, but than segment directions in only 1 of 24 (main reading 0.416 against 0.525; O1b with O3a and O4a 0.545 against 0.615).
fig8-linear-dynamics-simulation¶
In a simulated two-unit rotational linear system with 25 random initial states, the single-unit correlation between initial (preparatory) and later (peri-movement) activity is near zero (0.08), yet a PD in the PCA-based peri-movement space captures the initial state perfectly.
Criterion: x(t+1) = W x(t), W = [1, -2pi/360; 2pi/360, 1], 25 random initial states, 1000 simulations with initial states uniform in the unit disk and in a square, a released interval of 360 steps; reproduced when the median |r| (unit 1 initial rate against its mean over the interval) is below 0.3 and the PD in the 2-D PCA space of the released trajectories gives cross-validated performance above 0.99. The expected value 0.08 ± 0.3 is compared with the median r.
rng_sim = np.random.default_rng(6)
rows = []
for init in ("disk", "square"):
rs, perfs = [], []
for _ in range(1000):
x0, traj = mp.simulate_rotation(rng_sim, init=init)
later = traj[:, :, 0].mean(1)
rs.append(np.corrcoef(x0[:, 0], later)[0, 1])
T = traj.reshape(25, -1)
perfs.append(mp.loo_performance(x0[:, 0], mp.pca_scores(T, 2)))
rs, perfs = np.array(rs), np.array(perfs)
rows.append({"initial states": init, "median r": np.median(rs), "median |r|": np.median(np.abs(rs)),
"5th pct r": np.percentile(rs, 5), "95th pct r": np.percentile(rs, 95),
"min PD performance": perfs.min(), "median PD performance": np.median(perfs)})
sim_table = pd.DataFrame(rows)
print(sim_table.round(4).to_string(index=False))
initial states median r median |r| 5th pct r 95th pct r min PD performance median PD performance
disk 0.0044 0.1114 -0.2533 0.2718 1.0 1.0
square 0.0115 0.1473 -0.3312 0.3304 1.0 1.0
Result: reproduced. Over 1000 simulations the median correlation between unit 1's initial rate and its mean over the 360-step released interval is 0.004 (initial states in the unit disk) and 0.012 (in a square), within 0.08 ± 0.3; the median |r| is 0.111 and 0.147, below 0.3; the 2-D PCA-space PD gives leave-one-out performance of 1.0 in every simulation (minimum 1.0).
Summary¶
| Claim | Outcome |
|---|---|
| fig4d-mean-correlation-j-array | not reproduced |
| fig4abc-mean-correlation-other-datasets | blocked: data missing |
| fig4-broad-near-zero | not reproduced (assumption O1a) |
| fig4e-correlation-time-course | reproduced (assumption O1a; noise also passes) |
| suppfig1-population-averages-lose-tuning | reproduced (noise passes in 11 of 20 shuffles) |
| fig6d-perimovement-best-j-array | partially reproduced (assumption O1a) |
| fig6abc-perimovement-best-other-datasets | blocked: data missing |
| fig6-perimovement-best-each-dimensionality | not reproduced |
| fig6-maze-endpoint-speed-space-poor | not reproduced (assumption O1a) |
| fig6-dimensionality-trend | not reproduced (assumption O1a) |
| fig6e-summary-rank | blocked: data missing |
| prep-dimensionality-7-to-10 | not reproduced |
| fig7a-shuffled-space-poor | reproduced (weak: under O1a the intact space is below zero) |
| fig7b-window-start-slide | not reproduced (assumption O1b: under O1a the ratio is undefined) |
| ap-location-no-effect-j-array | reproduced (assumption: array identity as AP location) |
| ap-location-no-effect-other-datasets | blocked: data missing |
| fig7c-anterior-posterior-performance | blocked: other (upper limit at or below zero; normalisation undefined) |
| prep-activity-more-anterior | blocked: metadata insufficient (which array is anterior) |
| fig7d-pd-consistency | not reproduced |
| fig8-linear-dynamics-simulation | reproduced |
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 dandiset holds 27 maze conditions with 1 to 4
trials each after the delay rule (78 trials), where the paper's J-array dataset has 108 conditions and
about 2166 trials per neuron. Condition means are therefore noisy (median split-half reliability of
preparatory tuning 0.20): with condition labels shuffled, about as many units pass the > 5 spikes/s range
rule as with the true labels, and the fig4e and suppfig1 criteria are met by trial noise. Under the range
reading of tuning every PD space performs below zero, though above a shuffled-label null.
Several outcomes differ between the range and SD readings of "spikes/s of tuning"; the SD reading keeps
7 or 8 units, whose preparatory tuning is more reliable over trials (main reading, split-half r median
0.57 against 0.20 over all units, though the units are selected on the same trials). On them every space performs above zero, and the peri-movement space built from all
recorded units has the highest mean performance in 5 of the 6 O3/O4 settings (with O3b and O4b the
kinematic space is higher, 0.134 against 0.129), but no difference from the best other space is
significant (paired t p 0.47 to 0.98); built from those 7 or 8 units only, it is beaten by 6 to 8 spaces. The anterior-posterior claims use array identity, as the electrodes have
no positions, and the file does not say which array is PMd. Task-space details the Methods leave open (end of reach, segments, signed peak velocity)
are stated in maze_prep.py.
What was found about the dataset.
- Every unit's
electrodesentry points to the PMd electrode group, though unit ids carry two array prefixes (72 units with 1, 70 with 2) and the dataset describes recordings from M1 and PMd; no unit points to an M1 electrode, and nothing in the file says which prefix is PMd. - The electrodes table has no positions (x, y, z and impedance are NaN), so no unit has an anterior-posterior location.
- The dandiset holds one session of the J-array kind, 100 trials of 27 conditions; the single-electrode datasets of monkeys A, B and J, the EMG and 81 of the 108 conditions are not in it, and which trials and conditions were chosen is not said beyond "limited to 100 train trials and 100 test trials".
- The 100 test segments carry no condition, go cue or behaviour and only the 107 held-in units, so the paper's condition-wise analyses cannot use them; none duplicates a train trial.
- 22 of the 100 train trials have delays of 400 ms or less, which the paper excludes.
print(f"notebook ran in {(time.time() - T_START) / 60:.1f} min")
notebook ran in 1.4 min