Reproducing Even-Chen et al. (2019) from dandiset 000121¶
Paper. Even-Chen N, Sheffer B, Vyas S, Ryu SI, Shenoy KV (2019). Structure and variability of delay activity in premotor cortex. PLOS Computational Biology 15(2), e1006808. doi:10.1371/journal.pcbi.1006808. This attempt used the full text, which is open access.
Dandiset. DANDI:000121, version 0.220124.2156. The paper cites no DANDI version (it predates the dandiset; its data statement points to a Stanford zip file that now returns HTTP 403). This is the most recent published version. The earlier published version, 0.210815.0703, holds the same 12 files (same blobs and sizes) and differs only in metadata; the draft holds no files (see the data check).
What the dandiset holds. 12 sessions of delayed center-out reaching with two Utah arrays: all 7 sessions of monkey R and 5 of the 16 of monkey J (ring 1 of 4, vertical 2 of 4, 3-rings 2 of 4, horizontal none). Each file has a trials table, hand, cursor and eye positions, threshold crossings of the 192 electrodes as units, and LFP.
Code read. The paper names no analysis code. The conversion code,
catalystneuro/shenoy-lab-to-nwb at commit
319acb472372c9e4359ebacf13c8df3652a822b1 (center_out_task/), was read. It builds the units from
the rig's 1 ms threshold-crossing rasters of the two arrays (units 0 to 95 from spikeRaster, 96
to 191 from spikeRaster2), takes go_cue_time from the rig's timeTargetOn and
target_shown_time as that minus delayTime, assigns the units' electrode_group by the order of
a dictionary of groups, and stores positions with conversion NaN. It settles no analysis
choice.
Outside inputs. None.
No claim uses data from outside the dandiset. Random numbers come from generators seeded per section (the seeds are named at the top of each section).
| Claim | Status |
|---|---|
| session-configurations | tested |
| r-trial-counts | tested |
| j-sessions-complete | blocked: 11 of J's 16 sessions are not in the dandiset |
| delays-and-times | tested |
| channels-task-correlated | tested |
| pca-variance | tested (J on 5 of 16 sessions) |
| spatial-plane-variance | tested (J on 5 of 16 sessions, without horizontal) |
| plane-timecourse | tested |
| onoff-predicts-kinematics | tested |
| speed-axis-within-distance | tested (J on its 2 vertical sessions) |
| distance-after-speed | tested |
| speed-distance-projection | tested |
| single-trial-variance | tested |
| ellipse-orientation | tested |
| on-greater-than-off | tested |
| relative-deviation-decreases | tested |
| svm-ring-and-distance | tested (J ring on 1 of 4 sessions) |
| svm-3rings-by-ring | tested |
| arc-error-distance-independent | tested |
| distance-error-and-plane-projection | tested |
import json
import time
import urllib.request
import warnings
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from scipy import stats
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import KFold
import premotor
pd.set_option("display.width", 200)
pd.set_option("display.max_columns", 40)
pd.set_option("display.max_rows", 200)
warnings.filterwarnings("ignore", category=RuntimeWarning)
started = time.time()
sessions = premotor.load_all()
print(f"loaded {len(sessions)} sessions from dandiset {premotor.DANDISET_ID} version "
f"{premotor.DANDISET_VERSION} in {time.time() - started:.0f} s")
loaded 12 sessions from dandiset 000121 version 0.220124.2156 in 188 s
Data check¶
Each check of the metadata against the values is stated under its output.
API = "https://api.dandiarchive.org/api/dandisets/000121"
def get(url):
with urllib.request.urlopen(url) as r:
return json.load(r)
versions = [v["version"] for v in get(f"{API}/versions/")["results"]]
print("Versions:", versions)
listing = {}
for v in versions:
assets = get(f"{API}/versions/{v}/assets/?page_size=100")["results"]
listing[v] = sorted((a["path"], a.get("blob"), a["size"]) for a in assets)
print(f" {v}: {len(assets)} assets, {sum(a['size'] for a in assets) / 1e9:.1f} GB")
published = [v for v in versions if v != "draft"]
print("Published versions hold the same files (path, blob, size):",
listing[published[0]] == listing[published[1]])
meta = get(f"{API}/versions/{premotor.DANDISET_VERSION}/")
summary = meta["assetsSummary"]
print("assetsSummary of the pinned version: numberOfFiles", summary["numberOfFiles"],
"; numberOfBytes", f"{summary['numberOfBytes'] / 1e9:.1f} GB",
"; numberOfSubjects", summary["numberOfSubjects"])
print("Description:", meta["description"])
Versions: ['draft', '0.210815.0703', '0.220124.2156'] draft: 0 assets, 0.0 GB
0.210815.0703: 12 assets, 39.5 GB 0.220124.2156: 12 assets, 39.5 GB Published versions hold the same files (path, blob, size): True
assetsSummary of the pinned version: numberOfFiles 15 ; numberOfBytes 187.5 GB ; numberOfSubjects 2 Description: Extracellular physiology data of monkeys implanted with 96 channel Blackrock Utah arrays in the Motor cortex and Premotor cortex. Each file contains a sessions worth of data for one monkey (total 2) performing 1 of the four cursor movement task designs. The data contains hand, eye and cursor position data, LFP, sorted spikes and other task related trialized data.
Both published versions list the same 12 files with the same blobs; they differ only in metadata (schema version, contributor roles). The draft holds no assets. The archive's summary of the pinned version gives 15 files and 187.5 GB, while the version holds 12 files and 39.5 GB.
rows = []
for s in sessions:
t = s["trials"]
co = t[t.center_out]
delayed = co[co.successful & (co.delay >= premotor.MIN_DELAY_S)].delay
rows.append({
"monkey": s["monkey"], "date": s["date"], "trials": len(t), "successful": int(t.successful.sum()),
"center_out": len(co), "succ_center_out": int(co.successful.sum()),
"succ_co_no_delay": int((co.successful & (co.delay < premotor.MIN_DELAY_S)).sum()),
"delay_min_ms": round(1000 * delayed.min()), "delay_max_ms": round(1000 * delayed.max()),
"off_layout": int((co.successful & ~co.in_layout).sum()), "analysed": int(t.analysed.sum()),
"task_id": sorted(t.task_id.unique().tolist()), "version_id": sorted(t.version_id.unique().tolist()),
"targets": co[["tx", "ty"]].round(1).drop_duplicates().shape[0],
"hand_samples": s["n_hand_samples"], "hand_dt_ms": round(1000 * s["hand_dt_median"], 3),
"recorded_s": round(s["recorded_s"]),
"cursor_minus_hand_mm_p50_p99_max": tuple(np.round(s["meta"]["cursor_minus_hand_mm"], 2))})
overview = pd.DataFrame(rows)
print(overview.to_string(index=False))
print("\nDelays of successful center-out trials below", premotor.MIN_DELAY_S, "s, by value (ms):")
for m in "JR":
d = pd.concat([s["trials"] for s in sessions if s["monkey"] == m])
d = d[d.successful & d.center_out & (d.delay < premotor.MIN_DELAY_S)].delay
print(" ", m, (1000 * d).round().value_counts().sort_index().to_dict())
print("\nTarget positions (stated unit none; values):",
{s["date"]: (round(s["trials"].tx.min()), round(s["trials"].tx.max()), round(s["trials"].ty.min()),
round(s["trials"].ty.max())) for s in sessions})
monkey date trials successful center_out succ_center_out succ_co_no_delay delay_min_ms delay_max_ms off_layout analysed task_id version_id targets hand_samples hand_dt_ms recorded_s cursor_minus_hand_mm_p50_p99_max
J 2015-10-01 3426 3387 1709 1678 314 300 699 0 1364 [2.0] [21.0] 24 3985838 1.0 3986 (0.0, 1.86, 1439.74)
J 2015-10-02 4612 4580 2304 2275 428 300 700 0 1847 [2.0] [21.0] 24 5469603 1.0 5470 (0.0, 2.1, 1415.09)
J 2015-10-15 2654 2639 1327 1313 293 300 700 0 1020 [2.0] [21.0] 24 3082570 1.0 3083 (0.0, 2.83, 1338.29)
J 2016-01-27 4667 4571 2318 2253 443 300 700 0 1810 [2.0] [21.0] 48 5878154 1.0 5878 (0.0, 2.45, 1438.13)
J 2016-01-28 7020 6899 3499 3400 686 300 700 0 2714 [2.0] [21.0] 48 8106459 1.0 8106 (0.0, 2.59, 1461.43)
R 2017-01-15 5333 4407 2567 1840 0 400 999 0 1840 [2.0] [23.0] 24 6519979 1.0 6520 (0.0, 1.88, 1432.46)
R 2017-01-16 5040 4467 2466 2002 0 400 1000 0 2002 [2.0] [23.0] 24 6064247 1.0 6064 (0.0, 1.89, 1434.41)
R 2017-01-17 5398 4653 2678 1975 0 400 1000 69 1906 [2.0] [23.0] 30 6076796 1.0 6077 (0.0, 2.18, 1440.64)
R 2017-01-18 5188 4753 2572 2181 0 400 1000 0 2181 [2.0] [23.0] 24 6115923 1.0 6116 (0.0, 2.1, 1446.11)
R 2017-01-19 3701 3069 1703 1366 26 401 999 0 1340 [2.0] [23.0] 48 4914555 1.0 4915 (0.0, 1349.59, 1902.68)
R 2017-01-24 4303 3554 1948 1606 0 400 999 0 1606 [2.0] [23.0] 36 5619477 1.0 5619 (0.0, 2.08, 1427.09)
R 2017-01-25 5274 4598 2620 1978 0 401 1000 0 1978 [2.0] [23.0] 36 6154149 1.0 6154 (0.0, 2.35, 1444.99)
Delays of successful center-out trials below 0.25 s, by value (ms):
J {0.0: 1104, 1.0: 1060}
R {0.0: 26}
Target positions (stated unit none; values): {'2015-10-01': (0, 0, -120, 120), '2015-10-02': (0, 0, -120, 120), '2015-10-15': (-80, 80, -80, 80), '2016-01-27': (-120, 120, -120, 120), '2016-01-28': (-120, 120, -120, 120), '2017-01-15': (-120, 120, 0, 0), '2017-01-16': (-120, 120, 0, 0), '2017-01-17': (-80, 80, -120, 120), '2017-01-18': (0, 0, -120, 120), '2017-01-19': (-120, 120, -120, 120), '2017-01-24': (-80, 80, -80, 80), '2017-01-25': (-80, 80, -80, 80)}
task_id is 2 and version_id one value per monkey in every trial, so the files do not say which
target configuration a session used; it is inferred from the target positions (see
session-configurations). Target positions run to 120, the 12 cm of the outer ring, so they are
in mm. The hand and cursor positions are stated in meters with a conversion of NaN, but take
the same mm values (the target of 80 is reached at about 80, see the clock check below).
Cursor and hand agree to within 1.9 to 2.8 mm in 99% of samples in 11 sessions; in R's
2017-01-19 they differ by more than 1 m in over 1% of samples (up to 1903 mm, beyond any
screen position). Only the hand is used here. Hand samples are 1 ms apart and exist only within trials.
Delays of delayed trials run from 300 to 700 ms for J and from 400 to 1000 ms for R. Trials
without a delay have delays of 0 or 1 ms; J has them in every session, R only 26, all on
2017-01-19.
rows = []
for s in sessions:
u = s["units"]
pmd, lab = u.iloc[premotor.PMD_UNITS], u.iloc[premotor.FILE_PMD_UNITS]
rows.append({
"monkey": s["monkey"], "date": s["date"],
"units_0_95_label": "/".join(sorted(set(lab.group_label))),
"units_0_95_electrode_rows": f"{lab.electrode.min()}-{lab.electrode.max()} ({'/'.join(sorted(set(lab.electrode_group)))})",
"units_96_191_label": "/".join(sorted(set(pmd.group_label))),
"units_96_191_electrode_rows": f"{pmd.electrode.min()}-{pmd.electrode.max()} ({'/'.join(sorted(set(pmd.electrode_group)))})",
"lfp": s["meta"]["lfp_electrodes"],
"live_0_95": int((lab.n_spikes > 0).sum()), "live_96_191": int((pmd.n_spikes > 0).sum()),
"rate_0_95": round(lab.rate.median(), 1), "rate_96_191": round(pmd.rate.median(), 1),
"short_span_96_191": int((pmd["last"] - pmd["first"] < 0.9 * s["recorded_s"]).sum()),
"spikes_in_trials_96_191": round(pmd.in_trials.sum() / pmd.n_spikes.sum(), 4),
"trial_s/recorded_s/span_s": (round(s["trial_s"]), round(s["recorded_s"]), round(s["span_s"])),
"units_with_time_decreasing": [int(i) for i in np.flatnonzero(u.time_decreases)],
"unit0_span_s": (round(u["first"][0], 1), round(u["last"][0], 1)),
"unit1_n_spikes_vs_median": (int(u.n_spikes[1]), int(u.n_spikes.median()))})
print(pd.DataFrame(rows).to_string(index=False))
print("\nElectrode groups: 1 =", "M1", "; 2 = PMd (from general/extracellular_ephys)")
monkey date units_0_95_label units_0_95_electrode_rows units_96_191_label units_96_191_electrode_rows lfp live_0_95 live_96_191 rate_0_95 rate_96_191 short_span_96_191 spikes_in_trials_96_191 trial_s/recorded_s/span_s units_with_time_decreasing unit0_span_s unit1_n_spikes_vs_median
J 2015-10-01 2 0-95 (1) 1 96-191 (2) {'M1_1': (0, 95), 'M1_2': (0, 95), 'PMd_1': (96, 191), 'PMd_2': (96, 191), 'PMd_3': (96, 191)} 96 96 16.8 12.0 0 0.9993 (3982, 3986, 4736) [1] (0.0, 1323.0) (173016, 55736)
J 2015-10-02 2 0-95 (1) 1 96-191 (2) {'M1_1': (0, 95), 'M1_2': (0, 95), 'PMd_1': (96, 191), 'PMd_2': (96, 191), 'PMd_3': (96, 191)} 96 96 14.5 4.1 0 0.9994 (5465, 5470, 6479) [1] (0.1, 1543.9) (238073, 44293)
J 2015-10-15 2 0-95 (1) 1 96-191 (2) {'M1_1': (0, 95), 'M1_2': (0, 95), 'PMd_1': (96, 191), 'PMd_2': (96, 191)} 96 96 16.4 10.1 0 0.9993 (3080, 3083, 3664) [1] (0.1, 1422.8) (148937, 41504)
J 2016-01-27 2 0-95 (1) 1 96-191 (2) {'M1_1': (0, 95), 'M1_2': (0, 95), 'PMd_1': (96, 191), 'PMd_2': (96, 191), 'PMd_3': (96, 191)} 96 96 13.6 8.6 0 0.9994 (5873, 5878, 6900) [1] (0.1, 3119.5) (205715, 66638)
J 2016-01-28 2 0-95 (1) 1 96-191 (2) {'M1_1': (0, 95), 'M1_2': (0, 95), 'M1_3': (0, 95), 'PMd_1': (96, 191), 'PMd_2': (96, 191), 'PMd_3': (96, 191), 'PMd_4': (96, 191)} 96 96 14.7 11.0 0 0.9994 (8099, 8106, 9644) [1] (0.0, 2620.1) (192841, 97850)
R 2017-01-15 2 0-95 (1) 1 96-191 (2) {'M1_1': (0, 95), 'PMd_1': (96, 191)} 67 95 7.7 7.4 0 0.9995 (6515, 6520, 8299) [] (0.2, 8299.4) (34598, 48973)
R 2017-01-16 2 0-95 (1) 1 96-191 (2) {'M1_1': (0, 95), 'PMd_1': (96, 191)} 67 95 7.9 20.7 0 0.9994 (6059, 6064, 7157) [] (0.3, 7157.2) (30197, 97492)
R 2017-01-17 2 0-95 (1) 1 96-191 (2) {'M1_3': (0, 95), 'PMd_3': (96, 191)} 67 95 4.8 4.6 0 0.9994 (6071, 6077, 7039) [] (0.4, 7038.5) (29974, 28195)
R 2017-01-18 2 0-95 (1) 1 96-191 (2) {'M1_1': (0, 95), 'PMd_1': (96, 191)} 67 95 5.3 33.8 0 0.9994 (6111, 6116, 6781) [] (0.1, 6780.3) (32549, 125650)
R 2017-01-19 2 0-95 (1) 1 96-191 (2) {'M1_1': (0, 95), 'PMd_1': (96, 191)} 67 95 6.7 47.4 0 0.9995 (4911, 4915, 6699) [] (0.2, 6698.7) (28498, 136643)
R 2017-01-24 2 0-95 (1) 1 96-191 (2) {'M1_1': (0, 95), 'PMd_1': (96, 191)} 67 95 7.8 19.2 0 0.9995 (5615, 5619, 7578) [] (0.0, 7577.1) (31009, 79737)
R 2017-01-25 2 0-95 (1) 1 96-191 (2) {'M1_1': (0, 95), 'PMd_1': (96, 191)} 67 95 8.5 40.9 0 0.9994 (6149, 6154, 6633) [] (0.0, 6632.1) (37341, 152508)
Electrode groups: 1 = M1 ; 2 = PMd (from general/extracellular_ephys)
Units 0 to 95 carry electrode_group 2 (PMd), but their electrodes column points to electrode
rows 0 to 95, which belong to group 1 (M1); units 96 to 191 are labelled M1 and point to rows of
group 2 (PMd), the rows the PMd LFP uses. The conversion code assigns electrode_group by the
order of a dictionary of groups, not from the recording. The next cell asks which array carries
the delay-period target signal.
In every J file unit 1's spike times decrease once, from the end of the session back to its
start, unit 1 holds two to five times the median number of spikes, and unit 0 stops long before
the session ends: the boundary between units 0 and 1 in spike_times_index is misplaced, and
unit 1 holds the rest of unit 0's train followed by its own. Both are on the array of units 0 to
95, not the one used here as PMd. The loader sorts any unsorted train before counting.
Spikes, like hand samples, exist only within trials: 99.93 to 99.95% of the PMd spikes of every session fall inside a trial interval, although each file's span exceeds its summed trial time by 484 to 1963 s (the gaps between trials hold almost no spikes). The recorded time of the hand samples equals the summed trial time to within 7 s, so the day's mean rate, the spike count over the recorded time, is the rate over trial time. The per-session median PMd rates still vary widely, from 4.1 to 12.0 Hz in J and 4.6 to 47.4 Hz in R.
from sklearn.model_selection import StratifiedKFold # noqa: E402
ARRAY_SEED = 1
rows = []
for s in sessions:
t = s["trials"][s["trials"].analysed]
y = (np.round(np.degrees(np.arctan2(t.ty, t.tx))) * 1000 + np.round(np.hypot(t.tx, t.ty))).values
keep = pd.Series(y).map(pd.Series(y).value_counts()) >= 5
for name, units in (("units 0-95", premotor.FILE_PMD_UNITS), ("units 96-191", premotor.PMD_UNITS)):
row = {"monkey": s["monkey"], "date": s["date"], "array": name}
for window, counts in (("delay", s["delay_counts"]), ("move", s["move_counts"])):
X = counts[:, units].astype(float)
Z = np.c_[np.ones(len(t)), t.tx, t.ty]
n_sig = 0
for ch in range(X.shape[1]):
if X[:, ch].std() == 0:
continue
beta = np.linalg.lstsq(Z, X[:, ch], rcond=None)[0]
rss = ((X[:, ch] - Z @ beta) ** 2).sum()
tss = ((X[:, ch] - X[:, ch].mean()) ** 2).sum()
F = ((tss - rss) / 2) / (rss / (len(X) - 3))
n_sig += stats.f.sf(F, 2, len(X) - 3) < 0.05 / 96
_, acc = premotor.svm_cv(X[keep.values], y[keep.values], ARRAY_SEED, folds=5)
row[f"{window}_sig_channels"] = n_sig
row[f"{window}_accuracy"] = round(acc.mean(), 3)
row["chance"] = round(1 / len(np.unique(y[keep.values])), 3)
rows.append(row)
arrays = pd.DataFrame(rows)
print("Channels whose count regresses on target x, y (F-test, p < 0.05/96), and 5-fold SVM accuracy "
"for the target, in the 200 ms before the go cue (delay) and 200 to 400 ms after it (move):")
print(arrays.to_string(index=False))
ratio = arrays.assign(r=arrays.delay_accuracy / arrays.move_accuracy).pivot_table(
index="monkey", columns="array", values=["delay_sig_channels", "delay_accuracy", "r"], aggfunc="median")
print("\nMedians over sessions (r = delay accuracy / move accuracy):")
print(ratio.round(3).to_string())
Channels whose count regresses on target x, y (F-test, p < 0.05/96), and 5-fold SVM accuracy for the target, in the 200 ms before the go cue (delay) and 200 to 400 ms after it (move):
monkey date array delay_sig_channels delay_accuracy move_sig_channels move_accuracy chance
J 2015-10-01 units 0-95 34 0.108 58 0.157 0.042
J 2015-10-01 units 96-191 41 0.119 70 0.187 0.042
J 2015-10-02 units 0-95 40 0.105 52 0.153 0.042
J 2015-10-02 units 96-191 39 0.123 69 0.183 0.042
J 2015-10-15 units 0-95 48 0.126 75 0.305 0.042
J 2015-10-15 units 96-191 60 0.245 70 0.348 0.042
J 2016-01-27 units 0-95 51 0.061 84 0.188 0.021
J 2016-01-27 units 96-191 64 0.143 68 0.206 0.021
J 2016-01-28 units 0-95 59 0.095 87 0.302 0.021
J 2016-01-28 units 96-191 70 0.177 76 0.291 0.021
R 2017-01-15 units 0-95 0 0.049 16 0.069 0.042
R 2017-01-15 units 96-191 56 0.141 58 0.164 0.042
R 2017-01-16 units 0-95 1 0.048 11 0.062 0.042
R 2017-01-16 units 96-191 56 0.132 61 0.174 0.042
R 2017-01-17 units 0-95 0 0.054 0 0.071 0.042
R 2017-01-17 units 96-191 30 0.126 57 0.157 0.042
R 2017-01-18 units 0-95 6 0.048 0 0.057 0.042
R 2017-01-18 units 96-191 64 0.119 62 0.160 0.042
R 2017-01-19 units 0-95 2 0.020 0 0.015 0.021
R 2017-01-19 units 96-191 66 0.090 73 0.151 0.021
R 2017-01-24 units 0-95 1 0.031 1 0.030 0.028
R 2017-01-24 units 96-191 66 0.113 72 0.171 0.028
R 2017-01-25 units 0-95 5 0.026 0 0.024 0.028
R 2017-01-25 units 96-191 70 0.117 82 0.188 0.028
Medians over sessions (r = delay accuracy / move accuracy):
delay_accuracy delay_sig_channels r
array units 0-95 units 96-191 units 0-95 units 96-191 units 0-95 units 96-191
monkey
J 0.105 0.143 48.0 60.0 0.413 0.672
R 0.048 0.119 1.0 64.0 0.842 0.744
In R, units 0 to 95 carry almost no target information (at most 6 significant channels in the
delay window and 16 after the go cue, accuracy at most 0.03 above chance) and 29 of them never
fire: this fits the M1 array the paper calls too poor to use. Units 96 to 191 carry it in every
session (30 to 70 significant channels in the delay window). In J both arrays carry it; units 96
to 191 have the higher delay-window accuracy in all five sessions, more significant channels in
four of five, and a higher ratio of delay to movement accuracy (median 0.67 against 0.41), as
expected of PMd against M1. The data therefore agree
with the electrodes column and the LFP, not with the units' electrode_group: the analysis
takes units 96 to 191 as PMd. This is consistent with the data, not proven by them.
CLOCK_SEED = 2
rows = []
for s in sessions:
t = s["trials"][s["trials"].analysed]
y = (np.round(np.degrees(np.arctan2(t.ty, t.tx))) * 1000 + np.round(np.hypot(t.tx, t.ty))).values
keep = (pd.Series(y).map(pd.Series(y).value_counts()) >= 5).values
tc = s["tc_counts"].astype(float)
centers = premotor.TC_START_S + premotor.BIN_S * (np.arange(tc.shape[2]) + 0.5)
windows = {"-200 to 0 ms (before target)": (centers < 0),
"100 to 300 ms after target": (centers > 0.1) & (centers < 0.3)}
row = {"monkey": s["monkey"], "date": s["date"], "chance": round(1 / len(np.unique(y[keep])), 3)}
for name, w in windows.items():
X = tc[:, :, w].sum(axis=2)
row[name] = round(premotor.svm_cv(X[keep], y[keep], CLOCK_SEED, folds=5)[1].mean(), 3)
X = s["delay_counts"][:, premotor.PMD_UNITS].astype(float)
row["delay window"] = round(premotor.svm_cv(X[keep], y[keep], CLOCK_SEED, folds=5)[1].mean(), 3)
shifted = np.roll(y, 1)
keep2 = (pd.Series(shifted).map(pd.Series(shifted).value_counts()) >= 5).values
row["delay window, labels of the previous trial"] = round(
premotor.svm_cv(X[keep2], shifted[keep2], CLOCK_SEED, folds=5)[1].mean(), 3)
k = s["trials"][s["trials"].successful & s["trials"].center_out & s["trials"].max_speed.notna()]
row["onset before go cue (fraction)"] = round(float((k.reaction_time < 0).mean()), 3)
row["median reaction time (ms)"] = round(1000 * k.reaction_time.median())
row["hand at max speed, median along target (mm)"] = round(float(np.median(
(k.x_at_max * k.tx + k.y_at_max * k.ty) / np.hypot(k.tx, k.ty))))
rows.append(row)
clock = pd.DataFrame(rows)
print("5-fold SVM accuracy for the target from PMd counts, by window, and hand timing:")
print(clock.to_string(index=False))
5-fold SVM accuracy for the target from PMd counts, by window, and hand timing:
monkey date chance -200 to 0 ms (before target) 100 to 300 ms after target delay window delay window, labels of the previous trial onset before go cue (fraction) median reaction time (ms) hand at max speed, median along target (mm)
J 2015-10-01 0.042 0.041 0.114 0.119 0.062 0.035 349 31
J 2015-10-02 0.042 0.046 0.108 0.109 0.065 0.046 344 32
J 2015-10-15 0.042 0.045 0.240 0.249 0.060 0.000 331 35
J 2016-01-27 0.021 0.015 0.134 0.141 0.020 0.000 343 35
J 2016-01-28 0.021 0.021 0.187 0.184 0.024 0.000 306 36
R 2017-01-15 0.042 0.042 0.141 0.147 0.049 0.111 320 33
R 2017-01-16 0.042 0.033 0.131 0.127 0.053 0.095 345 32
R 2017-01-17 0.042 0.038 0.108 0.120 0.048 0.069 367 37
R 2017-01-18 0.042 0.043 0.115 0.125 0.053 0.071 371 37
R 2017-01-19 0.021 0.018 0.145 0.088 0.019 0.019 355 42
R 2017-01-24 0.028 0.025 0.123 0.121 0.031 0.011 379 41
R 2017-01-25 0.028 0.030 0.126 0.127 0.027 0.014 369 40
In every session the target cannot be decoded in the 200 ms before target onset (accuracy at chance), and can be 100 to 300 ms after it and in the 200 ms before the go cue; with each trial given the previous trial's target, accuracy in the delay window falls to near chance (at most 0.023 above it, against 0.067 to 0.207 above it with the true targets). The spike times therefore share the clock of the trial table, to the resolution of these windows. Median reaction times are 306 to 379 ms in every session, and the hand at maximum speed lies along the target direction (median 31 to 42 mm along it), so the hand shares the clock too. Movement onset (10% of peak speed) precedes the go cue in 7 to 11% of trials in R's four horizontal and vertical sessions and in at most 5% elsewhere; these trials are kept.
rows = []
for m in "JR":
ss = [s for s in sessions if s["monkey"] == m]
rates = np.array([s["units"].rate.values for s in ss])
counts = [len(s["trials"]) for s in ss]
c = np.corrcoef(rates)
np.fill_diagonal(c, np.nan)
print(f"{m}: trials per session {counts}; hand samples {[s['n_hand_samples'] for s in ss]}")
print(f" highest correlation between sessions of the 192 channel rates: {np.nanmax(c):.3f}; "
f"identical trial counts: {len(set(counts)) < len(counts)}")
first = [tuple(np.round(s["trials"][["go_cue_time", "tx", "ty"]].dropna().iloc[:5].values.ravel(), 3)) for s in ss]
print(" identical first five trials in two sessions:", len(set(first)) < len(first))
J: trials per session [3426, 4612, 2654, 4667, 7020]; hand samples [3985838, 5469603, 3082570, 5878154, 8106459] highest correlation between sessions of the 192 channel rates: 0.905; identical trial counts: False identical first five trials in two sessions: False R: trials per session [5333, 5040, 5398, 5188, 3701, 4303, 5274]; hand samples [6519979, 6064247, 6076796, 6115923, 4914555, 5619477, 6154149] highest correlation between sessions of the 192 channel rates: 0.933; identical trial counts: False identical first five trials in two sessions: False
No two sessions share trial or sample counts, their first trials, or their channel rates, so no recording was deposited twice. The analysis combines sessions of a monkey and treats each channel (electrode) as the same across sessions, as the paper does; this is an assumption the data cannot test, since threshold crossings carry no identity. Results that rest on it are also given per session where cheap.
for m in "JR":
d = pd.concat([s["trials"] for s in sessions if s["monkey"] == m])
k = d[d.successful & d.center_out]
print(f"{m}: successful center-out trials {len(k)}; with a hand glitch (|position| > "
f"{premotor.HAND_LIMIT_MM:.0f} mm) {int(k.glitch.sum())}; without kinematics {int(k.max_speed.isna().sum())}")
print(f" max speed (mm/s) range {k.max_speed.min():.0f} to {k.max_speed.max():.0f}, median {k.max_speed.median():.0f}; "
f"reaction time (ms) 5th to 95th percentile {1000 * k.reaction_time.quantile(0.05):.0f} to "
f"{1000 * k.reaction_time.quantile(0.95):.0f}")
a, _ = premotor.monkey_table(sessions, m)
size = a.groupby("condition").size()
print(f" analysed trials {len(a)} in {len(size)} conditions; trials per condition {size.min()} to {size.max()}")
J: successful center-out trials 10919; with a hand glitch (|position| > 250 mm) 67; without kinematics 67 max speed (mm/s) range 2 to 1117, median 389; reaction time (ms) 5th to 95th percentile 256 to 453 analysed trials 8755 in 96 conditions; trials per condition 36 to 142 R: successful center-out trials 12948; with a hand glitch (|position| > 250 mm) 25; without kinematics 25 max speed (mm/s) range 3 to 1398, median 462; reaction time (ms) 5th to 95th percentile -16 to 515 analysed trials 12853 in 132 conditions; trials per condition 26 to 181
session-configurations¶
The sessions are those the Methods list for each target configuration; for monkey R, ring on 2017-01-24 and 25, horizontal on 2017-01-17 and 18, vertical on 2017-01-15 and 16, and 3-rings on 2017-01-19; for J, ring, vertical and 3-rings on the dates listed.
rows = []
for s in sessions:
t = s["trials"]
co = t[t.center_out & t.successful]
off = co[~co.in_layout]
rows.append({"monkey": s["monkey"], "date": s["date"], "configuration in the file": s["config"],
"Methods": s["methods_config"], "match": s["config"] == s["methods_config"],
"radii (mm)": sorted(np.unique(np.round(np.hypot(co.tx, co.ty))).astype(int).tolist()),
"targets": co[["tx", "ty"]].round(1).drop_duplicates().shape[0],
"trials off the layout": len(off),
"their targets": sorted(set(zip(off.tx.round(1), off.ty.round(1))))})
configs = pd.DataFrame(rows)
print(configs.to_string(index=False))
print("\nSessions matching the Methods:", int(configs.match.sum()), "of", len(configs))
swapped = configs[~configs.match]
print("Mismatches that are horizontal and vertical exchanged:",
int((swapped["configuration in the file"].isin(["horizontal", "vertical"])
& swapped.Methods.isin(["horizontal", "vertical"])).sum()), "of", len(swapped))
for m in "JR":
listed = sorted(premotor.METHODS_CONFIG[m])
present = sorted(s["date"] for s in sessions if s["monkey"] == m)
print(f"\n{m}: the Methods list {len(listed)} sessions, the dandiset holds {len(present)}; missing:",
[d for d in listed if d not in present])
monkey date configuration in the file Methods match radii (mm) targets trials off the layout their targets
J 2015-10-01 vertical vertical True [10, 20, 30, 40, 50, 60, 70, 80, 90, 100, 110, 120] 24 0 []
J 2015-10-02 vertical vertical True [10, 20, 30, 40, 50, 60, 70, 80, 90, 100, 110, 120] 24 0 []
J 2015-10-15 ring ring True [80] 24 0 []
J 2016-01-27 3-rings 3-rings True [40, 80, 120] 48 0 []
J 2016-01-28 3-rings 3-rings True [40, 80, 120] 48 0 []
R 2017-01-15 horizontal vertical False [5, 16, 26, 37, 47, 57, 68, 78, 89, 99, 110, 120] 24 0 []
R 2017-01-16 horizontal vertical False [5, 16, 26, 37, 47, 57, 68, 78, 89, 99, 110, 120] 24 0 []
R 2017-01-17 vertical horizontal False [10, 20, 30, 40, 50, 60, 70, 80, 90, 100, 110, 120] 30 69 [(-80.0, 0.0), (-56.6, -56.6), (-56.6, 56.6), (56.6, -56.6), (56.6, 56.6), (80.0, 0.0)]
R 2017-01-18 vertical horizontal False [10, 20, 30, 40, 50, 60, 70, 80, 90, 100, 110, 120] 24 0 []
R 2017-01-19 3-rings 3-rings True [40, 80, 120] 48 0 []
R 2017-01-24 ring ring True [80] 36 0 []
R 2017-01-25 ring ring True [80] 36 0 []
Sessions matching the Methods: 8 of 12
Mismatches that are horizontal and vertical exchanged: 4 of 4
J: the Methods list 16 sessions, the dandiset holds 5; missing: ['2015-09-24', '2015-09-25', '2015-09-29', '2015-09-30', '2015-10-07', '2015-10-08', '2015-10-09', '2015-10-26', '2015-10-27', '2016-01-29', '2016-02-02']
R: the Methods list 7 sessions, the dandiset holds 7; missing: []
Result: partially reproduced. 8 of 12 sessions have the configuration the Methods give for their date. The other 4 are R's: 2017-01-15 and 16 hold horizontal targets where the Methods say vertical, and 2017-01-17 and 18 hold vertical targets where the Methods say horizontal, so all four mismatches are horizontal and vertical exchanged. 2017-01-17 also holds 69 successful trials to six ring targets outside its layout. The dandiset holds all 7 of R's sessions and 5 of J's 16; the 11 missing include all four of J's horizontal sessions.
r-trial-counts¶
Monkey R performed 3577 ring trials over two sessions, 3699 horizontal over two, 4040 vertical over two, and 912 3-rings trials in one session.
R_EXPECTED = {"ring": 3577, "horizontal": 3699, "vertical": 4040, "3-rings": 912}
rows = []
for s in sessions:
if s["monkey"] != "R":
continue
t = s["trials"]
co = t[t.center_out]
rows.append({"date": s["date"], "Methods": s["methods_config"], "file": s["config"],
"(a) successful, delayed": int((co.successful & (co.delay >= premotor.MIN_DELAY_S)).sum()),
"(b) delayed, any outcome": int((co.delay >= premotor.MIN_DELAY_S).sum()),
"(c) successful, with no-delay": int(co.successful.sum()),
"(b) also without a go cue": int((~(co.delay < premotor.MIN_DELAY_S)).sum()),
"(a) in the layout only": int((co.successful & (co.delay >= premotor.MIN_DELAY_S) & co.in_layout).sum())})
counts = pd.DataFrame(rows)
print(counts.to_string(index=False))
readings = ["(a) successful, delayed", "(b) delayed, any outcome", "(c) successful, with no-delay",
"(b) also without a go cue", "(a) in the layout only"]
for grouping in ("Methods", "file"):
g = counts.groupby(grouping)[readings].sum().reindex(list(R_EXPECTED))
g["paper"] = pd.Series(R_EXPECTED)
print(f"\nGrouped by the configuration of the {grouping}:")
print(g.to_string())
for r in readings:
ok = (g[r] - g.paper).abs() <= 0.1 * g.paper
print(f" {r}: within 10% for {int(ok.sum())} of 4 ({', '.join(g.index[ok])})")
date Methods file (a) successful, delayed (b) delayed, any outcome (c) successful, with no-delay (b) also without a go cue (a) in the layout only
2017-01-15 vertical horizontal 1840 1845 1840 2567 1840
2017-01-16 vertical horizontal 2002 2011 2002 2466 2002
2017-01-17 horizontal vertical 1975 1980 1975 2678 1906
2017-01-18 horizontal vertical 2181 2185 2181 2572 2181
2017-01-19 3-rings 3-rings 1340 1345 1366 1671 1340
2017-01-24 ring ring 1606 1613 1606 1948 1606
2017-01-25 ring ring 1978 1981 1978 2620 1978
Grouped by the configuration of the Methods:
(a) successful, delayed (b) delayed, any outcome (c) successful, with no-delay (b) also without a go cue (a) in the layout only paper
Methods
ring 3584 3594 3584 4568 3584 3577
horizontal 4156 4165 4156 5250 4087 3699
vertical 3842 3856 3842 5033 3842 4040
3-rings 1340 1345 1366 1671 1340 912
(a) successful, delayed: within 10% for 2 of 4 (ring, vertical)
(b) delayed, any outcome: within 10% for 2 of 4 (ring, vertical)
(c) successful, with no-delay: within 10% for 2 of 4 (ring, vertical)
(b) also without a go cue: within 10% for 0 of 4 ()
(a) in the layout only: within 10% for 2 of 4 (ring, vertical)
Grouped by the configuration of the file:
(a) successful, delayed (b) delayed, any outcome (c) successful, with no-delay (b) also without a go cue (a) in the layout only paper
file
ring 3584 3594 3584 4568 3584 3577
horizontal 3842 3856 3842 5033 3842 3699
vertical 4156 4165 4156 5250 4087 4040
3-rings 1340 1345 1366 1671 1340 912
(a) successful, delayed: within 10% for 3 of 4 (ring, horizontal, vertical)
(b) delayed, any outcome: within 10% for 3 of 4 (ring, horizontal, vertical)
(c) successful, with no-delay: within 10% for 3 of 4 (ring, horizontal, vertical)
(b) also without a go cue: within 10% for 0 of 4 ()
(a) in the layout only: within 10% for 3 of 4 (ring, horizontal, vertical)
Result: partially reproduced. Grouped by the Methods' dates, two of the four counts are within 10% under readings (a) to (c): ring 3584 against 3577 and vertical 3842 against 4040 under reading (a). Horizontal gives 4156 against 3699 and 3-rings 1340 against 912. Reading (b) drops failed trials without a go cue (delay undefined); counting them as well gives no match (ring 4568, horizontal 5250, vertical 5033, 3-rings 1671). Grouped by the geometry in the files, three of four are within 10% (horizontal 3842 against 3699, vertical 4156 against 4040), which fits the exchange of R's horizontal and vertical dates found above. The 3-rings session holds 1340 successful delayed trials where the paper used 912.
delays-and-times¶
Delays were uniform from 300 to 700 ms for J and 400 to 900 ms for R, 10% of trials had no delay; average reaction time was 346 ms (J) and 311 ms (R), average reach time 582 ms (J) and 482 ms (R), and reaction time was negatively correlated with delay (p < 1e-4 for both).
for m in "JR":
d = pd.concat([s["trials"].assign(date=s["date"]) for s in sessions if s["monkey"] == m])
k = d[d.successful & d.center_out]
delayed = k[k.delay >= premotor.MIN_DELAY_S]
kk = delayed[delayed.max_speed.notna() & delayed.onset.notna()]
rt = 1000 * kk.reaction_time
reach_a = 1000 * (kk.first_acquire - kk.onset)
reach_b = 1000 * (kk.first_acquire - kk.go_cue_time)
r, p = stats.pearsonr(kk.delay, kk.reaction_time)
print(f"{m}: delays of delayed trials {1000 * delayed.delay.min():.0f} to {1000 * delayed.delay.max():.0f} ms; "
f"deciles {np.round(1000 * delayed.delay.quantile(np.linspace(0, 1, 11))).astype(int).tolist()}")
print(f" no-delay fraction of successful center-out trials {np.mean(k.delay < premotor.MIN_DELAY_S):.3f} "
f"({int((k.delay < premotor.MIN_DELAY_S).sum())} of {len(k)})")
print(f" reaction time mean {rt.mean():.0f} ms (median {rt.median():.0f}); reach time (a) acquisition - onset "
f"mean {reach_a.mean():.0f} ms (median {reach_a.median():.0f}); (b) acquisition - go cue mean "
f"{reach_b.mean():.0f} ms (median {reach_b.median():.0f}); trials {len(kk)}, without kinematics "
f"{int(delayed.max_speed.isna().sum())}, without an onset {int((delayed.max_speed.notna() & delayed.onset.isna()).sum())}, "
f"without acquisition {int(kk.first_acquire.isna().sum())}")
print(f" reaction time against delay: r = {r:.3f}, p = {p:.1e}")
rt_trim = kk[(kk.reaction_time > 0.1) & (kk.reaction_time < 0.8)]
r2, p2 = stats.pearsonr(rt_trim.delay, rt_trim.reaction_time)
print(f" variant, reaction times 100 to 800 ms only ({len(rt_trim)} trials): mean {1000 * rt_trim.reaction_time.mean():.0f} ms, "
f"r = {r2:.3f}, p = {p2:.1e}")
v = kk[kk.onset_after_go.notna()]
rt_go = 1000 * (v.onset_after_go - v.go_cue_time)
rg, pg = stats.pearsonr(v.delay, v.onset_after_go - v.go_cue_time)
print(f" onset before the go cue in {np.mean(kk.reaction_time < 0):.3f} of trials; variant, onset searched only after "
f"the go cue: reaction time mean {rt_go.mean():.0f} ms (median {rt_go.median():.0f}), reach time (b) unchanged, "
f"trials without such an onset {len(kk) - len(v)} of {len(kk)}; r with delay {rg:.3f}, p = {pg:.1e}")
for key, label in (("onset_05", "5%"), ("onset_20", "20%")):
v = kk[kk[key].notna()]
rv, pv = stats.pearsonr(v.delay, v[key] - v.go_cue_time)
print(f" onset at {label} of peak speed: reaction time mean {1000 * (v[key] - v.go_cue_time).mean():.0f} ms, "
f"reach time (a) mean {1000 * (v.first_acquire - v[key]).mean():.0f} ms; r with delay {rv:.3f}, p = {pv:.1e}")
per = kk.groupby("date").apply(lambda x: pd.Series({
"rt_mean": 1000 * x.reaction_time.mean(), "reach_a": 1000 * (x.first_acquire - x.onset).mean(),
"reach_b": 1000 * (x.first_acquire - x.go_cue_time).mean(),
"r_delay": stats.pearsonr(x.delay, x.reaction_time)[0], "p": stats.pearsonr(x.delay, x.reaction_time)[1]}))
per[["rt_mean", "reach_a", "reach_b"]] = per[["rt_mean", "reach_a", "reach_b"]].round(0)
per["r_delay"] = per.r_delay.round(3)
print(per.to_string(formatters={"p": "{:.1e}".format}))
J: delays of delayed trials 300 to 700 ms; deciles [300, 338, 379, 420, 459, 500, 539, 581, 620, 660, 700]
no-delay fraction of successful center-out trials 0.198 (2164 of 10919)
reaction time mean 332 ms (median 329); reach time (a) acquisition - onset mean 255 ms (median 239); (b) acquisition - go cue mean 587 ms (median 575); trials 8700, without kinematics 54, without an onset 1, without acquisition 0
reaction time against delay: r = -0.013, p = 2.1e-01
variant, reaction times 100 to 800 ms only (8512 trials): mean 341 ms, r = -0.017, p = 1.2e-01
onset before the go cue in 0.014 of trials; variant, onset searched only after the go cue: reaction time mean 342 ms (median 330), reach time (b) unchanged, trials without such an onset 124 of 8700; r with delay -0.029, p = 7.9e-03
onset at 5% of peak speed: reaction time mean 303 ms, reach time (a) mean 286 ms; r with delay -0.034, p = 1.5e-03
onset at 20% of peak speed: reaction time mean 362 ms, reach time (a) mean 225 ms; r with delay -0.002, p = 8.5e-01
rt_mean reach_a reach_b r_delay p
date
2015-10-01 350.0 246.0 596.0 0.010 7.1e-01
2015-10-02 331.0 241.0 572.0 -0.031 1.8e-01
2015-10-15 337.0 271.0 608.0 0.003 9.3e-01
2016-01-27 354.0 277.0 631.0 -0.070 3.1e-03
2016-01-28 308.0 248.0 556.0 0.025 1.9e-01
R: delays of delayed trials 400 to 1000 ms; deciles [400, 454, 508, 564, 618, 677, 736, 795, 859, 924, 1000]
no-delay fraction of successful center-out trials 0.002 (26 of 12948)
reaction time mean 305 ms (median 361); reach time (a) acquisition - onset mean 181 ms (median 169); (b) acquisition - go cue mean 486 ms (median 519); trials 12906, without kinematics 0, without an onset 16, without acquisition 0
reaction time against delay: r = -0.264, p = 1.8e-205
variant, reaction times 100 to 800 ms only (11246 trials): mean 362 ms, r = -0.337, p = 1.2e-297
onset before the go cue in 0.059 of trials; variant, onset searched only after the go cue: reaction time mean 340 ms (median 367), reach time (b) unchanged, trials without such an onset 759 of 12906; r with delay -0.376, p = 0.0e+00
onset at 5% of peak speed: reaction time mean 272 ms, reach time (a) mean 221 ms; r with delay -0.252, p = 1.7e-184
onset at 20% of peak speed: reaction time mean 342 ms, reach time (a) mean 144 ms; r with delay -0.304, p = 1.6e-273
rt_mean reach_a reach_b r_delay p
date
2017-01-15 249.0 229.0 478.0 -0.167 6.0e-13
2017-01-16 271.0 196.0 467.0 -0.178 1.3e-15
2017-01-17 310.0 152.0 462.0 -0.199 3.8e-19
2017-01-18 306.0 143.0 449.0 -0.249 3.8e-32
2017-01-19 318.0 182.0 500.0 -0.473 1.1e-75
2017-01-24 353.0 192.0 545.0 -0.443 3.8e-78
2017-01-25 339.0 183.0 521.0 -0.530 8.6e-144
Result: partially reproduced. J's delays run from 300 to 700 ms, inside the stated range; R's run from 400 to 1000 ms, past the widened upper bound of 990 ms. Trials without a delay are 19.8% of J's successful center-out trials and 0.2% of R's (26 trials), not 10%. Mean reaction times are 332 ms (J) and 305 ms (R), within tolerance of 346 and 311 ms, and reach times read as acquisition minus go cue are 587 and 486 ms, within tolerance of 582 and 482 ms (read as acquisition minus movement onset they are 255 and 181 ms). R's mean reaction time includes the 5.9% of its trials whose onset precedes the go cue (J 1.4%); R's median is 361 ms and its mean over reaction times of 100 to 800 ms 362 ms, both outside 311 +- 31 ms. With onset searched only after the go cue the means are 342 ms (J) and 340 ms (R), the latter just inside the bound of 342 ms, with 124 (J) and 759 (R) trials left without an onset. Reaction time falls with delay in R (r = -0.264, p = 1.8e-205) but not in J (r = -0.013, p = 0.21); with onset at 5% of peak speed J gives r = -0.034, p = 1.5e-03, and with onset searched after the go cue r = -0.029, p = 7.9e-03, both above 1e-4.
channels-task-correlated¶
Of the 96 PMd electrodes, 74 for J and 80 for R had activity significantly correlated with the x and y location of the target (p < .05, Bonferroni over 96).
Result: reproduced. 72 of 96 channels for J and 80 of 96 for R regress significantly on target x and y, against the paper's 74 and 80. Within single sessions the counts are lower (J 39 to 70, R 30 to 70). On the units the file labels PMd the counts are 74 (J) and 40 (R).
pca-variance¶
The top six PCs of the condition-averaged delay activity explain cumulatively 58, 82, 88, 92, 94 and 95% of the variance for J (120 conditions) and 81, 89, 92, 94, 95 and 96% for R (132 conditions).
PAPER_CUM = {"J": [58, 82, 88, 92, 94, 95], "R": [81, 89, 92, 94, 95, 96]}
planes = {}
for m in "JR":
t, X = data[m]
planes[m] = premotor.fit_plane(t, X, premotor.N_PCS[m])
cum = 100 * np.cumsum(planes[m]["explained"])[:6]
within = np.abs(cum - PAPER_CUM[m]) <= 0.1 * np.array(PAPER_CUM[m])
print(f"{m}: {len(planes[m]['conditions'])} conditions; cumulative variance of PCs 1 to 6 (%): {np.round(cum, 1).tolist()}")
print(f" paper: {PAPER_CUM[m]}; within 10%: {within.tolist()}; outside by {np.round(np.maximum(np.abs(cum - PAPER_CUM[m]) - 0.1 * np.array(PAPER_CUM[m]), 0), 2).tolist()}")
variants = {
"day mean over the analysed delay windows": premotor.monkey_table(sessions, m, day_mean="trials"),
"units the file labels PMd": premotor.monkey_table(sessions, m, units=premotor.FILE_PMD_UNITS),
"no day-mean subtraction": (t, X + np.concatenate(
[s["units"].rate.values[premotor.PMD_UNITS][None, :].repeat((s["trials"].analysed).sum(), 0)
for s in sessions if s["monkey"] == m])),
}
for name, (tv, Xv) in variants.items():
pv = premotor.fit_plane(tv, Xv, premotor.N_PCS[m])
print(f" variant, {name}: {np.round(100 * np.cumsum(pv['explained'])[:6], 1).tolist()}")
J: 96 conditions; cumulative variance of PCs 1 to 6 (%): [59.8, 77.7, 84.0, 88.3, 90.9, 93.0]
paper: [58, 82, 88, 92, 94, 95]; within 10%: [True, True, True, True, True, True]; outside by [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
variant, day mean over the analysed delay windows: [62.7, 81.8, 86.7, 90.8, 93.0, 94.8] variant, units the file labels PMd: [28.4, 50.0, 64.2, 77.1, 84.7, 87.2] variant, no day-mean subtraction: [49.4, 68.8, 81.5, 90.8, 93.5, 95.0] R: 132 conditions; cumulative variance of PCs 1 to 6 (%): [72.8, 85.5, 90.1, 92.2, 93.7, 94.7] paper: [81, 89, 92, 94, 95, 96]; within 10%: [False, True, True, True, True, True]; outside by [0.14, 0.0, 0.0, 0.0, 0.0, 0.0]
variant, day mean over the analysed delay windows: [80.5, 87.1, 90.7, 92.2, 93.4, 94.2] variant, units the file labels PMd: [24.0, 36.6, 41.2, 45.0, 48.6, 51.9] variant, no day-mean subtraction: [80.6, 88.1, 92.0, 94.8, 96.4, 97.2]
Result: not reproduced. R has 132 conditions, as in the paper. For J (96 of its 120 conditions) all six cumulative values are within tolerance (59.8 to 93.0% against 58 to 95%). For R, PCs 2 to 6 are within tolerance (85.5 to 94.7% against 89 to 96%) and PC1 is 72.8% against 81 +- 8.1%, outside by 0.14 points. The criterion's partial case needs all of R's values within tolerance, so by its letter the result is not reproduced, though it misses by a small margin. The result depends on how the day's mean is taken: with the mean over the analysed delay windows R's PC1 is 80.5%, and without day-mean subtraction 80.6%. The main choice, the rate over the session's recorded time, is the rate over trial time, as the data check shows that spikes exist only within trials.
spatial-plane-variance¶
The spatial plane captures a large part of the variance of the condition-averaged activity for the ring, horizontal, vertical and 3-rings tasks; 66.8, 64.4, 47.7 and 68.2% for J and 33.7, 30.7, 31.0 and 30.1% for R (p < 1e-4).
PAPER_PLANE = {"J": {"ring": 66.8, "horizontal": 64.4, "vertical": 47.7, "3-rings": 68.2},
"R": {"ring": 33.7, "horizontal": 30.7, "vertical": 31.0, "3-rings": 30.1}}
PLANE_SEED = 3
def task_means(t, X, task_column):
"""Per task, the condition means centred within the task."""
out = {}
for task in sorted(t[task_column].unique()):
sel = (t[task_column] == task).values
M = pd.DataFrame(X[sel]).groupby(t.condition.values[sel]).mean().values
out[task] = M - M.mean(axis=0)
return out
def plane_shares(means, plane):
"""Per task, share of the variance of its condition means in the plane: (a) of all channels,
(b) of the top PCs."""
Q = premotor.plane_basis(plane)
return {task: {"a": 100 * ((Mc @ Q) ** 2).sum() / (Mc ** 2).sum(),
"b": 100 * ((Mc @ Q) ** 2).sum() / ((Mc @ plane["A"].T) ** 2).sum()}
for task, Mc in means.items()}
def refit_with_targets(X, plane, tx, ty):
"""The plane refit from the same PCs to other target positions."""
S = (X - plane["mu"]) @ plane["A"].T
ones = np.ones((len(S), 1))
bx = np.linalg.lstsq(np.hstack([ones, S]), tx, rcond=None)[0]
null = np.eye(S.shape[1]) - np.outer(bx[1:], bx[1:]) / (bx[1:] @ bx[1:])
by = np.linalg.lstsq(np.hstack([ones, S @ null]), ty, rcond=None)[0]
return {**plane, "wx": plane["A"].T @ bx[1:], "wy": plane["A"].T @ (null @ by[1:])}
plane_results = {}
for m in "JR":
t, X = data[m]
codes, names = pd.factorize(t.condition)
targets = t.groupby("condition")[["tx", "ty"]].first().loc[names].values
rng = np.random.default_rng(PLANE_SEED)
shuffled_planes = []
for _ in range(1000):
txy = targets[rng.permutation(len(targets))][codes]
shuffled_planes.append(refit_with_targets(X, planes[m], txy[:, 0], txy[:, 1]))
for task_column in ("methods_config", "config"):
means = task_means(t, X, task_column)
shares = plane_shares(means, planes[m])
null = pd.DataFrame([{k: v["a"] for k, v in plane_shares(means, sp).items()} for sp in shuffled_planes])
rows = []
for task, v in shares.items():
nul = null[task].values
rows.append({"task": task, "(a) of all channels %": round(v["a"], 1), "(b) of top PCs %": round(v["b"], 1),
"paper %": PAPER_PLANE[m].get(task), "shuffle null (a) median %": round(np.median(nul), 1),
"shuffle null (a) 99th pct %": round(np.quantile(nul, 0.99), 1),
"p": (1 + np.sum(nul >= v["a"])) / (1 + len(nul))})
res = pd.DataFrame(rows)
plane_results[(m, task_column)] = res
print(f"{m}, tasks by the {'Methods dates' if task_column == 'methods_config' else 'file geometry'}:")
print(res.to_string(index=False))
print(f" smallest p attainable with {len(shuffled_planes)} shuffles: {1 / (1 + len(shuffled_planes)):.1e}")
tv, Xv = premotor.monkey_table(sessions, m, day_mean="trials")
pv = premotor.fit_plane(tv, Xv, premotor.N_PCS[m])
print(" variant, day mean over the analysed delay windows, (a):",
{k: float(round(v["a"], 1)) for k, v in plane_shares(task_means(tv, Xv, "methods_config"), pv).items()})
J, tasks by the Methods dates:
task (a) of all channels % (b) of top PCs % paper % shuffle null (a) median % shuffle null (a) 99th pct % p
3-rings 44.3 46.9 68.2 20.7 54.6 0.045954
ring 54.4 58.7 66.8 24.8 51.7 0.004995
vertical 40.2 44.6 47.7 11.5 53.5 0.047952
J, tasks by the file geometry:
task (a) of all channels % (b) of top PCs % paper % shuffle null (a) median % shuffle null (a) 99th pct % p
3-rings 44.3 46.9 68.2 20.7 54.6 0.045954
ring 54.4 58.7 66.8 24.8 51.7 0.004995
vertical 40.2 44.6 47.7 11.5 53.5 0.047952
smallest p attainable with 1000 shuffles: 1.0e-03
variant, day mean over the analysed delay windows, (a): {'3-rings': 46.5, 'ring': 52.8, 'vertical': 42.9}
R, tasks by the Methods dates:
task (a) of all channels % (b) of top PCs % paper % shuffle null (a) median % shuffle null (a) 99th pct % p
3-rings 32.4 35.7 30.1 14.8 46.9 0.056943
horizontal 36.9 39.6 30.7 22.0 55.3 0.146853
ring 29.9 31.6 33.7 14.7 48.7 0.086913
vertical 29.5 32.5 31.0 14.9 52.1 0.161838
R, tasks by the file geometry:
task (a) of all channels % (b) of top PCs % paper % shuffle null (a) median % shuffle null (a) 99th pct % p
3-rings 32.4 35.7 30.1 14.8 46.9 0.056943
horizontal 29.5 32.5 30.7 14.9 52.1 0.161838
ring 29.9 31.6 33.7 14.7 48.7 0.086913
vertical 36.9 39.6 31.0 22.0 55.3 0.146853
smallest p attainable with 1000 shuffles: 1.0e-03
variant, day mean over the analysed delay windows, (a): {'3-rings': 24.6, 'horizontal': 34.5, 'ring': 22.9, 'vertical': 35.4}
Result: not reproduced. No pairing of readings puts more than three of the seven values within tolerance. For R the shares lie between 29.5 and 36.9% (reading a) against the paper's 30.1 to 33.7%, with ring at 29.9% against 33.7 +- 3.37%; for J (5 of 16 sessions, plane fit on 96 conditions) they are 54.4% (ring), 40.2% (vertical) and 44.3% (3-rings) against 66.8, 47.7 and 68.2%. Planes fit to targets shuffled across conditions capture median shares of 11.5 to 24.8%, and the observed shares exceed the shuffles with p from 0.005 to 0.16. The smallest p that 1000 shuffles can give is 1.0e-03, but the observed p values lie well above it, so the paper's p < 1e-4 does not hold against this null; the paper does not say what its null was.
plane-timecourse¶
After target onset the neural state in the spatial plane starts at the center, diverges after 50 ms (R 75 ms) and converges to the steady delay state after 300 ms (R 400 ms).
centers = premotor.TC_START_S + premotor.BIN_S * (np.arange(45) + 0.5)
bin_end = centers + premotor.BIN_S / 2
timecourse = {}
for m in "JR":
t, X = data[m]
plane = planes[m]
coords, valid = [], []
for s in sessions:
if s["monkey"] != m:
continue
a = s["trials"][s["trials"].analysed]
rates = s["tc_counts"] / premotor.BIN_S - s["units"].rate.values[premotor.PMD_UNITS][None, :, None]
kernel = np.ones(3) / 3
rates = np.apply_along_axis(lambda v: np.convolve(v, kernel, mode="same"), 2, rates)
xy = premotor.plane_coords(plane, np.transpose(rates, (0, 2, 1))) # trials x bins x 2
coords.append(xy)
# a bin is used while the go cue comes after it, including the bin the smoothing reaches
valid.append(a.delay.values[:, None] >= bin_end[None, :] + premotor.BIN_S)
coords, valid = np.concatenate(coords), np.concatenate(valid)
conv = pd.DataFrame(premotor.plane_coords(plane, X)).groupby(t.condition.values).mean()
d = np.full(len(centers), np.nan)
n_cond = np.zeros(len(centers), int)
for b in range(len(centers)):
per = []
for c, idx in t.groupby("condition").indices.items():
ok = idx[valid[idx, b]]
if len(ok) >= 5:
per.append(np.linalg.norm(coords[ok, b].mean(axis=0) - conv.loc[c].values))
if per:
d[b], n_cond[b] = np.mean(per), len(per)
d0 = np.nanmean(d[(centers > -0.1) & (centers < 0)])
after = centers > 0
dmin = np.nanmin(d[after])
D = (d - dmin) / (d0 - dmin)
diverge = centers[after][np.argmax(D[after] < 0.9)]
converge = centers[after][np.argmax(D[after] < 0.1)]
timecourse[m] = (d, D, n_cond)
print(f"{m}: d0 {d0:.1f} mm, dmin {dmin:.1f} mm; divergence (D < 0.9) at {1000 * diverge:.0f} ms, "
f"convergence (D < 0.1) at {1000 * converge:.0f} ms after target onset")
for lo, hi in ((0.8, 0.1), (0.95, 0.2), (0.9, 0.05)):
print(f" thresholds {lo} and {hi}: {1000 * centers[after][np.argmax(D[after] < lo)]:.0f} and "
f"{1000 * centers[after][np.argmax(D[after] < hi)]:.0f} ms")
print(" D(t) every 40 ms from -90 ms:", {int(round(1000 * c)): float(round(v, 2)) for c, v in zip(centers[5::2], D[5::2])})
print(f" conditions with at least 5 trials per bin: {n_cond.min()} to {n_cond.max()}")
pre = centers < 0
print(f" D over the bins before target onset (-200 to 0 ms): mean {np.nanmean(D[pre]):.2f}, SD {np.nanstd(D[pre]):.2f}, "
f"range {np.nanmin(D[pre]):.2f} to {np.nanmax(D[pre]):.2f}; pre-onset bins below 0.9: {int(np.sum(D[pre] < 0.9))} of {pre.sum()}")
# variant: at each bin, only trials whose 200 ms delay window starts after the bin (and its smoothing),
# with the converged state taken from those same trials
delays = np.concatenate([s["trials"][s["trials"].analysed].delay.values for s in sessions if s["monkey"] == m])
P_delay = premotor.plane_coords(plane, X)
d2 = np.full(len(centers), np.nan)
for b in range(len(centers)):
late = delays - premotor.DELAY_WINDOW_S >= bin_end[b] + premotor.BIN_S
per = []
for c, idx in t.groupby("condition").indices.items():
ok = idx[late[idx]]
if len(ok) >= 5:
per.append(np.linalg.norm(coords[ok, b].mean(axis=0) - P_delay[ok].mean(axis=0)))
if len(per) >= 0.5 * len(conv):
d2[b] = np.mean(per)
d20 = np.nanmean(d2[(centers > -0.1) & (centers < 0)])
d2min = np.nanmin(d2[after])
D2 = (d2 - d2min) / (d20 - d2min)
last = centers[np.isfinite(d2)].max()
hit = lambda thr: (f"{1000 * centers[after][np.nanargmax(D2[after] < thr)]:.0f} ms" if np.any(D2[after] < thr) else "never")
print(f" variant, converged state from trials whose delay window starts after the bin (bins up to {1000 * last:.0f} ms "
f"with at least half the conditions): divergence {hit(0.9)}, convergence {hit(0.1)}; minimum at "
f"{1000 * centers[after][np.nanargmin(D2[after])]:.0f} ms; pre-onset D SD {np.nanstd(D2[pre]):.2f}")
origin = coords[:, centers < 0][valid[:, centers < 0]]
print(f" state before target onset, mean over trials (mm): x {origin[:, 0].mean():.1f}, y {origin[:, 1].mean():.1f}")
fig, ax = plt.subplots(1, 2, figsize=(9, 3))
for i, m in enumerate("JR"):
ax[i].plot(1000 * centers, timecourse[m][0], color="tab:blue")
ax[i].axvline(0, color="grey", lw=0.5)
ax[i].set(title=f"monkey {m}", xlabel="time from target onset (ms)", ylabel="distance from delay state (mm)")
plt.tight_layout()
plt.show()
print("Plotted ranges (mm):", {m: (float(round(np.nanmin(timecourse[m][0]), 1)), float(round(np.nanmax(timecourse[m][0]), 1))) for m in "JR"})
J: d0 53.8 mm, dmin 6.8 mm; divergence (D < 0.9) at 70 ms, convergence (D < 0.1) at 270 ms after target onset
thresholds 0.8 and 0.1: 150 and 270 ms
thresholds 0.95 and 0.2: 70 and 230 ms
thresholds 0.9 and 0.05: 70 and 390 ms
D(t) every 40 ms from -90 ms: {-90: 1.03, -50: 1.02, -10: 0.96, 30: 1.0, 70: 0.89, 110: 0.85, 150: 0.68, 190: 0.4, 230: 0.19, 270: 0.08, 310: 0.12, 350: 0.1, 390: 0.03, 430: 0.0, 470: 0.04, 510: 0.12, 550: 0.17, 590: 0.2, 630: 0.24, 670: nan}
conditions with at least 5 trials per bin: 0 to 96
D over the bins before target onset (-200 to 0 ms): mean 1.01, SD 0.03, range 0.96 to 1.03; pre-onset bins below 0.9: 0 of 10
variant, converged state from trials whose delay window starts after the bin (bins up to 430 ms with at least half the conditions): divergence 70 ms, convergence 250 ms; minimum at 410 ms; pre-onset D SD 0.03
state before target onset, mean over trials (mm): x -2.1, y -30.4
R: d0 47.3 mm, dmin 10.2 mm; divergence (D < 0.9) at 90 ms, convergence (D < 0.1) at 350 ms after target onset
thresholds 0.8 and 0.1: 130 and 350 ms
thresholds 0.95 and 0.2: 10 and 310 ms
thresholds 0.9 and 0.05: 90 and 390 ms
D(t) every 40 ms from -90 ms: {-90: 0.87, -50: 1.1, -10: 0.97, 30: 0.9, 70: 0.91, 110: 0.85, 150: 0.74, 190: 0.63, 230: 0.45, 270: 0.29, 310: 0.16, 350: 0.08, 390: 0.03, 430: 0.04, 470: 0.03, 510: 0.02, 550: 0.04, 590: 0.02, 630: 0.01, 670: 0.01}
conditions with at least 5 trials per bin: 132 to 132
D over the bins before target onset (-200 to 0 ms): mean 0.94, SD 0.09, range 0.87 to 1.11; pre-onset bins below 0.9: 6 of 10
variant, converged state from trials whose delay window starts after the bin (bins up to 690 ms with at least half the conditions): divergence 10 ms, convergence 310 ms; minimum at 650 ms; pre-onset D SD 0.11
state before target onset, mean over trials (mm): x 17.6, y -13.1
Plotted ranges (mm): {'J': (6.8, 55.0), 'R': (10.2, 51.4)}
Result: partially reproduced. The distance from the delay state starts to fall (D below 0.9) 70 ms after target onset for J and 90 ms for R, against 50 and 75 ms, and reaches the delay state (D below 0.1) at 270 ms for J and 350 ms for R, against 300 and 400 ms. All four times are within tolerance, so by the letter of the criterion the result is reproduced; the outcome of partial reproduction is a judgment beyond the criterion, made because R's divergence cannot be told from baseline noise. J's times hold under the two checks added here: its D varies little before target onset (SD 0.03, range 0.96 to 1.03), and with the converged state taken only from trials whose delay window starts after the bin its times are 70 and 250 ms. R's D lies below 0.9 in 6 of the 10 bins before target onset (range 0.87 to 1.11, SD 0.09), and with the non-overlapping converged state R diverges at 10 ms and converges at 310 ms, both outside tolerance (75 +- 50 and 400 +- 50 ms). Both monkeys' times depend on the threshold: with other thresholds the divergence moves between 10 and 150 ms and the convergence between 230 and 390 ms, and J diverges at 150 ms with threshold 0.8 and converges at 390 ms with threshold 0.05, both outside tolerance. Before target onset the state sits at x -2.1, y -30.4 mm (J) and 17.6, -13.1 mm (R) of the plane's target coordinates, so the statement that it starts at the center is not borne out.
onoff-predicts-kinematics¶
On single trials, on-axis and off-axis deviations of the neural state in the spatial plane are positively correlated with the on-axis and off-axis deviations of the hand position at maximum speed (on-axis p = 1e-3 J, 4e-3 R; off-axis p = 4e-4 J, 2e-4 R).
devs = {}
for m in "JR":
t, X = data[m]
P = premotor.plane_coords(planes[m], X)
ok = t.max_speed.notna().values
tt = t[ok].reset_index(drop=True)
nd = premotor.deviations(P[ok], tt.condition.values)
kd = premotor.deviations(tt[["x_at_max", "y_at_max"]].values, tt.condition.values)
n_on, n_off = premotor.on_off(nd, tt.tx.values, tt.ty.values)
k_on, k_off = premotor.on_off(kd, tt.tx.values, tt.ty.values)
devs[m] = (tt, P[ok], n_on, n_off)
print(f"{m}: {len(tt)} trials ({int((~ok).sum())} without kinematics)")
for name, sel in (("(a) all tasks", np.ones(len(tt), bool)), ("(b) 3-rings", (tt.config == "3-rings").values)):
r1, p1 = stats.pearsonr(n_on[sel], k_on[sel])
r2, p2 = stats.pearsonr(n_off[sel], k_off[sel])
print(f" {name}: on-axis r = {r1:.3f}, p = {p1:.1e}; off-axis r = {r2:.3f}, p = {p2:.1e} (n = {sel.sum()})")
for d in sorted(tt.date.unique()):
sel = (tt.date == d).values
r1, p1 = stats.pearsonr(n_on[sel], k_on[sel])
r2, p2 = stats.pearsonr(n_off[sel], k_off[sel])
print(f" {d}: on r = {r1:.3f} (p = {p1:.1e}), off r = {r2:.3f} (p = {p2:.1e})")
sc = (tt.date + "|" + tt.condition).values
nd2 = premotor.deviations(P[ok], sc)
kd2 = premotor.deviations(tt[["x_at_max", "y_at_max"]].values, sc)
n_on2, n_off2 = premotor.on_off(nd2, tt.tx.values, tt.ty.values)
k_on2, k_off2 = premotor.on_off(kd2, tt.tx.values, tt.ty.values)
for name, sel in (("(a) all tasks", np.ones(len(tt), bool)), ("(b) 3-rings", (tt.config == "3-rings").values)):
r1, p1 = stats.pearsonr(n_on2[sel], k_on2[sel])
r2, p2 = stats.pearsonr(n_off2[sel], k_off2[sel])
print(f" variant, deviations from session x condition means, {name}: on-axis r = {r1:.3f}, p = {p1:.1e}; "
f"off-axis r = {r2:.3f}, p = {p2:.1e}")
rs1, ps1 = stats.spearmanr(n_on, k_on)
rs2, ps2 = stats.spearmanr(n_off, k_off)
print(f" variant, Spearman over all tasks: on {rs1:.3f} (p = {ps1:.1e}), off {rs2:.3f} (p = {ps2:.1e})")
J: 8701 trials (54 without kinematics) (a) all tasks: on-axis r = 0.091, p = 1.7e-17; off-axis r = 0.102, p = 1.3e-21 (n = 8701) (b) 3-rings: on-axis r = 0.101, p = 1.0e-11; off-axis r = 0.131, p = 1.4e-18 (n = 4471) 2015-10-01: on r = 0.094 (p = 4.7e-04), off r = 0.050 (p = 6.7e-02) 2015-10-02: on r = 0.070 (p = 2.7e-03), off r = 0.012 (p = 6.1e-01) 2015-10-15: on r = 0.102 (p = 1.1e-03), off r = 0.133 (p = 1.9e-05) 2016-01-27: on r = 0.070 (p = 3.2e-03), off r = 0.144 (p = 1.0e-09) 2016-01-28: on r = 0.125 (p = 6.1e-11), off r = 0.126 (p = 4.4e-11) variant, deviations from session x condition means, (a) all tasks: on-axis r = 0.092, p = 9.7e-18; off-axis r = 0.101, p = 3.1e-21 variant, deviations from session x condition means, (b) 3-rings: on-axis r = 0.100, p = 2.3e-11; off-axis r = 0.133, p = 4.2e-19 variant, Spearman over all tasks: on 0.078 (p = 4.3e-13), off 0.098 (p = 7.1e-20) R: 12853 trials (0 without kinematics) (a) all tasks: on-axis r = 0.022, p = 1.3e-02; off-axis r = 0.072, p = 2.2e-16 (n = 12853) (b) 3-rings: on-axis r = 0.051, p = 6.5e-02; off-axis r = 0.119, p = 1.2e-05 (n = 1340) 2017-01-15: on r = 0.050 (p = 3.4e-02), off r = 0.097 (p = 3.1e-05) 2017-01-16: on r = -0.044 (p = 4.9e-02), off r = 0.031 (p = 1.6e-01) 2017-01-17: on r = -0.003 (p = 8.8e-01), off r = 0.078 (p = 6.5e-04) 2017-01-18: on r = 0.029 (p = 1.8e-01), off r = 0.064 (p = 2.7e-03) 2017-01-19: on r = 0.051 (p = 6.5e-02), off r = 0.119 (p = 1.2e-05) 2017-01-24: on r = 0.012 (p = 6.3e-01), off r = 0.037 (p = 1.4e-01) 2017-01-25: on r = 0.028 (p = 2.2e-01), off r = 0.125 (p = 2.6e-08) variant, deviations from session x condition means, (a) all tasks: on-axis r = 0.020, p = 2.3e-02; off-axis r = 0.083, p = 3.0e-21 variant, deviations from session x condition means, (b) 3-rings: on-axis r = 0.051, p = 6.5e-02; off-axis r = 0.119, p = 1.2e-05 variant, Spearman over all tasks: on 0.033 (p = 1.5e-04), off 0.058 (p = 4.1e-11)
Result: reproduced. Over all tasks, neural and kinematic deviations correlate positively on axis (J r = 0.091, p = 1.7e-17; R r = 0.022, p = 1.3e-02) and off axis (J r = 0.102, p = 1.3e-21; R r = 0.072, p = 2.2e-16). With deviations from session by condition means, so that offsets between sessions cannot correlate, the values are nearly the same (on axis J r = 0.092, p = 9.7e-18, R r = 0.020, p = 2.3e-02; off axis J 0.101, R 0.083, p at most 3.1e-21). R's on-axis correlation is weak: in the 3-rings task alone it is not significant (r = 0.051, p = 6.5e-02), and within sessions it is positive with p < 0.05 only on 2017-01-15 and negative on 2017-01-16 (r = -0.044, p = 4.9e-02).
speed-axis-within-distance¶
The speed axis of delay activity varies with maximum speed for reaches to the same target distance, for 11 of 12 distances in J and 11 of 12 in R (p < 0.05), so it encodes speed independently of distance.
SPEED_SEED = 4
def speed_axis(Xc, speed, Q):
"""Regression of maximum speed on the activity with the plane projected out (unit-free axis)."""
Xp = Xc - (Xc @ Q) @ Q.T
fit = LinearRegression().fit(Xp, speed)
return fit, Xp
speed_fits = {}
for m in "JR":
t, X = data[m]
ok = t.max_speed.notna().values
tt, Xo = t[ok].reset_index(drop=True), X[ok]
Xc = Xo - Xo.mean(axis=0)
Q = premotor.plane_basis(planes[m])
lines = (tt.config.isin(["horizontal", "vertical"])).values
# the k-th distance of each line task (R's horizontal targets sit at 5.2 + 10.4 k mm)
tt["distance_rank"] = tt.groupby("config").distance.rank(method="dense").astype(int)
print(f"{m}: {len(tt)} trials with kinematics, {lines.sum()} horizontal or vertical "
f"({sorted(tt[lines].config.unique())}); distances (mm) by task "
f"{ {c: sorted(np.round(g.distance, 1).unique().tolist()) for c, g in tt[lines].groupby('config')} }")
for reading, fit_rows in (("(a) fit on all tasks", np.ones(len(tt), bool)), ("(b) fit on horizontal and vertical", lines)):
fit, Xp = speed_axis(Xc[fit_rows], tt.max_speed.values[fit_rows], Q)
proj = (Xc - (Xc @ Q) @ Q.T) @ fit.coef_
if reading.startswith("(a)"):
speed_fits[m] = (tt, Xc, fit.coef_)
# cross-validated projections (variant)
oof = np.full(len(tt), np.nan)
idx = np.flatnonzero(fit_rows)
for train, test in KFold(10, shuffle=True, random_state=SPEED_SEED).split(idx):
f2, _ = speed_axis(Xc[idx[train]], tt.max_speed.values[idx[train]], Q)
oof[idx[test]] = (Xc[idx[test]] - (Xc[idx[test]] @ Q) @ Q.T) @ f2.coef_
out = []
for dist in range(1, 13):
sel = lines & (tt.distance_rank == dist).values
r, p = stats.pearsonr(proj[sel], tt.max_speed[sel])
rc, pc = stats.pearsonr(oof[sel], tt.max_speed[sel]) if np.isfinite(oof[sel]).all() else (np.nan, np.nan)
out.append({"distance_rank": dist, "n": int(sel.sum()), "r": round(r, 3), "p": p,
"r_cv": round(rc, 3), "p_cv": pc})
out = pd.DataFrame(out)
n_sig = int(((out.r > 0) & (out.p < 0.05)).sum())
n_sig_cv = int(((out.r_cv > 0) & (out.p_cv < 0.05)).sum())
print(f" {reading}: significant positive at {n_sig} of {len(out)} distances; cross-validated {n_sig_cv} of {len(out)}")
print(out.to_string(index=False, formatters={"p": "{:.1e}".format, "p_cv": "{:.1e}".format}))
proj_a = (Xc - (Xc @ Q) @ Q.T) @ speed_fits[m][2]
split = []
for (task, date), g in tt[lines].groupby(["config", "date"]):
n_pos, ns = 0, []
for dist in range(1, 13):
idx = g.index[g.distance_rank == dist]
ns.append(len(idx))
r, p = stats.pearsonr(proj_a[idx], tt.max_speed[idx])
n_pos += (r > 0) & (p < 0.05)
split.append(f"{task} {date}: {n_pos} of 12 ({min(ns)} to {max(ns)} trials per distance)")
print(" reading (a), significant positive distances within each task and session:", "; ".join(split))
r_all, p_all = stats.pearsonr(tt.max_speed, tt.distance)
print(f" maximum speed against target distance over all trials: r = {r_all:.3f}, p = {p_all:.1e}")
J: 8701 trials with kinematics, 3211 horizontal or vertical (['vertical']); distances (mm) by task {'vertical': [10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0, 100.0, 110.0, 120.0]}
(a) fit on all tasks: significant positive at 6 of 12 distances; cross-validated 4 of 12
distance_rank n r p r_cv p_cv
1 267 0.039 5.3e-01 0.027 6.7e-01
2 264 0.129 3.6e-02 0.111 7.1e-02
3 270 0.191 1.7e-03 0.177 3.6e-03
4 274 0.178 3.2e-03 0.172 4.4e-03
5 256 0.187 2.7e-03 0.170 6.4e-03
6 266 0.132 3.1e-02 0.116 5.9e-02
7 277 0.084 1.6e-01 0.064 2.9e-01
8 269 0.030 6.2e-01 0.019 7.6e-01
9 265 0.216 3.9e-04 0.207 7.1e-04
10 270 -0.028 6.4e-01 -0.043 4.8e-01
11 266 0.087 1.6e-01 0.072 2.4e-01
12 267 0.053 3.9e-01 0.048 4.3e-01
(b) fit on horizontal and vertical: significant positive at 6 of 12 distances; cross-validated 5 of 12
distance_rank n r p r_cv p_cv
1 267 0.048 4.3e-01 0.013 8.3e-01
2 264 0.138 2.5e-02 0.109 7.6e-02
3 270 0.212 4.5e-04 0.160 8.4e-03
4 274 0.174 3.9e-03 0.128 3.4e-02
5 256 0.192 2.0e-03 0.160 1.0e-02
6 266 0.160 8.8e-03 0.121 4.9e-02
7 277 0.116 5.4e-02 0.074 2.2e-01
8 269 0.055 3.7e-01 0.029 6.4e-01
9 265 0.211 5.3e-04 0.176 4.1e-03
10 270 -0.019 7.6e-01 -0.047 4.4e-01
11 266 0.101 1.0e-01 0.043 4.9e-01
12 267 0.074 2.3e-01 0.046 4.6e-01
reading (a), significant positive distances within each task and session: vertical 2015-10-01: 1 of 12 (106 to 124 trials per distance); vertical 2015-10-02: 3 of 12 (143 to 164 trials per distance)
maximum speed against target distance over all trials: r = 0.813, p = 0.0e+00
R: 12853 trials with kinematics, 7929 horizontal or vertical (['horizontal', 'vertical']); distances (mm) by task {'horizontal': [5.2, 15.6, 26.1, 36.5, 47.0, 57.4, 67.8, 78.3, 88.7, 99.1, 109.6, 120.0], 'vertical': [10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0, 100.0, 110.0, 120.0]}
(a) fit on all tasks: significant positive at 11 of 12 distances; cross-validated 11 of 12
distance_rank n r p r_cv p_cv
1 658 -0.040 3.1e-01 -0.043 2.7e-01
2 659 0.159 4.0e-05 0.153 8.4e-05
3 659 0.280 2.5e-13 0.270 1.7e-12
4 660 0.358 2.1e-21 0.351 1.3e-20
5 657 0.346 5.7e-20 0.333 1.6e-18
6 659 0.384 1.4e-24 0.372 5.2e-23
7 659 0.313 1.8e-16 0.303 1.9e-15
8 682 0.258 7.6e-12 0.248 5.5e-11
9 659 0.261 9.9e-12 0.251 5.9e-11
10 659 0.227 3.9e-09 0.217 1.7e-08
11 659 0.259 1.5e-11 0.247 1.3e-10
12 659 0.197 3.5e-07 0.186 1.6e-06
(b) fit on horizontal and vertical: significant positive at 11 of 12 distances; cross-validated 11 of 12
distance_rank n r p r_cv p_cv
1 658 -0.031 4.3e-01 -0.028 4.8e-01
2 659 0.182 2.4e-06 0.172 8.7e-06
3 659 0.281 1.9e-13 0.266 4.0e-12
4 660 0.338 4.1e-19 0.325 1.1e-17
5 657 0.302 2.4e-15 0.283 1.5e-13
6 659 0.343 1.2e-19 0.328 5.1e-18
7 659 0.288 4.4e-14 0.274 8.1e-13
8 682 0.237 3.9e-10 0.214 1.7e-08
9 659 0.234 1.1e-09 0.214 2.9e-08
10 659 0.201 1.9e-07 0.184 1.9e-06
11 659 0.209 6.1e-08 0.187 1.3e-06
12 659 0.144 2.1e-04 0.126 1.2e-03
reading (a), significant positive distances within each task and session: horizontal 2017-01-15: 1 of 12 (152 to 154 trials per distance); horizontal 2017-01-16: 4 of 12 (166 to 168 trials per distance); vertical 2017-01-17: 10 of 12 (156 to 179 trials per distance); vertical 2017-01-18: 7 of 12 (181 to 182 trials per distance)
maximum speed against target distance over all trials: r = 0.764, p = 0.0e+00
Result: partially reproduced. For R the speed axis correlates positively with maximum speed at 11 of 12 distances under both readings, as in the paper, and also with cross-validated projections; the exception is the nearest distance. For J, whose horizontal sessions are not in the dandiset, it does at 6 of 12 distances under both readings (4 and 5 cross-validated), on 3211 vertical trials. Each distance rank pools both line tasks and both sessions of each; R's horizontal (5.2 + 10.4 k mm) and vertical (10 k mm) targets of the same rank differ by up to 5 mm, and rank 12 (both at 120 mm) is significant. Split by task and session, with 152 to 182 trials per distance (R) and 106 to 164 (J), R has 1, 4, 10 and 7 of 12 (horizontal 2017-01-15 and 16, vertical 01-17 and 18) and J 1 and 3, so the counts rest on pooling sessions.
distance-after-speed¶
After regressing out speed from the neural activity, a positive correlation with target distance remains (p < 1e-4 for J and R).
DIST_SEED = 5
distance_axes = {}
for m in "JR":
tt, Xc, ws = speed_fits[m]
s = tt.max_speed.values - tt.max_speed.mean()
Xr = Xc - np.outer(s, Xc.T @ s / (s @ s))
y = tt.distance.values
oof = np.empty(len(y))
for train, test in KFold(10, shuffle=True, random_state=DIST_SEED).split(Xr):
oof[test] = LinearRegression().fit(Xr[train], y[train]).predict(Xr[test])
r, p = stats.pearsonr(oof, y)
full = LinearRegression().fit(Xr, y)
ri, pi = stats.pearsonr(full.predict(Xr), y)
distance_axes[m] = full.coef_
print(f"{m}: cross-validated prediction of target distance from speed-residualised activity: r = {r:.3f}, p = {p:.1e} "
f"(n = {len(y)}); in-sample r = {ri:.3f}, p = {pi:.1e}")
print(f" residual correlation of every channel with speed (max |r|): "
f"{np.nanmax(np.abs([stats.pearsonr(Xr[:, c], s)[0] if Xr[:, c].std() > 0 else 0 for c in range(Xr.shape[1])])):.1e}")
for c in sorted(tt.config.unique()):
sel = (tt.config == c).values
if np.ptp(np.round(y[sel])) == 0:
continue # one distance (ring): no correlation to test
print(f" {c}: cross-validated r = {stats.pearsonr(oof[sel], y[sel])[0]:.3f}, p = {stats.pearsonr(oof[sel], y[sel])[1]:.1e}")
J: cross-validated prediction of target distance from speed-residualised activity: r = 0.123, p = 1.2e-30 (n = 8701); in-sample r = 0.182, p = 4.9e-66 residual correlation of every channel with speed (max |r|): 4.5e-16 3-rings: cross-validated r = 0.154, p = 5.1e-25 vertical: cross-validated r = 0.049, p = 5.6e-03
R: cross-validated prediction of target distance from speed-residualised activity: r = 0.233, p = 7.1e-158 (n = 12853); in-sample r = 0.257, p = 9.5e-193 residual correlation of every channel with speed (max |r|): 3.1e-16 3-rings: cross-validated r = 0.265, p = 5.2e-23 horizontal: cross-validated r = 0.355, p = 1.6e-114 vertical: cross-validated r = 0.214, p = 9.7e-44
Result: reproduced. With speed regressed out, cross-validated predictions of target distance correlate with distance (J r = 0.123, p = 1.2e-30; R r = 0.233, p = 7.1e-158). Within the tasks whose targets lie at several distances the correlation holds too: 3-rings J r = 0.154, R r = 0.265; vertical J 0.049 (p = 5.6e-03), R 0.214; horizontal R 0.355. The ring task, with one distance, is not tested.
speed-distance-projection¶
Projected onto the top PCs, the distance dimension is about half as large as the speed dimension.
for m in "JR":
A = planes[m]["A"]
_, _, ws = speed_fits[m]
wd = distance_axes[m]
ps, pd_ = np.linalg.norm(A @ (ws / np.linalg.norm(ws))), np.linalg.norm(A @ (wd / np.linalg.norm(wd)))
print(f"{m}: norm of the projection onto the top {A.shape[0]} PCs, speed axis {ps:.3f}, distance axis {pd_:.3f}; "
f"ratio {pd_ / ps:.2f}")
Q = premotor.plane_basis(planes[m])
wd_perp = wd - Q @ (Q.T @ wd)
print(f" variant, distance axis with the plane projected out: projection {np.linalg.norm(A @ (wd_perp / np.linalg.norm(wd_perp))):.3f}, "
f"ratio {np.linalg.norm(A @ (wd_perp / np.linalg.norm(wd_perp))) / ps:.2f}")
print(f" cosine between speed and distance axes: {ws @ wd / np.linalg.norm(ws) / np.linalg.norm(wd):.3f}")
J: norm of the projection onto the top 6 PCs, speed axis 0.324, distance axis 0.326; ratio 1.01 variant, distance axis with the plane projected out: projection 0.235, ratio 0.73 cosine between speed and distance axes: 0.124 R: norm of the projection onto the top 5 PCs, speed axis 0.314, distance axis 0.526; ratio 1.68 variant, distance axis with the plane projected out: projection 0.359, ratio 1.14 cosine between speed and distance axes: 0.427
Result: not reproduced. Projected onto the top PCs, the distance axis is as large as the speed axis for J (ratio 1.01) and larger for R (1.68), not about half. With the plane projected out of the distance axis the ratios are 0.73 and 1.14.
single-trial-variance¶
Across all single trials the speed dimension explains 2% (J) and 1% (R) of the variance of the neural state, and the spatial plane 20% (J) and 10% (R).
for m in "JR":
t, X = data[m]
Xc = X - X.mean(axis=0)
total = (Xc ** 2).sum()
_, _, ws = speed_fits[m]
u = ws / np.linalg.norm(ws)
Q = premotor.plane_basis(planes[m])
print(f"{m}: speed axis {100 * ((Xc @ u) ** 2).sum() / total:.2f}% of the single-trial variance; spatial plane "
f"{100 * ((Xc @ Q) ** 2).sum() / total:.1f}%; top {premotor.N_PCS[m]} PCs "
f"{100 * ((Xc @ planes[m]['A'].T) ** 2).sum() / total:.1f}% (trials {len(Xc)})")
tv, Xv = premotor.monkey_table(sessions, m, day_mean="trials")
Xvc = Xv - Xv.mean(axis=0)
pv = premotor.fit_plane(tv, Xv, premotor.N_PCS[m])
Qv = premotor.plane_basis(pv)
print(f" variant, day mean over the analysed delay windows: plane {100 * ((Xvc @ Qv) ** 2).sum() / (Xvc ** 2).sum():.1f}%, "
f"speed axis (same direction) {100 * ((Xvc @ u) ** 2).sum() / (Xvc ** 2).sum():.2f}%")
per = []
for d in sorted(t.date.unique()):
sel = (t.date == d).values
Xs = X[sel] - X[sel].mean(axis=0)
per.append(f"{d} {100 * ((Xs @ Q) ** 2).sum() / (Xs ** 2).sum():.1f}%")
print(" plane share within each session:", "; ".join(per))
J: speed axis 1.25% of the single-trial variance; spatial plane 15.7%; top 6 PCs 45.8% (trials 8755) variant, day mean over the analysed delay windows: plane 16.4%, speed axis (same direction) 1.17% plane share within each session: 2015-10-01 12.3%; 2015-10-02 11.2%; 2015-10-15 20.7%; 2016-01-27 14.8%; 2016-01-28 17.8% R: speed axis 0.86% of the single-trial variance; spatial plane 8.0%; top 5 PCs 63.9% (trials 12853)
variant, day mean over the analysed delay windows: plane 8.2%, speed axis (same direction) 0.76% plane share within each session: 2017-01-15 12.9%; 2017-01-16 6.5%; 2017-01-17 14.2%; 2017-01-18 8.7%; 2017-01-19 10.3%; 2017-01-24 10.6%; 2017-01-25 8.1%
Result: not reproduced. The speed axis explains 1.25% of the single-trial variance for J and 0.86% for R, against 2% and 1%; the plane explains 15.7% and 8.0%, against 20% and 10%. Only R's speed value is within tolerance. The values are lower than the paper's in all four cases and keep its order (J above R).
ellipse-orientation¶
The distributions of single-trial neural states for each condition are elliptical, and for most conditions their major axis points along the direction of the mean neural state from the origin.
for m in "JR":
t, X = data[m]
P = premotor.plane_coords(planes[m], X)
rows = []
for c, idx in t.groupby("condition").indices.items():
if len(idx) < 10:
continue
pts = P[idx]
mean = pts.mean(axis=0)
evals, evecs = np.linalg.eigh(np.cov(pts.T))
major = evecs[:, -1]
angle = np.degrees(np.arccos(min(1.0, abs(major @ mean) / np.linalg.norm(mean))))
rows.append({"condition": c, "config": t.config.iloc[idx[0]], "angle": angle,
"axis_ratio": np.sqrt(evals[-1] / evals[0]), "mean_radius": np.linalg.norm(mean)})
e = pd.DataFrame(rows)
w = stats.wilcoxon(e.angle - 45, alternative="less")
print(f"{m}: {len(e)} conditions; angle below 45 degrees in {int((e.angle < 45).sum())} ({np.mean(e.angle < 45):.2f}); "
f"median angle {e.angle.median():.1f} degrees; Wilcoxon against 45, p = {w.pvalue:.1e}; "
f"median ratio of major to minor axis {e.axis_ratio.median():.2f}")
print(" by configuration:")
print(e.groupby("config").angle.agg(conditions="count", median_angle="median",
below_45=lambda a: np.mean(a < 45)).round(2).to_string())
near = e.mean_radius < e.mean_radius.median()
print(f" conditions with the condition mean nearer the origin than the median: {np.mean(e.angle[near] < 45):.2f} below 45; "
f"farther: {np.mean(e.angle[~near] < 45):.2f}")
fig, ax = plt.subplots(1, 2, figsize=(9, 4.3))
for i, m in enumerate("JR"):
t, X = data[m]
P = premotor.plane_coords(planes[m], X)
sel = (t.config == "3-rings").values
for c, idx in t[sel].groupby("condition").indices.items():
g = P[np.flatnonzero(sel)[idx]]
mean = g.mean(axis=0)
evals, evecs = np.linalg.eigh(np.cov(g.T))
theta = np.linspace(0, 2 * np.pi, 60)
ell = mean + (evecs @ (np.sqrt(evals)[:, None] * np.stack([np.cos(theta), np.sin(theta)]))).T / 3
ax[i].plot(ell[:, 0], ell[:, 1], lw=0.8)
ax[i].set(title=f"monkey {m}, 3-rings, 1/3 SD ellipses", xlabel="x neural (mm)", ylabel="y neural (mm)", aspect="equal")
plt.tight_layout()
plt.show()
J: 96 conditions; angle below 45 degrees in 69 (0.72); median angle 21.0 degrees; Wilcoxon against 45, p = 4.1e-08; median ratio of major to minor axis 1.35
by configuration:
conditions median_angle below_45
config
3-rings 48 18.93 0.73
ring 24 30.46 0.62
vertical 24 21.13 0.79
conditions with the condition mean nearer the origin than the median: 0.75 below 45; farther: 0.69
R: 132 conditions; angle below 45 degrees in 99 (0.75); median angle 19.3 degrees; Wilcoxon against 45, p = 8.0e-10; median ratio of major to minor axis 1.56
by configuration:
conditions median_angle below_45
config
3-rings 48 31.06 0.65
horizontal 24 6.96 1.00
ring 36 13.42 0.78
vertical 24 38.21 0.67
conditions with the condition mean nearer the origin than the median: 0.67 below 45; farther: 0.83
Result: reproduced. The major axis lies within 45 degrees of the direction of the condition mean for 69 of 96 conditions in J (median 21.0 degrees, Wilcoxon p = 4.1e-08) and 99 of 132 in R (median 19.3 degrees, p = 8.0e-10). The median ratio of major to minor axis is 1.35 (J) and 1.56 (R).
on-greater-than-off¶
The mean absolute on-axis (distance) deviation of the neural state is larger than the off-axis (direction) deviation, in each task (J p = 9e-7 ring, 4e-7 horizontal, 1e-8 vertical, 3e-14 3-rings; R 8e-3, 4e-7, 2e-7 and 0.19).
PAPER_P = {"J": {"ring": 9e-7, "horizontal": 4e-7, "vertical": 1e-8, "3-rings": 3e-14},
"R": {"ring": 8e-3, "horizontal": 4e-7, "vertical": 2e-7, "3-rings": 0.19}}
onoff_rows = []
for m in "JR":
t, X = data[m]
P = premotor.plane_coords(planes[m], X)
dev = premotor.deviations(P, t.condition.values)
on, off = premotor.on_off(dev, t.tx.values, t.ty.values)
t = t.assign(on=np.abs(on), off=np.abs(off))
for task_column in ("methods_config", "config"):
for task in sorted(t[task_column].unique()):
g = t[t[task_column] == task]
per = g.groupby("condition")[["on", "off"]].mean()
pa = stats.ttest_rel(per.on, per.off).pvalue
pb = stats.ttest_ind(per.on, per.off, equal_var=False).pvalue
pc = stats.ttest_ind(g.on, g.off, equal_var=False).pvalue
onoff_rows.append({"monkey": m, "tasks by": "Methods" if task_column == "methods_config" else "file",
"task": task, "targets": len(per), "mean |on| mm": round(per.on.mean(), 2),
"mean |off| mm": round(per.off.mean(), 2), "targets on > off": int((per.on > per.off).sum()),
"(a) paired p": pa, "(b) Welch targets p": pb, "(c) Welch trials p": pc,
"paper p": PAPER_P[m][task]})
onoff = pd.DataFrame(onoff_rows)
print(onoff.to_string(index=False, formatters={k: "{:.1e}".format for k in ["(a) paired p", "(b) Welch targets p", "(c) Welch trials p", "paper p"]}))
main = onoff[(onoff["tasks by"] == "Methods") & (onoff["paper p"] < 0.05)]
for r in ["(a) paired p", "(b) Welch targets p", "(c) Welch trials p"]:
ok = (main["mean |on| mm"] > main["mean |off| mm"]) & (main[r] < 0.05)
print(f"{r}: on > off with p < 0.05 in {int(ok.sum())} of {len(main)} task and monkey pairs the paper finds significant")
monkey tasks by task targets mean |on| mm mean |off| mm targets on > off (a) paired p (b) Welch targets p (c) Welch trials p paper p
J Methods 3-rings 48 21.89 19.15 30 4.6e-03 4.3e-03 4.0e-14 3.0e-14
J Methods ring 24 21.65 21.15 13 6.9e-01 6.4e-01 3.1e-01 9.0e-07
J Methods vertical 24 17.98 15.00 19 4.3e-05 4.5e-02 1.0e-15 1.0e-08
J file 3-rings 48 21.89 19.15 30 4.6e-03 4.3e-03 4.0e-14 3.0e-14
J file ring 24 21.65 21.15 13 6.9e-01 6.4e-01 3.1e-01 9.0e-07
J file vertical 24 17.98 15.00 19 4.3e-05 4.5e-02 1.0e-15 1.0e-08
R Methods 3-rings 48 22.13 21.01 26 3.8e-01 5.2e-01 1.4e-01 1.9e-01
R Methods horizontal 24 17.29 14.55 22 3.0e-05 1.2e-01 3.1e-15 4.0e-07
R Methods ring 36 19.19 15.80 30 1.2e-04 7.4e-03 5.9e-20 8.0e-03
R Methods vertical 24 16.50 9.77 24 4.1e-07 5.7e-05 2.5e-119 2.0e-07
R file 3-rings 48 22.13 21.01 26 3.8e-01 5.2e-01 1.4e-01 1.9e-01
R file horizontal 24 16.50 9.77 24 4.1e-07 5.7e-05 2.5e-119 4.0e-07
R file ring 36 19.19 15.80 30 1.2e-04 7.4e-03 5.9e-20 8.0e-03
R file vertical 24 17.29 14.55 22 3.0e-05 1.2e-01 3.1e-15 2.0e-07
(a) paired p: on > off with p < 0.05 in 5 of 6 task and monkey pairs the paper finds significant
(b) Welch targets p: on > off with p < 0.05 in 4 of 6 task and monkey pairs the paper finds significant
(c) Welch trials p: on > off with p < 0.05 in 5 of 6 task and monkey pairs the paper finds significant
Result: partially reproduced. On-axis deviations exceed off-axis ones with p < 0.05 in 5 of the 6 task and monkey pairs under readings (a) and (c), and 4 of 6 under (b). The exception is J's ring task, which here is one session of four (13 of 24 targets with on > off, paired p = 0.69). R's 3-rings task shows no significant difference (paired p = 0.38), as in the paper.
relative-deviation-decreases¶
Relative deviations of the neural state (deviation as a fraction of reach distance), on and off axis, decrease with reach distance.
for m in "JR":
t, X = data[m]
P = premotor.plane_coords(planes[m], X)
on, off = premotor.on_off(premotor.deviations(P, t.condition.values), t.tx.values, t.ty.values)
dist = t.distance.values
for name, v in (("on-axis", on), ("off-axis", off)):
rel = np.abs(v) / dist
rho, p = stats.spearmanr(dist, rel)
rho_abs, p_abs = stats.spearmanr(dist, np.abs(v))
print(f"{m} {name}: Spearman of relative deviation with distance rho = {rho:.3f}, p = {p:.1e}; "
f"absolute deviation rho = {rho_abs:.3f}, p = {p_abs:.1e}")
binned = pd.DataFrame({"distance": np.round(dist / 10) * 10, "on": 100 * np.abs(on) / dist, "off": 100 * np.abs(off) / dist})
print(" mean relative deviation (% of distance) by distance (mm):")
print(binned.groupby("distance").mean().round(1).T.to_string())
J on-axis: Spearman of relative deviation with distance rho = -0.389, p = 0.0e+00; absolute deviation rho = 0.074, p = 3.3e-12 J off-axis: Spearman of relative deviation with distance rho = -0.378, p = 2.6e-295; absolute deviation rho = 0.103, p = 3.7e-22 mean relative deviation (% of distance) by distance (mm): distance 10.0 20.0 30.0 40.0 50.0 60.0 70.0 80.0 90.0 100.0 110.0 120.0 on 118.4 75.6 54.5 53.4 36.1 29.2 27.3 26.6 21.4 21.3 18.9 18.1 off 118.0 53.7 42.0 44.0 28.8 26.3 21.2 24.3 18.8 16.4 16.6 16.6 R on-axis: Spearman of relative deviation with distance rho = -0.404, p = 0.0e+00; absolute deviation rho = 0.073, p = 8.2e-17 R off-axis: Spearman of relative deviation with distance rho = -0.402, p = 0.0e+00; absolute deviation rho = 0.078, p = 5.5e-19 mean relative deviation (% of distance) by distance (mm): distance 10.0 20.0 30.0 40.0 50.0 60.0 70.0 80.0 90.0 100.0 110.0 120.0 on 169.5 80.9 64.7 48.4 37.1 30.9 25.1 24.0 18.5 18.5 16.0 16.9 off 131.9 58.3 42.7 40.3 27.7 21.6 17.6 19.8 14.1 13.0 11.3 13.5
Result: reproduced. Relative deviations fall with reach distance on axis (J rho = -0.389, R rho = -0.404) and off axis (J -0.378, R -0.402), all with p below 1e-290. The absolute deviations grow only slightly with distance (rho 0.073 to 0.103), so the fall comes mainly from dividing by distance, and the test would pass for any deviation that grows less than in proportion to distance; it discriminates little.
svm-ring-and-distance¶
Direction classification of the ring targets reaches 32.2% for J (24 targets) and 16% for R (36 targets); distance classification over 12 targets reaches 18.9% (J) and 21.4% (R).
SVM_SEED = 6
svm_out = {}
for m in "JR":
t, X = data[m]
ring = (t.config == "ring").values
y = np.round(np.degrees(np.arctan2(t.ty[ring], t.tx[ring]))).astype(int).values
pred, acc = premotor.svm_cv(X[ring], y, SVM_SEED)
print(f"{m} ring direction: {100 * acc.mean():.1f}% +- {100 * acc.std(ddof=1) / np.sqrt(len(acc)):.1f}% (SEM over 10 folds); "
f"chance {100 / len(np.unique(y)):.1f}%; trials {ring.sum()} in {len(np.unique(y))} targets, "
f"sessions {sorted(t[ring].date.unique())}")
for d in sorted(t[ring].date.unique()):
sel = ring & (t.date == d).values
ys = np.round(np.degrees(np.arctan2(t.ty[sel], t.tx[sel]))).astype(int).values
print(f" {d} alone: {100 * premotor.svm_cv(X[sel], ys, SVM_SEED)[1].mean():.1f}% ({sel.sum()} trials)")
rows, errors, errors_plane = [], [], []
P = premotor.plane_coords(planes[m], X)
for task in ("horizontal", "vertical"):
for side in (1, -1):
coord = t.tx if task == "horizontal" else t.ty
sel = ((t.config == task) & (np.sign(coord) == side)).values
if not sel.any():
continue
yd = np.round(t.distance[sel], 1).values
pred, acc = premotor.svm_cv(X[sel], (10 * yd).astype(int), SVM_SEED)
pred_p, acc_p = premotor.svm_cv(P[sel], (10 * yd).astype(int), SVM_SEED)
errors.append(np.abs(pred / 10 - yd) / 10) # cm
errors_plane.append(np.abs(pred_p / 10 - yd) / 10)
rows.append({"task": task, "side": side, "trials": int(sel.sum()), "accuracy %": 100 * acc.mean(),
"SEM %": 100 * acc.std(ddof=1) / np.sqrt(len(acc)), "plane-projected accuracy %": 100 * acc_p.mean()})
dist = pd.DataFrame(rows)
weighted = np.average(dist["accuracy %"], weights=dist.trials)
print(f"{m} distance (12 targets per side; chance 8.3%):")
print(dist.round(1).to_string(index=False))
print(f" accuracy weighted by trials: {weighted:.1f}%")
svm_out[m] = {"distance_errors": np.concatenate(errors), "distance_errors_plane": np.concatenate(errors_plane)}
Xl = premotor.monkey_table(sessions, m, units=premotor.FILE_PMD_UNITS)[1]
print(f" variant, units the file labels PMd: ring direction "
f"{100 * premotor.svm_cv(Xl[ring], y, SVM_SEED)[1].mean():.1f}%")
J ring direction: 24.7% +- 0.8% (SEM over 10 folds); chance 4.2%; trials 1020 in 24 targets, sessions ['2015-10-15']
2015-10-15 alone: 24.7% (1020 trials)
J distance (12 targets per side; chance 8.3%):
task side trials accuracy % SEM % plane-projected accuracy %
vertical 1 1603 19.1 0.9 17.9
vertical -1 1608 13.5 0.9 12.7
accuracy weighted by trials: 16.3%
variant, units the file labels PMd: ring direction 12.2%
R ring direction: 12.9% +- 0.4% (SEM over 10 folds); chance 2.8%; trials 3584 in 36 targets, sessions ['2017-01-24', '2017-01-25']
2017-01-24 alone: 11.1% (1606 trials)
2017-01-25 alone: 12.1% (1978 trials)
R distance (12 targets per side; chance 8.3%):
task side trials accuracy % SEM % plane-projected accuracy %
horizontal 1 1921 15.6 0.7 13.7
horizontal -1 1921 21.0 1.2 16.9
vertical 1 2043 18.2 0.8 16.8
vertical -1 2044 13.9 1.0 11.3
accuracy weighted by trials: 17.1%
variant, units the file labels PMd: ring direction 3.2%
Variant for this claim and svm-3rings-by-ring: hyperparameters chosen on the 3-rings task, as the paper did (its held-out 3-rings trials are not identified, so here the same 3-rings trials are used to choose them), then applied to the ring, distance and per-ring classifiers.
from sklearn.model_selection import GridSearchCV, StratifiedKFold # noqa: E402
from sklearn.pipeline import make_pipeline # noqa: E402
from sklearn.preprocessing import StandardScaler # noqa: E402
from sklearn.svm import SVC # noqa: E402
TUNE_SEED = 8
for m in "JR":
t, X = data[m]
sel = (t.config == "3-rings").values
y = (np.round(np.degrees(np.arctan2(t.ty[sel], t.tx[sel]))) * 1000 + np.round(t.distance[sel])).astype(int).values
grid = GridSearchCV(make_pipeline(StandardScaler(), SVC(kernel="rbf")),
{"svc__C": [0.3, 1, 3, 10], "svc__gamma": [0.1 / 96, 1 / 96, 3 / 96]},
cv=StratifiedKFold(3, shuffle=True, random_state=TUNE_SEED)).fit(X[sel], y)
C, gamma = grid.best_params_["svc__C"], grid.best_params_["svc__gamma"]
print(f"{m}: best on 3-rings (48 targets): C = {C}, gamma = {gamma * 96:.1f}/96, accuracy {100 * grid.best_score_:.1f}% "
f"(default C = 1, gamma = 1/96 after scaling: "
f"{100 * grid.cv_results_['mean_test_score'][list(grid.cv_results_['params']).index({'svc__C': 1, 'svc__gamma': 1 / 96})]:.1f}%)")
def tuned(Xs, ys):
accs = []
for tr, te in StratifiedKFold(10, shuffle=True, random_state=SVM_SEED).split(Xs, ys):
model = make_pipeline(StandardScaler(), SVC(kernel="rbf", C=C, gamma=gamma)).fit(Xs[tr], ys[tr])
accs.append(np.mean(model.predict(Xs[te]) == ys[te]))
return 100 * np.mean(accs)
ring = (t.config == "ring").values
yr = np.round(np.degrees(np.arctan2(t.ty[ring], t.tx[ring]))).astype(int).values
acc_d, n_d = [], []
for task in ("horizontal", "vertical"):
for side in (1, -1):
coord = t.tx if task == "horizontal" else t.ty
s_ = ((t.config == task) & (np.sign(coord) == side)).values
if s_.any():
acc_d.append(tuned(X[s_], (np.round(t.distance[s_], 1) * 10).astype(int).values))
n_d.append(s_.sum())
rings = [tuned(X[(t.config == "3-rings").values & (np.round(t.distance) == r).values],
np.round(np.degrees(np.arctan2(t.ty, t.tx)))[(t.config == "3-rings").values & (np.round(t.distance) == r).values].astype(int).values)
for r in (40, 80, 120)]
print(f" with these: ring direction {tuned(X[ring], yr):.1f}%, distance {np.average(acc_d, weights=n_d):.1f}%, "
f"3-rings by ring {[float(round(a, 1)) for a in rings]}%")
J: best on 3-rings (48 targets): C = 3, gamma = 0.1/96, accuracy 18.0% (default C = 1, gamma = 1/96 after scaling: 16.6%)
with these: ring direction 24.3%, distance 16.0%, 3-rings by ring [23.0, 35.6, 42.0]%
R: best on 3-rings (48 targets): C = 10, gamma = 0.1/96, accuracy 11.0% (default C = 1, gamma = 1/96 after scaling: 8.6%)
with these: ring direction 13.8%, distance 18.3%, 3-rings by ring [20.2, 22.9, 27.8]%
Result: not reproduced. All four accuracies are below the paper's: ring direction 24.7% (J, one session of four) and 12.9% (R) against 32.2 and 16%, and distance 16.3% (J, vertical only) and 17.1% (R) against 18.9 and 21.4%. All are well above chance (4.2, 2.8 and 8.3%). With hyperparameters chosen on the 3-rings trials the values are 24.3, 13.8, 16.0 and 18.3%.
svm-3rings-by-ring¶
Classifying the 16 directions of each ring of the 3-rings task separately, accuracy increases with ring distance; 23, 36 and 43% for J and 22, 28 and 29% for R.
RINGS_SEED = 7
PAPER_RINGS = {"J": [23, 36, 43], "R": [22, 28, 29]}
for m in "JR":
t, X = data[m]
accs, arc, ang = [], {}, {}
for radius in (40, 80, 120):
sel = ((t.config == "3-rings") & (np.round(t.distance) == radius)).values
y = np.round(np.degrees(np.arctan2(t.ty[sel], t.tx[sel]))).astype(int).values
pred, acc = premotor.svm_cv(X[sel], y, RINGS_SEED)
err = np.abs((pred - y + 180) % 360 - 180) # degrees
ang[radius] = err
arc[radius] = np.radians(err) * radius / 10 # cm
accs.append(100 * acc.mean())
print(f"{m} {radius // 10} cm ring: {100 * acc.mean():.1f}% +- {100 * acc.std(ddof=1) / np.sqrt(len(acc)):.1f}%, "
f"chance {100 / len(np.unique(y)):.2f}%, trials {sel.sum()}; mean angular error {err.mean():.1f} deg, "
f"mean arc error {arc[radius].mean():.2f} cm")
print(f" paper {PAPER_RINGS[m]}; within 10%: {[bool(abs(a - p) <= 0.1 * p) for a, p in zip(accs, PAPER_RINGS[m])]}; "
f"increasing: {accs[0] < accs[1] < accs[2]}")
svm_out[m].update(arc=arc, ang=ang)
J 4 cm ring: 23.4% +- 0.8%, chance 6.25%, trials 1519; mean angular error 41.0 deg, mean arc error 2.87 cm
J 8 cm ring: 34.2% +- 0.9%, chance 6.25%, trials 1493; mean angular error 24.3 deg, mean arc error 3.40 cm
J 12 cm ring: 42.3% +- 1.1%, chance 6.25%, trials 1512; mean angular error 18.3 deg, mean arc error 3.84 cm paper [23, 36, 43]; within 10%: [True, True, True]; increasing: True R 4 cm ring: 18.4% +- 1.4%, chance 6.25%, trials 450; mean angular error 52.0 deg, mean arc error 3.63 cm
R 8 cm ring: 21.9% +- 2.2%, chance 6.25%, trials 444; mean angular error 42.4 deg, mean arc error 5.92 cm R 12 cm ring: 22.0% +- 1.9%, chance 6.25%, trials 446; mean angular error 36.7 deg, mean arc error 7.69 cm paper [22, 28, 29]; within 10%: [False, False, False]; increasing: True
Result: partially reproduced. Accuracy rises with ring distance in both monkeys. J's values, 23.4, 34.2 and 42.3%, are within tolerance of 23, 36 and 43%; R's, 18.4, 21.9 and 22.0%, are below 22, 28 and 29%.
arc-error-distance-independent¶
Angular errors of direction classification decrease with ring distance, but errors measured as arc length are not significantly different between rings (F test, p > .05 for J and R).
for m in "JR":
arc, ang = svm_out[m]["arc"], svm_out[m]["ang"]
F, p = stats.f_oneway(*arc.values())
Fa, pa = stats.f_oneway(*ang.values())
print(f"{m}: mean angular error 4/8/12 cm {[float(round(ang[r].mean(), 1)) for r in (40, 80, 120)]} deg (ANOVA p = {pa:.1e}); "
f"mean arc error {[float(round(arc[r].mean(), 2)) for r in (40, 80, 120)]} cm; one-way ANOVA F = {F:.2f}, p = {p:.1e}")
limit = np.pi * 4
restricted = {r: a[a <= limit] for r, a in arc.items()}
Fr, pr = stats.f_oneway(*restricted.values())
print(f" variant, errors within the 4 cm ring's range (<= {limit:.2f} cm): means "
f"{[float(round(restricted[r].mean(), 2)) for r in (40, 80, 120)]} cm, F = {Fr:.2f}, p = {pr:.1e}; trials kept "
f"{[len(restricted[r]) for r in (40, 80, 120)]} of {[len(arc[r]) for r in (40, 80, 120)]}")
lev = stats.levene(*arc.values())
print(f" variant, Levene test of equal variances of arc errors: p = {lev.pvalue:.1e}; "
f"arc error 12 cm minus 4 cm ring {arc[120].mean() - arc[40].mean():.2f} cm")
J: mean angular error 4/8/12 cm [41.0, 24.3, 18.3] deg (ANOVA p = 1.7e-82); mean arc error [2.87, 3.4, 3.84] cm; one-way ANOVA F = 22.15, p = 2.7e-10 variant, errors within the 4 cm ring's range (<= 12.57 cm): means [2.87, 2.94, 3.05] cm, F = 1.49, p = 2.3e-01; trials kept [1519, 1451, 1443] of [1519, 1493, 1512] variant, Levene test of equal variances of arc errors: p = 4.6e-20; arc error 12 cm minus 4 cm ring 0.97 cm R: mean angular error 4/8/12 cm [52.0, 42.4, 36.7] deg (ANOVA p = 2.7e-07); mean arc error [3.63, 5.92, 7.69] cm; one-way ANOVA F = 52.16, p = 1.5e-22 variant, errors within the 4 cm ring's range (<= 12.57 cm): means [3.63, 4.39, 4.34] cm, F = 6.68, p = 1.3e-03; trials kept [450, 396, 354] of [450, 444, 446] variant, Levene test of equal variances of arc errors: p = 1.1e-14; arc error 12 cm minus 4 cm ring 4.06 cm
Result: not reproduced. Mean angular errors fall with ring distance (J 41.0, 24.3, 18.3 degrees; R 52.0, 42.4, 36.7), but arc-length errors grow with it (J 2.87, 3.40, 3.84 cm; R 3.63, 5.92, 7.69 cm), and the one-way ANOVA finds the rings different (J F = 22.15, p = 2.7e-10; R F = 52.16, p = 1.5e-22). Restricted to errors within the 4 cm ring's range, J's rings no longer differ (p = 0.23) and R's still do (p = 1.3e-03).
distance-error-and-plane-projection¶
Mean distance-classification error is slightly but significantly lower than the arc-length error of direction classification, and projecting the activity onto the spatial plane before classification raises it (J 2.64 to 3.03 cm; R given as 2.78 and 3.27 cm), to the level of the direction errors.
for m in "JR":
full, plane_err = svm_out[m]["distance_errors"], svm_out[m]["distance_errors_plane"]
w = stats.ttest_ind(plane_err, full, equal_var=False)
pr = stats.ttest_rel(plane_err, full)
print(f"{m}: distance error full {full.mean():.3f} cm, plane-projected {plane_err.mean():.3f} cm "
f"(Welch t = {w.statistic:.2f}, p = {w.pvalue:.1e}; paired over trials t = {pr.statistic:.2f}, p = {pr.pvalue:.1e}; n = {len(full)})")
for r in (40, 80, 120):
a = svm_out[m]["arc"][r]
print(f" direction arc error, {r // 10} cm ring {a.mean():.2f} cm; full distance error lower: "
f"p = {stats.ttest_ind(full, a, equal_var=False).pvalue:.1e}; plane-projected against arc: "
f"p = {stats.ttest_ind(plane_err, a, equal_var=False).pvalue:.1e}")
# variant: arc errors within the 4 cm ring's range, as in the Fig 4d inset
limit = np.pi * 4
for r in (40, 80, 120):
a = svm_out[m]["arc"][r]
a = a[a <= limit]
print(f" variant, arc error within the 4 cm ring's range (<= {limit:.2f} cm), {r // 10} cm ring {a.mean():.2f} cm "
f"(n = {len(a)}); full distance error lower: p = {stats.ttest_ind(full, a, equal_var=False).pvalue:.1e}; "
f"plane-projected against arc: p = {stats.ttest_ind(plane_err, a, equal_var=False).pvalue:.1e}")
J: distance error full 2.505 cm, plane-projected 2.801 cm (Welch t = 5.15, p = 2.7e-07; paired over trials t = 6.93, p = 4.9e-12; n = 3211) direction arc error, 4 cm ring 2.87 cm; full distance error lower: p = 3.1e-05; plane-projected against arc: p = 4.6e-01 direction arc error, 8 cm ring 3.40 cm; full distance error lower: p = 4.7e-16; plane-projected against arc: p = 7.1e-08 direction arc error, 12 cm ring 3.84 cm; full distance error lower: p = 3.0e-23; plane-projected against arc: p = 1.3e-14 variant, arc error within the 4 cm ring's range (<= 12.57 cm), 4 cm ring 2.87 cm (n = 1519); full distance error lower: p = 3.1e-05; plane-projected against arc: p = 4.6e-01 variant, arc error within the 4 cm ring's range (<= 12.57 cm), 8 cm ring 2.94 cm (n = 1451); full distance error lower: p = 3.0e-07; plane-projected against arc: p = 1.1e-01 variant, arc error within the 4 cm ring's range (<= 12.57 cm), 12 cm ring 3.05 cm (n = 1443); full distance error lower: p = 7.4e-10; plane-projected against arc: p = 5.3e-03 R: distance error full 2.503 cm, plane-projected 3.026 cm (Welch t = 13.46, p = 4.3e-41; paired over trials t = 16.86, p = 1.2e-62; n = 7929) direction arc error, 4 cm ring 3.63 cm; full distance error lower: p = 9.2e-13; plane-projected against arc: p = 1.0e-04 direction arc error, 8 cm ring 5.92 cm; full distance error lower: p = 3.8e-31; plane-projected against arc: p = 1.3e-23 direction arc error, 12 cm ring 7.69 cm; full distance error lower: p = 9.4e-36; plane-projected against arc: p = 4.6e-30 variant, arc error within the 4 cm ring's range (<= 12.57 cm), 4 cm ring 3.63 cm (n = 450); full distance error lower: p = 9.2e-13; plane-projected against arc: p = 1.0e-04 variant, arc error within the 4 cm ring's range (<= 12.57 cm), 8 cm ring 4.39 cm (n = 396); full distance error lower: p = 4.7e-21; plane-projected against arc: p = 3.2e-12 variant, arc error within the 4 cm ring's range (<= 12.57 cm), 12 cm ring 4.34 cm (n = 354); full distance error lower: p = 4.9e-23; plane-projected against arc: p = 3.2e-13
Result: reproduced. Under reading (a), mean distance errors are 2.505 cm (J) and 2.503 cm (R) with the full activity and 2.801 and 3.026 cm after projection onto the plane, all within tolerance; R's full value is inside its bound (2.502 cm) by 0.001 cm, so this match is fragile. Projection raises the error in both (Welch J p = 2.7e-07, R p = 4.3e-41; paired over trials J p = 4.9e-12, R p = 1.2e-62). The plane is fit on all trials, including each test fold, which favours the projected classifier, so the rise is if anything understated. The full-activity distance error is lower than the direction arc errors of every ring in both monkeys (p at most 3.1e-05). After projection it is no different from the arc error of J's 4 cm ring (p = 0.46) and differs from the others, so the claim that projection brings it to the level of the direction errors holds for 1 of 6 rings. With the arc errors restricted to the 4 cm ring's range, as in the paper's Figure 4d inset, the full-activity error is still lower than every ring's (p at most 3.1e-05), and the projected error is no different from J's 4 and 8 cm rings (p = 0.46 and 0.11) but differs from J's 12 cm ring (p = 5.3e-03) and from all of R's (p at most 1.0e-04), so that part holds for 2 of 6 rings.
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 |
|---|---|
| session-configurations | partially reproduced |
| r-trial-counts | partially reproduced |
| j-sessions-complete | blocked (data missing) |
| delays-and-times | partially reproduced |
| channels-task-correlated | reproduced |
| pca-variance | not reproduced |
| spatial-plane-variance | not reproduced |
| plane-timecourse | partially reproduced |
| onoff-predicts-kinematics | reproduced |
| speed-axis-within-distance | partially reproduced |
| distance-after-speed | reproduced |
| speed-distance-projection | not reproduced |
| single-trial-variance | not reproduced |
| ellipse-orientation | reproduced |
| on-greater-than-off | partially reproduced |
| relative-deviation-decreases | reproduced |
| svm-ring-and-distance | not reproduced |
| svm-3rings-by-ring | partially reproduced |
| arc-error-distance-independent | not reproduced |
| distance-error-and-plane-projection | reproduced |
How the analysis differs from the paper's¶
- Monkey J is represented by 5 of its 16 sessions and 96 of its 120 conditions; its PCA, plane and classifiers are fit on these alone. Monkey R's data are complete.
- The PMd array is taken as units 96 to 191, against the units'
electrode_grouplabels, from the data check; results on the labelled array are shown as variants where they bear on a claim. - The paper does not say how the day's average rate, movement onset, maximum speed, the SVM
hyperparameters, the folds or the null of its p-values were set; the choices here are listed in
the claims file and in
premotor.py. The day's mean rate is the rate over trial time, as spikes exist only within trials. - The spatial plane is fit on all trials, so the plane-projected distance classifier sees its test folds through the plane; the paper does not say how it handled this.
- R's horizontal and vertical tasks are named by the Methods' dates; the file geometry is shown as a variant. R's horizontal targets lie at 5.2 + 10.4 k mm, so its 12 distances are matched by rank.
What was found about the dandiset¶
- The units'
electrode_grouplabels are swapped: the array labelled PMd is the one whose units point to M1 electrode rows and which, in R, has 29 silent channels and almost no target information in the decoding (at most 6 significant channels per session in the delay window), though 40 of its 96 channels regress on target position when R's sessions are pooled. - In every J file the boundary between units 0 and 1 is misplaced, so unit 0 stops early and unit 1's spike times run past the session's end and restart.
- Positions are stated in meters with
conversionNaN but are in mm;task_idis the same in every trial, so the configuration must be inferred from the targets;is_successfulcarries a description of something else; electrode positions and impedances are empty. - The archive's summary gives 15 files and 187.5 GB for a version of 12 files and 39.5 GB.
- 11 of J's 16 sessions, including all its horizontal sessions, are not in the dandiset.
- The files place R's horizontal and vertical tasks on the dates the Methods give the other way round; R has almost no trials without a delay and delays up to 1000 ms; its 3-rings session holds 1340 successful delayed trials against the 912 the paper used.
- Spikes and hand samples exist only within trials; the gaps between trials (484 to 1963 s per session) hold almost no data.