Reproducing Guan, Aflalo et al. (2022) from dandiset 000147

Paper: Guan C, Aflalo T, Zhang CY, Amoruso E, Rosario ER, Pouratian N, Andersen RA (2022). Stability of motor representations after paralysis. eLife 11:e74478. doi:10.7554/eLife.74478. The paper is open access, and its full text, figure legends and methods were used (eLife, PMC9555862).

Data: DANDI:000147, version 0.221122.2256, streamed from the archive. The dandiset has two published versions, 0.220913.2243 and 0.221122.2256, with the same ten assets (same asset IDs and sizes). The later version is used because its metadata names this paper; the earlier one names the bioRxiv preprint. The authors' download script pins the earlier version, so the files are the ones the authors analysed.

What can and cannot be tested

The paper records from a Utah array in the posterior parietal cortex (PC-IP) of one tetraplegic participant, NS, over ten sessions on nine days, while she attempted to press the cued finger (thumb T, index I, middle M, ring R, pinky P) or did nothing (No-Go, X). The dandiset holds the main task of all ten sessions: sorted units, the cued finger, the online decoder's output, the index of the cue's grid position, the run and the trial number. It does not hold the calibration task, the no-movement control task, waveforms or sorting quality metrics. The fMRI and component model RDMs the paper compares against are not in the dandiset; the two fMRI models are in the authors' code repository and are read from there at the commit the paper archived.

Claim	Finding	Here
recording-counts	95.8 neurons and 408 trials per session	Tested
fig1b-online-accuracy	Online accuracy 86%	Tested
fig1b-errors-middle-ring-pinky	70% of online errors are middle-ring or ring-pinky	Tested
fig1s4a-offline-accuracy	Offline cross-validated accuracy 92%	Tested
fig1s4b-calibration-accuracy	Calibration task accuracy 96%	Blocked: the calibration task is not in the dandiset
fig1s5-tuned-fraction	89% of neurons tuned to a finger	Tested
fig1s6a-gaze-regression	Cue coordinates cannot be decoded by regression	Blocked: the files give a grid index without the grid's layout; run under two assumed layouts as a check
fig1s6b-gaze-classification	Cue location cannot be classified	Tested
fig1s6c-control-task	No finger decoding in the control task	Blocked: the control task is not in the dandiset
fig2d-rdm-structure	Thumb distinct; middle, ring and pinky close	Tested
fig2f-able-bodied-over-unstructured	RDMs match the able-bodied fMRI model better than the unstructured model	Tested, with the model from the authors' repository
fig2s5-mc-over-pcip	RDMs match the MC fMRI model better than the PC-IP model, near the noise ceiling	Tested, with the models from the authors' repository
fig2-typical-among-individuals	NS's RDM is typical among able-bodied participants	Blocked: the individual fMRI RDMs are not available
rdm-consistency	Session RDMs correlate at r = 0.90	Tested
fig3a-no-trend-to-unstructured	No trend toward the unstructured model over sessions	Tested
fig3b-confused-pairs-no-increase	Middle-ring and ring-pinky distances do not increase	Tested
fig3s1-all-pairs-small-decrease	All distances decrease slightly	Tested
fig4e-muscle-before-somatotopy	Muscle model emerges 170 ms before somatotopy	Blocked: the muscle and somatotopy model RDMs are not available
fig4s1-model-selection	Usage or muscle plus somatotopy fit best	Blocked: the component model RDMs are not available
fig4s4-well-isolated-units	Results hold for well-isolated units	Blocked: no L-ratio or waveforms in the dandiset

import os

os.environ.setdefault("OMP_NUM_THREADS", "1")
os.environ.setdefault("OPENBLAS_NUM_THREADS", "1")

import warnings

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import statsmodels.formula.api as smf
from scipy import stats
from sklearn.model_selection import StratifiedKFold, cross_val_predict
from statsmodels.stats.multitest import multipletests

import ppcfingers as pf

warnings.filterwarnings("ignore", category=FutureWarning)
plt.rcParams.update({"figure.dpi": 100, "axes.spines.top": False, "axes.spines.right": False})
pd.set_option("display.width", 160)
pd.set_option("display.max_columns", 30)

Load the data and check it against its description

pf.load_all streams the ten files and reduces each to spike counts per trial and unit, in the paper's window (0.5 to 1.5 s after the trial start) and in 50 ms bins from 1 s before to 3 s after it. The fMRI model RDMs are read from the authors' repository.

sessions = pf.load_all()
names = [s.name for s in sessions]
print(f"Dandiset {pf.DANDISET}, version {pf.VERSION}: {len(sessions)} files")
print("Subjects:", sorted({s.info["subject_id"] for s in sessions}))
print("Session descriptions:", sorted({s.info["session_description"] for s in sessions}))
print("session_id matches the file name in every file:", all(s.info["session_id"] == s.name for s in sessions))

overview = pd.DataFrame({
    "units": [s.info["n_units"] for s in sessions],
    "electrodes_with_units": [s.info["n_electrodes_with_units"] for s in sessions],
    "trials": [len(s.trials) for s in sessions],
    "runs": [s.trials.run.nunique() for s in sessions],
    "trials_per_run": [sorted(set(s.trials.groupby("run").size())) for s in sessions],
    "decoder_outputs": ["".join(sorted(set(s.trials.finger_pred))) for s in sessions],
    "cue_locations": [f"{s.trials.cue_location.min()}-{s.trials.cue_location.max()} ({s.trials.cue_location.nunique()})"
                      for s in sessions],
}, index=names)
print(overview.to_string())

print("\nTrials per cued condition:")
print(pd.DataFrame({s.name: s.trials.finger.value_counts() for s in sessions}).T[pf.CONDITIONS].to_string())

print("\nTiming and integrity, read from each file:")
keys = ["spike_time_range_s", "trial_time_range_s", "trial_durations_s", "n_obs_intervals",
        "obs_intervals_same_for_all_units", "trials_inside_obs_intervals", "runs_inside_epochs",
        "units_without_spikes", "units_with_unsorted_spike_times", "backward_jumps_in_spike_times",
        "trials_whose_range_meets_an_unsorted_span"]
print(pd.DataFrame({s.name: {k: s.info[k] for k in keys} for s in sessions}).T.to_string())

print("\nFiring rate per unit in the 0.5 to 1.5 s window (Hz), all trials:")
print(pd.DataFrame({s.name: {"median": s.info["window_rate_hz_median"], "range": s.info["window_rate_hz_range"]}
                    for s in sessions}).T.to_string())

print("\nTrials table column descriptions (the same in every file: "
      f"{all(s.info['trial_column_descriptions'] == sessions[0].info['trial_column_descriptions'] for s in sessions)}):")
for k, v in sessions[0].info["trial_column_descriptions"].items():
    print(f"  {k}: {v!r}")
print("Units table columns (the same in every file: "
      f"{all(s.info['units_columns'] == sessions[0].info['units_columns'] for s in sessions)}): {sessions[0].info['units_columns']}")

screens = pd.concat([s.trials.assign(session=s.name) for s in sessions]).groupby(["session", "run", "screen"])
screen_sizes = screens.size()
print(f"\nScreens (12 consecutive trials of a run): {len(screen_sizes)}, of {sorted(set(screen_sizes))} trials; "
      f"screens with all 12 cue locations once: {(screens.cue_location.nunique() == 12).sum()}; "
      f"with each condition twice: {screens.finger.apply(lambda f: (f.value_counts() == 2).all() and f.nunique() == 6).sum()}")
print("Screens per run:", sorted(set(pd.concat([s.trials.assign(session=s.name) for s in sessions])
                                     .groupby(["session", "run"]).screen.nunique())))

loaded = {name: pf.load_model_rdm(name) for name in ["Able-bodied (PC-IP)", "Able-bodied (MC)", "unstructured"]}
models = {name: vector for name, (vector, _) in loaded.items()}
print(f"\nModel RDMs from the authors' repository at commit {pf.MODEL_COMMIT[:10]}, pairs {pf.PAIR_NAMES}:")
for name, (vector, order) in loaded.items():
    print(f"  {name}: {np.round(vector, 3)}; condition order stored in the file: {''.join(order)}")

Dandiset 000147, version 0.221122.2256: 10 files
Subjects: ['P1']
Session descriptions: ['individual-finger press, main task']
session_id matches the file name in every file: True
             units  electrodes_with_units  trials  runs trials_per_run decoder_outputs cue_locations
2018-09-10     100                     74     384     8           [48]           IMPRT     1-12 (12)
2018-09-17a     93                     73     384     8           [48]          IMPRTX     1-12 (12)
2018-09-17b     88                     67     336     7           [48]          IMPRTX     1-12 (12)
2018-09-24     110                     75     480    10           [48]          IMPRTX     1-12 (12)
2018-09-26      95                     74     432     9           [48]          IMPRTX     1-12 (12)
2018-10-01      90                     68     432     9           [48]          IMPRTX     1-12 (12)
2018-10-12     100                     75     432     9           [48]          IMPRTX     1-12 (12)
2018-10-15      89                     69     432     9           [48]          IMPRTX     1-12 (12)
2018-10-17      99                     73     384     8           [48]          IMPRTX     1-12 (12)
2018-10-22      94                     72     384     8           [48]          IMPRTX     1-12 (12)

Trials per cued condition:
finger        T   I   M   R   P   X
2018-09-10   64  64  64  64  64  64
2018-09-17a  64  64  64  64  64  64
2018-09-17b  56  56  56  56  56  56
2018-09-24   80  80  80  80  80  80
2018-09-26   72  72  72  72  72  72
2018-10-01   72  72  72  72  72  72
2018-10-12   72  72  72  72  72  72
2018-10-15   72  72  72  72  72  72
2018-10-17   64  64  64  64  64  64
2018-10-22   64  64  64  64  64  64

Timing and integrity, read from each file:
             spike_time_range_s  trial_time_range_s trial_durations_s n_obs_intervals obs_intervals_same_for_all_units trials_inside_obs_intervals runs_inside_epochs units_without_spikes units_with_unsorted_spike_times backward_jumps_in_spike_times trials_whose_range_meets_an_unsorted_span
2018-09-10   (2284.01, 3913.72)  (2290.09, 3901.53)             [3.0]               8                             True                         384             8 of 8                    0                             100                           772                                         0
2018-09-17a   (589.82, 3136.46)   (596.27, 3122.68)             [3.0]               8                             True                         384             8 of 8                    0                              93                           707                                         0
2018-09-17b   (241.22, 1508.86)   (247.73, 1498.43)             [3.0]               7                             True                         336             7 of 7                    0                              87                           590                                         0
2018-09-24    (314.81, 2441.94)   (321.32, 2432.42)             [3.0]              10                             True                         480           10 of 10                    0                             110                          1016                                         0
2018-09-26    (241.72, 3904.21)   (247.77, 3894.67)             [3.0]               9                             True                         432             9 of 9                    0                              94                           838                                         0
2018-10-01      (83.6, 2747.63)     (90.8, 2730.44)             [3.0]               9                             True                         432             9 of 9                    0                              89                           714                                         0
2018-10-12     (158.07, 3068.5)   (161.35, 3057.81)             [3.0]               9                             True                         432             9 of 9                    0                             100                           767                                         0
2018-10-15     (569.8, 2801.85)   (576.24, 2793.53)             [3.0]               9                             True                         432             9 of 9                    0                              89                           747                                         0
2018-10-17      (64.2, 2716.71)    (70.11, 2660.19)             [3.0]               8                             True                         384             8 of 8                    0                              99                           781                                         0
2018-10-22     (99.63, 2692.08)   (105.49, 2680.49)             [3.0]               8                             True                         384             8 of 8                    0                              94                           693                                         0

Firing rate per unit in the 0.5 to 1.5 s window (Hz), all trials:
            median          range
2018-09-10    6.68  (0.253, 38.7)
2018-09-17a   3.79  (0.292, 24.0)
2018-09-17b   4.13  (0.333, 35.8)
2018-09-24    4.77  (0.113, 40.8)
2018-09-26    6.82  (0.315, 23.1)
2018-10-01     4.2  (0.345, 47.9)
2018-10-12    7.24   (0.12, 29.9)
2018-10-15    4.11  (0.947, 31.7)
2018-10-17    6.17  (0.643, 21.2)
2018-10-22    4.14  (0.529, 54.9)

Trials table column descriptions (the same in every file: True):
  Run: 'no description'
  TrialNumber: 'no description'
  start_time: 'Start time of epoch, in seconds'
  stop_time: 'Stop time of epoch, in seconds'
  finger: 'the finger cued during the trial'
  finger_pred: 'the finger predicted by the decoder during the trial'
  cue_location: 'screen location (index) where text cue was displayed'
Units table columns (the same in every file: True): ['spike_times', 'obs_intervals', 'electrodes', 'electrode_group']

Screens (12 consecutive trials of a run): 340, of [12] trials; screens with all 12 cue locations once: 340; with each condition twice: 340
Screens per run: [4]

Model RDMs from the authors' repository at commit 6915bf863a, pairs ['TI', 'TM', 'TR', 'TP', 'IM', 'IR', 'IP', 'MR', 'MP', 'RP']:
  Able-bodied (PC-IP): [0.241 0.256 0.248 0.246 0.032 0.106 0.122 0.064 0.149 0.089]; condition order stored in the file: TIMRP
  Able-bodied (MC): [0.374 0.448 0.513 0.467 0.214 0.353 0.377 0.2   0.294 0.18 ]; condition order stored in the file: TIMRP
  unstructured: [1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]; condition order stored in the file: TIMRP

Does the trial clock line up with the spikes?

The trials table gives a start time and a stop time 3 s later, with no description. The paper's window is defined from the cue. If the start time is the cue, the cued finger should be undecodable before it and become decodable a few hundred milliseconds after it, once the participant has looked at the cue and begun the movement. The cell below decodes the finger from 200 ms windows in each session separately (six classes, chance 1/6).

centers = np.round(np.arange(-0.9, 2.95, 0.1), 2)
curves = pd.DataFrame([pf.decoding_curve(s, centers) for s in sessions], index=names, columns=centers)

fig, ax = plt.subplots(figsize=(7, 3.2))
for name, row in curves.iterrows():
    ax.plot(centers, row.values, lw=1, alpha=0.7)
ax.axhline(1 / 6, color="k", ls="--", lw=1)
ax.axvline(0, color="0.6", lw=1)
ax.axvspan(*pf.WINDOW, color="0.92", zorder=0)
ax.set_xlabel("Window centre, s from trial start")
ax.set_ylabel("Finger accuracy (6 classes)")
plt.tight_layout()

pre = curves.loc[:, centers <= -0.1]
check = pd.DataFrame({
    "pre_start_mean": pre.mean(axis=1).round(3),
    "pre_start_max": pre.max(axis=1).round(3),
    "first_centre_above_0.4": [centers[np.argmax(r > 0.4)] for r in curves.values],
    "peak": curves.max(axis=1).round(3),
    "peak_centre": curves.idxmax(axis=1),
})
print(check.to_string())
print(f"\nPlotted accuracies range from {curves.values.min():.3f} to {curves.values.max():.3f}; none is off the axes.")

             pre_start_mean  pre_start_max  first_centre_above_0.4   peak  peak_centre
2018-09-10            0.176          0.214                     0.6  0.839          1.1
2018-09-17a           0.177          0.224                     0.6  0.833          1.1
2018-09-17b           0.152          0.202                     0.6  0.783          1.0
2018-09-24            0.182          0.210                     0.6  0.850          1.1
2018-09-26            0.173          0.204                     0.6  0.833          1.2
2018-10-01            0.192          0.225                     0.6  0.822          1.1
2018-10-12            0.189          0.248                     0.6  0.824          1.2
2018-10-15            0.182          0.215                     0.6  0.824          1.1
2018-10-17            0.168          0.201                     0.6  0.771          1.5
2018-10-22            0.192          0.219                     0.6  0.789          1.2

Plotted accuracies range from 0.119 to 0.850; none is off the axes.

Checks of the description against the data:

• Counts. Ten files from one subject, all described as the main task, with 336 to 480 trials in runs of 48 (eight of each condition) and 88 to 110 units. The six conditions are balanced in every session. These match the paper's design (Methods, "Finger flexion grid task").
• Subject. The files call the participant P1; the paper calls her NS. The dandiset's single subject and the ten session dates match the paper's sessions.
• Decoder output. In the first session (2018-09-10) the online decoder never output X, as the paper says. That session has 64 No-Go trials. The paper says 84, but its own total of 4016 scored trials is 4080 minus 64, so the data agree with the total and the 84 is a typo.
• Trial structure. Every run has 48 trials in four screens of 12 consecutive trials, and every one of the 340 screens cues each of the 12 locations once and each condition twice, as the methods describe (the methods say three to four screens per run; here every run has four).
• Clock. The trials table describes the start time only with the NWB default ("Start time of epoch, in seconds"), so the data have to show what it marks. In every session the finger is at chance before the trial start (mean accuracy 0.15 to 0.19 in windows ending at or before it, highest single window 0.25, against 1/6), first exceeds 0.4 in the window centred at 0.6 s, and peaks between 1.0 and 1.5 s, in all ten sessions. This is consistent with a start time at the cue onset and with spikes on the clock of the trials table, but it does not pin the start time to the cue: a shift of a few hundred milliseconds would look much the same.
• Observation intervals and epochs. Every unit has the same observation intervals, one per run; every trial's analysed range lies inside one, and every run lies inside its epoch.
• Spike times are not in ascending order within units. Within a unit, the spikes of the first few seconds of each run are stored before the rest of the previous run. No trial's analysed range meets a span where this happens, because the first trial of a run starts about 6 s into it, so sorting the spike times, as done here, gives the same counts. Tools that assume sorted spike times (for example a binary search over them) could still be misled near run boundaries.
• Cue location is described as the "screen location (index) where text cue was displayed", an index from 1 to 12 on the grid of three rows and four columns, with no map to the grid or to screen coordinates. The trials table's Run and TrialNumber columns read "no description".
• Units. The units table holds spike times, observation intervals, electrodes and electrode groups only: no waveforms and no sorting quality metrics.
• Model RDMs. Each of the three model files stores its condition order, and it is T, I, M, R, P in all three, the order used here.
• Electrodes. The location is given as PC-IP; coordinates and impedances are missing (NaN).

Claim recording-counts

"We recorded single and multineuron activity (95.8 ± standard deviation [SD] 6.7 neurons per session over 10 sessions)", and "408 ± SD 40.8 trials/session over 10 sessions".

units, trials = overview.units, overview.trials
print(f"Units per session: mean {units.mean():.1f}, SD {units.std(ddof=1):.1f}, range {units.min()} to {units.max()}, total {units.sum()}")
print(f"Trials per session: mean {trials.mean():.1f}, SD {trials.std(ddof=1):.1f}, range {trials.min()} to {trials.max()}, total {trials.sum()}")

Units per session: mean 95.8, SD 6.7, range 88 to 110, total 958
Trials per session: mean 408.0, SD 40.8, range 336 to 480, total 4080

Result: reproduced. The ten files hold 95.8 units per session (SD 6.7) and 408.0 trials per session (SD 40.8), 4080 trials in all, as the paper reports.

Claim fig1b-online-accuracy

"High classification accuracy during online control (86% ± SD 4% over 10 sessions; chance = 17%)", over 4016 trials. The No-Go trials of the first session are excluded, as the paper does, because the decoder could not output No-Go that day. The SD is across sessions, weighted by trials.

rng = np.random.default_rng(1471)
online = pd.concat([s.trials.assign(session=s.name) for s in sessions], ignore_index=True)
excluded = (online.session == names[0]) & (online.finger == "X")
scored = online[~excluded].copy()
scored["correct"] = scored.finger == scored.finger_pred

per_session = scored.groupby("session").agg(trials=("correct", "size"), accuracy=("correct", "mean"))
accuracy = scored.correct.mean()
print(per_session.round(3).to_string())
print(f"\nExcluded: {excluded.sum()} No-Go trials of {names[0]}. Scored: {len(scored)} trials.")
print(f"Online accuracy {accuracy:.4f}, weighted SD across sessions {pf.weighted_sd(per_session.accuracy, per_session.trials):.4f}")
print(f"Including the excluded trials: {(online.finger == online.finger_pred).mean():.4f}")

# Chance: the decoder's outputs shuffled across trials within each session.
def shuffled_accuracy():
    hits = sum((g.finger.to_numpy() == rng.permutation(g.finger_pred.to_numpy())).sum() for _, g in scored.groupby("session"))
    return hits / len(scored)
null = [shuffled_accuracy() for _ in range(500)]
print(f"Chance by shuffling the decoder's outputs (500 shuffles): mean {np.mean(null):.4f}, "
      f"95% range {np.percentile(null, 2.5):.4f} to {np.percentile(null, 97.5):.4f}; 1/6 = {1/6:.4f}")

confusion = pd.crosstab(scored.finger, scored.finger_pred).reindex(index=pf.CONDITIONS, columns=pf.CONDITIONS, fill_value=0)
fraction = confusion.div(confusion.sum(axis=1), axis=0)
fig, ax = plt.subplots(figsize=(3.8, 3.4))
im = ax.imshow(fraction.values, vmin=0, vmax=1, cmap="viridis")
ax.set_xticks(range(6), pf.CONDITIONS); ax.set_yticks(range(6), pf.CONDITIONS)
ax.set_xlabel("Decoded"); ax.set_ylabel("Cued"); plt.colorbar(im, fraction=0.046)
plt.tight_layout()
print("\nOnline confusion counts (rows cued, columns decoded):")
print(confusion.to_string())

             trials  accuracy
session                      
2018-09-10      320     0.838
2018-09-17a     384     0.805
2018-09-17b     336     0.851
2018-09-24      480     0.808
2018-09-26      432     0.852
2018-10-01      432     0.910
2018-10-12      432     0.880
2018-10-15      432     0.912
2018-10-17      384     0.859
2018-10-22      384     0.875

Excluded: 64 No-Go trials of 2018-09-10. Scored: 4016 trials.
Online accuracy 0.8596, weighted SD across sessions 0.0359
Including the excluded trials: 0.8461

Chance by shuffling the decoder's outputs (500 shuffles): mean 0.1694, 95% range 0.1580 to 0.1798; 1/6 = 0.1667

Online confusion counts (rows cued, columns decoded):
finger_pred    T    I    M    R    P    X
finger                                   
T            631   26    2    2    1   18
I              9  623   17    6    5   20
M              2   18  555   85    3   17
R              1    0   56  558   59    6
P              9    7    7   64  587    6
X             17   19   62   14    6  498

Result: reproduced. The decoder's recorded outputs match the cued finger on 85.96% of the 4016 scored trials, with a weighted SD across sessions of 3.6%, against the paper's 86% and 4%. Shuffling the outputs gives 16.9%, the paper's chance level of 17%.

Claim fig1b-errors-middle-ring-pinky

"The majority (70%) of the online classification errors were middle-ring or ring-pinky confusions." The paper does not say whether errors involving No-Go count in the denominator, so the fraction is shown for three denominators.

errors = scored[~scored.correct]
pair = errors.finger + errors.finger_pred
confused = pair.isin(["MR", "RM", "RP", "PR"])
involves_x = (errors.finger == "X") | (errors.finger_pred == "X")
variants = pd.DataFrame({
    "errors": [len(errors), int((errors.finger != "X").sum()), int((~involves_x).sum())],
    "middle_ring_or_ring_pinky": [int(confused.sum()), int(confused[errors.finger != "X"].sum()),
                                  int(confused[~involves_x].sum())],
}, index=["all errors", "errors on finger trials", "errors between two fingers"])
variants["fraction"] = (variants.middle_ring_or_ring_pinky / variants.errors).round(4)
print(variants.to_string())

by_session = errors[~involves_x].assign(confused=confused[~involves_x]).groupby("session").confused.agg(["size", "mean"])
print("\nBetween-finger errors by session:")
print(by_session.rename(columns={"size": "errors", "mean": "fraction"}).round(3).to_string())

                            errors  middle_ring_or_ring_pinky  fraction
all errors                     564                        264    0.4681
errors on finger trials        446                        264    0.5919
errors between two fingers     379                        264    0.6966

Between-finger errors by session:
             errors  fraction
session                      
2018-09-10       52     0.692
2018-09-17a      45     0.778
2018-09-17b      34     0.676
2018-09-24       58     0.638
2018-09-26       43     0.674
2018-10-01       26     0.769
2018-10-12       34     0.735
2018-10-15       17     0.824
2018-10-17       42     0.595
2018-10-22       28     0.714

Result: reproduced, for errors between two fingers. Of the 379 online errors in which both the cued and the decoded condition are fingers, 69.7% confuse middle with ring or ring with pinky, the paper's 70%. Counting errors that involve No-Go as well, the fraction is 46.8% of 564, and 59.2% of the 446 errors on finger trials. The paper's figure matches only the first denominator, which leaves out No-Go as Figure 2 does.

Claim fig1s4a-offline-accuracy

Offline cross-validated classification was 92% (SD 2%) over 4080 trials, with diagonal LDA on sorted firing rates in the 0.5 to 1.5 s window. The methods say leave-one-out cross-validation within each session. The authors' code (crossval_classify.yaml) uses eight stratified folds, drops units with a mean count below 0.5, z-scores, and fits LDA with identity covariance. Both are run here, with diagonal LDA on eight folds as a third variant.

rows = []
for s in sessions:
    X, y = s.counts.astype(float), s.trials.finger.to_numpy()
    authors = cross_val_predict(pf.authors_classifier(), X, y, cv=StratifiedKFold(8))
    rows.append({"session": s.name, "trials": len(y),
                 "diag_lda_leave_one_out": (pf.diag_lda_cv(X, y, None) == y).mean(),
                 "diag_lda_8_fold": (pf.diag_lda_cv(X, y, 8) == y).mean(),
                 "authors_pipeline_8_fold": (authors == y).mean()})
offline = pd.DataFrame(rows).set_index("session")
print(offline.round(3).to_string())
print()
for col in offline.columns[1:]:
    print(f"{col}: {np.average(offline[col], weights=offline.trials):.4f} "
          f"(weighted SD {pf.weighted_sd(offline[col], offline.trials):.4f}) over {offline.trials.sum()} trials")

             trials  diag_lda_leave_one_out  diag_lda_8_fold  authors_pipeline_8_fold
session                                                                              
2018-09-10      384                   0.945            0.940                    0.940
2018-09-17a     384                   0.951            0.938                    0.914
2018-09-17b     336                   0.917            0.926                    0.899
2018-09-24      480                   0.948            0.948                    0.929
2018-09-26      432                   0.944            0.944                    0.940
2018-10-01      432                   0.938            0.947                    0.935
2018-10-12      432                   0.917            0.917                    0.912
2018-10-15      432                   0.944            0.947                    0.944
2018-10-17      384                   0.898            0.898                    0.893
2018-10-22      384                   0.927            0.927                    0.927

diag_lda_leave_one_out: 0.9336 (weighted SD 0.0163) over 4080 trials
diag_lda_8_fold: 0.9338 (weighted SD 0.0153) over 4080 trials
authors_pipeline_8_fold: 0.9243 (weighted SD 0.0167) over 4080 trials

Result: reproduced. Diagonal LDA with leave-one-out cross-validation, as the methods describe, gives 93.4% (weighted SD 1.6%) over 4080 trials, and the classifier and eight folds of the authors' code give 92.4% (SD 1.7%), against the paper's 92% (SD 2%). The leave-one-out figure is 1.4 points above the paper's, within the tolerance of 2 points.

Claim fig1s5-tuned-fraction

"Most (89%) neurons were significantly tuned to individual finger-press movements (significance threshold: p < 0.05, FDR-corrected)": a two-tailed t-test of each finger's firing rate in the window against No-Go, neurons below 0.1 Hz excluded, neurons pooled across sessions. The paper does not say over which family the FDR correction runs. The main reading here corrects over all units and fingers pooled; the others are shown.

pvalues, kept, unit_session = [], 0, []
for s in sessions:
    p, keep = pf.tuning_pvalues(s)
    pvalues.append(p); kept += len(keep); unit_session += [s.name] * len(keep)
P = np.vstack(pvalues)
unit_session = np.array(unit_session)
total_units = int(overview.units.sum())
print(f"Units: {total_units}; below 0.1 Hz in the window and excluded: {total_units - kept}; tested: {kept}; "
      f"tests with an undefined p-value: {int(np.isnan(P).sum())}")
session_rate = np.concatenate([s.info["session_rate_hz"] for s in sessions])
print(f"Rate over the whole recording (all spikes over the observation intervals): min {session_rate.min():.3f} Hz, "
      f"units below 0.1 Hz by this measure: {(session_rate < pf.MIN_RATE_TUNING).sum()}")

def tuned(q):
    return (q < 0.05).any(axis=1)

pooled = multipletests(P.ravel(), method="fdr_bh")[1].reshape(P.shape)
per_session_q = np.empty_like(P)
for name in names:
    m = unit_session == name
    per_session_q[m] = multipletests(P[m].ravel(), method="fdr_bh")[1].reshape(P[m].shape)
per_unit_q = np.array([multipletests(p, method="fdr_bh")[1] for p in P])
variants = pd.Series({
    "FDR over all units and fingers pooled": tuned(pooled).mean(),
    "FDR within each session": tuned(per_session_q).mean(),
    "FDR within each unit (5 tests)": tuned(per_unit_q).mean(),
    "no correction": tuned(P).mean(),
})
print("\nFraction of units tuned to at least one finger:")
print(variants.round(4).to_string())
print("\nBy session (pooled FDR):")
print(pd.Series(tuned(pooled)).groupby(unit_session).mean().round(3).to_string())

Units: 958; below 0.1 Hz in the window and excluded: 0; tested: 958; tests with an undefined p-value: 0
Rate over the whole recording (all spikes over the observation intervals): min 0.147 Hz, units below 0.1 Hz by this measure: 0

Fraction of units tuned to at least one finger:
FDR over all units and fingers pooled    0.8789
FDR within each session                  0.8747
FDR within each unit (5 tests)           0.8392
no correction                            0.9165

By session (pooled FDR):
2018-09-10     0.940
2018-09-17a    0.828
2018-09-17b    0.864
2018-09-24     0.891
2018-09-26     0.926
2018-10-01     0.878
2018-10-12     0.850
2018-10-15     0.888
2018-10-17     0.869
2018-10-22     0.851

Result: reproduced. With FDR correction over all units and fingers pooled, 87.9% of the 958 units are tuned to at least one finger, against the paper's 89%. Correcting within each session gives 87.5%, within each unit 83.9%, and no correction 91.7%. The pooled and per-session readings are within 1.5 points of the paper. No unit fell below the 0.1 Hz cutoff, whether the rate is taken in the window over all trials (lowest 0.113 Hz) or over the whole recording (lowest 0.147 Hz).

Claim fig1s6a-gaze-regression (blocked)

"Cross-validated regression r2 values are close to 0 across sessions." The claim is blocked: the files give the cue location as an index from 1 to 12 without saying how it maps onto the grid of three rows and four columns. As a check of what the claim would give, the two coordinates are regressed on the spike counts in the 0.5 to 1.5 s window under the two plain orders of the index, along the rows and down the columns. Each session is fitted on its own with eight folds. The paper does not say whether the regression was regularised, so ordinary least squares and ridge regression (penalty chosen by inner cross-validation) are both shown. The outcome stays blocked because neither order is known to be the one used.

from sklearn.linear_model import LinearRegression, RidgeCV
from sklearn.model_selection import KFold
from sklearn.metrics import r2_score

rows = []
for s in sessions:
    X = s.counts.astype(float)
    for layout in ["row-major", "column-major"]:
        xy = pf.grid_coordinates(s.trials.cue_location.to_numpy(), layout)
        for model_name, model in [("OLS", LinearRegression()), ("ridge", RidgeCV(alphas=np.logspace(-1, 5, 13)))]:
            pred = cross_val_predict(model, X, xy, cv=KFold(8))
            rows.append({"session": s.name, "layout": layout, "model": model_name,
                         "r2_x": r2_score(xy[:, 0], pred[:, 0]), "r2_y": r2_score(xy[:, 1], pred[:, 1])})
gaze = pd.DataFrame(rows)
summary = gaze.groupby(["layout", "model"]).agg(
    r2_x_mean=("r2_x", "mean"), r2_x_max=("r2_x", "max"), r2_y_mean=("r2_y", "mean"), r2_y_max=("r2_y", "max"),
    sessions_both_at_or_below_0_05=("r2_x", lambda v: int(((v <= 0.05) & (gaze.loc[v.index, "r2_y"] <= 0.05)).sum())))
print(summary.round(4).to_string())
print("\nBy session, ridge:")
print(gaze[gaze.model == "ridge"].pivot(index="session", columns="layout", values=["r2_x", "r2_y"]).round(3).to_string())

                    r2_x_mean  r2_x_max  r2_y_mean  r2_y_max  sessions_both_at_or_below_0_05
layout       model                                                                          
column-major OLS      -0.2067   -0.1323    -0.3867   -0.2123                              10
             ridge     0.0436    0.0971    -0.0169    0.0037                               6
row-major    OLS      -0.1973   -0.1154    -0.2182   -0.1251                              10
             ridge     0.0404    0.0944     0.0441    0.1038                               5

By session, ridge:
                    r2_x                   r2_y          
layout      column-major row-major column-major row-major
session                                                  
2018-09-10         0.033     0.034       -0.024     0.033
2018-09-17a        0.042     0.044        0.002     0.064
2018-09-17b        0.054     0.012       -0.069     0.032
2018-09-24         0.006    -0.005        0.002     0.026
2018-09-26         0.039     0.077       -0.025     0.042
2018-10-01         0.006     0.024       -0.011    -0.013
2018-10-12         0.026     0.039       -0.009     0.018
2018-10-15         0.063     0.057       -0.005     0.043
2018-10-17         0.097     0.027       -0.033     0.093
2018-10-22         0.069     0.094        0.004     0.104

Result: blocked. Under either order of the index, ordinary least squares gives a cross-validated r2 below 0 for both coordinates in all ten sessions (session means -0.20 to -0.39), which would meet the criterion. With ridge regression the r2 reaches 0.097 for x and 0.104 for y, and both coordinates are at or below 0.05 in 5 (rows) or 6 (columns) of the ten sessions, which would not. The result depends on the regression as well as on the unknown layout, so the claim stays blocked.

Claim fig1s6b-gaze-classification

"A linear classifier (diagonal linear discriminant analysis [LDA]) could not classify the gaze location from neural activity during the attempted-movement period." The cue location (12 grid positions, chance 1/12) is classified from the same window and features as the fingers, with eight-fold cross-validation in each session. Two null distributions are built from 500 shuffles each: location labels shuffled across all trials, and shuffled only among trials of the same finger, which keeps any chance association between finger and location in the session.

rng = np.random.default_rng(1476)

def pooled_accuracy(label_sets):
    hits = sum((pf.diag_lda_cv(s.counts.astype(float), y, 8) == y).sum() for s, y in zip(sessions, label_sets))
    return hits / len(online)

locations = [s.trials.cue_location.to_numpy() for s in sessions]
fingers = [s.trials.finger.to_numpy() for s in sessions]
observed = pooled_accuracy(locations)
finger_accuracy = pooled_accuracy(fingers)
null_all = np.array([pooled_accuracy([rng.permutation(y) for y in locations]) for _ in range(500)])
null_within = np.array([pooled_accuracy([pf.shuffle_within(y, f, rng) for y, f in zip(locations, fingers)])
                        for _ in range(500)])
null_finger = np.array([pooled_accuracy([rng.permutation(y) for y in fingers]) for _ in range(100)])

authors_hits = sum((cross_val_predict(pf.authors_classifier(), s.counts.astype(float), y, cv=StratifiedKFold(8)) == y).sum()
                   for s, y in zip(sessions, locations))
residual_hits = 0
for s, y, f in zip(sessions, locations, fingers):
    X = s.counts.astype(float)
    for v in np.unique(f):
        X[f == v] -= X[f == v].mean(axis=0)
    residual_hits += (pf.diag_lda_cv(X, y, 8) == y).sum()

print(f"Cue location, diagonal LDA: {observed:.4f} (chance 1/12 = {1/12:.4f}), "
      f"binomial p = {stats.binomtest(round(observed * len(online)), len(online), 1/12, alternative='greater').pvalue:.2g}")
print(f"  null, shuffled across trials: 95% range {np.percentile(null_all, 2.5):.4f} to {np.percentile(null_all, 97.5):.4f}, max {null_all.max():.4f}")
print(f"  null, shuffled within finger: 95% range {np.percentile(null_within, 2.5):.4f} to {np.percentile(null_within, 97.5):.4f}, max {null_within.max():.4f}")
print(f"  shuffles at or above the observed accuracy: {(null_all >= observed).sum()} and {(null_within >= observed).sum()} of 500")
print(f"Cue location, the authors' pipeline: {authors_hits / len(online):.4f}")
print(f"Cue location, after removing each finger's mean activity: {residual_hits / len(online):.4f}")
print(f"Cued finger, same features and classifier: {finger_accuracy:.4f}; "
      f"shuffled 95% range {np.percentile(null_finger, 2.5):.4f} to {np.percentile(null_finger, 97.5):.4f}")

per_session = pd.DataFrame({
    "trials": [len(y) for y in locations],
    "location_accuracy": [(pf.diag_lda_cv(s.counts.astype(float), y, 8) == y).mean() for s, y in zip(sessions, locations)],
}, index=names)
per_session["binomial_p"] = [stats.binomtest(round(a * n), n, 1/12, alternative="greater").pvalue
                             for a, n in zip(per_session.location_accuracy, per_session.trials)]
print("\nBy session:")
print(per_session.round(4).to_string())

fig, ax = plt.subplots(figsize=(5, 2.8))
ax.hist(null_all, bins=30, alpha=0.6, label="shuffled across trials")
ax.hist(null_within, bins=30, alpha=0.6, label="shuffled within finger")
ax.axvline(observed, color="k", label="observed")
ax.set_xlabel("Cue-location accuracy (12 classes)"); ax.legend(frameon=False, fontsize=8)
plt.tight_layout()

# Trial structure. Each screen cues every location once, so location is tied to the screen's
# place in time only through the order within a screen. Is the order within a screen related to
# location, and does location stay decodable when no screen is split between training and test?
by_position = pd.crosstab(online.cue_location, pd.concat([s.trials.position_in_screen for s in sessions], ignore_index=True))
chi2 = stats.chi2_contingency(by_position)
print(f"\nCue location against position within the screen (12 x 12 counts, {by_position.values.min()} to {by_position.values.max()} per cell): "
      f"chi-square {chi2.statistic:.1f}, df {chi2.dof}, p = {chi2.pvalue:.3f}")
screen_ids = [s.trials.run.to_numpy() * 100 + s.trials.screen.to_numpy() for s in sessions]
run_ids = [s.trials.run.to_numpy() for s in sessions]
def grouped_accuracy(label_sets, groups):
    hits = sum((pf.diag_lda_cv(s.counts.astype(float), y, min(8, len(np.unique(g))), groups=g) == y).sum()
               for s, y, g in zip(sessions, label_sets, groups))
    return hits / len(online)
by_screen = grouped_accuracy(locations, screen_ids)
by_run = grouped_accuracy(locations, run_ids)
null_screen = np.array([grouped_accuracy([pf.shuffle_within(y, g, rng) for y, g in zip(locations, screen_ids)], screen_ids)
                        for _ in range(200)])
print(f"Cue location, eight folds that never split a screen: {by_screen:.4f}; that never split a run (leave one run out where a session has fewer than eight): {by_run:.4f}")
print(f"  null, locations shuffled within each screen, folds by screen (200 shuffles): 95% range "
      f"{np.percentile(null_screen, 2.5):.4f} to {np.percentile(null_screen, 97.5):.4f}, max {null_screen.max():.4f}; "
      f"at or above the observed: {(null_screen >= by_screen).sum()} of 200")

null_run = np.array([grouped_accuracy([pf.shuffle_within(y, g, rng) for y, g in zip(locations, run_ids)], run_ids)
                     for _ in range(200)])
print(f"  null, locations shuffled within each run, folds by run (200 shuffles): 95% range "
      f"{np.percentile(null_run, 2.5):.4f} to {np.percentile(null_run, 97.5):.4f}, max {null_run.max():.4f}; "
      f"at or above the observed: {(null_run >= by_run).sum()} of 200")

# Is the location signal a signal of time within the screen? Decode the position within the
# screen (12 classes) the same way, and decode location after removing each position's mean.
positions = [s.trials.position_in_screen.to_numpy() for s in sessions]
position_accuracy = pooled_accuracy(positions)
null_position = np.array([pooled_accuracy([pf.shuffle_within(y, g, rng) for y, g in zip(positions, screen_ids)])
                          for _ in range(200)])
position_removed_hits = 0
for s, y, q in zip(sessions, locations, positions):
    X = s.counts.astype(float)
    for v in np.unique(q):
        X[q == v] -= X[q == v].mean(axis=0)
    position_removed_hits += (pf.diag_lda_cv(X, y, 8) == y).sum()
print(f"\nPosition within the screen, diagonal LDA: {position_accuracy:.4f} (chance 1/12); null, positions shuffled "
      f"within each screen (200 shuffles): 95% range {np.percentile(null_position, 2.5):.4f} to "
      f"{np.percentile(null_position, 97.5):.4f}, max {null_position.max():.4f}; "
      f"at or above the observed: {(null_position >= position_accuracy).sum()} of 200")
print(f"Cue location, after removing each position's mean activity: {position_removed_hits / len(online):.4f}")

Cue location, diagonal LDA: 0.1154 (chance 1/12 = 0.0833), binomial p = 1e-12
  null, shuffled across trials: 95% range 0.0755 to 0.0944, max 0.0998
  null, shuffled within finger: 95% range 0.0805 to 0.0982, max 0.1022
  shuffles at or above the observed accuracy: 0 and 0 of 500
Cue location, the authors' pipeline: 0.1115
Cue location, after removing each finger's mean activity: 0.1174
Cued finger, same features and classifier: 0.9338; shuffled 95% range 0.1551 to 0.1809

By session:
             trials  location_accuracy  binomial_p
2018-09-10      384             0.0938      0.2542
2018-09-17a     384             0.1432      0.0001
2018-09-17b     336             0.1518      0.0000
2018-09-24      480             0.1000      0.1097
2018-09-26      432             0.1366      0.0001
2018-10-01      432             0.1019      0.0984
2018-10-12      432             0.1250      0.0020
2018-10-15      432             0.1042      0.0729
2018-10-17      384             0.1146      0.0203
2018-10-22      384             0.0911      0.3150

Cue location against position within the screen (12 x 12 counts, 13 to 44 per cell): chi-square 150.2, df 121, p = 0.037

Cue location, eight folds that never split a screen: 0.1110; that never split a run (leave one run out where a session has fewer than eight): 0.1130
  null, locations shuffled within each screen, folds by screen (200 shuffles): 95% range 0.0735 to 0.0941, max 0.0980; at or above the observed: 0 of 200

  null, locations shuffled within each run, folds by run (200 shuffles): 95% range 0.0745 to 0.0927, max 0.0978; at or above the observed: 0 of 200

Position within the screen, diagonal LDA: 0.1051 (chance 1/12); null, positions shuffled within each screen (200 shuffles): 95% range 0.0733 to 0.0924, max 0.0946; at or above the observed: 0 of 200
Cue location, after removing each position's mean activity: 0.1110

Result: not reproduced. The cue location is classified above chance: 11.5% against a chance level of 8.3% (binomial p = 1e-12), above all 500 shuffles of either kind (highest 10.0% and 10.2%). The authors' classifier gives 11.2%, and removing each finger's mean activity first gives 11.7%, so the location signal does not come from an association between finger and location. Five of the ten sessions are above chance at p < 0.05. With folds that never split a screen the accuracy is 11.1%, above all 200 shuffles of location within screens (highest 9.8%), and with folds that never split a run it is 11.3%, above all 200 shuffles of location within runs (highest 9.8%), so the signal does not come from slow drift between screens or runs. Location and position within the screen are weakly associated over all trials (chi-square p = 0.037), and the position within the screen can itself be decoded: 10.5%, above all 200 shuffles within screens (highest 9.5%). After removing each position's mean activity, location is still classified at 11.1%, above the highest of the shuffles of the unadjusted data (10.0%), so the location signal is not explained by position within the screen. The effect is small next to the finger signal in the same features (93.4%).

Claim fig2d-rdm-structure

The session-averaged RDM "shows that the thumb is distinct, while the middle, ring, and pinky fingers are close in neural space." Each session's RDM is the cross-validated Mahalanobis distance between fingers in the 0.5 to 1.5 s window, No-Go excluded, with five stratified folds.

rdms = np.stack([pf.crossnobis_rdm(s.counts, s.trials.finger.to_numpy()) for s in sessions])
mean_rdm = rdms.mean(axis=0)
print("Session-averaged RDM, max-scaled to 1:")
print(pf.square(mean_rdm / mean_rdm.max()).round(3).to_string())
print(f"\nRaw distances (per neuron): session means {rdms.mean(axis=1).min():.3f} to {rdms.mean(axis=1).max():.3f}; "
      f"all values {rdms.min():.3f} to {rdms.max():.3f}; negative values: {(rdms < 0).sum()}")

thumb = [i for i, p in enumerate(pf.PAIR_NAMES) if "T" in p]
mrp = [pf.PAIR_NAMES.index(p) for p in ["MR", "MP", "RP"]]
def holds(v):
    return v[mrp].max() < v[thumb].min()
print(f"\nLargest middle/ring/pinky distance {mean_rdm[mrp].max() / mean_rdm.max():.3f}, "
      f"smallest thumb distance {mean_rdm[thumb].min() / mean_rdm.max():.3f}: criterion holds: {holds(mean_rdm)}")
print(f"Sessions in which the criterion holds on their own RDM: {sum(holds(r) for r in rdms)} of {len(rdms)}")
print(f"Thumb's mean distance to the other fingers, by session: {np.round([r[thumb].mean() / r.max() for r in rdms], 2)}")

fig, ax = plt.subplots(figsize=(3.6, 3.2))
im = ax.imshow(pf.square(mean_rdm / mean_rdm.max()).values, vmin=0, vmax=1, cmap="viridis")
ax.set_xticks(range(5), pf.FINGERS); ax.set_yticks(range(5), pf.FINGERS); plt.colorbar(im, fraction=0.046)
ax.set_title("Mean RDM, max-scaled", fontsize=10)
plt.tight_layout()

Session-averaged RDM, max-scaled to 1:
       T      I      M      R      P
T  0.000  0.641  0.923  1.000  0.819
I  0.641  0.000  0.390  0.620  0.751
M  0.923  0.390  0.000  0.209  0.482
R  1.000  0.620  0.209  0.000  0.217
P  0.819  0.751  0.482  0.217  0.000

Raw distances (per neuron): session means 0.517 to 0.953; all values 0.154 to 1.728; negative values: 0

Largest middle/ring/pinky distance 0.482, smallest thumb distance 0.641: criterion holds: True
Sessions in which the criterion holds on their own RDM: 9 of 10
Thumb's mean distance to the other fingers, by session: [0.83 0.85 0.82 0.89 0.8  0.85 0.82 0.89 0.75 0.84]

Result: reproduced. In the session-averaged RDM (max-scaled), the distances among middle, ring and pinky are 0.21 to 0.48, and every distance involving the thumb is 0.64 or more. The pattern holds in 9 of the 10 sessions on their own.

Claim fig2f-able-bodied-over-unstructured

"The measured representational structure matched the able-bodied representational structure significantly over the unstructured model (p = 5.7 × 10–5, two-tailed t-test)", by whitened unbiased cosine similarity (WUC). Here WUC is rsatoolbox's cosine_cov, as in the authors' code, computed for each session's RDM. The paper's test is a t-test on 1000 bootstrap samples over sessions; here the sessions are resampled 1000 times, the t statistic is the mean difference over its bootstrap SD, and a paired t-test across sessions is shown as well.

rng = np.random.default_rng(1472)
similarity = pd.DataFrame({name: pf.wuc(vector, rdms) for name, vector in models.items()}, index=names)
ceiling_lower = np.array([pf.wuc(np.delete(rdms, i, axis=0).mean(axis=0), rdms[i])[0] for i in range(len(rdms))])
ceiling_upper = pf.wuc(rdms.mean(axis=0), rdms)
print(similarity.round(4).to_string())
print(f"\nMean WUC: " + ", ".join(f"{k} {v:.4f}" for k, v in similarity.mean().items()))
print(f"Noise ceiling: lower {ceiling_lower.mean():.4f}, upper {ceiling_upper.mean():.4f}")

boot_index = rng.integers(0, len(rdms), size=(1000, len(rdms)))
def compare(a, b):
    diff = (similarity[a] - similarity[b]).to_numpy()
    boot = diff[boot_index].mean(axis=1)
    t = diff.mean() / boot.std(ddof=1)
    return {"mean_difference": diff.mean(), "boot_ci_low": np.percentile(boot, 2.5), "boot_ci_high": np.percentile(boot, 97.5),
            "t_bootstrap": t, "p_bootstrap": 2 * stats.t.sf(abs(t), len(diff) - 1),
            "p_paired_t": stats.ttest_rel(similarity[a], similarity[b]).pvalue,
            "sessions_a_higher": int((diff > 0).sum())}
result = pd.Series(compare("Able-bodied (PC-IP)", "unstructured"))
print("\nAble-bodied (PC-IP) minus unstructured:")
print(result.to_string(float_format=lambda v: f"{v:.3g}"))

fig, ax = plt.subplots(figsize=(4.4, 3))
for i, name in enumerate(similarity.columns):
    ax.scatter(np.full(len(names), i) + rng.uniform(-0.08, 0.08, len(names)), similarity[name], s=14)
ax.axhspan(ceiling_lower.mean(), ceiling_upper.mean(), color="0.85", zorder=0)
ax.set_xticks(range(3), ["PC-IP fMRI", "MC fMRI", "unstructured"]); ax.set_ylabel("WUC per session")
plt.tight_layout()
print(f"\nPlotted WUC values range from {similarity.values.min():.3f} to {similarity.values.max():.3f}.")

             Able-bodied (PC-IP)  Able-bodied (MC)  unstructured
2018-09-10                0.9001            0.9754        0.7847
2018-09-17a               0.8923            0.9769        0.7942
2018-09-17b               0.8819            0.9693        0.7991
2018-09-24                0.9567            0.9597        0.7800
2018-09-26                0.9051            0.9802        0.7914
2018-10-01                0.9144            0.9575        0.7514
2018-10-12                0.8291            0.9647        0.7878
2018-10-15                0.8527            0.9776        0.8050
2018-10-17                0.8391            0.9639        0.7468
2018-10-22                0.8344            0.9604        0.7679

Mean WUC: Able-bodied (PC-IP) 0.8806, Able-bodied (MC) 0.9686, unstructured 0.7808
Noise ceiling: lower 0.9838, upper 0.9869

Able-bodied (PC-IP) minus unstructured:
mean_difference       0.0998
boot_ci_low           0.0744
boot_ci_high           0.127
t_bootstrap             7.34
p_bootstrap         4.35e-05
p_paired_t          5.86e-05
sessions_a_higher         10

Plotted WUC values range from 0.747 to 0.980.

Result: reproduced. The mean WUC is 0.881 with the able-bodied PC-IP model and 0.781 with the unstructured model. The PC-IP model is higher in all 10 sessions, the bootstrap 95% interval of the difference is 0.074 to 0.127, and the bootstrap t-test gives p = 4.4e-5, against the paper's 5.7e-5.

Claim fig2s5-mc-over-pcip

"Measured RDMs match the able-bodied MC fMRI RDM better than they match the able-bodied PC-IP fMRI RDM (p = 1.9 × 10–6)", and "the WUC similarity with the MC RDM was close to the noise ceiling."

result = pd.Series(compare("Able-bodied (MC)", "Able-bodied (PC-IP)"))
print("Able-bodied (MC) minus able-bodied (PC-IP):")
print(result.to_string(float_format=lambda v: f"{v:.3g}"))
mc = similarity["Able-bodied (MC)"]
print(f"\nMean WUC with MC {mc.mean():.4f}; lower noise ceiling {ceiling_lower.mean():.4f}; "
      f"difference {mc.mean() - ceiling_lower.mean():+.4f}")
print(f"Sessions where WUC with MC is at or above that session's lower ceiling: {(mc.to_numpy() >= ceiling_lower).sum()} of {len(mc)}")

Able-bodied (MC) minus able-bodied (PC-IP):
mean_difference        0.088
boot_ci_low           0.0603
boot_ci_high           0.112
t_bootstrap              6.8
p_bootstrap         7.92e-05
p_paired_t          9.96e-05
sessions_a_higher         10

Mean WUC with MC 0.9686; lower noise ceiling 0.9838; difference -0.0153
Sessions where WUC with MC is at or above that session's lower ceiling: 1 of 10

Result: reproduced. The mean WUC with the MC model is 0.969, higher than with the PC-IP model in all 10 sessions (difference 0.088, bootstrap 95% interval 0.060 to 0.112, p = 7.9e-5; the paper gives 1.9e-6). It is 0.015 below the lower noise ceiling of 0.984, inside the criterion's margin of 0.05, and at or above its own session's ceiling in 1 of 10 sessions.

Claim rdm-consistency

"Participant NS's finger RDMs were largely consistent across different recording sessions (average pairwise correlation, excluding the diagonal: r = 0.90 ± SD 0.04, min 0.83, max 0.99)." The Pearson correlation is taken between the ten distances of each pair of sessions.

r = np.corrcoef(rdms)
upper = r[np.triu_indices(len(rdms), 1)]
print(f"{len(upper)} session pairs: mean r {upper.mean():.4f}, SD {upper.std(ddof=1):.4f}, min {upper.min():.4f}, max {upper.max():.4f}")
spearman = stats.spearmanr(rdms.T).statistic[np.triu_indices(len(rdms), 1)]
print(f"Spearman, for comparison: mean {spearman.mean():.4f}, min {spearman.min():.4f}, max {spearman.max():.4f}")
print("\nMean r of each session with the others:")
print(pd.Series((r.sum(axis=1) - 1) / (len(rdms) - 1), index=names).round(3).to_string())

45 session pairs: mean r 0.9282, SD 0.0409, min 0.8172, max 0.9898
Spearman, for comparison: mean 0.9168, min 0.8061, max 1.0000

Mean r of each session with the others:
2018-09-10     0.953
2018-09-17a    0.952
2018-09-17b    0.936
2018-09-24     0.893
2018-09-26     0.939
2018-10-01     0.927
2018-10-12     0.921
2018-10-15     0.913
2018-10-17     0.935
2018-10-22     0.914

Result: reproduced. The 45 pairs of session RDMs correlate at a mean r of 0.928 (SD 0.041, min 0.817, max 0.990), against the paper's 0.90 (SD 0.04, min 0.83, max 0.99). The mean is 0.028 above the paper's, inside the tolerance of 0.03.

Claim fig3a-no-trend-to-unstructured

"We did not find conclusive evidence for a trend from the able-bodied model toward the unstructured model (linear-model session × model interaction: t(6) = 0.50, one-tailed t-test p = 0.32)." The authors' code for this panel uses the PC-IP fMRI model and the unstructured model. Here each session's WUC with the two models is regressed on session index, model and their interaction. A positive interaction means the unstructured model gains on the able-bodied one over sessions.

long = similarity[["Able-bodied (PC-IP)", "unstructured"]].reset_index(drop=True)
long["session"] = np.arange(len(long))
long = long.melt(id_vars="session", var_name="model", value_name="wuc")
long["unstructured"] = (long.model == "unstructured").astype(int)
fit = smf.ols("wuc ~ session * unstructured", long).fit()
t, df = fit.tvalues["session:unstructured"], fit.df_resid
print(fit.summary().tables[1])
print(f"\nInteraction t({df:.0f}) = {t:.3f}, one-tailed p (toward unstructured) = {stats.t.sf(t, df):.3f}")
for name in ["Able-bodied (PC-IP)", "unstructured", "Able-bodied (MC)"]:
    slope = stats.linregress(np.arange(len(names)), similarity[name])
    print(f"Slope of WUC with {name} over sessions: {slope.slope:+.4f} per session, two-tailed p = {slope.pvalue:.3f}")

# Other specifications, since the paper's t(6) is not explained: the difference between the two
# models' WUC regressed on session (8 residual df), with each able-bodied model.
for able in ["Able-bodied (PC-IP)", "Able-bodied (MC)"]:
    gap = stats.linregress(np.arange(len(names)), similarity["unstructured"] - similarity[able])
    t = gap.slope / gap.stderr
    print(f"Unstructured minus {able} on session: slope {gap.slope:+.4f}, t({len(names) - 2}) = {t:.3f}, "
          f"one-tailed p = {stats.t.sf(t, len(names) - 2):.3f}")
t6 = fit.tvalues["session:unstructured"]
print(f"The interaction t of the main fit read with 6 df: one-tailed p = {stats.t.sf(t6, 6):.3f}")

========================================================================================
                           coef    std err          t      P>|t|      [0.025      0.975]
----------------------------------------------------------------------------------------
Intercept                0.9210      0.016     58.547      0.000       0.888       0.954
session                 -0.0090      0.003     -3.050      0.008      -0.015      -0.003
unstructured            -0.1274      0.022     -5.726      0.000      -0.175      -0.080
session:unstructured     0.0061      0.004      1.473      0.160      -0.003       0.015
========================================================================================

Interaction t(16) = 1.473, one-tailed p (toward unstructured) = 0.080
Slope of WUC with Able-bodied (PC-IP) over sessions: -0.0090 per session, two-tailed p = 0.038
Slope of WUC with unstructured over sessions: -0.0028 per session, two-tailed p = 0.204
Slope of WUC with Able-bodied (MC) over sessions: -0.0012 per session, two-tailed p = 0.228
Unstructured minus Able-bodied (PC-IP) on session: slope +0.0061, t(8) = 1.296, one-tailed p = 0.116
Unstructured minus Able-bodied (MC) on session: slope -0.0017, t(8) = -1.035, one-tailed p = 0.835
The interaction t of the main fit read with 6 df: one-tailed p = 0.096

Result: reproduced. The session by model interaction is t(16) = 1.47, one-tailed p = 0.080, so there is no significant trend toward the unstructured model, as the paper reports (t(6) = 0.50, p = 0.32). The same t read with 6 degrees of freedom gives p = 0.096, and regressing the difference between the unstructured and PC-IP fits on session gives t(8) = 1.30, p = 0.116; with the MC model in place of PC-IP the slope is negative (p = 0.84). The statistic is three times the paper's, and the fit to the PC-IP model on its own does decline across sessions (slope -0.009 per session, p = 0.038), while the slopes for the unstructured and MC models are not significant. The outcome rests on specifications chosen here, since the paper's is not described, and a specification closer to the authors' could give a different p-value.

Claim fig3b-confused-pairs-no-increase

"The middle-ring and ring-pinky distances did not increase over the course of the experiment (across sessions: t(8) = −4.5, one-tailed t-test p > 0.99, BF = 0.03; across runs within sessions: t(82) = −0.45, one-tailed t-test p = 0.67, BF = 0.12)." Following the methods, an RDM is computed for each run of each session and its distances are regressed on finger pair, session index and run index. As in the paper, the degrees of freedom of each test are the number of RDMs it rests on minus two: sessions for the session term, runs for the run term. Three variants are shown: the session term from the ten session RDMs; the run term with a separate intercept per session; and both terms from per-run RDMs whose noise covariance is estimated from the whole session. The last variant is there because a run's five folds hold only one or two trials of each finger, so a fold's own noise covariance rests on few residuals and inflates the per-run distances. Bayes factors are not computed.

runs = pd.concat([pf.run_rdms(s).assign(session=i) for i, s in enumerate(sessions)], ignore_index=True)
runs_shared = pd.concat([pf.run_rdms(s, session_noise=True).assign(session=i) for i, s in enumerate(sessions)],
                        ignore_index=True)
session_level = pd.DataFrame([{"session": i, "pair": p, "distance": d}
                              for i, v in enumerate(rdms) for p, d in zip(pf.PAIR_NAMES, v)])
n_runs = runs.groupby(["session", "run"]).ngroups
print(f"Per-run RDMs: {n_runs} (runs per session {runs.groupby('session').run.nunique().tolist()}); "
      f"per-run distances {runs.distance.min():.3f} to {runs.distance.max():.3f}, "
      f"mean {runs.distance.mean():.3f}; session RDM mean {rdms.mean():.3f}")
print(f"Per-run distances with the session's noise covariance: {runs_shared.distance.min():.3f} to "
      f"{runs_shared.distance.max():.3f}, mean {runs_shared.distance.mean():.3f}")

def trend(pairs):
    sub = runs[runs.pair.isin(pairs)]
    full = smf.ols("distance ~ C(pair) + session + run", sub).fit()
    base = smf.ols("distance ~ C(pair)", sub).fit()
    out = {}
    for term, df in [("session", len(sessions) - 2), ("run", n_runs - 2)]:
        t = full.tvalues[term]
        one = smf.ols(f"distance ~ C(pair) + {term}", sub).fit()
        out[term] = {"coefficient": full.params[term], "t": t, "df": df,
                     "p_increase_one_tailed": stats.t.sf(t, df), "p_two_tailed": 2 * stats.t.sf(abs(t), df),
                     "cohen_f2": (one.rsquared - base.rsquared) / (1 - one.rsquared)}
    s_sub = session_level[session_level.pair.isin(pairs)]
    s_fit = smf.ols("distance ~ C(pair) + session", s_sub).fit()
    t = s_fit.tvalues["session"]
    out["session (session RDMs)"] = {"coefficient": s_fit.params["session"], "t": t, "df": len(sessions) - 2,
                                     "p_increase_one_tailed": stats.t.sf(t, len(sessions) - 2),
                                     "p_two_tailed": 2 * stats.t.sf(abs(t), len(sessions) - 2), "cohen_f2": np.nan}
    fe = smf.ols("distance ~ C(pair) + C(session) + run", sub).fit()
    t = fe.tvalues["run"]
    out["run (intercept per session)"] = {"coefficient": fe.params["run"], "t": t, "df": n_runs - 2,
                                          "p_increase_one_tailed": stats.t.sf(t, n_runs - 2),
                                          "p_two_tailed": 2 * stats.t.sf(abs(t), n_runs - 2), "cohen_f2": np.nan}
    sub_shared = runs_shared[runs_shared.pair.isin(pairs)]
    shared = smf.ols("distance ~ C(pair) + session + run", sub_shared).fit()
    base_shared = smf.ols("distance ~ C(pair)", sub_shared).fit()
    for term, df in [("session", len(sessions) - 2), ("run", n_runs - 2)]:
        t = shared.tvalues[term]
        one = smf.ols(f"distance ~ C(pair) + {term}", sub_shared).fit()
        out[f"{term} (session noise covariance)"] = {
            "coefficient": shared.params[term], "t": t, "df": df, "p_increase_one_tailed": stats.t.sf(t, df),
            "p_two_tailed": 2 * stats.t.sf(abs(t), df),
            "cohen_f2": (one.rsquared - base_shared.rsquared) / (1 - one.rsquared)}
    return pd.DataFrame(out).T

confused_trend = trend(["MR", "RP"])
print("\nMiddle-ring and ring-pinky distances:")
print(confused_trend.to_string(float_format=lambda v: f"{v:.4g}"))

fig, axes = plt.subplots(1, 2, figsize=(8, 3))
sub = runs[runs.pair.isin(["MR", "RP"])]
by_session = sub.groupby(["session", "pair"]).distance.mean().unstack()
by_run = sub.groupby(["run", "pair"]).distance.mean().unstack()
by_session.plot(ax=axes[0], marker="o"); axes[0].set_xlabel("Session index"); axes[0].set_ylabel("Mean per-run distance")
by_run.plot(ax=axes[1], marker="o"); axes[1].set_xlabel("Run index within session")
plt.tight_layout()
print(f"\nPlotted means range from {min(by_session.values.min(), by_run.values.min()):.3f} "
      f"to {max(by_session.values.max(), by_run.values.max()):.3f}.")

Per-run RDMs: 85 (runs per session [8, 8, 7, 10, 9, 9, 9, 9, 8, 8]); per-run distances 0.075 to 16.137, mean 5.060; session RDM mean 0.679
Per-run distances with the session's noise covariance: 0.033 to 2.128, mean 0.756

Middle-ring and ring-pinky distances:
                                    coefficient       t  df  p_increase_one_tailed  p_two_tailed  cohen_f2
session                                -0.07612  -3.557   8                 0.9963      0.007437   0.07653
run                                    -0.02524  -1.054  83                 0.8524        0.2951  0.006995
session (session RDMs)                -0.006456  -1.853   8                 0.9495        0.1009       NaN
run (intercept per session)            -0.03268   -1.36  83                 0.9113        0.1775       NaN
session (session noise covariance)    -0.001789 -0.6358   8                 0.7287        0.5427   0.00261
run (session noise covariance)        -0.004667  -1.481  83                 0.9288        0.1424   0.01339

Plotted means range from 0.933 to 2.877.

Result: reproduced. In every variant the session and run coefficients for the middle-ring and ring-pinky distances are negative, and the one-tailed p for an increase is 0.73 or more. In the main regression the session term is t(8) = -3.56 (paper -4.5) and the run term t(83) = -1.05 (paper -0.45).

Claim fig3s1-all-pairs-small-decrease

"Inter-finger distances instead decreased slightly (across sessions: t(8) = −4.0, two tailed t-test p = 0.004; across runs within sessions: t(82) = −2.4, two-tailed t-test p = 0.019), although the effect size was very small (across sessions: Cohen's f2 = 0.008; across runs within sessions: = 0.005)." The same regression as above, over all ten finger pairs. Cohen's f2 compares a model with finger pair and the one predictor against finger pair alone.

all_trend = trend(pf.PAIR_NAMES)
print("All ten finger pairs:")
print(all_trend.to_string(float_format=lambda v: f"{v:.4g}"))

All ten finger pairs:
                                    coefficient      t  df  p_increase_one_tailed  p_two_tailed  cohen_f2
session                                -0.04853 -2.276   8                 0.9738       0.05242  0.006417
run                                    -0.07453 -3.122  83                 0.9988      0.002476   0.01187
session (session RDMs)                -0.008941   -1.5   8                  0.914        0.1719       NaN
run (intercept per session)            -0.09857 -4.381  83                      1     3.423e-05       NaN
session (session noise covariance)     -0.01055  -3.76   8                 0.9972      0.005544   0.01719
run (session noise covariance)         -0.01033 -3.289  83                 0.9993      0.001476   0.01323

Result: reproduced. Over all ten pairs the session and run coefficients are negative in every variant, and Cohen's f2 is 0.006 for session and 0.012 for run in the main regression and 0.017 and 0.013 with the session's noise covariance (paper 0.008 and 0.005). In the main regression the session decrease is t(8) = -2.28, p = 0.052, and the run decrease t(83) = -3.12, p = 0.0025. With the noise covariance taken from the whole session, they are t(8) = -3.76, p = 0.0055, and t(83) = -3.29, p = 0.0015, close to the paper's t(8) = -4.0, p = 0.004. The per-run distances of the main regression are inflated (mean 5.06 against 0.76 with the session's covariance), so the variant with the session's covariance is the better estimate of the session term; the criterion, negative coefficients with small effect sizes, holds in both.

Summary

A failure to reproduce here may be a failure of this attempt, not of the paper: see the reasoning and the differences from the paper's method.

Claim	Outcome	Note
recording-counts	Reproduced	95.8 units and 408 trials per session
fig1b-online-accuracy	Reproduced	85.96% of 4016 trials
fig1b-errors-middle-ring-pinky	Reproduced	69.7% of errors between two fingers
fig1s4a-offline-accuracy	Reproduced	93.4% leave-one-out, 92.4% with the authors' pipeline
fig1s4b-calibration-accuracy	Blocked: data missing	Calibration task not in the dandiset
fig1s5-tuned-fraction	Reproduced	87.9% of 958 units
fig1s6a-gaze-regression	Blocked: metadata insufficient	No map from the cue index to the grid; under assumed layouts the result depends on the regression
fig1s6b-gaze-classification	Not reproduced	Location decoded at 11.5% against 8.3% chance, 11.1% with folds by screen or after removing position within the screen
fig1s6c-control-task	Blocked: data missing	Control task not in the dandiset
fig2d-rdm-structure	Reproduced	Thumb distinct; middle, ring and pinky close
fig2f-able-bodied-over-unstructured	Reproduced	WUC 0.881 against 0.781
fig2s5-mc-over-pcip	Reproduced	WUC 0.969 with MC, ceiling 0.984
fig2-typical-among-individuals	Blocked: data missing	Individual fMRI RDMs not available
rdm-consistency	Reproduced	Mean r 0.928
fig3a-no-trend-to-unstructured	Reproduced	Interaction one-tailed p = 0.080; 0.096 to 0.84 in other specifications
fig3b-confused-pairs-no-increase	Reproduced	Session and run terms negative
fig3s1-all-pairs-small-decrease	Reproduced	Negative, f2 at most 0.017
fig4e-muscle-before-somatotopy	Blocked: data missing	Muscle and somatotopy model RDMs not available
fig4s1-model-selection	Blocked: data missing	Component model RDMs not available
fig4s4-well-isolated-units	Blocked: data missing	No L-ratio or waveforms

How this analysis differs from the paper's:

• The fMRI model RDMs come from the authors' code repository, not from the dandiset.
• Offline classification is shown with leave-one-out cross-validation (the paper's methods) and with the authors' code (eight folds, a firing cutoff, z-scoring); they differ by one point.
• The noise covariance of the crossnobis distance uses OAS shrinkage, as the authors' code does; the paper names Ledoit-Wolf.
• The fraction of errors that are middle-ring or ring-pinky depends on whether errors involving No-Go count; the paper does not say, and only the count without them matches.
• The FDR family for tuning is not stated; the main reading pools all units and fingers.
• WUC p-values come from resampling sessions here; the paper used rsatoolbox's bootstrap t-test. Bayes factors are not computed.
• The paper's code for the session and run regressions is not in its repository. The degrees of freedom follow the paper's stated rule, and the per-run RDMs use the same five-fold estimate as the session RDMs, with a variant using the session's noise covariance. A run's folds hold one or two trials per finger, so the fold covariance is poorly estimated and inflates the per-run distances; the significance of the all-pairs session decrease depends on this choice.
• The paper's t(6) for the session by model interaction is not explained; three specifications are shown.
• The gaze regression is run under two assumed orders of the grid index, with and without a ridge penalty; the claim stays blocked.
• The 12-class cue-location classifier uses eight folds; the paper does not say how it cross-validated it.

What was found about the dandiset:

• It holds the main task of all ten sessions and matches the paper's counts exactly. The calibration task, the control task, waveforms and sorting quality metrics are not included, and the model RDMs the main comparisons rest on are not in it.
• Spike times are not in ascending order within units: the first seconds of each run are stored before the rest of the previous run. No analysed trial range is affected once they are sorted.
• The cue location is described as a screen location index, with no map to the grid or to screen coordinates. The trials table's Run and TrialNumber columns have no description, and the start time has only the NWB default description, though the decoding of the finger over time is consistent with the cue onset. The participant is P1 in the files and NS in the paper. Electrode coordinates are missing.
• The first session has 64 No-Go trials; the paper's text says 84, but its own total of 4016 scored trials implies 64.
• The two published versions hold the same files. The earlier one, which the authors' download script pins, names the preprint; the later one names the eLife paper.