Reproducing Churchland et al. (2012) from dandiset 000070¶
Paper. Churchland MM, Cunningham JP, Kaufman MT, Foster JD, Nuyujukian P, Ryu SI, Shenoy KV (2012). Neural population dynamics during reaching. Nature 487:51-56. doi:10.1038/nature11129. This attempt used the full text of the NIH author manuscript in PubMed Central (PMC3393826) and the Supplementary Information PDF from the publisher's site.
Dandiset. DANDI:000070, version 0.260528.0919. The paper cites no version (it predates the dandiset). Two versions are published, 0.251218.1714 and 0.260528.0919; they hold the same ten files (same blobs and sizes, printed below) and differ only in metadata (contributor affiliations, keywords, study target), none of which the analysis uses. The most recent published version is pinned, as the protocol asks when the paper cites none and no earlier version matches the paper better.
What the dandiset holds. Ten array sessions of monkeys Jenkins (J) and Nitschke (N) performing
the maze task, each with two 96-electrode arrays (PMd and M1). Two are the paper's array datasets:
J-array (2009-09-18) and N-array (2010-09-23). One more (sub-Nitschke_ses-20090910) is
consistent with the "N-9-10-2010" dataset of Supplementary Figure 5 (see the data check). The seven
single-electrode datasets, the EMG, the leech and the walking-monkey data are not in it. The units
are one per electrode (192 per session), not the paper's 146 and 218 sorted isolations.
Code read. The dandiset links its conversion code,
catalystneuro/shenoy-lab-to-nwb at commit
319acb4. The paper links no code; jPCA follows the Methods, the Supplementary Derivation and the
Python port of the lab's MATLAB codepack, bantin/jPCA at commit
3f8be9f (the lab's download page returned HTTP 403). See reach.py for where each parameter comes
from.
| Claim | Tested here |
|---|---|
| array-dataset-composition | yes |
| fig3ef-same-direction | yes |
| fig3-independent-of-curvature | yes |
| fig3-phase-from-prep-state | yes |
| var-captured-jpc-plane | yes |
| pc23-variance | yes |
| suppfig4-j-array-planes | yes |
| suppfig4-planes-range | yes |
| fig6a-angle-peak | yes |
| fig6b-skew-fit | yes |
| rotation-frequency | yes |
| suppfig3-shuffle-controls | yes |
| suppfig5-areas-variance | yes |
| suppfig5-m1-faster | yes |
| fig6-velocity-model | yes (simulated from the dandiset's hand kinematics) |
| fig6-complex-kinematic-model | yes (simulated from the dandiset's hand kinematics) |
| fig1-leech-and-walking | blocked: data missing |
| fig2-single-neuron-responses | blocked: data missing |
| fig3abcd-single-electrode-datasets | blocked: data missing |
| fig4-emg-no-rotation | blocked: data missing |
| fig5-generator-model | blocked: data missing |
| fig6cd-speed-independent-frequency | blocked: data missing |
import time
import warnings
import h5py
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import remfile
from dandi.dandiapi import DandiAPIClient
from scipy import stats
from scipy.optimize import minimize
import reach
pd.set_option("display.width", 200)
pd.set_option("display.max_columns", 30)
pd.set_option("display.max_rows", 200)
warnings.filterwarnings("ignore", category=RuntimeWarning)
started = time.time()
urls = reach.asset_urls()
sessions = {}
for name, path in reach.SESSIONS.items():
t = time.time()
sessions[name] = reach.load_session(path, urls[path], cursor=(name == "J-array"))
print(f"{name}: {path} loaded in {time.time() - t:.0f} s")
print(f"dandiset {reach.DANDISET_ID} version {reach.DANDISET_VERSION}")
J-array: sub-Jenkins/sub-Jenkins_ses-20090918_behavior+ecephys.nwb loaded in 14 s
N-array: sub-Nitschke/sub-Nitschke_ses-20100923_behavior+ecephys.nwb loaded in 40 s
N-0910: sub-Nitschke/sub-Nitschke_ses-20090910_behavior+ecephys.nwb loaded in 18 s dandiset 000070 version 0.260528.0919
Data check¶
Each check of the metadata against the values is stated under its output.
with DandiAPIClient() as client:
ds = client.get_dandiset(reach.DANDISET_ID)
versions = [v.identifier for v in ds.get_versions()]
per_version = {}
for v in versions:
d = client.get_dandiset(reach.DANDISET_ID, v)
per_version[v] = {a.path: (a.blob, a.size) for a in d.get_assets()}
meta = client.get_dandiset(reach.DANDISET_ID, reach.DANDISET_VERSION).get_raw_metadata()
print("Versions:", versions)
for v in versions:
same = per_version[v] == per_version[reach.DANDISET_VERSION]
print(f" {v}: {len(per_version[v])} files, same paths, blobs and sizes as the pinned version: {same}")
print("\nFiles at the pinned version:")
for p, (blob, size) in sorted(per_version[reach.DANDISET_VERSION].items()):
print(f" {p} {size / 1e9:.2f} GB")
summary_meta = meta["assetsSummary"]
print("\nassetsSummary: subjects", summary_meta["numberOfSubjects"], "; species",
[s["name"] for s in summary_meta["species"]], "; files", summary_meta["numberOfFiles"])
print("Description:", meta["description"])
print("Paper link:", [r["identifier"] for r in meta["relatedResource"] if r.get("relation") == "dcite:IsDescribedBy"])
Versions: ['draft', '0.251218.1714', '0.260528.0919'] draft: 10 files, same paths, blobs and sizes as the pinned version: True 0.251218.1714: 10 files, same paths, blobs and sizes as the pinned version: True 0.260528.0919: 10 files, same paths, blobs and sizes as the pinned version: True Files at the pinned version: sub-Jenkins/sub-Jenkins_ses-20090912_behavior+ecephys.nwb 2.71 GB sub-Jenkins/sub-Jenkins_ses-20090916_behavior+ecephys.nwb 3.97 GB sub-Jenkins/sub-Jenkins_ses-20090918_behavior+ecephys.nwb 4.49 GB sub-Jenkins/sub-Jenkins_ses-20090923_behavior+ecephys.nwb 3.77 GB sub-Nitschke/sub-Nitschke_ses-20090812_behavior+ecephys.nwb 7.38 GB sub-Nitschke/sub-Nitschke_ses-20090819_behavior+ecephys.nwb 6.55 GB sub-Nitschke/sub-Nitschke_ses-20090910_behavior+ecephys.nwb 6.60 GB sub-Nitschke/sub-Nitschke_ses-20090920_behavior+ecephys.nwb 5.99 GB sub-Nitschke/sub-Nitschke_ses-20090922_behavior+ecephys.nwb 4.45 GB sub-Nitschke/sub-Nitschke_ses-20100923_behavior+ecephys.nwb 7.43 GB assetsSummary: subjects 2 ; species ['Rhesus monkey', 'Macaca mulatta - Rhesus monkey'] ; files 10 Description: Monkeys recordings of Motor Cortex (M1) and dorsal Premotor Cortex (PMd) using two 96 channel high density Utah Arrays (Blackrock Microsystems) while performing reaching tasks with right hand. Paper link: ['https://doi.org/10.1038/nature11129']
The pinned version and the earlier published one hold the same ten files. The archive's summary counts two subjects, as the files do, and lists the same species twice under two names.
# Every file: subject, date, units, trials, conditions, and two faults found while loading
rows = []
for path in sorted(urls):
f = h5py.File(remfile.File(urls[path]), "r")
tr, u = f["intervals/trials"], f["units"]
starts = tr["start_time"][:]
index = u["spike_times_index"][:].astype(np.int64)
fixed, wraps = reach.fix_first_index(u["spike_times"], index)
cond = tr["maze_condition"][:] if "maze_condition" in tr else None
rows.append({
"file": path.split("/")[1].replace("_behavior+ecephys.nwb", ""),
"session_start": f["session_start_time"][()].decode()[:10],
"units": len(index), "trials": len(starts),
"maze_conditions": int(len(np.unique(cond[cond > 0]))) if cond is not None else None,
"has_correct_reach": "correct_reach" in tr,
"spike_times_layout": "chunked gzip" if u["spike_times"].chunks else "contiguous",
"index0_stored": int(index[0]), "index0_wraps_of_65536": wraps,
"clock_resets": int(np.sum(np.diff(starts) < 0)),
"trial_pairs_overlapping": int(np.sum(np.sort(starts)[1:] < tr["stop_time"][:][np.argsort(starts)][:-1])),
})
files = pd.DataFrame(rows)
print(files.to_string(index=False))
print("\nUnits table description in every file:",
{h5py.File(remfile.File(urls[p]), "r")["units"].attrs.get("description") for p in sorted(urls)[:1]})
file session_start units trials maze_conditions has_correct_reach spike_times_layout index0_stored index0_wraps_of_65536 clock_resets trial_pairs_overlapping sub-Jenkins_ses-20090912 2009-09-12 192 1588 27.0 True contiguous 63597 1 0 0 sub-Jenkins_ses-20090916 2009-09-16 192 2265 27.0 True contiguous 16207 7 0 0 sub-Jenkins_ses-20090918 2009-09-18 192 2636 108.0 True contiguous 61030 3 0 0 sub-Jenkins_ses-20090923 2009-09-23 192 2327 27.0 True contiguous 1018 4 0 0 sub-Nitschke_ses-20090812 2009-08-12 192 3038 27.0 True contiguous 26059 5 1 1064 sub-Nitschke_ses-20090819 2009-08-19 192 3126 27.0 True contiguous 24617 4 1 1064 sub-Nitschke_ses-20090910 2009-09-10 192 2911 108.0 True contiguous 31283 6 2 1415 sub-Nitschke_ses-20090920 2009-09-20 192 4044 NaN False contiguous 15025 5 1 1187 sub-Nitschke_ses-20090922 2009-09-22 192 2612 NaN False contiguous 30162 4 0 0 sub-Nitschke_ses-20100923 2010-09-23 192 3694 108.0 True chunked gzip 456806 0 3 2094
Units table description in every file: {'Autogenerated by NWBFile'}
Two faults of the files, found while loading and corrected in reach.py:
- In nine of the ten files (all but the 2010 file, which a later conversion wrote) the first entry
of
units/spike_times_indexis the true end of unit 0 modulo 65536. Unit 0 is cut short and the rest of its spikes open unit 1, whose times then run from late in the session back to its start. The true end is the one position, at the stored value plus a multiple of 65536, where the times jump back from the last trial to the first;fix_first_indexrestores it (shown for J-array below). This is a format error (the index of a ragged column does not mark the end of its row). - In five of the Nitschke files, including N-array (2010), the trial times go backwards at segment
boundaries (
clock_resets): the conversion code offsets only the first of the recording's segments, so later segments restart on the same clock and their trials overlap in time with earlier ones (trial_pairs_overlapping). Spike times per unit are therefore not sorted, and a spike cannot be placed in a trial by its time alone. The conversion wrote each unit's spikes trial by trial in table order, soassign_to_trialsplaces them by that order; every spike of every unit lands in its trial's interval (it raises otherwise). Position samples are likewise one block per trial.
The units are one per electrode, 192 per file. The paper's J-array and N-array datasets had 146 and 218 isolations; the dandiset does not hold those, and its Units table has no description of what its units are.
for name, s in sessions.items():
tr = s["trials"]
sel = reach.select_trials(tr)
cpc = tr[sel].maze_condition.value_counts()
print(f"{name}: {len(tr)} trials; successful {int((tr.task_success == 1).sum())}; flagged discard "
f"{int((tr.discard_trial == 1).sum())}; delay > 400 ms {int((tr.delay > 0.4).sum())}; "
f"consistent {int(tr.correct_reach.sum())}; novel mazes (condition 0) {int((tr.maze_condition == 0).sum())}")
print(f" selected {int(sel.sum())} trials in {len(cpc)} conditions, {cpc.min()} to {cpc.max()} per condition (median {cpc.median():.0f}); "
f"delay {tr.delay[sel].min():.3f} to {tr.delay[sel].max():.3f} s; "
f"movement onset minus go cue {np.percentile((tr.move_begins_time - tr.go_cue_time)[sel], 5):.3f} to "
f"{np.percentile((tr.move_begins_time - tr.go_cue_time)[sel], 95):.3f} s (5th to 95th percentile)")
print(f" condition 0 trials with novel_maze = 1: {int(((tr.maze_condition == 0) & (tr.novel_maze == 1)).sum())} of {int((tr.maze_condition == 0).sum())}")
print(f" session_start_time {s['meta']['session_start_time']}")
J-array: 2636 trials; successful 2570; flagged discard 0; delay > 400 ms 2278; consistent 2296; novel mazes (condition 0) 203 selected 1980 trials in 108 conditions, 14 to 22 per condition (median 18); delay 0.413 to 0.998 s; movement onset minus go cue 0.287 to 0.488 s (5th to 95th percentile) condition 0 trials with novel_maze = 1: 2 of 203 session_start_time 2009-09-18T00:00:00-07:00 N-array: 3694 trials; successful 3590; flagged discard 0; delay > 400 ms 3297; consistent 3135; novel mazes (condition 0) 419 selected 2791 trials in 108 conditions, 19 to 37 per condition (median 25); delay 0.414 to 0.999 s; movement onset minus go cue 0.316 to 0.464 s (5th to 95th percentile) condition 0 trials with novel_maze = 1: 150 of 419 session_start_time 2010-09-23T00:00:00-07:53 N-0910: 2911 trials; successful 2860; flagged discard 0; delay > 400 ms 2590; consistent 2021; novel mazes (condition 0) 828 selected 1795 trials in 108 conditions, 14 to 24 per condition (median 16); delay 0.414 to 0.998 s; movement onset minus go cue 0.301 to 0.432 s (5th to 95th percentile) condition 0 trials with novel_maze = 1: 129 of 828 session_start_time 2009-09-10T00:00:00-07:00
The Methods' selection (successful, consistent trials with delays above 400 ms) leaves 108 conditions in each session with at least 14 trials each. The delay is the go cue minus target onset; the go-cue-to-movement interval (reaction time) lies within 0.29 to 0.49 s between the 5th and 95th percentiles of every session, in the range of a reaching task. Condition 0 holds the novel mazes and is left out, as the paper analyses 108 conditions.
for name, s in sessions.items():
tr = s["trials"]
dur = (tr.stop_time - tr.start_time).sum()
n_spk = np.array([len(x) for x in s["spikes"]])
rate = n_spk / dur
print(f"{name}: {len(s['spikes'])} units ({pd.Series(s['unit_area']).value_counts().to_dict()}), "
f"electrode groups {pd.Series(s['unit_group']).value_counts().to_dict()}; electrode locations {s['meta']['electrode_locations']}")
print(f" one unit per electrode: {len(set(zip(s['unit_group'], range(len(s['spikes']))))) == len(s['spikes'])}; "
f"rate over trial time: median {np.median(rate):.1f}, range {rate.min():.1f} to {rate.max():.1f} spikes/s")
print(f" obs_intervals equal the trial intervals: {np.allclose(s['obs_first'], tr[['start_time', 'stop_time']].values)}; "
f"stored index wraps corrected: {s['index_wraps']}; spikes placed in their trials: {sum(o[-1] for o in s['spike_offsets'])} of {n_spk.sum()}")
s = sessions["J-array"]
st, fx = s["index_stored"], s["index_fixed"]
flat_counts = {"stored": (int(st[0]), int(st[1] - st[0])), "corrected": (int(fx[0]), int(fx[1] - fx[0]))}
print("\nJ-array spikes of unit 0 and unit 1, as stored and corrected:", flat_counts)
tr = s["trials"]
first_trial = lambda x: int(np.searchsorted(tr.start_time.values, x[0], "right") - 1)
last_trial = lambda x: int(np.searchsorted(tr.start_time.values, x[-1], "right") - 1)
print("Corrected: unit 0 spans trials", first_trial(s["spikes"][0]), "to", last_trial(s["spikes"][0]),
"; unit 1 spans trials", first_trial(s["spikes"][1]), "to", last_trial(s["spikes"][1]),
"; unit 1 times sorted:", bool(np.all(np.diff(s["spikes"][1]) >= 0)))
J-array: 192 units ({'PMd': 96, 'M1': 96}), electrode groups {'1': 96, '2': 96}; electrode locations ['M1 Motor Cortex', 'Pre-Motor Cortex, dorsal']
one unit per electrode: True; rate over trial time: median 32.8, range 15.6 to 178.0 spikes/s
obs_intervals equal the trial intervals: True; stored index wraps corrected: 3; spikes placed in their trials: 66676824 of 66676824
N-array: 192 units ({'PMd': 96, 'M1': 96}), electrode groups {'1': 96, '2': 96}; electrode locations ['M1 Motor Cortex', 'Pre-Motor Cortex, dorsal']
one unit per electrode: True; rate over trial time: median 35.9, range 13.5 to 75.3 spikes/s
obs_intervals equal the trial intervals: True; stored index wraps corrected: 0; spikes placed in their trials: 111572823 of 111572823
N-0910: 192 units ({'PMd': 96, 'M1': 96}), electrode groups {'1': 96, '2': 96}; electrode locations ['M1 Motor Cortex', 'Pre-Motor Cortex, dorsal']
one unit per electrode: True; rate over trial time: median 35.4, range 8.8 to 299.6 spikes/s
obs_intervals equal the trial intervals: True; stored index wraps corrected: 6; spikes placed in their trials: 94386885 of 94386885
J-array spikes of unit 0 and unit 1, as stored and corrected: {'stored': (61030, 461554), 'corrected': (257638, 264946)}
Corrected: unit 0 spans trials 0 to 2635 ; unit 1 spans trials 0 to 2635 ; unit 1 times sorted: True
Each unit sits on its own electrode, 96 on the PMd array (group 1) and 96 on the M1 array (group 2), and fires at 8.8 to 300 spikes/s over trial time (medians 33 to 36). These rates, and one unit per electrode, are those of threshold crossings or of all spikes on a channel, not of sorted single neurons; this is consistent with the data but the files do not say. Every spike of every unit is placed in its own trial. After the correction, unit 0 spans the whole session and unit 1's times are sorted.
for name, s in sessions.items():
tr = s["trials"]
sel = np.where(reach.select_trials(tr))[0]
xy = s["hand_xy"]
print(f"{name}: hand position stated in '{s['meta']['hand_unit']}' (conversion {s['meta']['hand_conversion']}); "
f"x {np.nanmin(xy[:, 0]):.1f} to {np.nanmax(xy[:, 0]):.1f}, y {np.nanmin(xy[:, 1]):.1f} to {np.nanmax(xy[:, 1]):.1f}; "
f"targets x {tr.hit_x.min():.0f} to {tr.hit_x.max():.0f}, y {tr.hit_y.min():.0f} to {tr.hit_y.max():.0f}")
rel, pos, vel, acc, jerk = reach.hand_kinematics(s, sel, tr.move_begins_time.values[sel])
speed = np.median(np.linalg.norm(vel, axis=2), axis=0) * 1000
print(" median hand speed (mm/s, values read as mm) at", {t: round(float(speed[t - reach.WINDOW_MS[0]]), 1) for t in (-200, -100, -50, 0, 50, 100, 200, 400)}, "ms from movement onset")
off = s["hand_offsets"]
ends = tr.move_ends_time.values[sel]
ok = ~np.isnan(ends)
idx = off[sel[ok]] + np.round((ends[ok] - s["hand_t"][off[sel[ok]]]) * 1000).astype(int)
d_hand = xy[idx] - tr[["hit_x", "hit_y"]].values[sel[ok]]
print(f" hand at movement end minus hit target: median x {np.median(d_hand[:, 0]):.1f}, y {np.median(d_hand[:, 1]):.1f}")
if s["cursor_xy"] is not None:
d_cur = s["cursor_xy"][idx] - tr[["hit_x", "hit_y"]].values[sel[ok]]
print(f" cursor at movement end minus hit target: median x {np.median(d_cur[:, 0]):.1f}, y {np.median(d_cur[:, 1]):.1f}; "
f"within the target half-width ({tr.target_size.iloc[0]}): {np.mean(np.all(np.abs(d_cur) <= tr.target_size.values[sel[ok], None], axis=1)):.2f}")
J-array: hand position stated in 'meters' (conversion nan); x -145.8 to 163.6, y -218.8 to 109.1; targets x -130 to 133, y -86 to 86
median hand speed (mm/s, values read as mm) at {-200: 1.2, -100: 1.5, -50: 15.6, 0: 144.8, 50: 517.5, 100: 730.2, 200: 510.8, 400: 56.1} ms from movement onset
hand at movement end minus hit target: median x 1.5, y -35.2
cursor at movement end minus hit target: median x 1.6, y -0.1; within the target half-width (13): 0.94
N-array: hand position stated in 'meters' (conversion nan); x -155.6 to 152.7, y -197.9 to 67.4; targets x -144 to 142, y -108 to 111
median hand speed (mm/s, values read as mm) at {-200: 1.2, -100: 2.3, -50: 14.6, 0: 163.5, 50: 439.4, 100: 623.2, 200: 534.4, 400: 103.4} ms from movement onset
hand at movement end minus hit target: median x 0.2, y -53.3
N-0910: hand position stated in 'meters' (conversion nan); x -165.4 to 151.3, y -206.3 to 69.4; targets x -144 to 145, y -108 to 111
median hand speed (mm/s, values read as mm) at {-200: 1.2, -100: 1.6, -50: 14.2, 0: 180.1, 50: 483.7, 100: 658.6, 200: 542.7, 400: 84.6} ms from movement onset
hand at movement end minus hit target: median x 0.7, y -52.2
Hand position is stated in metres, but spans about -220 to +165 on the same scale as the target positions (up to 145), and a reach peaks at about 600 to 730 per second: the values are in mm, the units of the screen coordinates. The unit attribute is wrong. The hand ends a median 35 mm below the hit target in J-array and 52 to 53 mm below it in the Nitschke sessions, while in J-array the cursor lands on it (94% within the target), so hand and screen coordinates differ by an offset in y that the files do not describe (the conversion code subtracts fixed offsets of 8 and 24 from hand y for the two monkeys). Hand speed is near zero 100 ms before the trial table's movement onset and rises from it in all three sessions, so the hand data share the trial table's clock.
# Clock check with a contrast: decode the reach direction (octant of the hit target) from spike
# counts around movement onset, against the 200 ms before target onset and against the labels of the
# next trial. Chance is about 1/8.
rng = np.random.default_rng(7001)
rows = []
for name, s in sessions.items():
tr = s["trials"]
sel = np.where(reach.select_trials(tr))[0]
y = (np.round(np.arctan2(tr.hit_y.values[sel], tr.hit_x.values[sel]) / (np.pi / 4)) % 8).astype(int)
mv, tg = tr.move_begins_time.values[sel], tr.target_presentation_time.values[sel]
X_move = reach.window_counts(s, sel, mv - 0.1, mv + 0.3)
X_pre = reach.window_counts(s, sel, tg - 0.2, tg)
seg = np.r_[0, np.cumsum(np.diff(tr.start_time.values) < 0)][sel]
row = {"session": name, "trials": len(sel), "largest class": np.bincount(y).max() / len(y),
"movement window": reach.decode_accuracy(X_move, y, rng),
"before target": reach.decode_accuracy(X_pre, y, rng),
"next trial's labels": reach.decode_accuracy(X_move[:-1], y[1:], rng)}
for g in np.unique(seg):
m = seg == g
row[f"clock segment {g} ({m.sum()} trials)"] = reach.decode_accuracy(X_move[m], y[m], rng)
rows.append(row)
print(pd.DataFrame(rows).set_index("session").T.round(3).to_string())
session J-array N-array N-0910 trials 1980.000 2791.000 1795.000 largest class 0.230 0.170 0.175 movement window 0.748 0.842 0.841 before target 0.133 0.142 0.130 next trial's labels 0.159 0.171 0.144 clock segment 0 (1980 trials) 0.748 NaN NaN clock segment 0 (1350 trials) NaN 0.834 NaN clock segment 1 (691 trials) NaN 0.865 NaN clock segment 2 (713 trials) NaN 0.849 NaN clock segment 3 (37 trials) NaN 0.514 NaN clock segment 0 (1018 trials) NaN NaN 0.827 clock segment 1 (489 trials) NaN NaN 0.836 clock segment 2 (288 trials) NaN NaN 0.826
In every session and in every clock segment of the 2010 and 20090910 files, the reach direction is decoded from the spikes around the trial table's movement onset far above chance, while the 200 ms before target onset and the next trial's labels give about chance. Spikes and trial events share a clock in every segment, and the placement of spikes by storage order puts them in the right trials.
# Which session is the "N-9-10-2010" dataset of Supplementary Figure 5? Compare the reached target
# of each maze condition (the most common hit target over its successful trials).
def condition_targets(s):
tr = s["trials"]
ok = (tr.task_success == 1) & (tr.maze_condition > 0)
return tr[ok].groupby("maze_condition")[["hit_x", "hit_y"]].agg(lambda v: v.mode().iloc[0])
targets = {name: condition_targets(s) for name, s in sessions.items()}
for a, b in [("N-0910", "N-array"), ("N-0910", "J-array"), ("N-array", "J-array")]:
common = targets[a].index.intersection(targets[b].index)
same = (targets[a].loc[common] == targets[b].loc[common]).all(axis=1)
print(f"{a} against {b}: {len(common)} conditions in both; same reached target {same.mean():.2f}")
print("Session dates in the files:", {n: s["meta"]["session_start_time"][:10] for n, s in sessions.items()})
N-0910 against N-array: 108 conditions in both; same reached target 1.00
N-0910 against J-array: 108 conditions in both; same reached target 0.00
N-array against J-array: 108 conditions in both; same reached target 0.00
Session dates in the files: {'J-array': '2009-09-18', 'N-array': '2010-09-23', 'N-0910': '2009-09-10'}
sub-Nitschke_ses-20090910 has 108 conditions, and for all of them the same reached target
as N-array (2010), while neither shares any with J-array. The other Nitschke sessions of 2009 have
27 conditions or none. The conversion code
that wrote it (maze_task/conversion_script.py) builds every session date as "2009" plus the folder
name, so a session recorded on 10 September 2010 would be dated 2009-09-10. This is consistent with
the file being the "N-9-10-2010" dataset of Supplementary Figure 5, which is otherwise absent, but it
is not shown; it is used here only as a variant.
Firing rates and jPCA¶
Shared by the claims below. Rates are averaged per condition over the selected trials, aligned on movement onset, and smoothed with a 20 ms Gaussian. jPCA uses -50 to +150 ms every 10 ms and 6 PCs.
prep = {}
for name, s in sessions.items():
tr = s["trials"]
sel = np.where(reach.select_trials(tr))[0]
align = tr.move_begins_time.values[sel]
cond = tr.maze_condition.values[sel]
labels, sums, n = reach.condition_counts(s, sel, align, cond)
prep[name] = {"sel": sel, "align": align, "cond": cond, "labels": labels, "sums": sums, "n": n,
"rates": reach.condition_rates(sums, n)}
res = {name: reach.jpca(p["rates"]) for name, p in prep.items()}
main = pd.DataFrame({name: reach.summary(r) for name, r in res.items()}).T
print(main.round(3).to_string())
print("Rates (spikes/s) over all conditions and times: ",
{n: (round(float(p["rates"].min()), 1), round(float(p["rates"].max()), 1)) for n, p in prep.items()})
plane_var pc1_var pc23_var pcs_var r2_ratio_6d r2_ratio_plane r2_skew r2_best angle_peak angle_median frac_ccw freq_hz
J-array 0.200 0.237 0.216 0.605 0.905 0.957 0.343 0.378 1.658 1.606 0.944 1.033
N-array 0.231 0.204 0.270 0.660 0.845 0.881 0.517 0.612 2.007 1.896 1.000 1.377
N-0910 0.188 0.171 0.241 0.591 0.874 0.932 0.526 0.602 1.833 1.792 0.981 1.524
Rates (spikes/s) over all conditions and times: {'J-array': (1.3, 303.1), 'N-array': (0.7, 266.7), 'N-0910': (0.1, 466.9)}
# Variants used in several claims: 24 ms smoothing, other windows, all successful trials (not only
# consistent ones), and 10 PCs.
variants = {}
for name in ("J-array", "N-array"):
p, s = prep[name], sessions[name]
v = {"main": reach.summary(res[name]),
"24 ms Gaussian": reach.summary(reach.jpca(reach.condition_rates(p["sums"], p["n"], 24.0))),
"window -100 to +100 ms": reach.summary(reach.jpca(p["rates"], times_ms=np.arange(-100, 101, 10))),
"window 0 to +200 ms": reach.summary(reach.jpca(p["rates"], times_ms=np.arange(0, 201, 10))),
"10 PCs": reach.summary(reach.jpca(p["rates"], n_pcs=10))}
tr = s["trials"]
sel_all = np.where(reach.select_trials(tr, consistent=False))[0]
lab, sums_all, n_all = reach.condition_counts(s, sel_all, tr.move_begins_time.values[sel_all], tr.maze_condition.values[sel_all])
v[f"all successful trials ({len(sel_all)})"] = reach.summary(reach.jpca(reach.condition_rates(sums_all, n_all)))
variants[name] = pd.DataFrame(v).T
print(name)
print(variants[name].round(3).to_string(), "\n")
J-array
plane_var pc1_var pc23_var pcs_var r2_ratio_6d r2_ratio_plane r2_skew r2_best angle_peak angle_median frac_ccw freq_hz
main 0.200 0.237 0.216 0.605 0.905 0.957 0.343 0.378 1.658 1.606 0.944 1.033
24 ms Gaussian 0.209 0.245 0.222 0.622 0.905 0.963 0.378 0.418 1.833 1.636 0.981 1.009
window -100 to +100 ms 0.175 0.222 0.217 0.598 0.886 0.963 0.443 0.500 1.658 1.650 0.972 1.446
window 0 to +200 ms 0.207 0.240 0.223 0.614 0.852 0.953 0.332 0.390 1.658 1.543 0.944 1.040
10 PCs 0.123 0.237 0.216 0.714 0.845 0.887 0.456 0.539 1.658 1.690 0.963 1.682
all successful trials (2045) 0.190 0.239 0.217 0.607 0.908 0.950 0.351 0.387 1.833 1.593 0.944 1.063
N-array
plane_var pc1_var pc23_var pcs_var r2_ratio_6d r2_ratio_plane r2_skew r2_best angle_peak angle_median frac_ccw freq_hz
main 0.231 0.204 0.270 0.660 0.845 0.881 0.517 0.612 2.007 1.896 1.000 1.377
24 ms Gaussian 0.235 0.209 0.276 0.673 0.844 0.881 0.542 0.642 2.007 1.907 1.000 1.360
window -100 to +100 ms 0.248 0.200 0.303 0.664 0.890 0.935 0.546 0.614 2.007 1.795 0.981 1.428
window 0 to +200 ms 0.208 0.199 0.235 0.645 0.819 0.887 0.492 0.600 1.833 1.850 0.981 1.393
10 PCs 0.162 0.204 0.270 0.753 0.843 0.913 0.561 0.666 1.833 1.822 0.972 1.639
all successful trials (2826) 0.231 0.206 0.269 0.662 0.842 0.879 0.515 0.611 2.007 1.892 1.000 1.372
# Shuffle controls (Supplementary Figures 2 and 3), used by suppfig3-shuffle-controls and
# fig3-phase-from-prep-state. Dividing point 50 ms after the go cue: the median go-cue time
# relative to movement onset plus 50 ms.
rng_shuffle = np.random.default_rng(7012)
shuffles = {}
for name in ("J-array", "N-array"):
p, tr = prep[name], sessions[name]["trials"]
go_rel = np.median((tr.go_cue_time - tr.move_begins_time).values[p["sel"]]) * 1000
t0 = int(round(go_rel + 50))
t0_index = t0 - reach.WINDOW_MS[0]
out = []
for kind, reps in ((1, 20), (2, 1), (3, 20)):
for r in range(reps):
sh = reach.jpca(reach.shuffle_rates(p["rates"], t0_index, kind, rng_shuffle))
out.append({"shuffle": kind, "run": r, **reach.summary(sh), "prep_r2": reach.prep_prediction_r2(sh["planes"][0])})
shuffles[name] = {"t0_ms": t0, "runs": pd.DataFrame(out)}
print(f"{name}: dividing point {t0} ms from movement onset")
J-array: dividing point -301 ms from movement onset
N-array: dividing point -320 ms from movement onset
array-dataset-composition¶
The two array datasets, monkey J-array recorded 2009-09-18 and monkey N-array recorded 2010-09-23, each present 108 reach conditions and comprise 146 and 218 isolations.
comp = pd.DataFrame({name: {"session_start": sessions[name]["meta"]["session_start_time"][:10],
"conditions": len(prep[name]["labels"]),
"units in the dandiset": len(sessions[name]["spikes"]),
"isolations in the paper": paper_n}
for name, paper_n in (("J-array", 146), ("N-array", 218))}).T
print(comp.to_string())
print("Tolerance on isolations (10%): J-array 131.4 to 160.6, N-array 196.2 to 239.8")
session_start conditions units in the dandiset isolations in the paper J-array 2009-09-18 108 192 146 N-array 2010-09-23 108 192 218 Tolerance on isolations (10%): J-array 131.4 to 160.6, N-array 196.2 to 239.8
Result: partially reproduced. Both sessions are in the dandiset with 108 conditions each. Each holds 192 units, one per electrode, against the paper's 146 (J-array) and 218 (N-array) isolations, outside the 10% tolerance in both. The sorted isolations the paper analysed are not in the dandiset, so every claim below is tested on the 192 units.
fig3ef-same-direction¶
In the first jPCA plane the neural state rotates in the same direction for all conditions of the J-array and N-array datasets, including reaches in opposite directions.
fig, axes = plt.subplots(1, 2, figsize=(10, 5))
for ax, name in zip(axes, ("J-array", "N-array")):
proj = res[name]["planes"][0]
col = plt.cm.RdYlGn_r((proj[:, 0, 0] - proj[:, 0, 0].min()) / np.ptp(proj[:, 0, 0]))
for c in range(proj.shape[0]):
ax.plot(proj[c, :, 0], proj[c, :, 1], color=col[c], lw=0.8)
ax.plot(proj[c, 0, 0], proj[c, 0, 1], "o", color=col[c], ms=3)
ax.set_title(f"{name}, -50 to +150 ms")
ax.set_xlabel("jPC1 (a.u.)"); ax.set_ylabel("jPC2 (a.u.)")
print(f"{name}: plotted jPC1 {proj[..., 0].min():.3f} to {proj[..., 0].max():.3f}, jPC2 {proj[..., 1].min():.3f} to {proj[..., 1].max():.3f} (all points inside the axes)")
plt.tight_layout(); plt.show()
for name in ("J-array", "N-array"):
rot = reach.net_rotation(res[name]["planes"][0])
print(f"{name}: {np.sum(rot > 0)} of {len(rot)} conditions rotate counter-clockwise ({np.mean(rot > 0):.3f}); "
f"net rotation median {np.degrees(np.median(rot)):.1f} deg, 5th to 95th percentile "
f"{np.degrees(np.percentile(rot, 5)):.1f} to {np.degrees(np.percentile(rot, 95)):.1f} deg")
print("\nFraction counter-clockwise under the variants:")
print(pd.DataFrame({n: variants[n].frac_ccw for n in variants}).round(3).to_string())
J-array: plotted jPC1 -0.870 to 1.264, jPC2 -0.863 to 1.083 (all points inside the axes) N-array: plotted jPC1 -1.246 to 1.460, jPC2 -0.784 to 1.342 (all points inside the axes)
J-array: 102 of 108 conditions rotate counter-clockwise (0.944); net rotation median 70.2 deg, 5th to 95th percentile 5.6 to 138.5 deg
N-array: 108 of 108 conditions rotate counter-clockwise (1.000); net rotation median 114.6 deg, 5th to 95th percentile 73.2 to 157.0 deg
Fraction counter-clockwise under the variants:
J-array N-array
10 PCs 0.963 0.972
24 ms Gaussian 0.981 1.000
all successful trials (2045) 0.944 NaN
all successful trials (2826) NaN 1.000
main 0.944 1.000
window -100 to +100 ms 0.972 0.981
window 0 to +200 ms 0.944 0.981
Result: reproduced. In the first jPCA plane 102 of 108 conditions (94.4%) rotate counter-clockwise in J-array and 108 of 108 in N-array, with median net rotations of 70.2 and 114.6 degrees over the 200 ms. Under the variants the fraction stays between 0.944 and 0.981 (J-array) and between 0.972 and 1.000 (N-array). jPCA orients the plane so that the net rotation is counter-clockwise, so this fraction is also high where no rotation was put in: 0.963 to 0.991 after shuffles 2 and 3, and 0.926 and 0.972 for the velocity model with no lead. It shows that the conditions turn the same way, not on its own that the rotation is specific to the data.
fig3-independent-of-curvature¶
The rotations do not relate directly to reach curvature; the neural state rotates the same way for clockwise-curving, counter-clockwise-curving and straight reaches.
kin = {name: reach.condition_kinematics(sessions[name], prep[name]["sel"], prep[name]["align"], prep[name]["cond"])
for name in ("J-array", "N-array")}
def logistic_slope(x, y):
if y.all() or not y.any():
return None
xs = (x - x.mean()) / x.std()
nll = lambda b: np.sum(np.logaddexp(0, b[0] + b[1] * xs) - y * (b[0] + b[1] * xs))
return minimize(nll, np.zeros(2)).x[1] / x.std()
for name in ("J-array", "N-array"):
turn = reach.heading_change(kin[name]["vel"])
rot = reach.net_rotation(res[name]["planes"][0])
cls = np.where(turn < -20, "clockwise", np.where(turn > 20, "counter-clockwise", "straight"))
tab = pd.DataFrame({"class": cls, "ccw": rot > 0}).groupby("class").ccw.agg(["size", "mean"])
print(f"{name}: heading change of the mean hand path {turn.min():.0f} to {turn.max():.0f} deg")
print(tab.rename(columns={"size": "conditions", "mean": "fraction rotating ccw"}).round(3).to_string())
rho, p = stats.spearmanr(turn, rot)
slope = logistic_slope(turn, (rot > 0).astype(float))
print(f" Spearman correlation of heading change with net neural rotation: rho {rho:.3f}, p {p:.3g}; "
f"logistic slope of rotation direction on heading change: {'not estimable (all conditions rotate one way)' if slope is None else f'{slope:.4f} per deg'}\n")
J-array: heading change of the mean hand path -213 to 225 deg
conditions fraction rotating ccw
class
clockwise 38 0.868
counter-clockwise 63 0.984
straight 7 1.000
Spearman correlation of heading change with net neural rotation: rho 0.227, p 0.018; logistic slope of rotation direction on heading change: 0.0101 per deg
N-array: heading change of the mean hand path -180 to 291 deg
conditions fraction rotating ccw
class
clockwise 44 1.0
counter-clockwise 53 1.0
straight 11 1.0
Spearman correlation of heading change with net neural rotation: rho -0.214, p 0.0264; logistic slope of rotation direction on heading change: not estimable (all conditions rotate one way)
# The hand paths by class (dotted: conditions whose neural state rotates clockwise), and the
# classes when the heading is taken over the unbroken run of samples around peak speed only, or
# with thresholds of 10 and 30 degrees.
fig, axes = plt.subplots(1, 2, figsize=(12, 6))
colours = {"clockwise": "tab:red", "counter-clockwise": "tab:blue", "straight": "0.3"}
rows = []
for ax, name in zip(axes, ("J-array", "N-array")):
pos, vel = kin[name]["pos"], kin[name]["vel"]
rot = reach.net_rotation(res[name]["planes"][0])
turn = reach.heading_change(vel)
cls = np.where(turn < -20, "clockwise", np.where(turn > 20, "counter-clockwise", "straight"))
for c in range(len(turn)):
xy = pos[c, -reach.WINDOW_MS[0]:] * 1000
ax.plot(xy[:, 0], xy[:, 1], color=colours[cls[c]], lw=1.2 if rot[c] < 0 else 0.6, ls="-" if rot[c] > 0 else ":")
ax.set_title(f"{name}: mean hand path, movement onset to +500 ms"); ax.set_xlabel("x (mm)"); ax.set_ylabel("y (mm)"); ax.set_aspect("equal")
print(f"{name}: plotted x {pos[:, -reach.WINDOW_MS[0]:, 0].min() * 1000:.0f} to {pos[:, -reach.WINDOW_MS[0]:, 0].max() * 1000:.0f}, "
f"y {pos[:, -reach.WINDOW_MS[0]:, 1].min() * 1000:.0f} to {pos[:, -reach.WINDOW_MS[0]:, 1].max() * 1000:.0f} mm")
print(" conditions rotating clockwise (heading change, net neural rotation, deg):",
[(int(kin[name]["labels"][c]), round(float(turn[c])), round(float(np.degrees(rot[c])))) for c in np.where(rot < 0)[0]])
for method, tc in (("all samples above 10% of peak speed", turn), ("unbroken run around peak speed", reach.heading_change(vel, contiguous=True))):
for thr in (10, 20, 30):
k = np.where(tc < -thr, "clockwise", np.where(tc > thr, "counter-clockwise", "straight"))
for c in ("clockwise", "counter-clockwise", "straight"):
rows.append({"session": name, "heading": method, "threshold deg": thr, "class": c,
"conditions": int((k == c).sum()), "fraction rotating ccw": float(np.mean(rot[k == c] > 0)) if (k == c).any() else np.nan})
plt.tight_layout(); plt.show()
print(pd.DataFrame(rows).pivot_table(index=["session", "heading", "threshold deg"], columns="class",
values=["conditions", "fraction rotating ccw"]).round(3).to_string())
J-array: plotted x -128 to 146, y -146 to 87 mm conditions rotating clockwise (heading change, net neural rotation, deg): [(2, 140, -139), (47, -162, -21), (48, -166, -1), (77, -180, -14), (78, -179, -19), (92, -164, -89)] N-array: plotted x -146 to 143, y -168 to 56 mm conditions rotating clockwise (heading change, net neural rotation, deg): []
conditions fraction rotating ccw
class clockwise counter-clockwise straight clockwise counter-clockwise straight
session heading threshold deg
J-array all samples above 10% of peak speed 10 41.0 65.0 2.0 0.878 0.985 1.0
20 38.0 63.0 7.0 0.868 0.984 1.0
30 38.0 60.0 10.0 0.868 0.983 1.0
unbroken run around peak speed 10 41.0 65.0 2.0 0.878 0.985 1.0
20 38.0 63.0 7.0 0.868 0.984 1.0
30 38.0 60.0 10.0 0.868 0.983 1.0
N-array all samples above 10% of peak speed 10 45.0 57.0 6.0 1.000 1.000 1.0
20 44.0 53.0 11.0 1.000 1.000 1.0
30 40.0 53.0 15.0 1.000 1.000 1.0
unbroken run around peak speed 10 44.0 57.0 7.0 1.000 1.000 1.0
20 43.0 53.0 12.0 1.000 1.000 1.0
30 39.0 53.0 16.0 1.000 1.000 1.0
Result: partially reproduced. In N-array every condition rotates counter-clockwise in all three classes of reach (44 clockwise, 53 counter-clockwise, 11 straight). In J-array 98.4% of the 63 counter-clockwise and all 7 straight reaches do, but 86.8% of the 38 clockwise reaches, below the 90% the criterion asks and above 75%. Rotation and curvature are weakly related, with opposite signs in the two datasets (Spearman rho 0.227 in J-array, -0.214 in N-array); in J-array the logistic slope of rotation direction on heading change is 0.0101 per degree. The hand paths drawn by class match the classes: the five clockwise-class conditions of J-array whose neural state turns clockwise are loops of 162 to 180 degrees, with net neural rotations of -1 to -21 degrees in four of them and -89 degrees in the fifth. With the heading taken over the unbroken run around peak speed, or with thresholds of 10 or 30 degrees, the J-array clockwise class has 38 or 41 conditions and 0.868 to 0.878 of them rotate counter-clockwise, so the outcome does not depend on how the classes are drawn.
fig3-phase-from-prep-state¶
The rotation phase follows from the preparatory state; knowing the initial state in the jPC plane, the subsequent states are reasonably predictable.
for name in ("J-array", "N-array"):
r2 = reach.prep_prediction_r2(res[name]["planes"][0])
sh3 = shuffles[name]["runs"].query("shuffle == 3").prep_r2
sh = {k: shuffles[name]["runs"].query("shuffle == @k").prep_r2.median() for k in (1, 2, 3)}
print(f"{name}: leave-one-condition-out R2 of the +150 ms state from the -50 ms state {r2:.3f}; "
f"shuffle 3 median {sh3.median():.3f} (range {sh3.min():.3f} to {sh3.max():.3f}); shuffle medians {({k: round(float(v), 3) for k, v in sh.items()})}")
for name in ("J-array", "N-array"):
print(f"{name} variants:", {k: round(reach.prep_prediction_r2(reach.jpca(prep[name]['rates'], times_ms=t)['planes'][0]), 3)
for k, t in (("-100 to +100", np.arange(-100, 101, 10)), ("0 to +200", np.arange(0, 201, 10)))})
J-array: leave-one-condition-out R2 of the +150 ms state from the -50 ms state 0.632; shuffle 3 median 0.671 (range 0.633 to 0.718); shuffle medians {1: 0.891, 2: 0.705, 3: 0.671}
N-array: leave-one-condition-out R2 of the +150 ms state from the -50 ms state 0.748; shuffle 3 median 0.656 (range 0.611 to 0.689); shuffle medians {1: 0.867, 2: 0.802, 3: 0.656}
J-array variants: {'-100 to +100': 0.692, '0 to +200': 0.636}
N-array variants: {'-100 to +100': 0.673, '0 to +200': 0.681}
# The R2 above per shuffle run (shuffle 1 against shuffle 3), and the same prediction with jPCA
# refitted leaving each condition out, together with a score of the phase alone: the last phase
# predicted as the first plus the other conditions' mean phase advance (`prep_prediction_cv`).
for name in ("J-array", "N-array"):
runs = shuffles[name]["runs"]
print(f"{name}, in-sample R2 per run (data {reach.prep_prediction_r2(res[name]['planes'][0]):.3f}):")
print(pd.DataFrame({f"shuffle {k}": runs.query("shuffle == @k").prep_r2.values for k in (1, 3)}).T.round(3).to_string(), "\n")
rng_cv = np.random.default_rng(7018)
rows = []
for name in ("J-array", "N-array"):
p = prep[name]
t0_index = shuffles[name]["t0_ms"] - reach.WINDOW_MS[0]
rows.append({"session": name, "rates": "data", "run": 0, **reach.prep_prediction_cv(p["rates"])})
for kind, reps in ((1, 10), (2, 1), (3, 10)):
for r in range(reps):
rows.append({"session": name, "rates": f"shuffle {kind}", "run": r,
**reach.prep_prediction_cv(reach.shuffle_rates(p["rates"], t0_index, kind, rng_cv))})
cv = pd.DataFrame(rows)
print(cv.groupby(["session", "rates"]).agg(runs=("r2", "size"), r2_median=("r2", "median"), r2_min=("r2", "min"), r2_max=("r2", "max"),
phase_median=("phase_score", "median"), phase_min=("phase_score", "min"), phase_max=("phase_score", "max"),
static_median=("static_phase_score", "median"), advance_deg_median=("advance_deg", "median")).round(3).to_string())
for name in ("J-array", "N-array"):
d, s3 = cv.query("session == @name and rates == 'data'").iloc[0], cv.query("session == @name and rates == 'shuffle 3'")
print(f"{name}: shuffle 3 runs at or above the data, cross-validated R2 {int((s3.r2 >= d.r2).sum())} of {len(s3)}, phase score {int((s3.phase_score >= d.phase_score).sum())} of {len(s3)}")
J-array, in-sample R2 per run (data 0.632):
0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19
shuffle 1 0.887 0.880 0.902 0.900 0.906 0.888 0.892 0.910 0.874 0.896 0.891 0.903 0.893 0.879 0.904 0.886 0.881 0.892 0.872 0.876
shuffle 3 0.677 0.674 0.676 0.718 0.642 0.669 0.673 0.675 0.667 0.646 0.701 0.649 0.665 0.695 0.670 0.672 0.692 0.633 0.646 0.657
N-array, in-sample R2 per run (data 0.748):
0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19
shuffle 1 0.861 0.886 0.868 0.866 0.881 0.864 0.847 0.889 0.849 0.863 0.886 0.886 0.871 0.871 0.872 0.850 0.873 0.860 0.835 0.817
shuffle 3 0.671 0.649 0.657 0.611 0.620 0.673 0.663 0.617 0.639 0.667 0.656 0.660 0.656 0.649 0.674 0.637 0.689 0.648 0.662 0.666
runs r2_median r2_min r2_max phase_median phase_min phase_max static_median advance_deg_median
session rates
J-array data 1 0.614 0.614 0.614 0.763 0.763 0.763 0.236 73.532
shuffle 1 10 0.723 0.628 0.793 0.779 0.673 0.857 0.782 3.294
shuffle 2 1 0.666 0.666 0.666 0.837 0.837 0.837 0.529 50.555
shuffle 3 10 0.630 0.596 0.701 0.811 0.697 0.847 0.334 64.858
N-array data 1 0.726 0.726 0.726 0.894 0.894 0.894 -0.356 113.091
shuffle 1 10 0.675 0.555 0.739 0.728 0.677 0.820 0.733 3.560
shuffle 2 1 0.787 0.787 0.787 0.887 0.887 0.887 0.110 80.754
shuffle 3 10 0.648 0.602 0.681 0.797 0.732 0.814 -0.153 99.295
J-array: shuffle 3 runs at or above the data, cross-validated R2 6 of 10, phase score 9 of 10
N-array: shuffle 3 runs at or above the data, cross-validated R2 0 of 10, phase score 0 of 10
Result: partially reproduced. One complex gain predicts the state at +150 ms from the state at -50 ms with cross-validated R2 0.632 (J-array) and 0.748 (N-array), both above 0.5. In N-array this is above every shuffle 3 run (median 0.656); in J-array it is below every shuffle 3 run (0.633 to 0.718, median 0.671). Shuffle 1, which has no rotation, scores higher than the data in every run (0.82 to 0.91), so this R2 rewards any consistent map from the first state to the last and is a weak test of the claim. With jPCA refitted leaving each condition out, the R2 is 0.615 (J-array) and 0.727 (N-array) against shuffle 3 medians of 0.628 and 0.649, and the pattern is the same: 8 of 10 shuffle 3 runs reach the J-array value and none the N-array value. Several shuffle 3 runs lie within 0.02 of the J-array value, so this count can move with small numerical differences between runs. The phase score, which asks whether the last phase follows from the first plus a shared advance, separates rotation from no rotation: the data score 0.763 (J-array) and 0.894 (N-array) against 0.236 and -0.356 for a prediction of no advance, while for shuffle 1 the two scores are equal (median 0.78 and 0.73). Shuffle 3 keeps a shared advance (median 65 and 99 degrees) and scores 0.811 and 0.797; 9 of 10 shuffle 3 runs reach the data's phase score in J-array and none in N-array, so the data exceed it in N-array and not in J-array, as for the criterion's measure.
var-captured-jpc-plane¶
The first jPCA plane captures an average of 28% of the total data variance (typically 20 to 35%).
pv = main.loc[["J-array", "N-array"], "plane_var"] * 100
print("Variance captured by the first jPCA plane (%):", pv.round(1).to_dict(), "; mean of the two", round(pv.mean(), 1))
print("Reading (a), mean within 25.2 to 30.8:", 25.2 <= pv.mean() <= 30.8, "; reading (b), each within 18 to 38.5:", bool(((pv >= 18) & (pv <= 38.5)).all()))
print("N-0910 (variant):", round(main.loc["N-0910", "plane_var"] * 100, 1))
print("\nUnder the variants (%):")
print((pd.DataFrame({n: variants[n].plane_var for n in variants}) * 100).round(1).to_string())
Variance captured by the first jPCA plane (%): {'J-array': 20.0, 'N-array': 23.1} ; mean of the two 21.6
Reading (a), mean within 25.2 to 30.8: False ; reading (b), each within 18 to 38.5: True
N-0910 (variant): 18.8
Under the variants (%):
J-array N-array
10 PCs 12.3 16.2
24 ms Gaussian 20.9 23.5
all successful trials (2045) 19.0 NaN
all successful trials (2826) NaN 23.1
main 20.0 23.1
window -100 to +100 ms 17.5 24.8
window 0 to +200 ms 20.7 20.8
Result: reproduced. The first jPCA plane captures 20.0% (J-array) and 23.1% (N-array) of the variance, both within the stated typical range widened to 18 to 38.5% (reading b). Their mean, 21.6%, is below the paper's average of 28 +- 2.8% (reading a). With 6 PCs the variants give 17.5 to 20.9% (J-array) and 20.8 to 24.8% (N-array). These are the dandiset's per-electrode units, not the paper's sorted isolations, so the fractions of variance compare different populations.
pc23-variance¶
The plane of PC2 and PC3 captures 29% of the total variance on average, about as much as the first jPCA plane.
t = main.loc[["J-array", "N-array"], ["pc1_var", "pc23_var", "plane_var", "pcs_var"]] * 100
t["jPC plane / PC2+PC3"] = t.plane_var / t.pc23_var
print(t.round(2).to_string())
print("Mean PC2+PC3 of the two datasets:", round(t.pc23_var.mean(), 1), "% (expected 29 +- 2.9)")
print("\nUnder the variants, PC2+PC3 (%):")
print((pd.DataFrame({n: variants[n].pc23_var for n in variants}) * 100).round(1).to_string())
pc1_var pc23_var plane_var pcs_var jPC plane / PC2+PC3
J-array 23.74 21.6 19.98 60.46 0.93
N-array 20.45 27.0 23.14 66.03 0.86
Mean PC2+PC3 of the two datasets: 24.3 % (expected 29 +- 2.9)
Under the variants, PC2+PC3 (%):
J-array N-array
10 PCs 21.6 27.0
24 ms Gaussian 22.2 27.6
all successful trials (2045) 21.7 NaN
all successful trials (2826) NaN 26.9
main 21.6 27.0
window -100 to +100 ms 21.7 30.3
window 0 to +200 ms 22.3 23.5
Result: partially reproduced. PC2 and PC3 together capture 21.6% (J-array) and 27.0% (N-array), mean 24.3%, below 29 +- 2.9% under both readings (N-array alone is within it). The jPCA plane captures 0.93 (J-array) and 0.86 (N-array) times the variance of PC2 and PC3, so the second part of the claim holds. These are the dandiset's per-electrode units, not the paper's sorted isolations, so the fractions of variance compare different populations.
suppfig4-j-array-planes and suppfig4-planes-range¶
In J-array the first three jPCA planes found within the top 10 PCs capture 22%, 19% and 12% of the variance (54% together), against 64% for the top six PCs.
Every dataset contains two, and usually three, planes with rotations, together capturing 50 to 70% of the variance.
planes10 = {}
for name in ("J-array", "N-array", "N-0910"):
r10 = reach.jpca(prep[name]["rates"], n_pcs=10, n_planes=3)
planes10[name] = {"plane 1 %": r10["plane_var"][0] * 100, "plane 2 %": r10["plane_var"][1] * 100,
"plane 3 %": r10["plane_var"][2] * 100, "three planes %": r10["plane_var"].sum() * 100,
"top 6 PCs %": res[name]["pc_var"].sum() * 100, "top 10 PCs %": r10["pc_var"].sum() * 100,
**{f"plane {k + 1} ccw": float(np.mean(reach.net_rotation(r10["planes"][k]) > 0)) for k in range(3)},
**{f"plane {k + 1} Hz": r10["eigenvalues_rad_s"][k] / (2 * np.pi) for k in range(3)}}
print(pd.DataFrame(planes10).T.round(3).to_string())
print("\nExpected for J-array: 22 +- 2.2, 19 +- 1.9, 12 +- 1.2, total 54 +- 5.4, top six PCs 64 +- 6.4; range for every dataset 45 to 77% (N-0910 is a variant)")
plane 1 % plane 2 % plane 3 % three planes % top 6 PCs % top 10 PCs % plane 1 ccw plane 2 ccw plane 3 ccw plane 1 Hz plane 2 Hz plane 3 Hz J-array 12.333 14.824 12.807 39.964 60.455 71.379 0.963 0.991 0.926 1.682 1.253 0.713 N-array 16.236 18.079 15.756 50.071 66.026 75.273 0.972 0.954 0.963 1.639 1.301 0.584 N-0910 15.059 16.461 9.356 40.876 59.091 70.002 1.000 0.991 0.926 1.797 1.520 0.915 Expected for J-array: 22 +- 2.2, 19 +- 1.9, 12 +- 1.2, total 54 +- 5.4, top six PCs 64 +- 6.4; range for every dataset 45 to 77% (N-0910 is a variant)
Result: not reproduced. In J-array the first three jPCA planes from the top 10 PCs capture 12.3%, 14.8% and 12.8% (40.0% together), against 22%, 19% and 12% (54%) in the paper; only the third plane is within tolerance. The top six PCs capture 60.5%, within 64 +- 6.4%. All three planes rotate consistently, with 92.6% to 99.1% of conditions counter-clockwise. These are the dandiset's per-electrode units, not the paper's sorted isolations, so the fractions of variance compare different populations.
Result: partially reproduced. The three planes capture 50.1% of the variance in N-array, inside the widened range of 45 to 77%, and 40.0% in J-array, below it. In both datasets the first two planes have 95.4% or more of the conditions rotating in the plane's direction. In ses-20090910 the three planes capture 40.9%. These are the dandiset's per-electrode units, not the paper's sorted isolations, so the fractions of variance compare different populations.
fig6a-angle-peak¶
The distribution of the angle between the neural state and its derivative in the first jPCA plane peaks near pi/2.
fig, ax = plt.subplots(figsize=(7, 3.5))
for name in ("J-array", "N-array"):
ang = reach.state_angles(res[name]["planes"][0])
h, e = np.histogram(ang, bins=36, range=(-np.pi, np.pi))
ax.plot((e[:-1] + e[1:]) / 2 / np.pi, h, label=name)
print(f"{name}: {ang.size} angles; peak {reach.angle_peak(ang) / np.pi:.3f} pi; median {np.median(ang) / np.pi:.3f} pi; "
f"fraction within pi/4 of pi/2 {np.mean(np.abs(ang - np.pi / 2) < np.pi / 4):.3f}; peak with 24 and 48 bins "
f"{reach.angle_peak(ang, 24) / np.pi:.3f} and {reach.angle_peak(ang, 48) / np.pi:.3f} pi; largest bin count {h.max()}")
ax.set_xlabel("angle between x and dx/dt (units of pi)"); ax.set_ylabel("count"); ax.legend(); plt.show()
print("\nPeak (units of pi) under the variants:")
print((pd.DataFrame({n: variants[n].angle_peak for n in variants}) / np.pi).round(3).to_string())
J-array: 2160 angles; peak 0.528 pi; median 0.511 pi; fraction within pi/4 of pi/2 0.624; peak with 24 and 48 bins 0.542 and 0.562 pi; largest bin count 198 N-array: 2160 angles; peak 0.639 pi; median 0.603 pi; fraction within pi/4 of pi/2 0.819; peak with 24 and 48 bins 0.625 and 0.646 pi; largest bin count 310
Peak (units of pi) under the variants:
J-array N-array
10 PCs 0.528 0.583
24 ms Gaussian 0.583 0.639
all successful trials (2045) 0.583 NaN
all successful trials (2826) NaN 0.639
main 0.528 0.639
window -100 to +100 ms 0.528 0.639
window 0 to +200 ms 0.528 0.583
Result: reproduced. The angle between the state and its derivative peaks at 0.528 pi in J-array and 0.639 pi in N-array (medians 0.511 and 0.603 pi), within 0.3 to 0.7 pi. Under the variants the peak stays between 0.528 and 0.583 pi (J-array) and between 0.583 and 0.639 pi (N-array), and with 24 or 48 bins it lies between 0.542 and 0.646 pi.
fig6b-skew-fit¶
For the neural data the rotational dynamics matrix fits the data nearly as well as the best unconstrained linear dynamics matrix, in 6-D and in the first jPCA plane.
print(main.loc[["J-array", "N-array", "N-0910"], ["r2_skew", "r2_best", "r2_ratio_6d", "r2_ratio_plane"]].round(3).to_string())
print("\n6-D ratio under the variants:")
print(pd.DataFrame({n: variants[n].r2_ratio_6d for n in variants}).round(3).to_string())
print("\nPlane ratio under the variants:")
print(pd.DataFrame({n: variants[n].r2_ratio_plane for n in variants}).round(3).to_string())
r2_skew r2_best r2_ratio_6d r2_ratio_plane
J-array 0.343 0.378 0.905 0.957
N-array 0.517 0.612 0.845 0.881
N-0910 0.526 0.602 0.874 0.932
6-D ratio under the variants:
J-array N-array
10 PCs 0.845 0.843
24 ms Gaussian 0.905 0.844
all successful trials (2045) 0.908 NaN
all successful trials (2826) NaN 0.842
main 0.905 0.845
window -100 to +100 ms 0.886 0.890
window 0 to +200 ms 0.852 0.819
Plane ratio under the variants:
J-array N-array
10 PCs 0.887 0.913
24 ms Gaussian 0.963 0.881
all successful trials (2045) 0.950 NaN
all successful trials (2826) NaN 0.879
main 0.957 0.881
window -100 to +100 ms 0.963 0.935
window 0 to +200 ms 0.953 0.887
Result: reproduced. R2(M_skew) / R2(M) is 0.905 (J-array) and 0.845 (N-array) in 6-D and 0.957 and 0.881 in the first plane, above the thresholds of 0.7 and 0.75. Under the variants the ratios stay at 0.819 or more in 6-D and 0.879 or more in the plane. The absolute R2 of M_skew is 0.343 (J-array) and 0.517 (N-array).
rotation-frequency¶
The rotations are brief and fast, up to about 2.5 Hz; the largest eigenvalues of M_skew in the monkey A and B datasets are 8.8 to 13.5 radians/s.
for name in ("J-array", "N-array", "N-0910"):
ev = res[name]["eigenvalues_rad_s"]
print(f"{name}: eigenvalue magnitudes of M_skew (rad/s) {np.round(ev, 2).tolist()}; first plane {ev[0] / (2 * np.pi):.2f} Hz")
print("\nFirst-plane frequency (Hz) under the variants:")
print(pd.DataFrame({n: variants[n].freq_hz for n in variants}).round(2).to_string())
J-array: eigenvalue magnitudes of M_skew (rad/s) [6.49, 4.65, 2.22]; first plane 1.03 Hz
N-array: eigenvalue magnitudes of M_skew (rad/s) [8.65, 4.67, 0.92]; first plane 1.38 Hz
N-0910: eigenvalue magnitudes of M_skew (rad/s) [9.58, 6.84, 1.92]; first plane 1.52 Hz
First-plane frequency (Hz) under the variants:
J-array N-array
10 PCs 1.68 1.64
24 ms Gaussian 1.01 1.36
all successful trials (2045) 1.06 NaN
all successful trials (2826) NaN 1.37
main 1.03 1.38
window -100 to +100 ms 1.45 1.43
window 0 to +200 ms 1.04 1.39
# Is J-array's low frequency a slower rotation, or a rotation that starts after the window does?
# The first plane of the main fit, projected from -300 to +300 ms; its mean angular speed over the
# conditions at each time; and the fit with 200 ms windows starting at other times.
times_long = np.arange(-300, 301, 10)
fig, axes = plt.subplots(3, 2, figsize=(11, 8), sharex=True)
for col, name in enumerate(("J-array", "N-array")):
r = res[name]
proj = reach.project(prep[name]["rates"], {**r, "transform": {**r["transform"], "times_ms": times_long}})
x, dx = proj[:, :-1], np.diff(proj, axis=1)
ang_speed = np.mean((x[..., 0] * dx[..., 1] - x[..., 1] * dx[..., 0]) / np.sum(x ** 2, axis=2), axis=0) / reach.DT_S / (2 * np.pi)
for k in range(2):
axes[k, col].plot(times_long, proj[:, :, k].T, lw=0.4, color="0.3", alpha=0.5)
axes[k, col].axvspan(-50, 150, color="tab:orange", alpha=0.15)
axes[k, col].set_ylabel(f"jPC{k + 1} (a.u.)")
axes[0, col].set_title(f"{name}, first plane of the -50 to +150 ms fit")
axes[2, col].plot(times_long[:-1] + 5, ang_speed); axes[2, col].axhline(0, color="0.6", lw=0.5)
axes[2, col].axvspan(-50, 150, color="tab:orange", alpha=0.15)
axes[2, col].set_ylabel("mean angular speed (Hz)"); axes[2, col].set_xlabel("ms from movement onset")
print(f"{name}: jPC1 {proj[..., 0].min():.3f} to {proj[..., 0].max():.3f}, jPC2 {proj[..., 1].min():.3f} to {proj[..., 1].max():.3f}; "
"mean angular speed (Hz) at", {int(t + 5): round(float(v), 2) for t, v in zip(times_long[:-1], ang_speed) if t % 50 == 0})
plt.tight_layout(); plt.show()
rows = []
for name in ("J-array", "N-array"):
for start in (-150, -100, -50, 0, 50):
sm = reach.summary(reach.jpca(prep[name]["rates"], times_ms=np.arange(start, start + 201, 10)))
rows.append({"session": name, "window": f"{start} to {start + 200} ms", "freq_hz": sm["freq_hz"],
"plane_var": sm["plane_var"], "frac_ccw": sm["frac_ccw"], "r2_ratio_6d": sm["r2_ratio_6d"]})
print(pd.DataFrame(rows).round(3).to_string(index=False))
# Helper check: a steady 2 Hz rotation in 192 units, present throughout or starting at movement onset
rng_syn = np.random.default_rng(7019)
t_all = np.arange(*reach.WINDOW_MS)
phase0, mix = rng_syn.uniform(0, 2 * np.pi, 108), rng_syn.normal(size=(2, 192))
for label, start in (("throughout", reach.WINDOW_MS[0]), ("from movement onset", 0)):
ang = phase0[:, None] + 2 * np.pi * 2.0 * np.clip(t_all - start, 0, None) / 1000
syn_rates = 20 + 10 * np.stack([np.cos(ang), np.sin(ang)], axis=-1) @ mix + rng_syn.normal(0, 0.5, (108, len(t_all), 192))
f = {w: round(reach.jpca(syn_rates, times_ms=np.arange(w, w + 201, 10))["freq_hz"], 2) for w in (-50, 0)}
print(f"Synthetic 2 Hz rotation {label}: frequency with the window from -50 ms {f[-50]} Hz, from 0 ms {f[0]} Hz")
del syn_rates
J-array: jPC1 -0.870 to 1.264, jPC2 -0.876 to 1.083; mean angular speed (Hz) at {-295: -0.09, -245: 0.63, -195: -0.64, -145: 0.49, -95: 0.99, -45: 1.32, 5: 1.26, 55: 0.94, 105: 0.74, 155: 0.92, 205: 0.58, 255: 0.02}
N-array: jPC1 -1.246 to 1.460, jPC2 -0.784 to 1.342; mean angular speed (Hz) at {-295: -0.69, -245: -0.68, -195: 2.04, -145: 0.92, -95: 1.52, -45: 1.38, 5: 1.28, 55: 1.94, 105: 1.74, 155: 1.38, 205: 0.87, 255: 0.51}
session window freq_hz plane_var frac_ccw r2_ratio_6d J-array -150 to 50 ms 1.305 0.203 0.954 0.910 J-array -100 to 100 ms 1.446 0.175 0.972 0.886 J-array -50 to 150 ms 1.033 0.200 0.944 0.905 J-array 0 to 200 ms 1.040 0.207 0.944 0.852 J-array 50 to 250 ms 1.126 0.191 0.972 0.849 N-array -150 to 50 ms 1.411 0.270 1.000 0.879 N-array -100 to 100 ms 1.428 0.248 0.981 0.890 N-array -50 to 150 ms 1.377 0.231 1.000 0.845 N-array 0 to 200 ms 1.393 0.208 0.981 0.819 N-array 50 to 250 ms 1.427 0.161 0.963 0.850
Synthetic 2 Hz rotation throughout: frequency with the window from -50 ms 1.99 Hz, from 0 ms 1.99 Hz
Synthetic 2 Hz rotation from movement onset: frequency with the window from -50 ms 1.49 Hz, from 0 ms 1.99 Hz
Result: partially reproduced. The first plane rotates at 1.38 Hz in N-array (8.65 radians/s), inside 1.26 to 2.75 Hz, and at 1.03 Hz in J-array (6.49 radians/s), below it. No variant exceeds 1.68 Hz, so the upper bound of about 2.5 Hz holds; Among the variants above, J-array reaches the range only with the -100 to +100 ms window (1.45 Hz) or 10 PCs (1.68 Hz); the -150 to +50 ms window below also does (1.31 Hz). J-array's lower frequency is not the effect of a window that starts before a rotation does (on a synthetic 2 Hz rotation that starts at movement onset, the -50 to +150 ms window reads 1.49 Hz and a window from 0 ms reads 1.99 Hz): in J-array windows starting at 0 or +50 ms give 1.04 and 1.13 Hz. Its mean angular speed in the plane is about 1.3 Hz from -60 to +20 ms and falls to 0.7 to 0.9 Hz by +100 ms, and windows starting at -150 or -100 ms give 1.31 and 1.45 Hz, so J-array rotates faster before movement onset than over the window used here. The paper chose a window per dataset; a J-array window starting before -50 ms would reach the range.
suppfig3-shuffle-controls¶
Three shuffle controls that keep preparatory activity and invert or reassign the peri-movement activity largely remove the rotational structure found by jPCA.
for name in ("J-array", "N-array"):
runs = shuffles[name]["runs"]
orig = main.loc[name]
g = runs.groupby("shuffle").agg(runs=("r2_skew", "size"), r2_skew_median=("r2_skew", "median"), r2_skew_max=("r2_skew", "max"),
ratio_6d_median=("r2_ratio_6d", "median"), frac_ccw_median=("frac_ccw", "median"),
plane_var_median=("plane_var", "median"))
g["all runs below original"] = runs.groupby("shuffle").r2_skew.max() < orig.r2_skew
g["median / original"] = g.r2_skew_median / orig.r2_skew
print(f"{name}: original R2(M_skew) {orig.r2_skew:.3f}, 6-D ratio {orig.r2_ratio_6d:.3f}, fraction ccw {orig.frac_ccw:.3f}; dividing point {shuffles[name]['t0_ms']} ms")
print(g.round(3).to_string(), "\n")
fig, axes = plt.subplots(1, 4, figsize=(14, 3.5))
p = prep["J-array"]
t0_index = shuffles["J-array"]["t0_ms"] - reach.WINDOW_MS[0]
demo_rng = np.random.default_rng(7013)
for ax, kind in zip(axes, (0, 1, 2, 3)):
rates = p["rates"] if kind == 0 else reach.shuffle_rates(p["rates"], t0_index, kind, demo_rng)
proj = reach.jpca(rates)["planes"][0]
for c in range(proj.shape[0]):
ax.plot(proj[c, :, 0], proj[c, :, 1], lw=0.6)
ax.set_title("J-array original" if kind == 0 else f"J-array shuffle {kind}")
print(f"panel {kind}: jPC1 {proj[..., 0].min():.3f} to {proj[..., 0].max():.3f}, jPC2 {proj[..., 1].min():.3f} to {proj[..., 1].max():.3f}")
plt.tight_layout(); plt.show()
J-array: original R2(M_skew) 0.343, 6-D ratio 0.905, fraction ccw 0.944; dividing point -301 ms
runs r2_skew_median r2_skew_max ratio_6d_median frac_ccw_median plane_var_median all runs below original median / original
shuffle
1 20 0.003 0.007 0.089 0.565 0.061 True 0.010
2 1 0.245 0.245 0.845 0.991 0.151 True 0.714
3 20 0.279 0.302 0.911 0.963 0.156 True 0.814
N-array: original R2(M_skew) 0.517, 6-D ratio 0.845, fraction ccw 1.000; dividing point -320 ms
runs r2_skew_median r2_skew_max ratio_6d_median frac_ccw_median plane_var_median all runs below original median / original
shuffle
1 20 0.004 0.009 0.090 0.602 0.062 True 0.008
2 1 0.485 0.485 0.907 0.991 0.202 True 0.938
3 20 0.484 0.524 0.885 0.963 0.206 False 0.937
panel 0: jPC1 -0.870 to 1.264, jPC2 -0.863 to 1.083
panel 1: jPC1 -0.825 to 1.026, jPC2 -1.127 to 1.165
panel 2: jPC1 -0.852 to 0.887, jPC2 -0.774 to 1.029
panel 3: jPC1 -1.213 to 0.852, jPC2 -1.191 to 0.823
Result: not reproduced. Shuffle 1 removes the rotational structure: R2(M_skew) falls to a median of 1.0% (J-array) and 0.8% (N-array) of the original, and the fraction of conditions rotating counter-clockwise to 0.565 and 0.602. Shuffles 2 and 3 do not meet the criterion. In J-array they leave 71% and 81% of the original R2(M_skew) (shuffle 2 is under the 75% threshold, shuffle 3 is not); in N-array they leave 94% each, and one run of shuffle 3 exceeds the original. After shuffles 2 and 3 the median fraction rotating counter-clockwise is 0.963 to 0.991 in both datasets. The dividing point here is one time per dataset relative to movement onset (-301 and -320 ms), where the paper used 50 ms after the go cue of each condition's own alignment. As the supplement says, shuffles 2 and 3 are not expected to remove all rotational structure; they are meant to break the relation of the rotation to the preparatory state, which R2(M_skew) does not measure. The criterion fixed R2(M_skew) before the analysis, and the outcome follows it.
suppfig5-areas-variance and suppfig5-m1-faster¶
Analysed separately, the PMd and M1 arrays each show clear rotations; the jPCA plane captures 29% (PMd) and 29% (M1) in J-array and 32% and 28% in N-array.
The angular velocity of the rotations is slightly higher for M1 than for PMd in each dataset.
rows = []
for name in ("J-array", "N-array", "N-0910"):
for area in ("PMd", "M1"):
m = sessions[name]["unit_area"] == area
for label, rates in (("main", prep[name]["rates"]), ("24 ms Gaussian", reach.condition_rates(prep[name]["sums"], prep[name]["n"], 24.0))):
for wlabel, times in (("-50 to +150", reach.ANALYSIS_TIMES), ("-100 to +100", np.arange(-100, 101, 10))):
r = reach.jpca(rates[:, :, m], times_ms=times)
sm = reach.summary(r)
rows.append({"session": name, "area": area, "units": int(m.sum()), "smoothing": label, "window": wlabel,
"plane var %": sm["plane_var"] * 100, "frac ccw": sm["frac_ccw"], "freq Hz": sm["freq_hz"],
"r2 ratio 6d": sm["r2_ratio_6d"]})
areas = pd.DataFrame(rows)
print(areas.round(3).to_string(index=False))
mainrows = areas[(areas.smoothing == "main") & (areas.window == "-50 to +150")].set_index(["session", "area"])
print("\nExpected (%): J-array PMd 29, M1 29; N-array PMd 32, M1 28; N-9-10-2010 PMd 26, M1 24 (variant)")
for name in ("J-array", "N-array", "N-0910"):
f = mainrows.loc[name, "freq Hz"]
print(f"{name}: M1 {f['M1']:.2f} Hz, PMd {f['PMd']:.2f} Hz, M1 higher: {f['M1'] > f['PMd']}")
v = areas.pivot_table(index=["session", "smoothing", "window"], columns="area", values="freq Hz")
v["M1 higher"] = v.M1 > v.PMd
print("\nM1 against PMd frequency under the variants:")
print(v.round(2).to_string())
session area units smoothing window plane var % frac ccw freq Hz r2 ratio 6d
J-array PMd 96 main -50 to +150 23.726 0.935 0.830 0.848
J-array PMd 96 main -100 to +100 25.125 0.954 0.777 0.836
J-array PMd 96 24 ms Gaussian -50 to +150 24.281 0.944 0.829 0.846
J-array PMd 96 24 ms Gaussian -100 to +100 25.538 0.954 0.781 0.835
J-array M1 96 main -50 to +150 22.788 0.991 1.822 0.785
J-array M1 96 main -100 to +100 23.615 0.991 1.836 0.807
J-array M1 96 24 ms Gaussian -50 to +150 23.302 0.991 1.774 0.786
J-array M1 96 24 ms Gaussian -100 to +100 24.193 0.981 1.801 0.810
N-array PMd 96 main -50 to +150 18.066 0.917 1.018 0.766
N-array PMd 96 main -100 to +100 18.563 0.954 1.050 0.845
N-array PMd 96 24 ms Gaussian -50 to +150 18.321 0.917 1.010 0.767
N-array PMd 96 24 ms Gaussian -100 to +100 18.836 0.954 1.033 0.841
N-array M1 96 main -50 to +150 24.966 0.972 1.541 0.800
N-array M1 96 main -100 to +100 27.635 0.972 1.538 0.805
N-array M1 96 24 ms Gaussian -50 to +150 25.275 0.972 1.525 0.800
N-array M1 96 24 ms Gaussian -100 to +100 28.020 0.972 1.508 0.803
N-0910 PMd 96 main -50 to +150 15.388 0.833 0.892 0.625
N-0910 PMd 96 main -100 to +100 16.181 0.926 1.145 0.818
N-0910 PMd 96 24 ms Gaussian -50 to +150 15.794 0.833 0.892 0.629
N-0910 PMd 96 24 ms Gaussian -100 to +100 16.564 0.917 1.131 0.821
N-0910 M1 96 main -50 to +150 21.572 0.981 1.701 0.807
N-0910 M1 96 main -100 to +100 22.520 0.935 1.744 0.846
N-0910 M1 96 24 ms Gaussian -50 to +150 21.876 0.981 1.680 0.808
N-0910 M1 96 24 ms Gaussian -100 to +100 22.916 0.944 1.713 0.842
Expected (%): J-array PMd 29, M1 29; N-array PMd 32, M1 28; N-9-10-2010 PMd 26, M1 24 (variant)
J-array: M1 1.82 Hz, PMd 0.83 Hz, M1 higher: True
N-array: M1 1.54 Hz, PMd 1.02 Hz, M1 higher: True
N-0910: M1 1.70 Hz, PMd 0.89 Hz, M1 higher: True
M1 against PMd frequency under the variants:
area M1 PMd M1 higher
session smoothing window
J-array 24 ms Gaussian -100 to +100 1.80 0.78 True
-50 to +150 1.77 0.83 True
main -100 to +100 1.84 0.78 True
-50 to +150 1.82 0.83 True
N-0910 24 ms Gaussian -100 to +100 1.71 1.13 True
-50 to +150 1.68 0.89 True
main -100 to +100 1.74 1.14 True
-50 to +150 1.70 0.89 True
N-array 24 ms Gaussian -100 to +100 1.51 1.03 True
-50 to +150 1.53 1.01 True
main -100 to +100 1.54 1.05 True
-50 to +150 1.54 1.02 True
Result: partially reproduced. The jPCA plane captures 23.7% (PMd) and 22.8% (M1) in J-array and 18.1% and 25.0% in N-array, against 29, 29, 32 and 28%; none is within the 10% tolerance (N-array M1 misses its lower bound of 25.2% by 0.2 points). All four lie within the widened typical range of 18 to 38.5%, and 91.7% to 99.1% of conditions rotate counter-clockwise in each array. In ses-20090910 the values are 15.4% (PMd) and 21.6% (M1), against 26% and 24% for the paper's N-9-10-2010. These are the dandiset's per-electrode units, not the paper's sorted isolations, so the fractions of variance compare different populations.
Result: reproduced. The first plane rotates faster for M1 than for PMd in J-array (1.82 against 0.83 Hz) and N-array (1.54 against 1.02 Hz), and in ses-20090910 (1.70 against 0.89 Hz), under every variant shown. The differences are larger than the paper's "slightly": in J-array the M1 frequency is more than twice the PMd one.
fig6-velocity-model and fig6-complex-kinematic-model¶
A simulated population tuned to hand velocity and speed, with preparatory activity tuned to the endpoint, shows no robust or consistent rotations; M_skew fits it poorly and its angle distribution peaks slightly above zero.
A simulated population tuned to hand position, velocity, acceleration and jerk produces multiphasic responses but no consistent rotations.
One simulated unit per recorded unit, from each session's condition-mean hand kinematics, with
spikes drawn per trial at the recorded trial counts (reach.simulate_population; the choices the
Methods leave open are listed in its docstring). A lead of 0 ms is a variant.
rng_models = np.random.default_rng(7016)
rows, sims = [], {}
for name in ("J-array", "N-array"):
p = prep[name]
for kind in ("velocity", "complex"):
for lead in (100, 0):
sim = reach.simulate_population(kind, kin[name], p["rates"], p["n"], rng_models, lead_ms=lead)
r = reach.jpca(sim)
sm = reach.summary(r)
sims[(name, kind, lead)] = r
rows.append({"session": name, "model": kind, "lead ms": lead, **sm,
"sim rate range": f"{sim.min():.1f} to {sim.max():.1f}"})
rows.append({"session": name, "model": "neural data", "lead ms": None, **main.loc[name].to_dict(), "sim rate range": ""})
models = pd.DataFrame(rows)
models["angle_peak / pi"] = models.angle_peak / np.pi
print(models[["session", "model", "lead ms", "r2_ratio_6d", "r2_ratio_plane", "angle_peak / pi", "frac_ccw", "plane_var", "freq_hz", "sim rate range"]].round(3).to_string(index=False))
fig, axes = plt.subplots(1, 2, figsize=(10, 4.5))
for ax, kind in zip(axes, ("velocity", "complex")):
proj = sims[("J-array", kind, 100)]["planes"][0]
for c in range(proj.shape[0]):
ax.plot(proj[c, :, 0], proj[c, :, 1], lw=0.6)
ax.set_title(f"J-array, {kind} model")
print(f"{kind}: jPC1 {proj[..., 0].min():.3f} to {proj[..., 0].max():.3f}, jPC2 {proj[..., 1].min():.3f} to {proj[..., 1].max():.3f}")
plt.tight_layout(); plt.show()
session model lead ms r2_ratio_6d r2_ratio_plane angle_peak / pi frac_ccw plane_var freq_hz sim rate range J-array velocity 100.0 0.361 0.765 0.806 0.630 0.193 0.836 0.0 to 318.2 J-array velocity 0.0 0.362 0.564 0.194 0.926 0.490 0.730 0.0 to 390.9 J-array complex 100.0 0.082 0.159 0.861 0.565 0.190 0.586 0.0 to 413.1 J-array complex 0.0 0.289 0.785 0.806 0.602 0.538 0.662 0.5 to 375.2 J-array neural data NaN 0.905 0.957 0.528 0.944 0.200 1.033 N-array velocity 100.0 0.255 0.746 0.361 0.778 0.217 0.815 0.0 to 322.2 N-array velocity 0.0 0.337 0.545 0.417 0.972 0.205 1.062 0.0 to 282.2 N-array complex 100.0 0.132 0.207 0.750 0.528 0.275 0.466 0.0 to 246.5 N-array complex 0.0 0.207 0.583 0.917 0.750 0.409 0.657 2.6 to 248.8 N-array neural data NaN 0.845 0.881 0.639 1.000 0.231 1.377 velocity: jPC1 -2.204 to 2.431, jPC2 -2.183 to 1.414 complex: jPC1 -0.496 to 0.643, jPC2 -0.347 to 0.328
# What drives the models' angle peak: the preparatory envelope, the lead and the seed. The main
# models above decay their preparatory tuning with a 100 ms time constant from movement onset;
# here it is held on to the end of the window (step) or left out, and the main model is drawn
# with four more seeds.
syn = np.exp(-np.arange(0, 200, 10) / 100.0)[None, :, None] * np.random.default_rng(7020).normal(size=(108, 1, 2))
print(f"Helper check: angle peak of a pure decay toward the origin {reach.angle_peak(reach.state_angles(syn)) / np.pi:.3f} pi\n")
rng_var = np.random.default_rng(7017)
rows, angles = [], {}
for name in ("J-array", "N-array"):
p = prep[name]
for kind in ("velocity", "complex"):
runs = [("decay 100 ms", {}, lead, s) for lead in (100, 0) for s in range(1, 5)]
runs += [(env, kw, lead, 0) for env, kw in (("step", {"prep_decay_ms": None}),
("no preparatory tuning", {"prep_onset_ms": reach.WINDOW_MS[1]}))
for lead in (100, 0)]
for env, kw, lead, seed in runs:
r = reach.jpca(reach.simulate_population(kind, kin[name], p["rates"], p["n"], rng_var, lead_ms=lead, **kw))
sm = reach.summary(r)
angles[(name, kind, env, lead, seed)] = reach.state_angles(r["planes"][0])
rows.append({"session": name, "model": kind, "preparatory envelope": env, "lead ms": lead, "seed": seed,
"angle_peak / pi": sm["angle_peak"] / np.pi, "angle_median / pi": sm["angle_median"] / np.pi,
"r2_ratio_6d": sm["r2_ratio_6d"], "frac_ccw": sm["frac_ccw"]})
mvar = pd.DataFrame(rows)
print(mvar.round(3).to_string(index=False))
print("\nFraction of runs whose angle peak lies within pi/4 of zero, by model, envelope and lead:")
print(mvar.assign(below=mvar["angle_peak / pi"].abs() < 0.25).groupby(["model", "preparatory envelope", "lead ms", "session"]).below.mean().unstack().to_string())
fig, axes = plt.subplots(2, 2, figsize=(12, 6), sharex=True)
centres = np.linspace(-1, 1, 37)[:-1] + 1 / 36
for (i, name), (j, kind) in ((a, b) for a in enumerate(("J-array", "N-array")) for b in enumerate(("velocity", "complex"))):
ax = axes[i, j]
for env, colour in (("decay 100 ms", "tab:blue"), ("step", "tab:orange"), ("no preparatory tuning", "tab:green")):
for lead, ls in ((100, "-"), (0, ":")):
h, _ = np.histogram(angles[(name, kind, env, lead, 1 if env == "decay 100 ms" else 0)], bins=36, range=(-np.pi, np.pi))
ax.plot(centres, h, color=colour, ls=ls, label=f"{env}, lead {lead} ms")
ax.set_title(f"{name}, {kind} model"); ax.set_ylabel("count")
axes[1, 0].set_xlabel("angle between x and dx/dt (units of pi)"); axes[1, 1].set_xlabel("angle between x and dx/dt (units of pi)")
axes[0, 0].legend(fontsize=7); plt.tight_layout(); plt.show()
Helper check: angle peak of a pure decay toward the origin 0.972 pi
session model preparatory envelope lead ms seed angle_peak / pi angle_median / pi r2_ratio_6d frac_ccw
J-array velocity decay 100 ms 100 1 0.694 0.312 0.376 0.630
J-array velocity decay 100 ms 100 2 0.861 0.291 0.342 0.639
J-array velocity decay 100 ms 100 3 0.861 0.312 0.330 0.657
J-array velocity decay 100 ms 100 4 0.806 0.324 0.389 0.639
J-array velocity decay 100 ms 0 1 0.250 0.279 0.385 0.954
J-array velocity decay 100 ms 0 2 0.250 0.287 0.395 0.954
J-array velocity decay 100 ms 0 3 0.194 0.298 0.375 0.935
J-array velocity decay 100 ms 0 4 0.361 0.285 0.375 0.926
J-array velocity step 100 0 0.583 0.313 0.551 0.546
J-array velocity step 0 0 0.194 0.217 0.327 0.907
J-array velocity no preparatory tuning 100 0 0.472 0.297 0.482 0.556
J-array velocity no preparatory tuning 0 0 0.028 0.032 0.088 0.574
J-array complex decay 100 ms 100 1 0.806 0.316 0.087 0.593
J-array complex decay 100 ms 100 2 0.806 0.235 0.094 0.630
J-array complex decay 100 ms 100 3 0.806 0.234 0.080 0.657
J-array complex decay 100 ms 100 4 0.861 0.268 0.092 0.611
J-array complex decay 100 ms 0 1 0.861 0.126 0.260 0.611
J-array complex decay 100 ms 0 2 0.806 0.124 0.270 0.620
J-array complex decay 100 ms 0 3 -0.917 0.127 0.255 0.611
J-array complex decay 100 ms 0 4 -0.917 0.128 0.255 0.620
J-array complex step 100 0 0.750 0.285 0.124 0.741
J-array complex step 0 0 -0.972 0.091 0.469 0.546
J-array complex no preparatory tuning 100 0 0.694 0.153 0.061 0.611
J-array complex no preparatory tuning 0 0 -0.083 0.080 0.608 0.556
N-array velocity decay 100 ms 100 1 0.361 0.370 0.269 0.861
N-array velocity decay 100 ms 100 2 0.361 0.377 0.274 0.824
N-array velocity decay 100 ms 100 3 0.361 0.370 0.282 0.815
N-array velocity decay 100 ms 100 4 0.528 0.371 0.240 0.833
N-array velocity decay 100 ms 0 1 0.417 0.366 0.325 0.963
N-array velocity decay 100 ms 0 2 0.417 0.343 0.315 0.963
N-array velocity decay 100 ms 0 3 0.417 0.356 0.309 0.963
N-array velocity decay 100 ms 0 4 0.417 0.333 0.316 0.954
N-array velocity step 100 0 0.472 0.266 0.189 0.583
N-array velocity step 0 0 0.361 0.350 0.381 0.963
N-array velocity no preparatory tuning 100 0 0.583 0.223 0.298 0.667
N-array velocity no preparatory tuning 0 0 0.028 0.078 0.076 0.593
N-array complex decay 100 ms 100 1 0.861 0.076 0.148 0.574
N-array complex decay 100 ms 100 2 0.806 0.090 0.166 0.546
N-array complex decay 100 ms 100 3 0.806 -0.006 0.122 0.519
N-array complex decay 100 ms 100 4 0.694 0.154 0.137 0.574
N-array complex decay 100 ms 0 1 0.194 0.240 0.233 0.713
N-array complex decay 100 ms 0 2 0.250 0.261 0.223 0.694
N-array complex decay 100 ms 0 3 0.417 0.249 0.226 0.676
N-array complex decay 100 ms 0 4 0.306 0.239 0.219 0.713
N-array complex step 100 0 0.750 0.170 0.188 0.787
N-array complex step 0 0 0.972 0.168 0.202 0.546
N-array complex no preparatory tuning 100 0 0.750 0.194 0.146 0.556
N-array complex no preparatory tuning 0 0 0.694 0.130 0.405 0.583
Fraction of runs whose angle peak lies within pi/4 of zero, by model, envelope and lead:
session J-array N-array
model preparatory envelope lead ms
complex decay 100 ms 0 0.00 0.25
100 0.00 0.00
no preparatory tuning 0 1.00 0.00
100 0.00 0.00
step 0 0.00 0.00
100 0.00 0.00
velocity decay 100 ms 0 0.25 0.00
100 0.00 0.00
no preparatory tuning 0 1.00 1.00
100 0.00 0.00
step 0 1.00 0.00
100 0.00 0.00
Result: not reproduced. The simulated velocity populations fit M_skew poorly, with 6-D ratios of 0.361 (J-array) and 0.255 (N-array) against 0.905 and 0.845 for the data, and fewer of their conditions rotate counter-clockwise (0.630 and 0.778). Their angle distributions peak at 0.806 pi and 0.361 pi, not below pi/4 as the criterion asks, so the criterion fails in both datasets. The angle peak depends on choices the Methods leave open. A pure decay toward the origin peaks near pi (0.972 pi), so the decaying preparatory tuning pushes the peak up: with the lead of 100 ms the peak stays between 0.36 and 0.86 pi over five seeds and under every envelope, and with no lead it lies between 0.19 and 0.42 pi with preparatory tuning (decaying or held). Only with no preparatory tuning and no lead does it peak near zero in both datasets (0.028 pi), as the paper reports for its velocity model. The poor fit of M_skew holds under every variant (6-D ratio 0.08 to 0.55).
Result: not reproduced. The complex kinematic populations fit M_skew poorly, with 6-D ratios of 0.082 (J-array) and 0.132 (N-array), 0.289 and 0.207 with no lead, and 0.565 and 0.528 of their conditions rotate counter-clockwise. Their angle peaks are at 0.861 pi and 0.750 pi (0.806 and 0.917 pi with no lead), not below pi/4, so the criterion fails in both datasets. Under the other envelopes, leads and seeds the angle distribution is broad with its peak between 0.19 and 0.97 pi or near -pi; within pi/4 of zero it peaks in J-array only with no preparatory tuning and no lead (-0.083 pi), and in N-array only in one of five seeds with no lead (0.194 pi). The fit of M_skew separates the model from the data as the paper reports (6-D ratio 0.06 to 0.61 across the variants); under none of the choices tried here does the angle distribution peak near zero in both datasets, as in Figure 6a.
Blocked claims¶
These claims rest on data the dandiset does not hold, and that the paper did not deposit elsewhere.
- fig1-leech-and-walking: the leech imaging (made by K. Briggman and W. Kristan) and the walking monkey's wireless recording are not in the dandiset.
- fig2-single-neuron-responses and fig3abcd-single-electrode-datasets: the 469 single-electrode isolations of monkeys A, B, J (J1 to J4) and N are not in the dandiset, which holds only array sessions.
- fig4-emg-no-rotation and fig5-generator-model: EMG was not recorded for the array datasets (Methods), and the EMG of the single-electrode datasets is not in the dandiset.
- fig6cd-speed-independent-frequency: the speed task of monkeys A and B is not in the dandiset.
print(f"Total run time {(time.time() - started) / 60:.1f} min")
Total run time 16.3 min
Summary¶
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.
| Claim | Outcome |
|---|---|
| array-dataset-composition | partially reproduced |
| fig3ef-same-direction | reproduced |
| fig3-independent-of-curvature | partially reproduced |
| fig3-phase-from-prep-state | partially reproduced |
| var-captured-jpc-plane | reproduced |
| pc23-variance | partially reproduced |
| suppfig4-j-array-planes | not reproduced |
| suppfig4-planes-range | partially reproduced |
| fig6a-angle-peak | reproduced |
| fig6b-skew-fit | reproduced |
| rotation-frequency | partially reproduced |
| suppfig3-shuffle-controls | not reproduced |
| suppfig5-areas-variance | partially reproduced |
| suppfig5-m1-faster | reproduced |
| fig6-velocity-model | not reproduced |
| fig6-complex-kinematic-model | not reproduced |
| fig1-leech-and-walking, fig2-single-neuron-responses, fig3abcd-single-electrode-datasets, fig4-emg-no-rotation, fig5-generator-model, fig6cd-speed-independent-frequency | blocked: data missing |
How the analysis differs from the paper's. The units are the dandiset's 192 per session, one per electrode, where the paper used 146 and 218 sorted isolations. The rates are not PC-smoothed, and the Gaussian is 20 ms (24 ms is a variant). The jPCA window is -50 to +150 ms from movement onset for every dataset, the codepack example's, where the paper chose a start per dataset. The parameters the Methods leave open come from the Python port of the jPCA codepack, as the lab's own download page could not be fetched. The shuffle controls use one dividing time per dataset relative to movement onset. The kinematic models are simulated here from the dandiset's hand kinematics, with a lead, a preparatory time course and modulation depths chosen here. The paper's averages over nine datasets are compared with the two array datasets the dandiset holds.
What was found about the dandiset.
- It holds the paper's two array datasets and eight further array sessions the paper does not describe. The seven single-electrode datasets, the EMG, the leech and the walking-monkey data are not in it, nor are the sorted isolations of the array datasets: the Units table has one unit per electrode and no description of what a unit is.
- In nine of the ten files the first entry of
units/spike_times_indexwrapped at 65536, so unit 0 is cut short and its remaining spikes open unit 1. - In five Nitschke files, N-array among them, the clock restarts at recording segments after the first, so trials overlap in time and spike times are not sorted; spikes can be placed in trials only by their storage order.
- Hand position is stated in metres but is in mm, and is offset in y from the screen coordinates by 35 to 53 mm, which the files do not describe.
- The conversion code dates every session of its first conversion in 2009. The file
ses-20090910has the 108 conditions and reached targets of N-array (2010), and is consistent with the N-9-10-2010 dataset of Supplementary Figure 5, but its date cannot confirm this.