Reproducing Perich, Gallego & Miller (2018) from dandiset 000688¶
Paper. Perich MG, Gallego JA, Miller LE (2018). A neural population mechanism for rapid learning. Neuron 100, 964-976. doi:10.1016/j.neuron.2018.09.030. The claims file was written from the accepted author manuscript in PubMed Central (PMC6250582), with its eight supplement files (Figures S1 to S7 and their legends); this attempt read the same text.
Dandiset. DANDI:000688, version 0.250122.1735. The paper cites no DANDI version (it predates the dandiset). This is the only published version; the draft has the same 111 files and was modified two days later, so the published version is pinned.
What the dandiset holds. 111 sessions of four monkeys (C, J, M, T) in center-out (CO) and random-target (RT) reaching, with sorted units, cursor kinematics and a trial table. The trial tables have no column for the epoch (Baseline, Force or Rotation, Washout) or the perturbation, and none of the files is described as a curl field or visuomotor rotation session. The dual-array (M1 and PMd) center-out sessions of Monkeys C (2016) and M (February and March 2014) have unit counts that match the paper's ranges, except one Monkey C session with too few PMd units, but each holds about 200 reward trials with small, steady reach errors: the size of one Baseline epoch, without the Force or Rotation and Washout epochs. The data check below shows this. Only the claims that use Baseline trials alone can be tested.
| Claim | Tested here |
|---|---|
| methods-neuron-counts | tested |
| fig3b-pmd-higher-dimensional | tested (Baseline trials, as the method uses) |
| figs4d-dimensionality-robust-to-size | tested |
| figs4c-null-potent-index-centered | tested (subspaces from Baseline, as the method uses) |
| fig6c-null-leads-potent | tested for Baseline; the end-of-Force part cannot be |
| the 23 other claims (behavior, covariance across epochs, GLMs, Washout, Mahalanobis distance, VR) | blocked, data missing |
import os
os.environ.setdefault("OMP_NUM_THREADS", "1")
os.environ.setdefault("OPENBLAS_NUM_THREADS", "1")
import itertools
import time
import warnings
from concurrent.futures import ThreadPoolExecutor
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from scipy import stats
from sklearn.discriminant_analysis import LinearDiscriminantAnalysis
from sklearn.model_selection import StratifiedKFold, cross_val_score
import perich
pd.set_option("display.width", 220)
pd.set_option("display.max_columns", 40)
pd.set_option("display.max_rows", 200)
warnings.filterwarnings("ignore", category=RuntimeWarning)
started = time.time()
urls = perich.asset_urls()
with ThreadPoolExecutor(16) as ex:
meta = dict(zip(urls, ex.map(perich.read_metadata, urls.values())))
survey = pd.DataFrame([{"session": k, "monkey": k[4], "task": k.split("_")[1].split("-")[1],
"date": k.split("-")[-1], **v} for k, v in meta.items()]).sort_values("session")
print(f"read the metadata of {len(survey)} files of dandiset {perich.DANDISET_ID} version "
f"{perich.DANDISET_VERSION} in {time.time() - started:.0f} s")
SESSIONS = perich.CANDIDATE_SESSIONS
sessions = perich.load_sessions(SESSIONS, urls)
by_id = {s["sid"]: s for s in sessions}
print(f"loaded spikes, kinematics and trials of {len(sessions)} sessions in {time.time() - started:.0f} s")
read the metadata of 111 files of dandiset 000688 version 0.250122.1735 in 88 s
loaded spikes, kinematics and trials of 21 sessions in 108 s
Data check¶
Each check of the metadata against the values is stated under its output.
print("Files per monkey and task:")
print(pd.crosstab(survey.monkey, survey.task).to_string())
print("\nTrial table columns (sets of columns and how many files have each):")
print(survey.trial_columns.apply(", ".join).value_counts().to_string())
print("\nExperiment descriptions:", survey.experiment.value_counts().to_dict())
print("\nSession descriptions (text, number of files):")
print(survey.session_description.value_counts().to_string())
warned = survey[survey.session_description.str.contains("Warning")]
print("Files whose session description holds a warning:", list(warned.session))
survey["arrays"] = np.where((survey.n_M1 > 0) & (survey.n_PMd > 0), "M1+PMd", np.where(survey.n_PMd > 0, "PMd", "M1"))
print("\nArrays per monkey (sessions):")
print(pd.crosstab(survey.monkey, survey.arrays).to_string())
survey["rewards"] = survey.results.apply(lambda r: r.get("R", 0))
fp = survey.groupby(["n_trials", "n_M1", "n_PMd", "n_spikes"]).session.apply(list)
print("\nFiles sharing trial count, unit counts per area and total spike count:", [v for v in fp if len(v) > 1] or "none")
print("\nAll files, with units per area and trial outcomes:")
print(survey[["session", "task", "arrays", "n_M1", "n_PMd", "n_trials", "rewards", "last_stop"]].to_string(index=False))
Files per monkey and task:
task CO RT
monkey
C 53 15
J 3 0
M 22 6
T 6 6
Trial table columns (sets of columns and how many files have each):
trial_columns
go_cue_time, result, start_time, stop_time, target_corners, target_dir, target_id, target_on_time 84
go_cue_time_array, num_attempted, num_targets, result, start_time, stop_time, target_dir, target_id, target_size 27
Experiment descriptions: {'Center-out delayed reaching task': 84, 'Random target reaching task': 27}
Session descriptions (text, number of files):
session_description
Monkey C performing center-out reaching task. 46
Monkey M performing center-out reaching task. 22
Monkey C performing random target reaching task. 15
Monkey C performing center-out reaching task. Warning: This session with Monkey 7
Monkey M performing random target reaching task. 6
Monkey T performing center-out reaching task. 6
Monkey T performing random target reaching task. 6
Monkey J performing center-out reaching task. 3
Files whose session description holds a warning: ['sub-C_ses-CO-20160909', 'sub-C_ses-CO-20160912', 'sub-C_ses-CO-20160914', 'sub-C_ses-CO-20160915', 'sub-C_ses-CO-20160919', 'sub-C_ses-CO-20160921', 'sub-C_ses-CO-20160923']
Arrays per monkey (sessions):
arrays M1 M1+PMd
monkey
C 54 14
J 3 0
M 2 26
T 0 12
Files sharing trial count, unit counts per area and total spike count: none
All files, with units per area and trial outcomes:
session task arrays n_M1 n_PMd n_trials rewards last_stop
sub-C_ses-CO-20131003 CO M1 71 0 193 159 659.508
sub-C_ses-CO-20131022 CO M1 41 0 182 152 601.594
sub-C_ses-CO-20131023 CO M1 61 0 224 193 714.455
sub-C_ses-CO-20131031 CO M1 49 0 247 218 1161.706
sub-C_ses-CO-20131101 CO M1 41 0 280 246 885.990
sub-C_ses-CO-20131203 CO M1 48 0 208 166 1268.851
sub-C_ses-CO-20131204 CO M1 52 0 178 163 635.943
sub-C_ses-CO-20131219 CO M1 65 0 212 180 708.550
sub-C_ses-CO-20131220 CO M1 67 0 233 198 762.641
sub-C_ses-CO-20150309 CO M1 74 0 788 606 3461.347
sub-C_ses-CO-20150311 CO M1 88 0 1067 912 3628.740
sub-C_ses-CO-20150312 CO M1 91 0 1058 907 3416.915
sub-C_ses-CO-20150313 CO M1 86 0 1179 1039 3708.057
sub-C_ses-CO-20150319 CO M1 72 0 1155 1026 3743.715
sub-C_ses-CO-20150629 CO M1 49 0 295 179 1201.780
sub-C_ses-CO-20150630 CO M1 44 0 243 178 905.540
sub-C_ses-CO-20150701 CO M1 49 0 282 196 1029.207
sub-C_ses-CO-20150703 CO M1 52 0 229 185 984.451
sub-C_ses-CO-20150706 CO M1 46 0 221 184 921.608
sub-C_ses-CO-20150707 CO M1 42 0 222 189 900.677
sub-C_ses-CO-20150708 CO M1 64 0 214 190 865.748
sub-C_ses-CO-20150709 CO M1 61 0 218 196 900.515
sub-C_ses-CO-20150710 CO M1 55 0 231 203 946.780
sub-C_ses-CO-20150713 CO M1 58 0 218 185 912.872
sub-C_ses-CO-20150714 CO M1 66 0 227 199 901.526
sub-C_ses-CO-20150715 CO M1 60 0 212 192 842.948
sub-C_ses-CO-20150716 CO M1 61 0 207 191 814.617
sub-C_ses-CO-20151103 CO M1 38 0 320 250 1200.216
sub-C_ses-CO-20151104 CO M1 59 0 440 320 1641.194
sub-C_ses-CO-20151106 CO M1 60 0 303 258 1297.013
sub-C_ses-CO-20151109 CO M1 65 0 210 175 895.942
sub-C_ses-CO-20151110 CO M1 61 0 219 199 872.423
sub-C_ses-CO-20151112 CO M1 42 0 202 188 815.704
sub-C_ses-CO-20151113 CO M1 53 0 200 188 792.862
sub-C_ses-CO-20151116 CO M1 55 0 201 169 793.646
sub-C_ses-CO-20151117 CO M1 58 0 221 186 838.780
sub-C_ses-CO-20151119 CO M1 71 0 226 208 901.422
sub-C_ses-CO-20151120 CO M1 55 0 198 169 803.936
sub-C_ses-CO-20151201 CO M1 53 0 219 189 851.867
sub-C_ses-CO-20160909 CO M1+PMd 55 191 335 155 1051.200
sub-C_ses-CO-20160912 CO M1+PMd 79 162 326 201 1121.617
sub-C_ses-CO-20160914 CO M1+PMd 95 258 322 195 1074.391
sub-C_ses-CO-20160915 CO M1+PMd 76 233 291 182 959.728
sub-C_ses-CO-20160919 CO M1+PMd 60 235 276 186 963.667
sub-C_ses-CO-20160921 CO M1+PMd 74 232 271 176 904.827
sub-C_ses-CO-20160923 CO M1+PMd 76 223 278 203 903.367
sub-C_ses-CO-20160929 CO M1+PMd 74 114 253 209 899.306
sub-C_ses-CO-20161005 CO M1+PMd 82 167 252 202 898.315
sub-C_ses-CO-20161006 CO M1+PMd 63 192 251 210 902.003
sub-C_ses-CO-20161007 CO M1+PMd 70 137 214 171 783.001
sub-C_ses-CO-20161011 CO M1+PMd 71 171 239 199 901.538
sub-C_ses-CO-20161013 CO M1+PMd 87 138 271 239 1033.276
sub-C_ses-CO-20161021 CO M1+PMd 84 211 331 286 1219.682
sub-C_ses-RT-20131009 RT M1 74 0 187 155 1102.790
sub-C_ses-RT-20131010 RT M1 74 0 193 171 811.358
sub-C_ses-RT-20131011 RT M1 88 0 175 162 999.753
sub-C_ses-RT-20131028 RT M1 77 0 174 151 748.503
sub-C_ses-RT-20131029 RT M1 67 0 194 179 827.826
sub-C_ses-RT-20131209 RT M1 65 0 205 178 949.782
sub-C_ses-RT-20131210 RT M1 63 0 203 181 854.338
sub-C_ses-RT-20131212 RT M1 61 0 198 182 866.181
sub-C_ses-RT-20131213 RT M1 67 0 198 179 911.995
sub-C_ses-RT-20131217 RT M1 57 0 197 177 1013.845
sub-C_ses-RT-20131218 RT M1 62 0 188 178 801.297
sub-C_ses-RT-20150316 RT M1 75 0 654 587 3533.975
sub-C_ses-RT-20150317 RT M1 78 0 694 623 3555.968
sub-C_ses-RT-20150318 RT M1 75 0 921 847 3806.741
sub-C_ses-RT-20150320 RT M1 67 0 841 772 3509.538
sub-J_ses-CO-20160405 CO M1 38 0 322 203 1276.733
sub-J_ses-CO-20160406 CO M1 18 0 296 195 1003.163
sub-J_ses-CO-20160407 CO M1 19 0 331 199 1076.518
sub-M_ses-CO-20140203 CO M1+PMd 34 80 245 182 1361.671
sub-M_ses-CO-20140217 CO M1+PMd 44 104 255 208 1282.556
sub-M_ses-CO-20140218 CO M1+PMd 38 121 259 225 1234.635
sub-M_ses-CO-20140303 CO M1+PMd 52 66 233 209 1049.480
sub-M_ses-CO-20140304 CO M1+PMd 39 76 235 203 1310.543
sub-M_ses-CO-20140306 CO M1+PMd 43 86 264 219 1254.860
sub-M_ses-CO-20140307 CO M1+PMd 26 66 254 216 1191.687
sub-M_ses-CO-20140626 CO M1 26 0 1198 671 5027.153
sub-M_ses-CO-20140627 CO M1 25 0 1250 609 5177.220
sub-M_ses-CO-20140929 CO M1+PMd 21 124 805 733 4369.554
sub-M_ses-CO-20141203 CO M1+PMd 10 48 223 194 1509.702
sub-M_ses-CO-20150511 CO M1+PMd 24 44 621 536 4328.158
sub-M_ses-CO-20150512 CO M1+PMd 30 60 563 443 6277.970
sub-M_ses-CO-20150610 CO M1+PMd 18 31 208 190 949.504
sub-M_ses-CO-20150611 CO M1+PMd 16 39 215 196 1089.981
sub-M_ses-CO-20150612 CO M1+PMd 15 22 206 175 1185.509
sub-M_ses-CO-20150615 CO M1+PMd 17 24 205 174 1163.632
sub-M_ses-CO-20150616 CO M1+PMd 17 37 211 183 1072.404
sub-M_ses-CO-20150617 CO M1+PMd 17 40 223 175 1199.071
sub-M_ses-CO-20150623 CO M1+PMd 20 32 207 183 929.340
sub-M_ses-CO-20150625 CO M1+PMd 27 32 219 174 1038.991
sub-M_ses-CO-20150626 CO M1+PMd 22 40 208 184 901.485
sub-M_ses-RT-20140114 RT M1+PMd 59 106 252 215 1313.837
sub-M_ses-RT-20140115 RT M1+PMd 70 94 246 216 1294.975
sub-M_ses-RT-20140116 RT M1+PMd 67 94 198 174 1027.706
sub-M_ses-RT-20140214 RT M1+PMd 22 94 201 178 1121.662
sub-M_ses-RT-20140221 RT M1+PMd 32 95 209 187 951.610
sub-M_ses-RT-20140224 RT M1+PMd 32 89 182 165 814.843
sub-T_ses-CO-20130819 CO M1+PMd 9 46 179 144 897.052
sub-T_ses-CO-20130821 CO M1+PMd 6 53 197 164 899.331
sub-T_ses-CO-20130823 CO M1+PMd 4 41 256 213 1091.611
sub-T_ses-CO-20130903 CO M1+PMd 8 45 216 168 941.450
sub-T_ses-CO-20130905 CO M1+PMd 10 55 212 175 904.132
sub-T_ses-CO-20130909 CO M1+PMd 12 47 280 218 1226.097
sub-T_ses-RT-20130820 RT M1+PMd 10 44 154 144 897.745
sub-T_ses-RT-20130822 RT M1+PMd 8 56 247 224 1179.117
sub-T_ses-RT-20130830 RT M1+PMd 3 33 204 180 910.215
sub-T_ses-RT-20130904 RT M1+PMd 4 69 162 149 798.905
sub-T_ses-RT-20130906 RT M1+PMd 12 41 216 196 1063.453
sub-T_ses-RT-20130910 RT M1+PMd 6 63 229 210 994.168
No trial table has a column for the epoch or the perturbation. The experiment descriptions say only "center-out" or "random target", and the session descriptions only which monkey performed which task, except in seven Monkey C sessions (2016-09-09 to 2016-09-23), where the text ends in "Warning: This session with Monkey": the warning is cut off in the file and in the archive's metadata, so what it warns of is unknown. These are the seven sessions whose cursor positions do not match the stored targets (below). The paper's sessions are center-out sessions of Monkeys C and M with arrays in both M1 and PMd (STAR Methods, Neural recordings). Monkey C has such sessions only from 2016-09-09 to 2016-10-21 (14). Monkey M has 7 from 2014-02-03 to 2014-03-07 and 13 later ones (2014-09-29, 2014-12-03, two in May 2015 and nine in June 2015), with 10 to 30 M1 and 22 to 124 PMd units and 174 to 733 rewards. Of the later ones, only 2015-05-12 (30 M1, 60 PMd) is inside the paper's Monkey M ranges widened by 10%; the others have too few M1 units (2014-09-29 has 21) or too few PMd units. The 7 February and March sessions match the paper's count of Monkey M sessions and its bounds (66 and 121 PMd, 26 M1), so they are taken as the paper's, and 2015-05-12, with 443 rewards, is left out. That choice rests on the fit of the whole block, not on a rule applied to each session: C 2016-10-13 and 2016-10-21 are also outside the Baseline reward range and are kept among the Monkey C candidates (the second session set leaves them out). The 21 sessions of the two February to March 2014 and September to October 2016 blocks are the candidates analysed here. No two files share their trial count, unit counts and spike count, so no recording is deposited twice.
PAPER = {"C": {"PMd": (137, 256), "M1": (55, 93)}, "M": {"PMd": (66, 121), "M1": (26, 51)}}
BASELINE_REWARDS = {"CF": (170, 225), "VR": (154, 217)} # STAR Methods, Behavioral task
cand = []
for s in sessions:
tr = s["trials"]
cand.append({"session": s["sid"], "monkey": s["monkey"], "n_M1": int((s["area"] == "M1").sum()),
"n_PMd": int((s["area"] == "PMd").sum()), "trials": len(tr),
"rewards": int((tr.result == "R").sum()),
"minutes": round((tr.stop.iloc[-1] - tr.start.iloc[0]) / 60, 1)})
cand = pd.DataFrame(cand)
lo = min(v[0] for v in BASELINE_REWARDS.values()); hi = max(v[1] for v in BASELINE_REWARDS.values())
cand["in_baseline_range"] = cand.rewards.between(lo, hi)
print(cand.to_string(index=False))
print(f"\nPaper: Baseline 170-225 rewards (CF) or 154-217 successful trials (VR); Force 201-337, Rotation 219-316; "
f"Washout 153-404 (CF) or 162-348 (VR). A whole CF session therefore has at least {170+201+153} rewards.")
print("Reward trials per candidate session:", cand.rewards.min(), "to", cand.rewards.max(),
f"; within the Baseline range {lo}-{hi}:", int(cand.in_baseline_range.sum()), "of", len(cand))
SETS = {"all 21": list(cand.session), "19 within Baseline reward range": list(cand.session[cand.in_baseline_range])}
for k, v in SETS.items():
print(f"session set '{k}':", len(v), "sessions;", pd.Series([x[4] for x in v]).value_counts().to_dict())
session monkey n_M1 n_PMd trials rewards minutes in_baseline_range
sub-C_ses-CO-20160909 C 55 191 335 155 17.4 True
sub-C_ses-CO-20160912 C 79 162 326 201 18.6 True
sub-C_ses-CO-20160914 C 95 258 322 195 17.9 True
sub-C_ses-CO-20160915 C 76 233 291 182 15.9 True
sub-C_ses-CO-20160919 C 60 235 276 186 16.0 True
sub-C_ses-CO-20160921 C 74 232 271 176 15.0 True
sub-C_ses-CO-20160923 C 76 223 278 203 15.0 True
sub-C_ses-CO-20160929 C 74 114 253 209 14.9 True
sub-C_ses-CO-20161005 C 82 167 252 202 15.0 True
sub-C_ses-CO-20161006 C 63 192 251 210 15.0 True
sub-C_ses-CO-20161007 C 70 137 214 171 13.0 True
sub-C_ses-CO-20161011 C 71 171 239 199 15.0 True
sub-C_ses-CO-20161013 C 87 138 271 239 17.1 False
sub-C_ses-CO-20161021 C 84 211 331 286 20.3 False
sub-M_ses-CO-20140203 M 34 80 245 182 22.5 True
sub-M_ses-CO-20140217 M 44 104 255 208 20.7 True
sub-M_ses-CO-20140218 M 38 121 259 225 20.3 True
sub-M_ses-CO-20140303 M 52 66 233 209 17.3 True
sub-M_ses-CO-20140304 M 39 76 235 203 21.8 True
sub-M_ses-CO-20140306 M 43 86 264 219 20.5 True
sub-M_ses-CO-20140307 M 26 66 254 216 19.8 True
Paper: Baseline 170-225 rewards (CF) or 154-217 successful trials (VR); Force 201-337, Rotation 219-316; Washout 153-404 (CF) or 162-348 (VR). A whole CF session therefore has at least 524 rewards.
Reward trials per candidate session: 155 to 286 ; within the Baseline range 154-225: 19 of 21
session set 'all 21': 21 sessions; {'C': 14, 'M': 7}
session set '19 within Baseline reward range': 19 sessions; {'C': 12, 'M': 7}
Every candidate session has 154 to 286 reward trials, while a whole curl field session of the paper has at least 524 (Baseline, Force and Washout together). 19 of the 21 fall inside the paper's Baseline range (154 to 225); C 2016-10-13 (239) and C 2016-10-21 (286) do not. The files therefore hold one epoch each, consistent with the Baseline epoch only; the kinematics below check that no block of perturbed reaches is present. The paper has 16 sessions (Monkey C five CF and four VR, Monkey M four CF and three VR). The 7 Monkey M sessions match its count; Monkey C has 14 candidates for 9, and nothing in the files says which are the paper's or which were CF or VR sessions. Two session sets are used throughout: all 21 candidates, and the 19 within the Baseline reward range.
rows = []
for s in sessions:
tr = s["trials"]; r = tr[tr.result == "R"]
t, pos = s["pos_t"], s["pos"]
end = pos[np.clip(np.searchsorted(t, r.stop.values), 0, len(t) - 1)]
errs = perich.initial_errors(s)
errs["dev"] = errs.error - errs.groupby("target_id").error.transform("mean")
rho_s = stats.spearmanr(np.arange(len(errs)), errs.dev)
# end of movement: last sample before the stop at which speed exceeds 20% of the trial's peak,
# compared with the assumed acquisition time (stop - 0.5 s hold); finite speeds only
with np.errstate(all="ignore"):
speed = np.hypot(*np.gradient(pos, t, axis=0).T)
move_end = []
for go, stop, acq in zip(r.go_cue, r.stop, r.acquire):
a, b = np.searchsorted(t, [go, stop])
sp = speed[a:b]; ok = np.isfinite(sp)
if ok.sum() < 5:
continue
fast = np.flatnonzero(ok & (sp > 0.2 * sp[ok].max()))
move_end.append(t[a + fast[-1]] - acq)
run = errs.dev.abs().rolling(8).mean()
rho = stats.spearmanr(np.arange(len(errs)), errs.dev.abs())
rows.append({"session": s["sid"], "pos_unit": s["pos_unit"], "pos_hz": round(1 / np.median(np.diff(t))),
"repeated_timestamps": int((np.diff(t) <= 0).sum()),
"target_radius_cm": round(np.nanmedian(r.target_radius), 2),
"cursor_radius_at_stop": round(np.median(np.hypot(*end.T)), 2),
"peak_speed": round(errs.peak_speed.median(), 1),
"cursor_enters_target": round(r.cursor_enters_target.notna().mean(), 2),
"stop_minus_go_s": round((r.stop - r.go_cue).median(), 2),
"delay_s": f"{(r.go_cue - r.target_on).min():.2f}-{(r.go_cue - r.target_on).max():.2f}",
"move_end_minus_acquire_s": round(np.median(move_end), 2),
"signed_dev_first8": round(errs.dev[:8].mean(), 1), "signed_dev_last8": round(errs.dev[-8:].mean(), 1),
"rho_trial_signed": round(rho_s.statistic, 2), "p_rho_signed": round(rho_s.pvalue, 3),
"abs_dev_first8": round(errs.dev.abs()[:8].mean(), 1),
"abs_dev_last8": round(errs.dev.abs()[-8:].mean(), 1),
"max_8trial_mean": round(run.max(), 1), "median_abs_dev": round(errs.dev.abs().median(), 1),
"rho_trial_absdev": round(rho.statistic, 2), "p_rho": round(rho.pvalue, 3)})
kin = pd.DataFrame(rows)
print("Kinematics of reward trials (positions in the file's unit; errors in degrees, 150 ms after movement onset at 20% of peak speed,")
print("minus the session's mean per target):")
print(kin.to_string(index=False))
print("\nEnd of movement minus assumed acquisition (s), median over sessions:", kin.move_end_minus_acquire_s.median(),
"; range", kin.move_end_minus_acquire_s.min(), "to", kin.move_end_minus_acquire_s.max())
Kinematics of reward trials (positions in the file's unit; errors in degrees, 150 ms after movement onset at 20% of peak speed,
minus the session's mean per target):
session pos_unit pos_hz repeated_timestamps target_radius_cm cursor_radius_at_stop peak_speed cursor_enters_target stop_minus_go_s delay_s move_end_minus_acquire_s signed_dev_first8 signed_dev_last8 rho_trial_signed p_rho_signed abs_dev_first8 abs_dev_last8 max_8trial_mean median_abs_dev rho_trial_absdev p_rho
sub-C_ses-CO-20160909 cm 100 0 8.0 4.31 13.5 0.00 1.21 0.50-2.34 0.15 -3.2 5.4 0.21 0.010 4.1 5.9 10.4 5.1 -0.08 0.336
sub-C_ses-CO-20160912 cm 100 0 8.0 4.22 13.1 0.00 1.23 0.51-1.50 0.14 3.7 -1.7 -0.09 0.187 6.4 5.3 10.1 5.3 0.04 0.595
sub-C_ses-CO-20160914 cm 100 0 8.0 4.22 13.8 0.00 1.22 0.50-1.50 0.14 -6.4 -5.3 0.05 0.469 6.8 5.6 13.3 5.4 -0.01 0.922
sub-C_ses-CO-20160915 cm 100 0 8.0 4.22 14.5 0.00 1.18 0.50-1.50 0.11 2.3 0.8 -0.06 0.390 4.5 4.1 13.6 4.3 -0.08 0.261
sub-C_ses-CO-20160919 cm 100 0 8.0 4.28 14.0 0.00 1.19 0.51-1.50 0.13 -2.9 3.1 -0.02 0.808 9.5 6.5 9.5 5.4 0.00 0.983
sub-C_ses-CO-20160921 cm 100 0 8.0 4.26 13.3 0.00 1.19 0.50-1.49 0.13 0.2 -5.1 -0.14 0.056 6.7 9.4 11.5 5.0 -0.07 0.349
sub-C_ses-CO-20160923 cm 100 8531 8.0 4.21 13.4 0.00 1.18 0.50-1.50 0.13 12.5 -4.6 0.12 0.081 25.3 41.0 68.5 31.7 0.12 0.101
sub-C_ses-CO-20160929 cm 100 0 8.0 8.00 25.3 0.87 1.16 0.51-1.49 0.11 0.8 2.3 0.04 0.524 4.4 6.4 10.4 4.6 0.07 0.328
sub-C_ses-CO-20161005 cm 100 0 8.0 7.91 26.0 0.89 1.17 0.50-1.50 0.12 -0.8 5.2 0.22 0.001 9.0 6.8 11.9 4.4 0.00 0.972
sub-C_ses-CO-20161006 cm 100 0 8.0 7.91 25.8 0.89 1.16 0.50-1.49 0.10 -0.5 -2.3 0.08 0.246 3.9 6.5 9.6 4.7 0.03 0.680
sub-C_ses-CO-20161007 cm 100 0 8.0 7.94 25.1 0.87 1.18 0.50-1.50 0.10 4.2 0.2 -0.10 0.200 4.7 4.9 11.7 3.6 -0.12 0.126
sub-C_ses-CO-20161011 cm 100 5014 8.0 7.91 23.8 0.88 1.19 0.51-3.05 0.11 3.4 1.2 0.02 0.784 23.5 15.1 73.5 15.4 -0.21 0.003
sub-C_ses-CO-20161013 cm 100 0 8.0 7.91 25.3 0.89 1.18 0.50-1.50 0.09 1.9 -3.9 -0.05 0.418 4.7 6.1 10.1 4.4 -0.08 0.196
sub-C_ses-CO-20161021 cm 100 0 8.0 7.98 25.8 0.90 1.14 0.50-1.50 0.09 -0.1 -5.0 -0.08 0.177 4.7 5.9 9.1 4.1 0.04 0.542
sub-M_ses-CO-20140203 cm 100 0 8.0 7.81 21.0 0.01 1.36 0.51-1.49 0.09 -1.2 1.0 -0.02 0.818 9.1 5.6 9.8 4.6 -0.05 0.479
sub-M_ses-CO-20140217 cm 100 0 8.0 7.87 23.7 0.89 1.23 0.50-1.50 0.09 2.0 -0.9 -0.14 0.044 4.3 2.9 8.4 3.3 0.00 0.992
sub-M_ses-CO-20140218 cm 100 0 8.0 7.85 22.7 0.01 1.27 0.50-1.49 0.09 2.6 1.0 0.05 0.420 8.8 4.8 8.8 3.3 -0.09 0.198
sub-M_ses-CO-20140303 cm 100 0 8.0 7.89 23.4 0.89 1.24 0.50-1.50 0.09 10.2 3.3 0.16 0.019 14.0 6.4 14.0 3.4 0.02 0.821
sub-M_ses-CO-20140304 cm 100 0 8.0 7.79 20.6 0.89 1.34 0.50-1.50 0.09 2.6 0.3 0.02 0.756 4.7 4.4 9.6 4.4 -0.05 0.447
sub-M_ses-CO-20140306 cm 100 0 8.0 7.82 20.1 0.89 1.37 0.50-15.66 0.10 4.0 -1.1 NaN NaN 8.4 7.2 25.4 4.9 NaN NaN
sub-M_ses-CO-20140307 cm 100 0 8.0 7.76 23.0 0.88 1.27 0.50-1.50 0.09 3.8 -0.3 -0.13 0.055 6.2 4.8 8.7 4.1 0.00 0.960
End of movement minus assumed acquisition (s), median over sessions: 0.1 ; range 0.09 to 0.15
Cursor positions are stated in cm and sampled at 100 Hz. The instructed delay spans about 0.5 to 1.5 s, as the Methods say, with a few longer delays in C 2016-09-09, C 2016-10-11 and M 2014-03-06 (up to 15.7 s). The reach errors show no block of perturbed reaches of the size a curl field or rotation produces. The mean absolute deviation of the first eight reward trials is 3.9 to 14.0 degrees, against 2.9 to 9.4 for the last eight, outside two sessions. It is higher at the start by 3 to 8 degrees in M 2014-03-03 (14.0 against 6.4, signed 10.2 against 3.3, with a significant trend in the signed deviation, p = 0.019), M 2014-02-03 (9.1 against 5.6), M 2014-02-18 (8.8 against 4.8) and C 2016-09-19 (9.5 against 6.5), and large throughout in C 2016-09-23 (25.3 and 41.0) and C 2016-10-11 (23.5 and 15.1). The last two come with repeated cursor timestamps, which make the speed, and so the movement onset, undefined on some trials (the count is in the table). Their signed deviations are 12.5 and 3.4 degrees over the first eight trials and -4.6 and 1.2 over the last eight, with no significant trend in the signed deviation (Spearman rho 0.12, p = 0.081, and 0.02, p = 0.784); in C 2016-10-11 the absolute deviation decreases significantly over the session (rho -0.21, p = 0.003). Four sessions have a significant trend in the signed deviation (C 2016-09-09, C 2016-10-05, M 2014-02-17, M 2014-03-03, rho -0.14 to 0.22), with first and last eight-trial means within 10.2 degrees; rho is undefined for M 2014-03-06. None shows the large offset that decays over tens of trials that the first trials of a perturbation epoch would.
In the seven Monkey C sessions with the cut-off warning (2016-09-09 to 2016-09-23) the cursor ends
rewarded reaches about 4.2 cm from the center with a median peak speed near 14 cm/s, while the
target corners place every target at 8 cm and the later sessions end at about 7.9 cm with peak
speeds near 25 cm/s. Either the positions of those sessions are scaled differently or the target
corners are wrong; the paper's statement that Monkey C's first three sessions used 6 cm targets
does not fit either 4.2 or 8 cm. In the seven early Monkey C sessions the cursor never enters the square that
target_corners gives on a rewarded trial, and in M 2014-02-03 and M 2014-02-18 it enters it on 1%
of them (0.01 in the table). The analyses here use no positions; target acquisition is taken as the
trial's stop time minus the 0.5 s target hold (helper docstring). As a check of that choice, the
hand's speed last exceeds 20% of its peak a median of 0.09 to 0.15 s after this acquisition time
(per session), so the reach ends just after it, as when the cursor enters the target while still
slowing. This is consistent with a stop time at the end of the hold, but does not prove it.
def decode(sess, align, t0, t1, rng):
tr = sess["trials"]; r = tr[(tr.result == "R") & tr[align].notna()]
X = np.stack([np.searchsorted(sp, r[align].values + t1) - np.searchsorted(sp, r[align].values + t0)
for sp in sess["spikes"]], axis=1)
y = r.target_id.values
cv = StratifiedKFold(5, shuffle=True, random_state=0)
lda = LinearDiscriminantAnalysis(solver="lsqr", shrinkage="auto")
acc = cross_val_score(lda, X, y, cv=cv).mean()
acc_shuf = cross_val_score(lda, X, rng.permutation(y), cv=cv).mean()
return acc, acc_shuf
rng = np.random.default_rng(1)
WINDOWS = [("target_on", -0.3, 0.0), ("target_on", 0.1, 0.4), ("go_cue", 0.0, 0.3), ("go_cue", 0.3, 0.6),
("target_on", -1.3, -1.0)]
rows = []
for s in sessions:
row = {"session": s["sid"]}
for a, t0, t1 in WINDOWS:
acc, sh = decode(s, a, t0, t1, rng)
row[f"{a}[{t0},{t1}]"] = round(acc, 2)
row[f"shuffled {a}[{t0},{t1}]"] = round(sh, 2)
rows.append(row)
dec = pd.DataFrame(rows)
print("Cross-validated decoding of the target (8 targets, chance 0.125) from all units' spike counts:")
print(dec[["session"] + [c for c in dec.columns if c != "session" and not c.startswith("shuffled")]].to_string(index=False))
print("\nShuffled labels, mean over sessions per window:", dec[[c for c in dec.columns if c.startswith("shuffled")]].mean().round(3).to_dict())
/home/runner/work/_temp/executed/.venv/lib/python3.13/site-packages/sklearn/model_selection/_split.py:812: UserWarning: The least populated class in y has only 2 members, which is less than n_splits=5. warnings.warn( /home/runner/work/_temp/executed/.venv/lib/python3.13/site-packages/sklearn/covariance/_shrunk_covariance.py:349: UserWarning: Only one sample available. You may want to reshape your data array warnings.warn( /home/runner/work/_temp/executed/.venv/lib/python3.13/site-packages/sklearn/covariance/_empirical_covariance.py:100: UserWarning: Only one sample available. You may want to reshape your data array warnings.warn( /home/runner/work/_temp/executed/.venv/lib/python3.13/site-packages/sklearn/covariance/_shrunk_covariance.py:349: UserWarning: Only one sample available. You may want to reshape your data array warnings.warn( /home/runner/work/_temp/executed/.venv/lib/python3.13/site-packages/sklearn/covariance/_empirical_covariance.py:100: UserWarning: Only one sample available. You may want to reshape your data array warnings.warn( /home/runner/work/_temp/executed/.venv/lib/python3.13/site-packages/sklearn/model_selection/_split.py:812: UserWarning: The least populated class in y has only 2 members, which is less than n_splits=5. warnings.warn( /home/runner/work/_temp/executed/.venv/lib/python3.13/site-packages/sklearn/covariance/_shrunk_covariance.py:349: UserWarning: Only one sample available. You may want to reshape your data array warnings.warn( /home/runner/work/_temp/executed/.venv/lib/python3.13/site-packages/sklearn/covariance/_empirical_covariance.py:100: UserWarning: Only one sample available. You may want to reshape your data array warnings.warn( /home/runner/work/_temp/executed/.venv/lib/python3.13/site-packages/sklearn/covariance/_shrunk_covariance.py:349: UserWarning: Only one sample available. You may want to reshape your data array warnings.warn( /home/runner/work/_temp/executed/.venv/lib/python3.13/site-packages/sklearn/covariance/_empirical_covariance.py:100: UserWarning: Only one sample available. You may want to reshape your data array warnings.warn( /home/runner/work/_temp/executed/.venv/lib/python3.13/site-packages/sklearn/model_selection/_split.py:812: UserWarning: The least populated class in y has only 2 members, which is less than n_splits=5. warnings.warn( /home/runner/work/_temp/executed/.venv/lib/python3.13/site-packages/sklearn/covariance/_shrunk_covariance.py:349: UserWarning: Only one sample available. You may want to reshape your data array warnings.warn( /home/runner/work/_temp/executed/.venv/lib/python3.13/site-packages/sklearn/covariance/_empirical_covariance.py:100: UserWarning: Only one sample available. You may want to reshape your data array warnings.warn( /home/runner/work/_temp/executed/.venv/lib/python3.13/site-packages/sklearn/covariance/_shrunk_covariance.py:349: UserWarning: Only one sample available. You may want to reshape your data array warnings.warn( /home/runner/work/_temp/executed/.venv/lib/python3.13/site-packages/sklearn/covariance/_empirical_covariance.py:100: UserWarning: Only one sample available. You may want to reshape your data array warnings.warn( /home/runner/work/_temp/executed/.venv/lib/python3.13/site-packages/sklearn/model_selection/_split.py:812: UserWarning: The least populated class in y has only 2 members, which is less than n_splits=5. warnings.warn( /home/runner/work/_temp/executed/.venv/lib/python3.13/site-packages/sklearn/covariance/_shrunk_covariance.py:349: UserWarning: Only one sample available. You may want to reshape your data array warnings.warn( /home/runner/work/_temp/executed/.venv/lib/python3.13/site-packages/sklearn/covariance/_empirical_covariance.py:100: UserWarning: Only one sample available. You may want to reshape your data array warnings.warn( /home/runner/work/_temp/executed/.venv/lib/python3.13/site-packages/sklearn/covariance/_shrunk_covariance.py:349: UserWarning: Only one sample available. You may want to reshape your data array warnings.warn( /home/runner/work/_temp/executed/.venv/lib/python3.13/site-packages/sklearn/covariance/_empirical_covariance.py:100: UserWarning: Only one sample available. You may want to reshape your data array warnings.warn( /home/runner/work/_temp/executed/.venv/lib/python3.13/site-packages/sklearn/model_selection/_split.py:812: UserWarning: The least populated class in y has only 2 members, which is less than n_splits=5. warnings.warn( /home/runner/work/_temp/executed/.venv/lib/python3.13/site-packages/sklearn/covariance/_shrunk_covariance.py:349: UserWarning: Only one sample available. You may want to reshape your data array warnings.warn( /home/runner/work/_temp/executed/.venv/lib/python3.13/site-packages/sklearn/covariance/_empirical_covariance.py:100: UserWarning: Only one sample available. You may want to reshape your data array warnings.warn( /home/runner/work/_temp/executed/.venv/lib/python3.13/site-packages/sklearn/covariance/_shrunk_covariance.py:349: UserWarning: Only one sample available. You may want to reshape your data array warnings.warn( /home/runner/work/_temp/executed/.venv/lib/python3.13/site-packages/sklearn/covariance/_empirical_covariance.py:100: UserWarning: Only one sample available. You may want to reshape your data array warnings.warn(
/home/runner/work/_temp/executed/.venv/lib/python3.13/site-packages/sklearn/model_selection/_split.py:812: UserWarning: The least populated class in y has only 2 members, which is less than n_splits=5. warnings.warn( /home/runner/work/_temp/executed/.venv/lib/python3.13/site-packages/sklearn/covariance/_shrunk_covariance.py:349: UserWarning: Only one sample available. You may want to reshape your data array warnings.warn( /home/runner/work/_temp/executed/.venv/lib/python3.13/site-packages/sklearn/covariance/_empirical_covariance.py:100: UserWarning: Only one sample available. You may want to reshape your data array warnings.warn( /home/runner/work/_temp/executed/.venv/lib/python3.13/site-packages/sklearn/covariance/_shrunk_covariance.py:349: UserWarning: Only one sample available. You may want to reshape your data array warnings.warn( /home/runner/work/_temp/executed/.venv/lib/python3.13/site-packages/sklearn/covariance/_empirical_covariance.py:100: UserWarning: Only one sample available. You may want to reshape your data array warnings.warn( /home/runner/work/_temp/executed/.venv/lib/python3.13/site-packages/sklearn/model_selection/_split.py:812: UserWarning: The least populated class in y has only 2 members, which is less than n_splits=5. warnings.warn( /home/runner/work/_temp/executed/.venv/lib/python3.13/site-packages/sklearn/covariance/_shrunk_covariance.py:349: UserWarning: Only one sample available. You may want to reshape your data array warnings.warn( /home/runner/work/_temp/executed/.venv/lib/python3.13/site-packages/sklearn/covariance/_empirical_covariance.py:100: UserWarning: Only one sample available. You may want to reshape your data array warnings.warn( /home/runner/work/_temp/executed/.venv/lib/python3.13/site-packages/sklearn/covariance/_shrunk_covariance.py:349: UserWarning: Only one sample available. You may want to reshape your data array warnings.warn( /home/runner/work/_temp/executed/.venv/lib/python3.13/site-packages/sklearn/covariance/_empirical_covariance.py:100: UserWarning: Only one sample available. You may want to reshape your data array warnings.warn( /home/runner/work/_temp/executed/.venv/lib/python3.13/site-packages/sklearn/model_selection/_split.py:812: UserWarning: The least populated class in y has only 2 members, which is less than n_splits=5. warnings.warn( /home/runner/work/_temp/executed/.venv/lib/python3.13/site-packages/sklearn/covariance/_shrunk_covariance.py:349: UserWarning: Only one sample available. You may want to reshape your data array warnings.warn( /home/runner/work/_temp/executed/.venv/lib/python3.13/site-packages/sklearn/covariance/_empirical_covariance.py:100: UserWarning: Only one sample available. You may want to reshape your data array warnings.warn( /home/runner/work/_temp/executed/.venv/lib/python3.13/site-packages/sklearn/covariance/_shrunk_covariance.py:349: UserWarning: Only one sample available. You may want to reshape your data array warnings.warn( /home/runner/work/_temp/executed/.venv/lib/python3.13/site-packages/sklearn/covariance/_empirical_covariance.py:100: UserWarning: Only one sample available. You may want to reshape your data array warnings.warn( /home/runner/work/_temp/executed/.venv/lib/python3.13/site-packages/sklearn/model_selection/_split.py:812: UserWarning: The least populated class in y has only 2 members, which is less than n_splits=5. warnings.warn( /home/runner/work/_temp/executed/.venv/lib/python3.13/site-packages/sklearn/covariance/_shrunk_covariance.py:349: UserWarning: Only one sample available. You may want to reshape your data array warnings.warn( /home/runner/work/_temp/executed/.venv/lib/python3.13/site-packages/sklearn/covariance/_empirical_covariance.py:100: UserWarning: Only one sample available. You may want to reshape your data array warnings.warn( /home/runner/work/_temp/executed/.venv/lib/python3.13/site-packages/sklearn/covariance/_shrunk_covariance.py:349: UserWarning: Only one sample available. You may want to reshape your data array warnings.warn( /home/runner/work/_temp/executed/.venv/lib/python3.13/site-packages/sklearn/covariance/_empirical_covariance.py:100: UserWarning: Only one sample available. You may want to reshape your data array warnings.warn( /home/runner/work/_temp/executed/.venv/lib/python3.13/site-packages/sklearn/model_selection/_split.py:812: UserWarning: The least populated class in y has only 2 members, which is less than n_splits=5. warnings.warn( /home/runner/work/_temp/executed/.venv/lib/python3.13/site-packages/sklearn/covariance/_shrunk_covariance.py:349: UserWarning: Only one sample available. You may want to reshape your data array warnings.warn( /home/runner/work/_temp/executed/.venv/lib/python3.13/site-packages/sklearn/covariance/_empirical_covariance.py:100: UserWarning: Only one sample available. You may want to reshape your data array warnings.warn( /home/runner/work/_temp/executed/.venv/lib/python3.13/site-packages/sklearn/covariance/_shrunk_covariance.py:349: UserWarning: Only one sample available. You may want to reshape your data array warnings.warn( /home/runner/work/_temp/executed/.venv/lib/python3.13/site-packages/sklearn/covariance/_empirical_covariance.py:100: UserWarning: Only one sample available. You may want to reshape your data array warnings.warn(
Cross-validated decoding of the target (8 targets, chance 0.125) from all units' spike counts:
session target_on[-0.3,0.0] target_on[0.1,0.4] go_cue[0.0,0.3] go_cue[0.3,0.6] target_on[-1.3,-1.0]
sub-C_ses-CO-20160909 0.10 0.75 0.86 0.97 0.14
sub-C_ses-CO-20160912 0.12 0.73 0.94 1.00 0.15
sub-C_ses-CO-20160914 0.17 0.74 0.89 0.99 0.17
sub-C_ses-CO-20160915 0.08 0.79 0.95 1.00 0.12
sub-C_ses-CO-20160919 0.15 0.73 0.93 0.99 0.10
sub-C_ses-CO-20160921 0.15 0.81 0.93 1.00 0.12
sub-C_ses-CO-20160923 0.12 0.78 0.95 1.00 0.11
sub-C_ses-CO-20160929 0.09 0.74 0.95 1.00 0.10
sub-C_ses-CO-20161005 0.12 0.67 0.97 0.99 0.12
sub-C_ses-CO-20161006 0.13 0.69 0.93 0.99 0.13
sub-C_ses-CO-20161007 0.12 0.55 0.91 0.98 0.11
sub-C_ses-CO-20161011 0.09 0.64 0.94 0.99 0.12
sub-C_ses-CO-20161013 0.11 0.61 0.94 1.00 0.10
sub-C_ses-CO-20161021 0.07 0.71 0.95 1.00 0.11
sub-M_ses-CO-20140203 0.09 0.66 0.80 0.93 0.18
sub-M_ses-CO-20140217 0.13 0.75 0.91 0.96 0.12
sub-M_ses-CO-20140218 0.14 0.70 0.86 0.96 0.12
sub-M_ses-CO-20140303 0.12 0.74 0.93 0.98 0.12
sub-M_ses-CO-20140304 0.15 0.73 0.85 0.96 0.11
sub-M_ses-CO-20140306 0.16 0.55 0.86 0.96 0.09
sub-M_ses-CO-20140307 0.15 0.70 0.88 0.96 0.12
Shuffled labels, mean over sessions per window: {'shuffled target_on[-0.3,0.0]': 0.131, 'shuffled target_on[0.1,0.4]': 0.128, 'shuffled go_cue[0.0,0.3]': 0.123, 'shuffled go_cue[0.3,0.6]': 0.126, 'shuffled target_on[-1.3,-1.0]': 0.133}
In every session the target cannot be decoded before target onset (windows ending at onset) and can be decoded after it, and best around the go cue, while shuffled labels stay at chance. The spikes therefore share the clock of the trial table.
methods-neuron-counts¶
Per session, 137 to 256 PMd and 55 to 93 M1 neurons were isolated in Monkey C and 66 to 121 PMd and 26 to 51 M1 neurons in Monkey M; neurons above 5 Hz, pooled over all CF and VR sessions, number 918 in M1 and 2221 in PMd.
Criterion: units as sorted in the data, per area and session; each count within its monkey's range widened by 10% of each bound; the 5 Hz rule is the trial-averaged mean rate over all trials of the session. The files hold only Baseline trials, so the rate is taken over the trials there: the reward trials from start to stop (main), and every trial of the file (variant). The paper's 16 sessions cannot be identified, so the two session sets of the data check are both shown.
rows = []
for s in sessions:
tr = s["trials"]
rate_r = perich.mean_rates(s["spikes"], tr[tr.result == "R"])
rate_all = perich.mean_rates(s["spikes"], tr)
for area in ("M1", "PMd"):
m = s["area"] == area
rows.append({"session": s["sid"], "monkey": s["monkey"], "area": area, "n": int(m.sum()),
"above5_reward": int((rate_r[m] > 5).sum()), "above5_all": int((rate_all[m] > 5).sum()),
"above1_reward": int((rate_r[m] > 1).sum())})
counts = pd.DataFrame(rows)
wide = counts.pivot(index="session", columns="area", values=["n", "above5_reward", "above5_all"])
print(wide.to_string())
def within(n, lo, hi):
return (lo - 0.1 * lo) <= n <= (hi + 0.1 * hi)
summary = []
for name, sids in SETS.items():
c = counts[counts.session.isin(sids)]
for mk in ("C", "M"):
for area in ("PMd", "M1"):
v = c[(c.monkey == mk) & (c.area == area)]
lo, hi = PAPER[mk][area]
out = v[~v.n.apply(lambda n: within(n, lo, hi))]
summary.append({"set": name, "monkey": mk, "area": area, "paper": f"{lo}-{hi}",
"data": f"{v.n.min()}-{v.n.max()}", "sessions": len(v),
"outside widened range": ", ".join(f"{a} ({b})" for a, b in zip(out.session, out.n)) or "none"})
for area, exp in (("M1", 918), ("PMd", 2221)):
c2 = c[c.area == area]
summary.append({"set": name, "monkey": "pooled", "area": area, "paper": f"{exp} above 5 Hz",
"data": f"{c2.above5_reward.sum()} (reward trials), {c2.above5_all.sum()} (all trials)",
"sessions": c2.session.nunique(),
"outside widened range": ("within 10%" if abs(c2.above5_reward.sum() - exp) <= 0.1 * exp else "outside 10%")
+ " (this set, not 16 sessions; see next cell)"})
print()
print(pd.DataFrame(summary).to_string(index=False))
n above5_reward above5_all
area M1 PMd M1 PMd M1 PMd
session
sub-C_ses-CO-20160909 55 191 28 80 28 81
sub-C_ses-CO-20160912 79 162 52 100 51 99
sub-C_ses-CO-20160914 95 258 61 145 60 145
sub-C_ses-CO-20160915 76 233 61 154 59 153
sub-C_ses-CO-20160919 60 235 40 154 40 155
sub-C_ses-CO-20160921 74 232 45 134 45 134
sub-C_ses-CO-20160923 76 223 44 134 44 134
sub-C_ses-CO-20160929 74 114 44 55 44 54
sub-C_ses-CO-20161005 82 167 58 102 58 106
sub-C_ses-CO-20161006 63 192 47 102 48 102
sub-C_ses-CO-20161007 70 137 47 79 47 81
sub-C_ses-CO-20161011 71 171 50 70 51 71
sub-C_ses-CO-20161013 87 138 57 57 57 57
sub-C_ses-CO-20161021 84 211 63 108 63 108
sub-M_ses-CO-20140203 34 80 14 6 15 6
sub-M_ses-CO-20140217 44 104 16 8 16 9
sub-M_ses-CO-20140218 38 121 15 15 15 16
sub-M_ses-CO-20140303 52 66 28 8 28 8
sub-M_ses-CO-20140304 39 76 28 18 28 18
sub-M_ses-CO-20140306 43 86 22 21 22 22
sub-M_ses-CO-20140307 26 66 15 15 15 16
set monkey area paper data sessions outside widened range
all 21 C PMd 137-256 114-258 14 sub-C_ses-CO-20160929 (114)
all 21 C M1 55-93 55-95 14 none
all 21 M PMd 66-121 66-121 7 none
all 21 M M1 26-51 26-52 7 none
all 21 pooled M1 918 above 5 Hz 835 (reward trials), 834 (all trials) 21 within 10% (this set, not 16 sessions; see next cell)
all 21 pooled PMd 2221 above 5 Hz 1565 (reward trials), 1575 (all trials) 21 outside 10% (this set, not 16 sessions; see next cell)
19 within Baseline reward range C PMd 137-256 114-258 12 sub-C_ses-CO-20160929 (114)
19 within Baseline reward range C M1 55-93 55-95 12 none
19 within Baseline reward range M PMd 66-121 66-121 7 none
19 within Baseline reward range M M1 26-51 26-52 7 none
19 within Baseline reward range pooled M1 918 above 5 Hz 715 (reward trials), 714 (all trials) 19 outside 10% (this set, not 16 sessions; see next cell)
19 within Baseline reward range pooled PMd 2221 above 5 Hz 1400 (reward trials), 1410 (all trials) 19 outside 10% (this set, not 16 sessions; see next cell)
# How far the pooled counts depend on which 9 of the Monkey C candidates are the paper's: every choice
# of 9 Monkey C sessions from each set, with the 7 Monkey M sessions.
for name, sids in SETS.items():
cs = [x for x in sids if x[4] == "C"]; ms = [x for x in sids if x[4] == "M"]
a5 = counts.set_index(["session", "area"]).above5_reward
res = []
for combo in itertools.combinations(cs, 9):
ss_ = list(combo) + ms
ok = all(within(int(counts[(counts.session == x) & (counts.area == ar)].n.iloc[0]), *PAPER[x[4]][ar])
for x in ss_ for ar in ("M1", "PMd"))
res.append((sum(a5[(x, "M1")] for x in ss_), sum(a5[(x, "PMd")] for x in ss_), ok))
res = pd.DataFrame(res, columns=["M1", "PMd", "ranges_hold"])
both = res[(abs(res.M1 - 918) <= 91.8) & (abs(res.PMd - 2221) <= 222.1)]
print(f"{name}: {len(res)} choices of 9 C sessions; pooled M1 {res.M1.min()}-{res.M1.max()}, PMd {res.PMd.min()}-{res.PMd.max()}; "
f"ranges hold in {int(res.ranges_hold.sum())}; both pooled counts within 10% in {len(both)}, "
f"of which ranges also hold in {int(both.ranges_hold.sum())}")
all 21: 2002 choices of 9 C sessions; pooled M1 535-634, PMd 844-1224; ranges hold in 715; both pooled counts within 10% in 0, of which ranges also hold in 0 19 within Baseline reward range: 220 choices of 9 C sessions; pooled M1 535-603, PMd 947-1196; ranges hold in 55; both pooled counts within 10% in 0, of which ranges also hold in 0
Result: partially reproduced, resting on an assumption about which sessions are the paper's. Monkey M's seven sessions give 66 to 121 PMd and 26 to 52 M1 units, the paper's ranges (52 is within 10% of 51). Thirteen of the 14 Monkey C candidates lie inside the widened ranges (PMd 137 to 258, M1 55 to 95); C 2016-09-29 has 114 PMd units, below 137 - 13.7. The ranges hold for 715 of the 2002 ways to choose 9 Monkey C sessions from all 21 candidates, those that leave out C 2016-09-29. The pooled counts above 5 Hz do not hold under any choice of 9 Monkey C sessions with the 7 Monkey M sessions: at most 634 M1 and 1224 PMd units (paper 918 and 2221). Even all 21 sessions give 835 M1 (within 10% of 918) and 1565 PMd (outside). Rates over every trial of the file give the same counts within 10 units. Monkey M has few PMd units above 5 Hz (6 to 21 per session). The pooled verdicts in the first table are for each whole set, not for 16 sessions; the second cell gives the test on 16 sessions.
fig3b-pmd-higher-dimensional¶
In every session the PMd manifold had more dimensions than the M1 manifold (median 14 against 8.5).
Method from STAR Methods, Dimensionality reduction, paragraph 2 (helper dimensionality): Baseline
reward trials, rates as M2, trial-averaged per target; noise from differences between random pairs
of trials to the same target, 1000 repetitions, the 99% limit of each noise eigenvalue. Readings:
threshold rule (a) or (b) (O9), window (a) target onset to acquisition or (b) the whole trial (O3),
kernel width as a standard deviation or a full width at half maximum (O4), and the scale of the
noise, which the paper does not give (raw differences, the text's wording; difference / sqrt 2;
difference / sqrt(2K), the noise of an average of K trials). All units of each area are used.
The claims file expected the kernel reading (O4) not to change an outcome; here it does, and the
results are shown under both. The paper's 16 sessions are 9 of the 14 Monkey C candidates and the 7
Monkey M sessions; the medians are also computed for every such choice.
t0 = time.time()
def dims_for(i):
s = sessions[i]
rng = np.random.default_rng([3, i]) # one generator per session, independent of thread order
out_fits, out_rows = {}, []
for window in ("a", "b"):
for kern in ("sd", "fwhm"):
fit = perich.latent_fit(s, window, perich.KERNEL_SD_S[kern])
out_fits[(s["sid"], window, kern)] = fit
tids = fit["trials"].target_id.to_numpy()
for area in ("M1", "PMd"):
d = perich.dimensionality(fit[area]["rates"], tids, rng)
for scale, (da, db) in d.items():
out_rows.append({"session": s["sid"], "monkey": s["monkey"], "window": window, "kernel": kern,
"scale": scale, "area": area, "rule_a": da, "rule_b": db})
return out_fits, out_rows
fits, rows = {}, []
with ThreadPoolExecutor(4) as ex:
for f_, r_ in ex.map(dims_for, range(len(sessions))):
fits.update(f_); rows += r_
dims = pd.DataFrame(rows)
print(f"{len(dims)} dimensionality estimates in {time.time() - t0:.0f} s")
long = dims.melt(id_vars=["session", "monkey", "window", "kernel", "scale", "area"], value_vars=["rule_a", "rule_b"],
var_name="rule", value_name="dim")
piv = long.pivot_table(index=["rule", "scale", "window", "kernel", "session"], columns="area", values="dim").reset_index()
rows = []
for name, sids in SETS.items():
for key, g in piv[piv.session.isin(sids)].groupby(["rule", "scale", "window", "kernel"]):
rows.append({"set": name, "rule": key[0], "scale": key[1], "window": key[2], "kernel": key[3],
"M1_median": g.M1.median(), "PMd_median": g.PMd.median(),
"M1_range": f"{g.M1.min()}-{g.M1.max()}", "PMd_range": f"{g.PMd.min()}-{g.PMd.max()}",
"PMd>M1": f"{int((g.PMd > g.M1).sum())}/{len(g)}",
"medians_in_tol": bool(abs(g.M1.median() - 8.5) <= 0.85 and abs(g.PMd.median() - 14) <= 1.4)})
dimsum = pd.DataFrame(rows)
print(dimsum.to_string(index=False))
504 dimensionality estimates in 51 s
set rule scale window kernel M1_median PMd_median M1_range PMd_range PMd>M1 medians_in_tol
all 21 rule_a raw a fwhm 0.0 0.0 0.0-0.0 0.0-0.0 0/21 False
all 21 rule_a raw a sd 0.0 0.0 0.0-0.0 0.0-0.0 0/21 False
all 21 rule_a raw b fwhm 0.0 0.0 0.0-0.0 0.0-0.0 0/21 False
all 21 rule_a raw b sd 0.0 0.0 0.0-0.0 0.0-0.0 0/21 False
all 21 rule_a sqrt2 a fwhm 0.0 0.0 0.0-1.0 0.0-1.0 1/21 False
all 21 rule_a sqrt2 a sd 0.0 0.0 0.0-1.0 0.0-1.0 1/21 False
all 21 rule_a sqrt2 b fwhm 0.0 0.0 0.0-0.0 0.0-0.0 0/21 False
all 21 rule_a sqrt2 b sd 0.0 0.0 0.0-0.0 0.0-0.0 0/21 False
all 21 rule_a sqrt2K a fwhm 47.0 78.0 18.0-69.0 46.0-94.0 21/21 False
all 21 rule_a sqrt2K a sd 38.0 53.0 21.0-56.0 45.0-64.0 19/21 False
all 21 rule_a sqrt2K b fwhm 29.0 49.0 9.0-52.0 14.0-99.0 19/21 False
all 21 rule_a sqrt2K b sd 26.0 46.0 10.0-48.0 17.0-71.0 18/21 False
all 21 rule_b raw a fwhm 0.0 0.0 0.0-0.0 0.0-0.0 0/21 False
all 21 rule_b raw a sd 0.0 0.0 0.0-0.0 0.0-0.0 0/21 False
all 21 rule_b raw b fwhm 0.0 0.0 0.0-0.0 0.0-0.0 0/21 False
all 21 rule_b raw b sd 0.0 0.0 0.0-0.0 0.0-0.0 0/21 False
all 21 rule_b sqrt2 a fwhm 0.0 0.0 0.0-1.0 0.0-1.0 1/21 False
all 21 rule_b sqrt2 a sd 0.0 0.0 0.0-1.0 0.0-1.0 1/21 False
all 21 rule_b sqrt2 b fwhm 0.0 0.0 0.0-0.0 0.0-0.0 0/21 False
all 21 rule_b sqrt2 b sd 0.0 0.0 0.0-0.0 0.0-0.0 0/21 False
all 21 rule_b sqrt2K a fwhm 9.0 13.0 6.0-14.0 8.0-19.0 21/21 True
all 21 rule_b sqrt2K a sd 7.0 10.0 5.0-10.0 8.0-13.0 21/21 False
all 21 rule_b sqrt2K b fwhm 6.0 8.0 4.0-22.0 7.0-29.0 20/21 False
all 21 rule_b sqrt2K b sd 6.0 7.0 4.0-13.0 6.0-16.0 20/21 False
19 within Baseline reward range rule_a raw a fwhm 0.0 0.0 0.0-0.0 0.0-0.0 0/19 False
19 within Baseline reward range rule_a raw a sd 0.0 0.0 0.0-0.0 0.0-0.0 0/19 False
19 within Baseline reward range rule_a raw b fwhm 0.0 0.0 0.0-0.0 0.0-0.0 0/19 False
19 within Baseline reward range rule_a raw b sd 0.0 0.0 0.0-0.0 0.0-0.0 0/19 False
19 within Baseline reward range rule_a sqrt2 a fwhm 0.0 0.0 0.0-1.0 0.0-1.0 1/19 False
19 within Baseline reward range rule_a sqrt2 a sd 0.0 0.0 0.0-1.0 0.0-1.0 1/19 False
19 within Baseline reward range rule_a sqrt2 b fwhm 0.0 0.0 0.0-0.0 0.0-0.0 0/19 False
19 within Baseline reward range rule_a sqrt2 b sd 0.0 0.0 0.0-0.0 0.0-0.0 0/19 False
19 within Baseline reward range rule_a sqrt2K a fwhm 46.0 78.0 18.0-66.0 46.0-94.0 19/19 False
19 within Baseline reward range rule_a sqrt2K a sd 37.0 53.0 21.0-50.0 45.0-64.0 18/19 False
19 within Baseline reward range rule_a sqrt2K b fwhm 29.0 45.0 9.0-52.0 14.0-99.0 17/19 False
19 within Baseline reward range rule_a sqrt2K b sd 26.0 44.0 10.0-42.0 17.0-71.0 16/19 False
19 within Baseline reward range rule_b raw a fwhm 0.0 0.0 0.0-0.0 0.0-0.0 0/19 False
19 within Baseline reward range rule_b raw a sd 0.0 0.0 0.0-0.0 0.0-0.0 0/19 False
19 within Baseline reward range rule_b raw b fwhm 0.0 0.0 0.0-0.0 0.0-0.0 0/19 False
19 within Baseline reward range rule_b raw b sd 0.0 0.0 0.0-0.0 0.0-0.0 0/19 False
19 within Baseline reward range rule_b sqrt2 a fwhm 0.0 0.0 0.0-1.0 0.0-1.0 1/19 False
19 within Baseline reward range rule_b sqrt2 a sd 0.0 0.0 0.0-1.0 0.0-1.0 1/19 False
19 within Baseline reward range rule_b sqrt2 b fwhm 0.0 0.0 0.0-0.0 0.0-0.0 0/19 False
19 within Baseline reward range rule_b sqrt2 b sd 0.0 0.0 0.0-0.0 0.0-0.0 0/19 False
19 within Baseline reward range rule_b sqrt2K a fwhm 9.0 13.0 6.0-14.0 8.0-19.0 19/19 True
19 within Baseline reward range rule_b sqrt2K a sd 7.0 11.0 5.0-10.0 8.0-13.0 19/19 False
19 within Baseline reward range rule_b sqrt2K b fwhm 6.0 8.0 4.0-22.0 7.0-29.0 18/19 False
19 within Baseline reward range rule_b sqrt2K b sd 6.0 8.0 4.0-13.0 6.0-16.0 18/19 False
main = piv[(piv.rule == "rule_b") & (piv.scale == "sqrt2K")]
print("Rule (b), noise of the trial average: per session, window (a) and (b), standard-deviation kernel")
show = main[main.kernel == "sd"].pivot(index="session", columns="window", values=["M1", "PMd"])
print(show.to_string())
print("\nRule (b), noise of the trial average: per session, window (a) and (b), full-width kernel")
print(main[main.kernel == "fwhm"].pivot(index="session", columns="window", values=["M1", "PMd"]).to_string())
# The paper used 9 Monkey C sessions and the 7 Monkey M sessions: every choice of 9 of the 14 C candidates
c_sids = [x for x in SETS["all 21"] if x[4] == "C"]; m_sids = [x for x in SETS["all 21"] if x[4] == "M"]
CHOICES = [list(c) + m_sids for c in itertools.combinations(c_sids, 9)]
rows = []
for (w, k), g in main.groupby(["window", "kernel"]):
g = g.set_index("session")
med = np.array([[g.loc[ss, "M1"].median(), g.loc[ss, "PMd"].median()] for ss in CHOICES])
every = np.array([bool((g.loc[ss, "PMd"] > g.loc[ss, "M1"]).all()) for ss in CHOICES])
tol = (np.abs(med[:, 0] - 8.5) <= 0.85) & (np.abs(med[:, 1] - 14) <= 1.4)
rows.append({"window": w, "kernel": k, "choices": len(CHOICES), "M1_median_range": f"{med[:, 0].min()}-{med[:, 0].max()}",
"PMd_median_range": f"{med[:, 1].min()}-{med[:, 1].max()}", "PMd>M1_in_every_session": int(every.sum()),
"medians_in_tol": int(tol.sum()), "both": int((tol & every).sum())})
print("\nRule (b), noise of the trial average, over the", len(CHOICES), "choices of 9 Monkey C sessions with the 7 Monkey M sessions:")
print(pd.DataFrame(rows).to_string(index=False))
for mk in ("C", "M"):
g = main[(main.kernel == "sd") & (main.session.str[4] == mk)]
for w in ("a", "b"):
gw = g[g.window == w]
print(f"Monkey {mk}, window {w}: M1 median {gw.M1.median()}, PMd median {gw.PMd.median()}, PMd > M1 in {int((gw.PMd > gw.M1).sum())}/{len(gw)}")
fig, axes = plt.subplots(1, 2, figsize=(9, 3), sharey=True)
for ax, w in zip(axes, ("a", "b")):
g = main[(main.kernel == "sd") & (main.window == w)]
bins = np.arange(0, max(g.PMd.max(), g.M1.max()) + 2) - 0.5
ax.hist(g.M1, bins, alpha=0.6, label="M1", color="k"); ax.hist(g.PMd, bins, alpha=0.6, label="PMd", color="r")
ax.set_title(f"window ({w}), rule (b), noise of the average"); ax.set_xlabel("dimensionality"); ax.legend()
axes[0].set_ylabel("sessions")
plt.tight_layout(); plt.show()
print("All values are inside the axes.")
Rule (b), noise of the trial average: per session, window (a) and (b), standard-deviation kernel
M1 PMd
window a b a b
session
sub-C_ses-CO-20160909 5.0 5.0 12.0 9.0
sub-C_ses-CO-20160912 6.0 6.0 10.0 7.0
sub-C_ses-CO-20160914 6.0 5.0 12.0 8.0
sub-C_ses-CO-20160915 6.0 4.0 11.0 7.0
sub-C_ses-CO-20160919 7.0 5.0 10.0 7.0
sub-C_ses-CO-20160921 7.0 5.0 11.0 8.0
sub-C_ses-CO-20160923 7.0 5.0 11.0 7.0
sub-C_ses-CO-20160929 8.0 6.0 10.0 7.0
sub-C_ses-CO-20161005 8.0 6.0 10.0 7.0
sub-C_ses-CO-20161006 8.0 6.0 10.0 7.0
sub-C_ses-CO-20161007 8.0 7.0 11.0 8.0
sub-C_ses-CO-20161011 8.0 6.0 10.0 6.0
sub-C_ses-CO-20161013 7.0 6.0 10.0 7.0
sub-C_ses-CO-20161021 8.0 5.0 10.0 6.0
sub-M_ses-CO-20140203 6.0 5.0 8.0 7.0
sub-M_ses-CO-20140217 7.0 7.0 11.0 9.0
sub-M_ses-CO-20140218 6.0 5.0 12.0 11.0
sub-M_ses-CO-20140303 8.0 7.0 9.0 9.0
sub-M_ses-CO-20140304 7.0 6.0 11.0 8.0
sub-M_ses-CO-20140306 10.0 13.0 13.0 16.0
sub-M_ses-CO-20140307 6.0 4.0 10.0 9.0
Rule (b), noise of the trial average: per session, window (a) and (b), full-width kernel
M1 PMd
window a b a b
session
sub-C_ses-CO-20160909 6.0 5.0 16.0 12.0
sub-C_ses-CO-20160912 7.0 6.0 14.0 8.0
sub-C_ses-CO-20160914 8.0 6.0 15.0 10.0
sub-C_ses-CO-20160915 7.0 4.0 14.0 8.0
sub-C_ses-CO-20160919 9.0 5.0 12.0 7.0
sub-C_ses-CO-20160921 9.0 6.0 15.0 10.0
sub-C_ses-CO-20160923 9.0 5.0 13.0 8.0
sub-C_ses-CO-20160929 10.0 7.0 12.0 8.0
sub-C_ses-CO-20161005 9.0 7.0 12.0 7.0
sub-C_ses-CO-20161006 9.0 6.0 13.0 9.0
sub-C_ses-CO-20161007 9.0 7.0 12.0 8.0
sub-C_ses-CO-20161011 9.0 6.0 13.0 7.0
sub-C_ses-CO-20161013 9.0 6.0 13.0 8.0
sub-C_ses-CO-20161021 10.0 6.0 13.0 7.0
sub-M_ses-CO-20140203 6.0 5.0 8.0 8.0
sub-M_ses-CO-20140217 8.0 7.0 12.0 9.0
sub-M_ses-CO-20140218 7.0 5.0 14.0 11.0
sub-M_ses-CO-20140303 9.0 7.0 10.0 10.0
sub-M_ses-CO-20140304 8.0 5.0 12.0 8.0
sub-M_ses-CO-20140306 14.0 22.0 19.0 29.0
sub-M_ses-CO-20140307 7.0 4.0 12.0 10.0
Rule (b), noise of the trial average, over the 2002 choices of 9 Monkey C sessions with the 7 Monkey M sessions:
window kernel choices M1_median_range PMd_median_range PMd>M1_in_every_session medians_in_tol both
a fwhm 2002 8.0-9.0 12.0-13.0 2002 910 910
a sd 2002 7.0-7.5 10.0-11.0 2002 0 0
b fwhm 2002 5.5-6.0 8.0-9.0 715 0 0
b sd 2002 5.0-6.0 7.0-8.0 715 0 0
Monkey C, window a: M1 median 7.0, PMd median 10.0, PMd > M1 in 14/14
Monkey C, window b: M1 median 5.5, PMd median 7.0, PMd > M1 in 13/14
Monkey M, window a: M1 median 7.0, PMd median 11.0, PMd > M1 in 7/7
Monkey M, window b: M1 median 6.0, PMd median 9.0, PMd > M1 in 7/7
All values are inside the axes.
# Check of `dimensionality` on synthetic data of known rank: 8 targets, 20 bins, 60 neurons, a rank-5
# trial-averaged signal and independent single-trial noise, 25 trials per target.
rng = np.random.default_rng(31)
rows = []
for noise_sd in (0.5, 1.0, 2.0):
for rep_ in range(3):
sig = rng.normal(size=(8 * 20, 5)) @ rng.normal(size=(5, 60))
sig = sig.reshape(8, 20, 60)
trials, tids = [], []
for tg in range(8):
for _ in range(25):
trials.append(sig[tg] + noise_sd * rng.normal(size=(20, 60))); tids.append(tg)
d = perich.dimensionality(trials, tids, rng)
rows.append({"noise_sd": noise_sd, "rep": rep_, **{f"{k}_rule_{r}": v[j] for k, v in d.items() for j, r in enumerate("ab")}})
print("Estimated dimensionality of rank-5 synthetic data, by scale of the noise and rule:")
print(pd.DataFrame(rows).to_string(index=False))
Estimated dimensionality of rank-5 synthetic data, by scale of the noise and rule:
noise_sd rep raw_rule_a raw_rule_b sqrt2_rule_a sqrt2_rule_b sqrt2K_rule_a sqrt2K_rule_b
0.5 0 5 5 5 5 28 5
0.5 1 5 5 5 5 21 5
0.5 2 5 5 5 5 22 5
1.0 0 5 5 5 5 33 5
1.0 1 5 5 5 5 16 5
1.0 2 5 5 5 5 23 5
2.0 0 5 4 5 5 45 5
2.0 1 5 5 5 5 26 5
2.0 2 5 5 5 5 26 5
Result: partially reproduced, under the kernel reading adopted for every claim and resting on an assumption about the scale of the noise. With rule (b), the noise of the trial average, window (a) and a kernel of 100 ms standard deviation, the medians over all 21 sessions are 7 for M1 and 10 for PMd (paper 8.5 and 14, both outside tolerance) and PMd exceeds M1 in 21 of 21 sessions (19 of 19 in the second session set, medians 7 and 11). Over the 2002 choices of 16 sessions (9 Monkey C with the 7 Monkey M), the medians are 7 to 7.5 (M1) and 10 to 11 (PMd), outside tolerance in every choice, and PMd is above M1 in every session in all 2002, so the outcome does not depend on which 16 sessions are the paper's. The closest session is M 2014-03-03, 9 (PMd) against 8 (M1). By monkey, C has medians 7 and 10 and M 7 and 11. No other reading under this kernel meets the criterion: window (b) gives PMd above M1 in 20 of 21, and rule (a) in 19 and 18 of 21. With a kernel of 100 ms full width at half maximum the medians would be 9 and 13 over 21 sessions, within tolerance (in 910 of the 2002 choices of 16), which would make the claim reproduced; that reading is not the one adopted (see "One kernel reading for every claim"). The outcome was reproduced, on the full-width reading, until the third verification round agreed to this change. Rule (a) gives medians of 26 to 47 for M1 and 46 to 78 for PMd. With the noise taken as the raw difference of two trials, or that difference over the square root of two, no component exceeds the noise and every estimate is 0. On synthetic data of rank 5, rule (b) with the noise of the trial average returns 5 in all nine runs, while rule (a) with the same scale returns 16 to 45, which is why rule (a) gives large values here; the raw and root-two scales also return 5 (4 once) on these synthetic data, whose signal is strong against the single-trial noise, so the synthetic check does not rule them out in general, only shows that the adopted scale with rule (b) recovers a known rank.
figs4d-dimensionality-robust-to-size¶
With random subsamples of equal, fixed numbers of neurons, estimated dimensionality fell for small samples but PMd stayed higher-dimensional than M1 for every neuron count tested.
As fig3b-pmd-higher-dimensional on random subsamples of 40, 50, 60 and 70 neurons per area and all neurons, mean over sessions; sessions with fewer neurons than a size are left out of that size. Run with the noise of the trial average, under both kernel readings (the full-width kernel is the one whose fig3b medians match the paper) and both windows; both rules are reported. 100 subsamples per size (the claims file's standard choice), with 100 noise repetitions per subsample to fit the time limit (the paper's 1000 are kept for "all"); with 100 repetitions the 99th percentile is close to the largest of them, which can only raise the noise limit slightly.
t0 = time.time()
SIZES = (40, 50, 60, 70)
N_SUB, SUB_REPS = 100, 100
def subsample_for(i):
s = sessions[i]
rng = np.random.default_rng([4, i])
out = []
for window, kern in itertools.product(("a", "b"), ("sd", "fwhm")):
fit = fits[(s["sid"], window, kern)]
tids = fit["trials"].target_id.to_numpy()
for area in ("M1", "PMd"):
n_units = len(fit[area]["units"])
for size in SIZES:
if n_units < size:
continue
for k in range(N_SUB):
sub = rng.choice(n_units, size, replace=False)
d = perich.dimensionality([r[:, sub] for r in fit[area]["rates"]], tids, rng, reps=SUB_REPS)
out.append({"session": s["sid"], "window": window, "kernel": kern, "area": area, "size": size,
"rule_a": d["sqrt2K"][0], "rule_b": d["sqrt2K"][1]})
return out
rows = []
with ThreadPoolExecutor(4) as ex:
for r_ in ex.map(subsample_for, range(len(sessions))):
rows += r_
sub = pd.DataFrame(rows)
print(f"{len(sub)} subsample estimates in {time.time() - t0:.0f} s")
full = dims[dims.scale == "sqrt2K"].assign(size="all")[["session", "window", "kernel", "area", "size", "rule_a", "rule_b"]]
per_sess = pd.concat([sub.groupby(["session", "window", "kernel", "area", "size"])[["rule_a", "rule_b"]].mean().reset_index(), full])
per_sess["size"] = per_sess["size"].astype(str)
for name, sids in SETS.items():
g = per_sess[per_sess.session.isin(sids)]
for rule in ("rule_b", "rule_a"):
t = g.groupby(["kernel", "window", "size", "area"])[rule].agg(["mean", "sem", "count"]).unstack("area")
t[("diff", "PMd-M1")] = t[("mean", "PMd")] - t[("mean", "M1")]
print(f"\nset '{name}', {rule}: mean and s.e.m. over sessions, and the number of sessions")
print(t.round(2).to_string())
55200 subsample estimates in 373 s
set 'all 21', rule_b: mean and s.e.m. over sessions, and the number of sessions
mean sem count diff
area M1 PMd M1 PMd M1 PMd PMd-M1
kernel window size
fwhm a 40 8.03 9.90 0.42 0.28 17 21 1.88
50 8.07 10.70 0.27 0.32 15 21 2.63
60 8.53 11.33 0.25 0.35 13 21 2.80
70 8.65 11.85 0.29 0.43 11 19 3.20
all 8.52 13.05 0.38 0.48 21 21 4.52
b 40 6.74 8.07 0.89 0.62 17 21 1.33
50 5.94 8.56 0.20 0.74 15 21 2.62
60 6.00 8.97 0.22 0.84 13 21 2.97
70 6.13 9.16 0.25 1.01 11 19 3.03
all 6.52 9.62 0.80 1.02 21 21 3.10
sd a 40 6.77 8.39 0.27 0.18 17 21 1.62
50 6.87 8.90 0.24 0.20 15 21 2.03
60 7.12 9.35 0.23 0.21 13 21 2.23
70 7.20 9.66 0.26 0.24 11 19 2.46
all 7.10 10.57 0.25 0.24 21 21 3.48
b 40 5.82 7.12 0.49 0.37 17 21 1.30
50 5.45 7.50 0.19 0.41 15 21 2.05
60 5.48 7.73 0.19 0.44 13 21 2.26
70 5.59 7.80 0.24 0.50 11 19 2.21
all 5.90 8.10 0.40 0.47 21 21 2.19
set 'all 21', rule_a: mean and s.e.m. over sessions, and the number of sessions
mean sem count diff
area M1 PMd M1 PMd M1 PMd PMd-M1
kernel window size
fwhm a 40 32.52 34.08 0.75 0.60 17 21 1.56
50 38.90 41.90 0.91 0.76 15 21 3.00
60 44.53 48.96 1.31 0.97 13 21 4.43
70 50.65 54.89 2.16 1.18 11 19 4.24
all 45.48 75.71 2.99 3.19 21 21 30.24
b 40 23.49 23.36 1.25 1.10 17 21 -0.13
50 26.93 27.98 1.23 1.40 15 21 1.05
60 29.20 31.75 1.76 1.80 13 21 2.56
70 33.31 35.89 1.93 2.40 11 19 2.58
all 28.29 52.81 2.39 5.07 21 21 24.52
sd a 40 30.97 33.46 0.67 0.43 17 21 2.49
50 35.41 39.76 0.73 0.61 15 21 4.35
60 38.62 44.51 1.08 0.76 13 21 5.89
70 41.15 47.48 1.46 1.00 11 19 6.33
all 37.57 53.86 1.93 1.07 21 21 16.29
b 40 22.31 23.18 1.26 0.98 17 21 0.87
50 24.83 27.01 1.18 1.31 15 21 2.18
60 26.66 30.25 1.68 1.56 13 21 3.59
70 30.67 33.39 1.98 2.01 11 19 2.72
all 27.62 45.00 2.18 3.18 21 21 17.38
set '19 within Baseline reward range', rule_b: mean and s.e.m. over sessions, and the number of sessions
mean sem count diff
area M1 PMd M1 PMd M1 PMd PMd-M1
kernel window size
fwhm a 40 8.02 9.93 0.47 0.31 15 19 1.91
50 7.98 10.75 0.29 0.35 13 19 2.77
60 8.45 11.37 0.27 0.39 11 19 2.92
70 8.53 11.91 0.32 0.48 9 17 3.38
all 8.42 13.05 0.41 0.53 19 19 4.63
b 40 6.88 8.21 1.00 0.68 15 19 1.34
50 5.95 8.73 0.24 0.80 13 19 2.78
60 6.01 9.17 0.26 0.92 11 19 3.16
70 6.16 9.39 0.31 1.12 9 17 3.24
all 6.58 9.84 0.89 1.11 19 19 3.26
sd a 40 6.76 8.42 0.31 0.20 15 19 1.66
50 6.82 8.93 0.26 0.22 13 19 2.11
60 7.08 9.41 0.26 0.23 11 19 2.32
70 7.13 9.71 0.30 0.27 9 17 2.58
all 7.05 10.63 0.27 0.27 19 19 3.58
b 40 5.90 7.25 0.56 0.39 15 19 1.34
50 5.47 7.64 0.22 0.44 13 19 2.17
60 5.48 7.88 0.22 0.48 11 19 2.39
70 5.61 7.95 0.28 0.55 9 17 2.35
all 5.95 8.26 0.44 0.51 19 19 2.32
set '19 within Baseline reward range', rule_a: mean and s.e.m. over sessions, and the number of sessions
mean sem count diff
area M1 PMd M1 PMd M1 PMd PMd-M1
kernel window size
fwhm a 40 32.37 34.07 0.82 0.66 15 19 1.70
50 38.65 41.89 0.98 0.84 13 19 3.24
60 43.89 49.00 1.36 1.07 11 19 5.11
70 49.91 54.93 2.34 1.31 9 17 5.02
all 44.05 74.84 3.04 3.45 19 19 30.79
b 40 23.36 23.07 1.41 1.19 15 19 -0.30
50 26.63 27.62 1.39 1.52 13 19 0.99
60 28.52 31.27 1.97 1.95 11 19 2.75
70 32.62 35.30 2.23 2.65 9 17 2.68
all 27.21 51.21 2.50 5.48 19 19 24.00
sd a 40 30.71 33.39 0.71 0.48 15 19 2.68
50 34.95 39.65 0.70 0.67 13 19 4.70
60 37.96 44.47 1.03 0.84 11 19 6.51
70 40.17 47.41 1.31 1.11 9 17 7.25
all 36.37 53.89 1.85 1.18 19 19 17.53
b 40 22.09 22.90 1.42 1.06 15 19 0.82
50 24.19 26.59 1.25 1.41 13 19 2.40
60 25.61 29.70 1.76 1.67 11 19 4.09
70 29.49 32.66 2.12 2.18 9 17 3.17
all 26.16 43.37 2.08 3.29 19 19 17.21
# Same comparison on the sessions that have both areas at each size (paired), rule (b)
for name, sids in SETS.items():
g = per_sess[per_sess.session.isin(sids)]
for k, w in itertools.product(("sd", "fwhm"), ("a", "b")):
for size in [str(x) for x in SIZES] + ["all"]:
p = g[(g.kernel == k) & (g.window == w) & (g["size"] == size)].pivot(index="session", columns="area", values="rule_b").dropna()
if len(p):
print(f"set '{name}', {k} kernel, window {w}, size {size}: {len(p)} sessions with both areas, "
f"PMd > M1 in {int((p.PMd > p.M1).sum())}, mean PMd {p.PMd.mean():.2f}, M1 {p.M1.mean():.2f}")
set 'all 21', sd kernel, window a, size 40: 17 sessions with both areas, PMd > M1 in 17, mean PMd 8.36, M1 6.77 set 'all 21', sd kernel, window a, size 50: 15 sessions with both areas, PMd > M1 in 15, mean PMd 8.64, M1 6.87 set 'all 21', sd kernel, window a, size 60: 13 sessions with both areas, PMd > M1 in 13, mean PMd 9.02, M1 7.12 set 'all 21', sd kernel, window a, size 70: 11 sessions with both areas, PMd > M1 in 11, mean PMd 9.35, M1 7.20 set 'all 21', sd kernel, window a, size all: 21 sessions with both areas, PMd > M1 in 21, mean PMd 10.57, M1 7.10 set 'all 21', sd kernel, window b, size 40: 17 sessions with both areas, PMd > M1 in 17, mean PMd 7.01, M1 5.82 set 'all 21', sd kernel, window b, size 50: 15 sessions with both areas, PMd > M1 in 15, mean PMd 6.89, M1 5.45 set 'all 21', sd kernel, window b, size 60: 13 sessions with both areas, PMd > M1 in 13, mean PMd 6.83, M1 5.48 set 'all 21', sd kernel, window b, size 70: 11 sessions with both areas, PMd > M1 in 11, mean PMd 6.99, M1 5.59 set 'all 21', sd kernel, window b, size all: 21 sessions with both areas, PMd > M1 in 20, mean PMd 8.10, M1 5.90 set 'all 21', fwhm kernel, window a, size 40: 17 sessions with both areas, PMd > M1 in 17, mean PMd 10.04, M1 8.03 set 'all 21', fwhm kernel, window a, size 50: 15 sessions with both areas, PMd > M1 in 15, mean PMd 10.57, M1 8.07 set 'all 21', fwhm kernel, window a, size 60: 13 sessions with both areas, PMd > M1 in 13, mean PMd 11.13, M1 8.53 set 'all 21', fwhm kernel, window a, size 70: 11 sessions with both areas, PMd > M1 in 11, mean PMd 11.61, M1 8.65 set 'all 21', fwhm kernel, window a, size all: 21 sessions with both areas, PMd > M1 in 21, mean PMd 13.05, M1 8.52 set 'all 21', fwhm kernel, window b, size 40: 17 sessions with both areas, PMd > M1 in 16, mean PMd 8.09, M1 6.74 set 'all 21', fwhm kernel, window b, size 50: 15 sessions with both areas, PMd > M1 in 15, mean PMd 7.69, M1 5.94 set 'all 21', fwhm kernel, window b, size 60: 13 sessions with both areas, PMd > M1 in 13, mean PMd 7.65, M1 6.00 set 'all 21', fwhm kernel, window b, size 70: 11 sessions with both areas, PMd > M1 in 11, mean PMd 7.82, M1 6.13 set 'all 21', fwhm kernel, window b, size all: 21 sessions with both areas, PMd > M1 in 20, mean PMd 9.62, M1 6.52 set '19 within Baseline reward range', sd kernel, window a, size 40: 15 sessions with both areas, PMd > M1 in 15, mean PMd 8.40, M1 6.76 set '19 within Baseline reward range', sd kernel, window a, size 50: 13 sessions with both areas, PMd > M1 in 13, mean PMd 8.64, M1 6.82 set '19 within Baseline reward range', sd kernel, window a, size 60: 11 sessions with both areas, PMd > M1 in 11, mean PMd 9.05, M1 7.08 set '19 within Baseline reward range', sd kernel, window a, size 70: 9 sessions with both areas, PMd > M1 in 9, mean PMd 9.38, M1 7.13 set '19 within Baseline reward range', sd kernel, window a, size all: 19 sessions with both areas, PMd > M1 in 19, mean PMd 10.63, M1 7.05 set '19 within Baseline reward range', sd kernel, window b, size 40: 15 sessions with both areas, PMd > M1 in 15, mean PMd 7.15, M1 5.90 set '19 within Baseline reward range', sd kernel, window b, size 50: 13 sessions with both areas, PMd > M1 in 13, mean PMd 6.99, M1 5.47 set '19 within Baseline reward range', sd kernel, window b, size 60: 11 sessions with both areas, PMd > M1 in 11, mean PMd 6.91, M1 5.48 set '19 within Baseline reward range', sd kernel, window b, size 70: 9 sessions with both areas, PMd > M1 in 9, mean PMd 7.09, M1 5.61 set '19 within Baseline reward range', sd kernel, window b, size all: 19 sessions with both areas, PMd > M1 in 18, mean PMd 8.26, M1 5.95 set '19 within Baseline reward range', fwhm kernel, window a, size 40: 15 sessions with both areas, PMd > M1 in 15, mean PMd 10.09, M1 8.02 set '19 within Baseline reward range', fwhm kernel, window a, size 50: 13 sessions with both areas, PMd > M1 in 13, mean PMd 10.62, M1 7.98 set '19 within Baseline reward range', fwhm kernel, window a, size 60: 11 sessions with both areas, PMd > M1 in 11, mean PMd 11.16, M1 8.45 set '19 within Baseline reward range', fwhm kernel, window a, size 70: 9 sessions with both areas, PMd > M1 in 9, mean PMd 11.67, M1 8.53 set '19 within Baseline reward range', fwhm kernel, window a, size all: 19 sessions with both areas, PMd > M1 in 19, mean PMd 13.05, M1 8.42 set '19 within Baseline reward range', fwhm kernel, window b, size 40: 15 sessions with both areas, PMd > M1 in 14, mean PMd 8.28, M1 6.88 set '19 within Baseline reward range', fwhm kernel, window b, size 50: 13 sessions with both areas, PMd > M1 in 13, mean PMd 7.81, M1 5.95 set '19 within Baseline reward range', fwhm kernel, window b, size 60: 11 sessions with both areas, PMd > M1 in 11, mean PMd 7.76, M1 6.01 set '19 within Baseline reward range', fwhm kernel, window b, size 70: 9 sessions with both areas, PMd > M1 in 9, mean PMd 7.95, M1 6.16 set '19 within Baseline reward range', fwhm kernel, window b, size all: 19 sessions with both areas, PMd > M1 in 18, mean PMd 9.84, M1 6.58
Result: reproduced, resting on the same assumption as fig3b-pmd-higher-dimensional. Under rule (b) the mean PMd dimensionality exceeds the mean M1 dimensionality at every size (40, 50, 60, 70 and all neurons) under both kernels, both windows and both session sets, by 1.30 to 4.63. With the full-width kernel and window (a), the reading whose fig3b medians match the paper, the difference is 1.88 at 40 neurons and 4.52 with all neurons; with the standard-deviation kernel and window (a), 1.62 and 3.48. Among sessions with both areas at a size, PMd is higher in every session except one at "all" in window (b) under both kernels (20 of 21) and one at 40 neurons with the full-width kernel in window (b) (16 of 17). Rule (a) gives PMd above M1 at every size except 40 neurons with the full-width kernel in window (b) (differences -0.13 and -0.30 in the two sets). The PMd mean falls at small sizes in every reading (window a, full-width kernel, 13.05 with all neurons against 9.90 at 40; standard-deviation kernel 10.57 against 8.39). The M1 mean does not fall consistently: in window (a) with the full-width kernel it is 8.52 with all neurons, 8.53 and 8.65 at 60 and 70, and 8.03 at 40; with the standard-deviation kernel 7.10, 7.12 and 7.20, and 6.77 at 40; in window (b) with the full-width kernel it is higher at 40 (6.74) than with all neurons (6.52). The statement's fall for small samples is therefore shown for PMd only; it is outside the criterion. The means at 40 neurons cover 17 sessions for M1, against 21 for all neurons. The M1 means at 50 to 70 neurons come from fewer sessions (11 to 15), as Monkey M has 26 to 52 M1 units.
figs4c-null-potent-index-centered¶
Single neurons contributed to both output-potent and output-null subspaces; the null/potent index was centered on zero with virtually no neurons near -1 or 1.
Index per PMd neuron from the product of the PCA loading matrix (16 PMd components) and V (M3, O7), as (sum of absolute potent weights - sum of absolute null weights) / (sum of both); the sum of squares is the second reading. Reproduced when the mean over neurons lies between -0.1 and 0.1 and fewer than 1% of neurons have an absolute index above 0.9. All PMd units of the candidate sessions are pooled. The control shuffles the potent weights across neurons against the null weights, 1000 times per session (STAR Methods, "randomly shuffling the weights across neurons").
rng = np.random.default_rng(5)
rows, pooled = [], {}
for (sid, window, kern), fit in fits.items():
for power in (1, 2):
idx = perich.null_potent_index(fit["PMd"]["axes"], fit["potent"], fit["null"], power)
pooled.setdefault((window, kern, power), {})[sid] = idx
summ = []
for name, sids in SETS.items():
for (window, kern, power), d in pooled.items():
v = np.concatenate([d[x] for x in sids])
summ.append({"set": name, "window": window, "kernel": kern, "magnitude": "abs" if power == 1 else "squares",
"neurons": len(v), "mean": round(v.mean(), 3), "sd": round(v.std(), 3),
"median": round(np.median(v), 3), "frac_abs>0.9": round((np.abs(v) > 0.9).mean(), 4),
"min": round(v.min(), 3), "max": round(v.max(), 3),
"meets": bool(-0.1 <= v.mean() <= 0.1 and (np.abs(v) > 0.9).mean() < 0.01)})
print(pd.DataFrame(summ).to_string(index=False))
v = np.concatenate(list(pooled[("a", "fwhm", 2)].values()))
print(f"\nWindow (a), full-width kernel, squares, all 21 sessions: mean {v.mean():.6f}")
for mk in ("C", "M"):
for key in (("a", "sd", 1), ("a", "fwhm", 1), ("a", "fwhm", 2)):
v = np.concatenate([x for sid, x in pooled[key].items() if sid[4] == mk])
print(f"Monkey {mk}, window {key[0]}, {key[1]} kernel, {'abs' if key[2] == 1 else 'squares'}: {len(v)} neurons, "
f"mean {v.mean():.3f}, fraction |index| > 0.9 {(np.abs(v) > 0.9).mean():.4f}")
# Shuffle control for the main reading (window a, sd kernel, sums of absolute values)
sh_means, sh_frac, sh_sd = [], [], []
weights = [(np.abs(fits[(s["sid"], "a", "sd")]["PMd"]["axes"] @ fits[(s["sid"], "a", "sd")]["potent"]).sum(1),
np.abs(fits[(s["sid"], "a", "sd")]["PMd"]["axes"] @ fits[(s["sid"], "a", "sd")]["null"]).sum(1)) for s in sessions]
for k in range(1000):
vals = []
for wp, wn in weights:
wp = rng.permutation(wp)
vals.append((wp - wn) / (wp + wn))
v = np.concatenate(vals)
sh_means.append(v.mean()); sh_frac.append((np.abs(v) > 0.9).mean()); sh_sd.append(v.std())
print(f"\nShuffled control (window a, sd kernel, absolute values, 21 sessions, 1000 shuffles): mean {np.mean(sh_means):.3f} "
f"(2.5-97.5%: {np.percentile(sh_means, 2.5):.3f} to {np.percentile(sh_means, 97.5):.3f}), sd {np.mean(sh_sd):.3f}, "
f"fraction |index| > 0.9: {np.mean(sh_frac):.4f}")
per = pd.DataFrame([{"session": sid, "mean": round(v.mean(), 3), "sd": round(v.std(), 3), "n": len(v)}
for sid, v in pooled[("a", "sd", 1)].items()])
print("\nPer session, main reading:"); print(per.to_string(index=False))
v = np.concatenate(list(pooled[("a", "sd", 1)].values()))
plt.figure(figsize=(4, 3)); plt.hist(v, np.linspace(-1, 1, 41), color="k")
plt.xlabel("null/potent index"); plt.ylabel("PMd neurons"); plt.title("window (a), sd kernel, all 21 sessions")
plt.tight_layout(); plt.show()
print(f"Index range {v.min():.3f} to {v.max():.3f}; all points inside the axes.")
set window kernel magnitude neurons mean sd median frac_abs>0.9 min max meets
all 21 a sd abs 3263 0.043 0.217 0.050 0.0000 -0.661 0.683 True
all 21 a sd squares 3263 0.079 0.376 0.089 0.0015 -0.946 0.923 True
all 21 a fwhm abs 3263 0.056 0.222 0.057 0.0000 -0.701 0.691 True
all 21 a fwhm squares 3263 0.100 0.384 0.117 0.0009 -0.923 0.914 False
all 21 b sd abs 3263 0.017 0.209 0.021 0.0000 -0.632 0.659 True
all 21 b sd squares 3263 0.032 0.366 0.041 0.0015 -0.918 0.912 True
all 21 b fwhm abs 3263 0.035 0.216 0.029 0.0000 -0.667 0.741 True
all 21 b fwhm squares 3263 0.059 0.374 0.057 0.0021 -0.941 0.945 True
19 within Baseline reward range a sd abs 2914 0.041 0.217 0.046 0.0000 -0.661 0.683 True
19 within Baseline reward range a sd squares 2914 0.075 0.377 0.083 0.0017 -0.946 0.923 True
19 within Baseline reward range a fwhm abs 2914 0.054 0.223 0.054 0.0000 -0.701 0.691 True
19 within Baseline reward range a fwhm squares 2914 0.097 0.385 0.111 0.0010 -0.923 0.914 True
19 within Baseline reward range b sd abs 2914 0.018 0.209 0.021 0.0000 -0.632 0.659 True
19 within Baseline reward range b sd squares 2914 0.035 0.365 0.048 0.0017 -0.918 0.912 True
19 within Baseline reward range b fwhm abs 2914 0.035 0.217 0.029 0.0000 -0.667 0.741 True
19 within Baseline reward range b fwhm squares 2914 0.061 0.376 0.060 0.0024 -0.941 0.945 True
Window (a), full-width kernel, squares, all 21 sessions: mean 0.100085
Monkey C, window a, sd kernel, abs: 2664 neurons, mean 0.028, fraction |index| > 0.9 0.0000
Monkey C, window a, fwhm kernel, abs: 2664 neurons, mean 0.045, fraction |index| > 0.9 0.0000
Monkey C, window a, fwhm kernel, squares: 2664 neurons, mean 0.081, fraction |index| > 0.9 0.0008
Monkey M, window a, sd kernel, abs: 599 neurons, mean 0.109, fraction |index| > 0.9 0.0000
Monkey M, window a, fwhm kernel, abs: 599 neurons, mean 0.103, fraction |index| > 0.9 0.0000
Monkey M, window a, fwhm kernel, squares: 599 neurons, mean 0.186, fraction |index| > 0.9 0.0017
Shuffled control (window a, sd kernel, absolute values, 21 sessions, 1000 shuffles): mean 0.036 (2.5-97.5%: 0.033 to 0.038), sd 0.483, fraction |index| > 0.9: 0.0318
Per session, main reading:
session mean sd n
sub-C_ses-CO-20160909 0.015 0.183 191
sub-C_ses-CO-20160912 0.011 0.214 162
sub-C_ses-CO-20160914 -0.007 0.219 258
sub-C_ses-CO-20160915 0.004 0.203 233
sub-C_ses-CO-20160919 0.052 0.231 235
sub-C_ses-CO-20160921 0.028 0.210 232
sub-C_ses-CO-20160923 0.007 0.205 223
sub-C_ses-CO-20160929 0.062 0.232 114
sub-C_ses-CO-20161005 0.015 0.208 167
sub-C_ses-CO-20161006 0.044 0.249 192
sub-C_ses-CO-20161007 0.023 0.209 137
sub-C_ses-CO-20161011 0.047 0.208 171
sub-C_ses-CO-20161013 0.081 0.220 138
sub-C_ses-CO-20161021 0.051 0.212 211
sub-M_ses-CO-20140203 0.143 0.226 80
sub-M_ses-CO-20140217 0.112 0.198 104
sub-M_ses-CO-20140218 0.135 0.207 121
sub-M_ses-CO-20140303 0.095 0.212 66
sub-M_ses-CO-20140304 0.088 0.201 76
sub-M_ses-CO-20140306 0.081 0.200 86
sub-M_ses-CO-20140307 0.092 0.226 66
Index range -0.661 to 0.683; all points inside the axes.
Result: reproduced, resting on the reading of "magnitude" as the sum of absolute values. Over 3263 PMd units of the 21 sessions, the mean index is 0.017 to 0.056 across the four window and kernel readings, and no unit has an absolute index above 0.9 (range -0.701 to 0.741). With sums of squares, three of four readings also meet the criterion; window (a) with the full-width kernel gives a mean of 0.100085, just above the bound of 0.1 (0.09% of units above 0.9). The shuffled control has a similar mean (0.036) and a wider spread (standard deviation 0.483 against 0.217, 3.2% of units above 0.9). The index leans positive in Monkey M: its session means are 0.081 to 0.143 under the main reading, and its 599 pooled units have a mean of 0.109 (above the bound of 0.1 if Monkey M were judged alone), against 0.028 for the 2664 Monkey C units (session means -0.007 to 0.081).
fig6c-null-leads-potent¶
In Baseline, output-null activity preceded output-potent activity by 35 ms on average across the 16 sessions, a lag significantly below zero (t-test, p = 0.02), and output-null also led at the end of Force.
Method from STAR Methods, Calculating the lag (helper cca_lag): Baseline trials, the 700 ms after
target presentation, potent activity shifted from -350 to 350 ms in 10 ms steps against null
activity, CCA at each lag, the lag of the maximum first canonical correlation per session; trials
concatenated within a session. Readings of the time resolution: (a) 10 ms bins smoothed as M2,
(b) 50 ms latent activity interpolated to 10 ms. The subspaces come from the PCA and W of the
window readings (a) and (b) (O3), under both kernel readings (O4). One-sample two-sided t-test of
the session lags against zero. A maximum within one step of the edge of the search (|lag| of 340 or
350 ms) is flagged, and the test is shown with such sessions kept and dropped. As the paper's test
is over 16 sessions, it is also run on every choice of 9 of the 14 Monkey C candidates with the 7
Monkey M sessions. The end-of-Force part cannot be tested: no Force trials are in the dandiset.
t0 = time.time()
def lags_for(s):
out_rows, out_curves = [], {}
for window in ("a", "b"):
for kern in ("sd", "fwhm"):
fit = fits[(s["sid"], window, kern)]
onsets = fit["trials"].target_on.to_numpy()
for res in ("a", "b"):
lat = perich.lag_grid_latents(s, fit["PMd"]["units"], onsets, fit["PMd"]["axes"], res,
perich.KERNEL_SD_S[kern])
lag, lags, ccs = perich.cca_lag([l @ fit["null"] for l in lat], [l @ fit["potent"] for l in lat])
out_curves[(s["sid"], window, kern, res)] = ccs
out_rows.append({"session": s["sid"], "monkey": s["monkey"], "window": window, "kernel": kern,
"resolution": res, "lag_ms": round(lag * 1000), "peak_cc": round(ccs.max(), 3),
"cc_at_0": round(ccs[np.argmin(np.abs(lags))], 3)})
return out_rows, out_curves
rows, curves = [], {}
with ThreadPoolExecutor(4) as ex:
for r_, c_ in ex.map(lags_for, sessions):
rows += r_; curves.update(c_)
lagdf = pd.DataFrame(rows)
# A maximum at the last or second-to-last lag of the search (|lag| >= 340 ms) is an edge maximum: the
# correlation may still rise outside the range, so the lag is not located. Shown with and without them.
lagdf["edge"] = lagdf.lag_ms.abs() >= 340
print(f"{len(lagdf)} lag estimates in {time.time() - t0:.0f} s")
print("Edge maxima:", lagdf[lagdf.edge][["session", "window", "kernel", "resolution", "lag_ms"]].to_dict("records"))
summ = []
for (name, sids), drop_edge in itertools.product(SETS.items(), (False, True)):
d_ = lagdf[lagdf.session.isin(sids) & ~(lagdf.edge & drop_edge)]
for key, g in d_.groupby(["window", "kernel", "resolution"]):
t = stats.ttest_1samp(g.lag_ms, 0)
summ.append({"set": name, "edge_maxima": "dropped" if drop_edge else "kept", "window": key[0], "kernel": key[1], "resolution": key[2], "sessions": len(g),
"mean_lag_ms": round(g.lag_ms.mean(), 1), "sem": round(g.lag_ms.sem(), 1),
"median": g.lag_ms.median(), "negative": int((g.lag_ms < 0).sum()), "zero": int((g.lag_ms == 0).sum()),
"t": round(t.statistic, 2), "p": round(t.pvalue, 4),
"within_-35+-3.5": bool(abs(g.lag_ms.mean() + 35) <= 3.5)})
print(pd.DataFrame(summ).to_string(index=False))
168 lag estimates in 96 s
Edge maxima: [{'session': 'sub-C_ses-CO-20161006', 'window': 'a', 'kernel': 'sd', 'resolution': 'b', 'lag_ms': -350}, {'session': 'sub-C_ses-CO-20161011', 'window': 'a', 'kernel': 'sd', 'resolution': 'b', 'lag_ms': -340}]
set edge_maxima window kernel resolution sessions mean_lag_ms sem median negative zero t p within_-35+-3.5
all 21 kept a fwhm a 21 -3.3 3.7 0.0 8 7 -0.91 0.3750 False
all 21 kept a fwhm b 21 -5.7 5.8 0.0 9 5 -0.99 0.3360 False
all 21 kept a sd a 21 -18.1 7.9 -30.0 13 3 -2.29 0.0327 False
all 21 kept a sd b 21 -51.4 22.8 -20.0 15 1 -2.26 0.0353 False
all 21 kept b fwhm a 21 -4.8 3.4 -10.0 11 5 -1.39 0.1796 False
all 21 kept b fwhm b 21 -7.1 5.6 -10.0 11 5 -1.27 0.2177 False
all 21 kept b sd a 21 -5.2 6.7 -20.0 13 2 -0.78 0.4467 False
all 21 kept b sd b 21 -2.4 10.3 -10.0 12 3 -0.23 0.8192 False
all 21 dropped a fwhm a 21 -3.3 3.7 0.0 8 7 -0.91 0.3750 False
all 21 dropped a fwhm b 21 -5.7 5.8 0.0 9 5 -0.99 0.3360 False
all 21 dropped a sd a 21 -18.1 7.9 -30.0 13 3 -2.29 0.0327 False
all 21 dropped a sd b 19 -20.5 8.9 -20.0 13 1 -2.30 0.0338 False
all 21 dropped b fwhm a 21 -4.8 3.4 -10.0 11 5 -1.39 0.1796 False
all 21 dropped b fwhm b 21 -7.1 5.6 -10.0 11 5 -1.27 0.2177 False
all 21 dropped b sd a 21 -5.2 6.7 -20.0 13 2 -0.78 0.4467 False
all 21 dropped b sd b 21 -2.4 10.3 -10.0 12 3 -0.23 0.8192 False
19 within Baseline reward range kept a fwhm a 19 -2.1 3.5 0.0 7 7 -0.61 0.5515 False
19 within Baseline reward range kept a fwhm b 19 -4.2 5.0 0.0 8 5 -0.84 0.4138 False
19 within Baseline reward range kept a sd a 19 -18.9 8.6 -30.0 12 3 -2.21 0.0404 False
19 within Baseline reward range kept a sd b 19 -54.7 25.0 -20.0 14 1 -2.19 0.0420 False
19 within Baseline reward range kept b fwhm a 19 -5.3 3.0 -10.0 10 5 -1.76 0.0962 False
19 within Baseline reward range kept b fwhm b 19 -7.9 4.9 -10.0 10 5 -1.62 0.1225 False
19 within Baseline reward range kept b sd a 19 -1.6 6.9 -10.0 11 2 -0.23 0.8225 False
19 within Baseline reward range kept b sd b 19 -11.6 8.8 -20.0 12 3 -1.32 0.2029 False
19 within Baseline reward range dropped a fwhm a 19 -2.1 3.5 0.0 7 7 -0.61 0.5515 False
19 within Baseline reward range dropped a fwhm b 19 -4.2 5.0 0.0 8 5 -0.84 0.4138 False
19 within Baseline reward range dropped a sd a 19 -18.9 8.6 -30.0 12 3 -2.21 0.0404 False
19 within Baseline reward range dropped a sd b 17 -20.6 9.7 -20.0 12 1 -2.13 0.0494 False
19 within Baseline reward range dropped b fwhm a 19 -5.3 3.0 -10.0 10 5 -1.76 0.0962 False
19 within Baseline reward range dropped b fwhm b 19 -7.9 4.9 -10.0 10 5 -1.62 0.1225 False
19 within Baseline reward range dropped b sd a 19 -1.6 6.9 -10.0 11 2 -0.23 0.8225 False
19 within Baseline reward range dropped b sd b 19 -11.6 8.8 -20.0 12 3 -1.32 0.2029 False
print("Per session lags (ms), standard-deviation kernel:")
print(lagdf[lagdf.kernel == "sd"].pivot(index="session", columns=["window", "resolution"], values="lag_ms").to_string())
for mk in ("C", "M"):
g = lagdf[(lagdf.monkey == mk) & (lagdf.kernel == "sd")]
print(f"Monkey {mk}:", g.groupby(["window", "resolution"]).lag_ms.agg(["mean", "count"]).round(1).to_dict("index"))
fig, axes = plt.subplots(2, 2, figsize=(9, 6), sharex=True)
for ax, (w, res) in zip(axes.ravel(), itertools.product(("a", "b"), ("a", "b"))):
for s in sessions:
ax.plot(perich.LAGS_S * 1000, curves[(s["sid"], w, "sd", res)], lw=0.8)
ax.axvline(0, color="k", lw=0.5); ax.set_title(f"window ({w}), resolution ({res}), sd kernel, 21 sessions")
allc = np.concatenate([curves[(s["sid"], w, "sd", res)] for s in sessions])
print(f"window ({w}), resolution ({res}): plotted correlations range {allc.min():.3f} to {allc.max():.3f}; all inside the axes.")
for ax in axes[1]:
ax.set_xlabel("lag of potent relative to null (ms; negative: null leads)")
for ax in axes[:, 0]:
ax.set_ylabel("first canonical correlation")
plt.tight_layout(); plt.show()
Per session lags (ms), standard-deviation kernel:
window a b
resolution a b a b
session
sub-C_ses-CO-20160909 10 10 20 20
sub-C_ses-CO-20160912 -30 -40 60 -80
sub-C_ses-CO-20160914 -50 -50 -10 -10
sub-C_ses-CO-20160915 -10 -20 -30 0
sub-C_ses-CO-20160919 80 70 60 80
sub-C_ses-CO-20160921 20 30 40 50
sub-C_ses-CO-20160923 -70 -90 -30 -70
sub-C_ses-CO-20160929 -30 -40 20 30
sub-C_ses-CO-20161005 -70 -90 -20 -20
sub-C_ses-CO-20161006 -40 -350 20 0
sub-C_ses-CO-20161007 -30 -50 -30 -40
sub-C_ses-CO-20161011 -80 -340 0 -30
sub-C_ses-CO-20161013 -30 -50 -40 60
sub-C_ses-CO-20161021 10 10 -40 110
sub-M_ses-CO-20140203 0 10 -20 -30
sub-M_ses-CO-20140217 0 -10 0 0
sub-M_ses-CO-20140218 -30 -30 -10 -10
sub-M_ses-CO-20140303 0 -10 -20 -20
sub-M_ses-CO-20140304 -10 -20 -20 -30
sub-M_ses-CO-20140306 10 0 -20 -20
sub-M_ses-CO-20140307 -30 -20 -40 -40
Monkey C: {('a', 'a'): {'mean': -22.9, 'count': 14}, ('a', 'b'): {'mean': -71.4, 'count': 14}, ('b', 'a'): {'mean': 1.4, 'count': 14}, ('b', 'b'): {'mean': 7.1, 'count': 14}}
Monkey M: {('a', 'a'): {'mean': -8.6, 'count': 7}, ('a', 'b'): {'mean': -11.4, 'count': 7}, ('b', 'a'): {'mean': -18.6, 'count': 7}, ('b', 'b'): {'mean': -21.4, 'count': 7}}
window (a), resolution (a): plotted correlations range 0.411 to 0.837; all inside the axes.
window (a), resolution (b): plotted correlations range 0.418 to 0.825; all inside the axes.
window (b), resolution (a): plotted correlations range 0.454 to 0.896; all inside the axes.
window (b), resolution (b): plotted correlations range 0.457 to 0.858; all inside the axes.
# The paper's t-test is over 16 sessions (9 Monkey C, 7 Monkey M): the same test on every choice of
# 9 of the 14 Monkey C candidates with the 7 Monkey M sessions, per reading, edge maxima kept or dropped.
rows = []
sig_any = {False: np.zeros(len(CHOICES), bool), True: np.zeros(len(CHOICES), bool)}
for (w, k, res), g in lagdf.groupby(["window", "kernel", "resolution"]):
g = g.set_index("session")
for drop_edge in (False, True):
ps, means = [], []
for ss in CHOICES:
v = g.loc[ss]
v = v[~v.edge] if drop_edge else v
ps.append(stats.ttest_1samp(v.lag_ms, 0).pvalue); means.append(v.lag_ms.mean())
ps, means = np.array(ps), np.array(means)
sig = (ps < 0.05) & (means < 0)
sig_any[drop_edge] |= sig
rows.append({"window": w, "kernel": k, "resolution": res, "edge_maxima": "dropped" if drop_edge else "kept",
"mean_lag_range_ms": f"{means.min():.1f} to {means.max():.1f}", "median_p": round(np.median(ps), 3),
"significant_negative": f"{int(sig.sum())}/{len(CHOICES)}", "fraction": round(sig.mean(), 2),
"within_-35+-3.5_and_significant": int((sig & (np.abs(means + 35) <= 3.5)).sum())})
print(pd.DataFrame(rows).to_string(index=False))
for drop_edge, v in sig_any.items():
print(f"edge maxima {'dropped' if drop_edge else 'kept'}: significantly negative under at least one of the 8 readings "
f"in {int(v.sum())} of {len(CHOICES)} choices ({v.mean():.2f})")
window kernel resolution edge_maxima mean_lag_range_ms median_p significant_negative fraction within_-35+-3.5_and_significant
a fwhm a kept -9.4 to 3.8 0.502 11/2002 0.01 0
a fwhm a dropped -9.4 to 3.8 0.502 11/2002 0.01 0
a fwhm b kept -14.4 to 5.6 0.501 6/2002 0.00 0
a fwhm b dropped -14.4 to 5.6 0.501 6/2002 0.00 0
a sd a kept -30.6 to -4.4 0.095 763/2002 0.38 0
a sd a dropped -30.6 to -4.4 0.095 763/2002 0.38 0
a sd b kept -73.8 to -10.0 0.076 478/2002 0.24 0
a sd b dropped -35.0 to -4.3 0.071 816/2002 0.41 12
b fwhm a kept -11.2 to 0.6 0.178 279/2002 0.14 0
b fwhm a dropped -11.2 to 0.6 0.178 279/2002 0.14 0
b fwhm b kept -16.9 to 3.1 0.256 181/2002 0.09 0
b fwhm b dropped -16.9 to 3.1 0.256 181/2002 0.09 0
b sd a kept -19.4 to 3.8 0.340 157/2002 0.08 0
b sd a dropped -19.4 to 3.8 0.340 157/2002 0.08 0
b sd b kept -23.8 to 11.9 0.612 39/2002 0.02 0
b sd b dropped -23.8 to 11.9 0.612 39/2002 0.02 0
edge maxima kept: significantly negative under at least one of the 8 readings in 1013 of 2002 choices (0.51)
edge maxima dropped: significantly negative under at least one of the 8 readings in 1052 of 2002 choices (0.53)
Result: partially reproduced, under one reading, for Baseline only, and resting on the use of all 21 candidate sessions. The mean Baseline lag over the 21 sessions is negative (output-null first) under all 8 readings (-2.4 to -51.4 ms). It is significantly below zero only with window (a) and the 100 ms standard-deviation kernel: -18.1 ms (t = -2.29, p = 0.033) at 10 ms resolution and -51.4 ms (p = 0.035) with interpolated 50 ms activity. In the latter, two sessions peak within one step of the edge of the search (C 2016-10-06 at -350 ms and C 2016-10-11 at -340 ms); without them the mean is -20.5 ms over 19 sessions (p = 0.034). The full-width kernel and window (b) give p from 0.18 to 0.82. No reading over 21 sessions is within 3.5 ms of -35 ms. On the 2002 choices of 16 sessions the result is fragile: the main reading (window a, standard-deviation kernel, 10 ms) is significantly negative in 763 (38%), the 50 ms reading in 478 (24%; 816, 41%, with edge maxima dropped), and at least one of the 8 readings in 1013 (51%; 1052 with edge maxima dropped). Twelve choices, all with edge maxima dropped at 50 ms resolution, give a significant mean within tolerance of -35 ms. In the second session set the main reading is significant (p = 0.040). By monkey, window (a) and 10 ms resolution give -22.9 ms for C and -8.6 ms for M. The end-of-Force part cannot be tested. The outcome rests on the standard-deviation reading of the kernel, which is the reading adopted for every claim (below); with the full-width reading no reading is significant (p 0.18 to 0.38), which would make it not reproduced.
One kernel reading for every claim¶
The claims file names two readings of the kernel's "width 100 ms" (O4) and expected them not to change an outcome. They do, in opposite directions for fig3b-pmd-higher-dimensional and fig6c-null-leads-potent, so a reading chosen per claim would take the favourable one for each. The cell below gives, for each kernel reading applied to all claims, the numbers each criterion is judged on (rule b and the noise of the trial average for the dimensionality; all 21 sessions).
# Per kernel reading: the fig3b main reading, the figs4c readings, the figs4d comparison and the fig6c t-tests
rows = []
for kern in ("sd", "fwhm"):
d3 = dimsum[(dimsum.set == "all 21") & (dimsum.rule == "rule_b") & (dimsum.scale == "sqrt2K") & (dimsum.kernel == kern)]
f3 = "; ".join(f"window {r.window}: medians {r.M1_median}/{r.PMd_median}, PMd>M1 {r['PMd>M1']}, in tolerance {r.medians_in_tol}"
for _, r in d3.iterrows())
meets4c = []
for (w, k, power), d in pooled.items():
if k == kern:
v = np.concatenate(list(d.values()))
meets4c.append(bool(-0.1 <= v.mean() <= 0.1 and (np.abs(v) > 0.9).mean() < 0.01))
g4 = per_sess[per_sess.kernel == kern].groupby(["window", "size", "area"]).rule_b.mean().unstack("area")
lag_sig = []
for (w, res), g in lagdf[lagdf.kernel == kern].groupby(["window", "resolution"]):
t = stats.ttest_1samp(g.lag_ms, 0)
lag_sig.append(f"window {w} res {res}: {g.lag_ms.mean():.1f} ms, p {t.pvalue:.3f}")
rows.append({"kernel": kern, "fig3b": f3, "figs4c readings meeting": f"{sum(meets4c)}/{len(meets4c)}",
"figs4d PMd>M1 at every size and window": bool((g4.PMd > g4.M1).all()), "fig6c": "; ".join(lag_sig)})
for r in rows:
print(f"Kernel reading '{r['kernel']}':")
for k, v in r.items():
if k != "kernel":
print(f" {k}: {v}")
Kernel reading 'sd': fig3b: window a: medians 7.0/10.0, PMd>M1 21/21, in tolerance False; window b: medians 6.0/7.0, PMd>M1 20/21, in tolerance False figs4c readings meeting: 4/4 figs4d PMd>M1 at every size and window: True fig6c: window a res a: -18.1 ms, p 0.033; window a res b: -51.4 ms, p 0.035; window b res a: -5.2 ms, p 0.447; window b res b: -2.4 ms, p 0.819 Kernel reading 'fwhm': fig3b: window a: medians 9.0/13.0, PMd>M1 21/21, in tolerance True; window b: medians 6.0/8.0, PMd>M1 20/21, in tolerance False figs4c readings meeting: 3/4 figs4d PMd>M1 at every size and window: True fig6c: window a res a: -3.3 ms, p 0.375; window a res b: -5.7 ms, p 0.336; window b res a: -4.8 ms, p 0.180; window b res b: -7.1 ms, p 0.218
Kernel reading adopted for every claim: a standard deviation of 100 ms. STAR Methods,
Dimensionality reduction, says only "a non-causal Gaussian kernel of width 100 ms". The protocol
lets the author adopt the reading the paper best supports, as an assumption, when readings disagree
(the claims file's header instead says such a claim is blocked; the protocol's later rule is
followed here, and the conflict is recorded). The standard deviation is better supported: it is the
usual way a Gaussian smoothing kernel's width is stated, and the same group's later paper on
recordings that include Monkey C (Gallego, Perich et al. 2020, Nat. Neurosci., in paper/) states
its kernel width the same way, "a Gaussian kernel (s.d.: 50 ms)", applied there to square-root
transformed multiunit threshold crossings in 30 ms bins. Applied to all claims, the
standard-deviation reading gives figs4d, figs4c and fig6c their recorded outcomes, and fig3b
medians of 7 and 10 with PMd above M1 in every session, which the criterion counts as partially
reproduced. Under the full-width reading fig3b would be reproduced and no fig6c reading would be
significant (not reproduced), so neither reading favours every claim.
Claims not tested¶
The dandiset holds no Force, Rotation or Washout trials, and no label of which sessions were curl
field or visuomotor rotation sessions (data check). Every claim that needs those trials is blocked
with data_missing: fig2d-cf-behavioral-adaptation, fig7bc-vr-behavior-similar,
figs2-covariance-preserved-example, figs3a-covariance-like-baseline-split,
figs3b-covariance-preserved-vr, figs5e-glm-cell-counts, figs5e-pot-null-similar-training,
fig3e-local-rpr2-unchanged, fig3g-local-models-generalize, fig4d-null-m1-fails-early,
fig4c-null-tracks-behavior, figs6-pc-m1-fails, figs6de-per-monkey, figs5f-rpr2-distributions,
figs7a-unnormalized-change, figs7e-dimensionality-range, figs7f-factor-analysis,
figs7cd-single-neuron-models, fig5-washout-null-m1, fig6e-planning-change-in-null,
fig7de-vr-all-generalize, fig7f-vr-washout-unchanged and cf-vr-dissociation.
Summary¶
| Claim | Outcome |
|---|---|
| methods-neuron-counts | partially reproduced (assumed: which sessions are the paper's) |
| fig3b-pmd-higher-dimensional | partially reproduced (assumed: scale of the noise, standard-deviation kernel): PMd above M1 in every session, medians 7 and 10 outside tolerance |
| figs4d-dimensionality-robust-to-size | reproduced (assumed: scale of the noise) |
| figs4c-null-potent-index-centered | reproduced (assumed: magnitude as absolute values) |
| fig6c-null-leads-potent | partially reproduced, Baseline only (standard-deviation kernel, all 21 candidate sessions; significant in 38% of 16-session choices) |
| 23 claims on Force, Rotation and Washout trials | blocked, data missing |
A failure to reproduce here may be a failure of this attempt, not of the paper: see the reasoning and the differences from the paper's method.
How the analysis differs from the paper's. The paper's 16 sessions cannot be identified; all 21 dual-array center-out sessions of Monkeys C and M are used, with the 19 whose reward counts fit the paper's Baseline range as a second set, and every choice of 16 sessions (9 of the 14 Monkey C candidates with the 7 Monkey M sessions) for the session-level tests. All of each file's reward trials are taken as Baseline. Target acquisition is the trial's stop time minus the 0.5 s hold. The scale of the noise in the dimensionality estimate is not given; the noise of the trial average is assumed, as the text's raw differences give no dimensions at all. The kernel width (standard deviation or full width at half maximum) changes the fig3b medians and the fig6c significance, which the claims file did not expect, and no single reading gives the recorded outcome of both: the full-width reading would make fig3b reproduced and fig6c not reproduced. The standard-deviation reading is adopted for every claim, as the better supported; under it fig3b is partially reproduced and fig6c partially reproduced. Figure S4D uses 100 subsamples per size with 100 noise repetitions per subsample.
What was found about the dataset. It holds only the Baseline epoch of the sessions: about 200 reward trials each, with no perturbation block. The Force or Rotation and Washout epochs, and any label of a session as curl field or visuomotor rotation, are missing, so 23 of the 28 claims cannot be tested. Nothing marks which sessions are the paper's. In seven Monkey C sessions the session description ends in a cut-off warning, and the cursor ends reaches at about 4.2 cm while the targets are stored at 8 cm; in two Monkey C sessions the cursor timestamps repeat thousands of times; in several sessions the cursor never enters the stored target square. The pooled counts of units above 5 Hz are far below the paper's.
print(f"notebook ran in {(time.time() - started) / 60:.1f} min")
notebook ran in 11.3 min