Reproducing Athalye et al. (2023) from dandiset 000404¶
Paper. Athalye VR, Khanna P, Gowda S, Orsborn AL, Costa RM, Carmena JM (2023). Invariant neural dynamics drive commands to control different movements. Current Biology 33, 2962-2976. doi:10.1016/j.cub.2023.06.027. The claims file (claims/000404/10.1016-j.cub.2023.06.027.yaml) was written from the accepted author manuscript in PubMed Central (PMC10527529) with its supplement; this attempt read the same text.
Dandiset. DANDI:000404, version 0.230605.2024. This is the version the paper's Data availability cites (doi:10.48324/dandi.000404/0.230605.2024), and the only published version. The draft holds the same 13 files (compared below).
What the dandiset holds. 13 NWB files, 9 sessions of Monkey G and 4 of Monkey J: per task frame (60 Hz for G, 200 Hz for J) the spike counts of the decoder's units, the decoder state, cursor, target and obstacle series, the task entry, a flag for decoder updates, and a trials table (go cue to reward). It holds no decoder matrices, although the paper says decoder parameters were deposited. The decoder gain is recovered here from the data by least squares (exactly for Monkey G's Kalman filter), and the command is computed from it. The target series of Monkey J's files is wrong (the same target for every trial), so targets are read from the cursor position at reward.
Paper's code. https://github.com/pkhanna104/bmi_dynamics_code at commit 416187366f84f2b7367ac54205aa321480e1acd1, as the claims file cites it. One file of it is read here: resim_ppf/jeev_obs_positions_from_amy.mat, the positions of Monkey J's obstacle-task targets, used to recover the code's obstacle trial types (an outside input, listed in the attempt record).
| Claim | Tested here |
|---|---|
| conditions-per-session | tested |
| figs1b-conditions-per-command | tested |
| figs1d-command-subtrajectories-differ | tested |
| figs1e-next-commands-differ | tested |
| fig3d-population-distance | tested |
| fig3b-single-neuron-distance | tested |
| fig3e-fractions-significant | tested |
| fig4c-dynamics-predict-activity-given-command | tested |
| fig4de-condition-specific-activity-predicted | tested |
| fig4f-condition-specific-component | tested |
| fig4g-majority-predicted | tested |
| figs3-low-dimensional-decaying-dynamics | tested |
| fig5c-next-activity | tested |
| fig5d-four-dimensions | tested |
| fig5e-next-command | tested |
| fig5f-condition-specific-next-command | tested |
| fig5h-turn-predicted | tested |
| figs5b-dynamics-beyond-behavior | tested |
| figs5f-nonlinear-no-better | blocked, not runnable: the rsLDS package (ssm) cannot be installed in the entry's locked environment (see the attempt record) |
| figs3f-residual-early-late | tested |
| figs4b-generalize-command-sets | tested |
| figs4g-generalize-condition-classes | tested |
| figs6cd-pairwise-distances-predicted | tested |
| figs6g-distance-tracks-behavior | tested |
| figs6k-next-command-pairs | tested |
| fig6c-dynamics-reduce-input | tested in ofc_simulations.ipynb |
| fig6d-decoder-null-dynamics-no-saving | tested in ofc_simulations.ipynb |
| fig6g-ofc-uses-distinct-patterns | tested in ofc_simulations.ipynb |
| monkey-j-linear-command | blocked, data missing: the files hold neither the PPF's update nor its parameters, which the criterion needs; the estimated readings are shown |
In [1]:
import os
os.environ["OPENBLAS_NUM_THREADS"] = "1"
os.environ["OMP_NUM_THREADS"] = "1"
os.environ["MKL_NUM_THREADS"] = "1"
import time
import warnings
from concurrent.futures import ProcessPoolExecutor, ThreadPoolExecutor
import numpy as np
import pandas as pd
import requests
import scipy.stats
import statsmodels.formula.api as smf
import analyses as an
import bmi
warnings.filterwarnings("ignore")
pd.set_option("display.width", 200)
pd.set_option("display.max_columns", 40)
pd.set_option("display.max_rows", 300)
NSHUFF = 1000 # Methods, "Behavior-preserving shuffle of activity": 1000 shuffled datasets
started = time.time()
In [2]:
sessions = bmi.session_list()
amy = bmi.load_amy()
def load(x):
m, s, path, _, _ = x
d = bmi.load_raw(path)
dec = bmi.fit_decoder(d, m)
B = bmi.bin_session(d, m, dec["K"], amy)
B = bmi.prepare(B, dec, m)
return (m, s), d, dec, B
with ThreadPoolExecutor(4) as ex:
loaded = list(ex.map(load, sessions))
RAW = {k: d for k, d, _, _ in loaded}
DEC = {k: dec for k, _, dec, _ in loaded}
BB = {k: B for k, _, _, B in loaded}
KEYS = sorted(BB)
print(f"loaded {len(KEYS)} files of dandiset {bmi.DANDISET_ID} version {bmi.DANDISET_VERSION} in {time.time() - started:.0f} s")
loaded 13 files of dandiset 000404 version 0.230605.2024 in 48 s
Data check¶
In [3]:
api = "https://api.dandiarchive.org/api/dandisets/000404/versions/"
vers = requests.get(api).json()["results"]
print("Versions:", [(v["version"], v["asset_count"], v["size"]) for v in vers])
def asset_table(v):
a = requests.get(api + v + "/assets/?page_size=100").json()["results"]
return {x["path"]: (x["size"], x.get("blob") or x.get("zarr")) for x in a}
pub = asset_table(bmi.DANDISET_VERSION)
draft = requests.get(api + "draft/assets/?page_size=100").json()["results"]
print("Published files:", len(pub), "; draft files:", len(draft),
"; same paths and sizes:", sorted((x["path"], x["size"]) for x in draft) == sorted((p, s[0]) for p, s in pub.items()))
rows = []
for k in KEYS:
d, B, dec = RAW[k], BB[k], DEC[k]
T = B["trials"]
rows.append(dict(monkey=k[0], session=k[1], identifier=d["identifier"], subject=d["subject"], rate_hz=d["rate"],
frames=len(d["spike_counts"]), minutes=round(len(d["spike_counts"]) / d["rate"] / 60, 1),
units=d["spike_counts"].shape[1], task_entries=len(np.unique(d["te_num"])),
trials=len(d["start"]), trials_analysed=int((T.cond != "dropped").sum()),
co_trials=int(((T.task == 0) & (T.cond != "dropped")).sum()),
obs_trials=int(((T.task == 1) & (T.cond != "dropped")).sum()),
bins=len(B["cmd"]), conditions=B["trials"].query("cond != 'dropped'").cond.nunique()))
info = pd.DataFrame(rows)
print(info.to_string(index=False))
print("\nUnits per session, Monkey G: min", info[info.monkey == "G"].units.min(), "max", info[info.monkey == "G"].units.max(),
"median", info[info.monkey == "G"].units.median(), "(Methods: 35-151, median 48); Monkey J:",
sorted(info[info.monkey == "J"].units.unique()), "(Methods: 20 or 21)")
print("Sessions: G", (info.monkey == "G").sum(), "J", (info.monkey == "J").sum(), "(paper: 9 and 4)")
d0 = RAW[("G", 0)]
dj = RAW[("J", 0)]
print("\nStated units, Monkey G:", d0["meta"])
print("Stated units, Monkey J:", dj["meta"])
print("Descriptions, Monkey J:", {k: v for k, v in dj["descr"].items() if k in ("cursor", "decoder_state", "target_state", "obs_details")})
for k in [("G", 0), ("J", 0)]:
d = RAW[k]
print(f"{k}: cursor x range {d['cursor'][:, 0].min():.3f} to {d['cursor'][:, 0].max():.3f}, y {d['cursor'][:, 1].min():.3f} to {d['cursor'][:, 1].max():.3f}; "
f"decoder velocity range {d['decoder_state'][:, 2:4].min():.3f} to {d['decoder_state'][:, 2:4].max():.3f}; "
f"update frames {int(d['update_bmi'].sum())} of {len(d['update_bmi'])}")
Versions: [('draft', 13, 1046740856), ('0.230605.2024', 13, 1046740856)]
Published files: 13 ; draft files: 13 ; same paths and sizes: True
monkey session identifier subject rate_hz frames minutes units task_entries trials trials_analysed co_trials obs_trials bins conditions
G 0 MonkeyG_session0 monk_g 60.0 137749 38.3 44 3 287 287 64 223 11971 17
G 1 MonkeyG_session1 monk_g 60.0 98575 27.4 35 2 193 193 64 129 7252 18
G 2 MonkeyG_session2 monk_g 60.0 237585 66.0 148 3 527 527 64 463 16505 19
G 3 MonkeyG_session3 monk_g 60.0 201709 56.0 53 4 474 474 64 410 15870 19
G 4 MonkeyG_session4 monk_g 60.0 141347 39.3 151 3 344 344 64 280 12383 20
G 5 MonkeyG_session5 monk_g 60.0 152485 42.4 42 4 398 398 127 271 11696 18
G 6 MonkeyG_session6 monk_g 60.0 68778 19.1 50 3 156 156 63 93 4536 18
G 7 MonkeyG_session7 monk_g 60.0 46309 12.9 48 2 164 164 64 100 5692 16
G 8 MonkeyG_session8 monk_g 60.0 49950 13.9 48 2 139 139 64 75 5412 17
J 0 MonkeyJ_session0 monk_j 200.0 852000 71.0 20 2 370 359 199 160 8848 16
J 1 MonkeyJ_session1 monk_j 200.0 544997 45.4 20 2 168 158 16 142 5110 16
J 2 MonkeyJ_session2 monk_j 200.0 312000 26.0 20 2 223 218 90 128 5525 16
J 3 MonkeyJ_session3 monk_j 200.0 432000 36.0 21 2 181 176 80 96 6015 17
Units per session, Monkey G: min 35 max 151 median 48.0 (Methods: 35-151, median 48); Monkey J: [np.int64(20), np.int64(21)] (Methods: 20 or 21)
Sessions: G 9 J 4 (paper: 9 and 4)
Stated units, Monkey G: {'cursor': 'cm', 'decoder_state': 'cm, cm/sec', 'obs_details': 'text', 'obstacle_position': 'cm', 'spike_counts': 'cnt', 'target_state': 'cm', 'te_num': 'int', 'update_bmi': 'binary indicator'}
Stated units, Monkey J: {'cursor': 'm', 'decoder_state': 'm, m/sec', 'obs_details': 'text', 'obstacle_position': 'cm', 'spike_counts': 'cnt', 'target_state': 'cm', 'te_num': 'int', 'update_bmi': 'binary indicator'}
Descriptions, Monkey J: {'cursor': '2D cursor position (x, y) (meters)', 'decoder_state': 'decoder state (2d-pos, 2d-vel, offset) (meters, m/sec)', 'obs_details': 'obstacle shape (position (x,y), major/minor radius, and orientation of obstacle) (cm)', 'target_state': 'target state (2d-pos, 2d-vel, offset) -- location of target (cm)'}
('G', 0): cursor x range -25.000 to 10.502, y -12.233 to 14.000; decoder velocity range -7.382 to 15.425; update frames 22962 of 137749
('J', 0): cursor x range -0.082 to 0.158, y 0.000 to 0.258; decoder velocity range -0.106 to 0.186; update frames 852000 of 852000
Monkey G's files are in cm and update the decoder every 6th frame (10 Hz at 60 Hz), as stated. Monkey J's cursor and decoder state are in m and update every 5 ms frame; Monkey J's workspace spans about 0.24 m, as the Methods give (24 cm). The session counts and unit counts match the Methods. The draft holds the same files as the published version.
In [4]:
# The decoder: least squares of the velocity after each update on the state before it and the update's counts.
dec_rows = []
for k in KEYS:
dec = DEC[k]
ff = [v["frac_fit"] for v in dec["per_te"].values()]
dec_rows.append(dict(monkey=k[0], session=k[1], updates=dec["n_updates"], frac_updates_fit=round(dec["frac_kept"], 4),
vaf_command_fitted_updates=round(dec["vaf_command_kept"], 6), vaf_command_all=round(dec["vaf_command"], 4),
vaf_velocity_all=round(dec["vaf_velocity"], 5),
min_frac_fit_per_task_entry=round(min(ff), 4)))
print(pd.DataFrame(dec_rows).to_string(index=False))
# the gain applied to the binned counts reproduces the decoder's velocity update at the update frames
for k in [("G", 0), ("J", 0)]:
B = BB[k]
print(k, "commands per bin, magnitude percentiles 5/50/95:", np.round(np.percentile(np.linalg.norm(B["cmd"], axis=1), [5, 50, 95]), 3),
"(G cm/s per update; J cm/s summed over 20 frames)")
monkey session updates frac_updates_fit vaf_command_fitted_updates vaf_command_all vaf_velocity_all min_frac_fit_per_task_entry
G 0 22961 0.9994 1.000000 0.9994 0.99977 0.9993
G 1 16427 0.9997 1.000000 1.0000 0.99998 0.9993
G 2 39915 0.9918 1.000000 0.9461 0.97412 0.9833
G 3 33905 0.9905 1.000000 0.9790 0.99725 0.9858
G 4 23689 0.9933 1.000000 0.9283 0.98881 0.9864
G 5 25420 0.9988 1.000000 0.9969 0.99939 0.9969
G 6 11494 0.9951 1.000000 0.9929 0.99874 0.9924
G 7 7717 0.9978 1.000000 0.9994 0.99992 0.9971
G 8 8333 0.9981 1.000000 0.9967 0.99867 0.9973
J 0 851999 1.0000 0.989157 0.9892 0.99985 0.2672
J 1 544996 1.0000 0.988897 0.9889 0.99985 0.4726
J 2 311999 1.0000 0.995668 0.9957 0.99995 0.6384
J 3 431999 1.0000 0.992453 0.9925 0.99992 0.6932
('G', 0) commands per bin, magnitude percentiles 5/50/95: [0.368 1.369 3.223] (G cm/s per update; J cm/s summed over 20 frames)
('J', 0) commands per bin, magnitude percentiles 5/50/95: [1.045 3.77 8.227] (G cm/s per update; J cm/s summed over 20 frames)
The files hold no decoder matrices. For Monkey G the fit is exact to rounding on 98 to 100% of updates (variance accounted for 0.99999 or more on those): those are the steady-state Kalman filter, so the count coefficients are its gain. One gain fits at least 98% of the updates of every task entry of a session (smallest fraction in the table), so each session used one decoder for both tasks, as the Methods say. The updates left out are the first of each block, where the filter had not yet converged, and a few runs late in some blocks. For Monkey J the linear map accounts for 98.9 to 99.6% of the variance of the PPF's velocity update (claim monkey-j-linear-command).
In [5]:
# Clocks. (1) Spike counts against decoder state: the gain fit at frame offsets; only offset 0 is exact.
def lag_vaf(k, lag):
d = RAW[k]
S = d["decoder_state"]
Y = d["spike_counts"].astype(float)
if k[0] == "G":
up = np.where(d["update_bmi"] == 1)[0]
up = up[(up >= 10) & (up < len(S) - 10)]
cs = np.vstack([np.zeros((1, Y.shape[1])), np.cumsum(Y, 0)])
yb = cs[up + 1 + lag] - cs[up - 5 + lag]
else:
up = np.arange(10, len(S) - 10)
yb = Y[up + lag]
X = np.hstack([S[up - 1], yb])
b, *_ = np.linalg.lstsq(X, S[up][:, 2:4], rcond=None)
r = S[up][:, 2:4] - X @ b
u = S[up][:, 2:4] - S[up - 1] @ b[:5]
return 1 - (r ** 2).sum() / ((u - u.mean(0)) ** 2).sum()
for k in [("G", 0), ("G", 4), ("J", 0)]:
print(k, "variance of the update accounted for at count offsets -2..2 frames:",
[round(lag_vaf(k, l), 4) for l in (-2, -1, 0, 1, 2)])
# (2) Trials table against the cursor: distance of the cursor from the trial's start target at go and
# from the trial's end target at reward, with the table's times and shifted by 1 s.
rows = []
for k in KEYS:
d, B = RAW[k], BB[k]
rate = d["rate"]
T = B["trials"]
T = T[T.cond != "dropped"]
sc = 100.0 if k[0] == "J" else 1.0
for sh in (-1.0, 0.0, 1.0):
go = np.clip(T.go.values + int(sh * rate), 0, len(d["cursor"]) - 1)
rw = np.clip(T.reward.values + int(sh * rate), 0, len(d["cursor"]) - 1)
end_true = d["cursor"][T.reward.values]
start_true = d["cursor"][T.go.values]
rows.append(dict(monkey=k[0], session=k[1], shift_s=sh,
go_dist_cm=np.median(np.linalg.norm(d["cursor"][go] - start_true, axis=1)) * sc if sh else 0.0,
reward_dist_cm=np.median(np.linalg.norm(d["cursor"][rw] - end_true, axis=1)) * sc if sh else 0.0,
reward_radius_cm=np.median(np.linalg.norm(end_true - (bmi.J_CENTER if k[0] == "J" else 0), axis=1)) * sc))
clk = pd.DataFrame(rows)
print(clk.pivot_table(index=["monkey", "session"], columns="shift_s", values=["go_dist_cm", "reward_dist_cm"]).round(2).to_string())
# speed of the cursor at go and reward against 1 s before and after
spd = []
for k in KEYS:
d, B = RAW[k], BB[k]
T = B["trials"][B["trials"].cond != "dropped"]
v = np.linalg.norm(d["decoder_state"][:, 2:4], axis=1) * (100.0 if k[0] == "J" else 1.0)
r = int(d["rate"])
spd.append(dict(monkey=k[0], session=k[1], speed_at_go=np.median(v[T.go.values]), speed_go_plus_0p5s=np.median(v[T.go.values + r // 2]),
speed_at_reward=np.median(v[T.reward.values]), speed_reward_minus_0p5s=np.median(v[T.reward.values - r // 2])))
print(pd.DataFrame(spd).round(2).to_string(index=False))
('G', 0) variance of the update accounted for at count offsets -2..2 frames: [np.float64(0.7716), np.float64(0.8614), np.float64(0.9995), np.float64(0.8427), np.float64(0.7265)]
('G', 4) variance of the update accounted for at count offsets -2..2 frames: [np.float64(0.6581), np.float64(0.7802), np.float64(0.9465), np.float64(0.778), np.float64(0.6404)]
('J', 0) variance of the update accounted for at count offsets -2..2 frames: [np.float64(0.0399), np.float64(0.042), np.float64(0.9892), np.float64(0.0357), np.float64(0.0384)]
go_dist_cm reward_dist_cm
shift_s -1.0 0.0 1.0 -1.0 0.0 1.0
monkey session
G 0 8.61 0.0 2.63 2.59 0.0 9.59
1 8.47 0.0 3.31 2.89 0.0 9.84
2 8.92 0.0 4.69 4.38 0.0 10.14
3 9.12 0.0 3.87 4.10 0.0 10.47
4 8.70 0.0 3.82 3.86 0.0 10.22
5 9.14 0.0 4.95 4.56 0.0 10.73
6 8.61 0.0 4.99 4.88 0.0 10.52
7 8.75 0.0 2.91 3.56 0.0 10.06
8 8.71 0.0 3.08 2.79 0.0 9.88
J 0 3.83 0.0 4.59 4.01 0.0 3.44
1 3.59 0.0 3.96 3.79 0.0 2.94
2 3.78 0.0 4.10 4.58 0.0 2.97
3 2.95 0.0 3.38 3.47 0.0 2.76
monkey session speed_at_go speed_go_plus_0p5s speed_at_reward speed_reward_minus_0p5s
G 0 1.12 2.57 1.82 2.41
G 1 1.55 3.82 2.08 2.79
G 2 2.77 5.11 3.31 3.89
G 3 1.73 2.36 2.45 2.64
G 4 1.57 2.27 2.19 2.51
G 5 2.56 3.68 3.25 3.68
G 6 2.52 3.26 2.98 3.53
G 7 1.14 1.59 1.78 2.16
G 8 1.81 3.74 1.96 2.58
J 0 3.36 5.74 2.99 4.60
J 1 3.29 5.11 3.23 4.11
J 2 3.20 5.10 4.30 5.17
J 3 2.81 4.35 3.23 3.81
The spike counts and the decoder state share one clock: the gain fit is exact only with the counts at offset 0 (the counts of frames i-5 to i for an update at frame i, as the paper's code bins them). The trials table shares the cursor's clock: the cursor is slow at go and during the target hold before reward and fast half a second into the movement, and it moves away from its go and reward positions when the times are shifted by 1 s. This is consistent with the table marking go and reward, though it does not prove the exact frames.
In [6]:
# Target and obstacle series. Monkey G: is the target series at reward where the cursor ends?
for k in [("G", 0), ("G", 1)]:
d, B = RAW[k], BB[k]
T = B["trials"]
tsr = d["target_state"][T.reward.values][:, :2]
print(k, "target series at reward, unique values:", np.unique(tsr.round(2), axis=0).shape[0],
"; median distance from the cursor at reward", round(float(np.median(np.linalg.norm(tsr - d["cursor"][T.reward.values], axis=1))), 2), "cm")
obs = d["obstacle_position"][T.go.values + 2]
oo = T.task.values == 1
ang_obs = np.degrees(np.arctan2(obs[oo, 1], obs[oo, 0])) % 360
print(" obstacle trials: obstacle direction minus target direction from the cursor at reward (deg), unique:",
np.unique(np.round((ang_obs - 45 * T.target.values[oo]) % 360)), "; obstacle distance from the center (cm):",
np.unique(np.linalg.norm(obs[oo], axis=1).round(2)))
for k in [k for k in KEYS if k[0] == "J"]:
d, B = RAW[k], BB[k]
T = B["trials"]
tsr = d["target_state"][T.reward.values][:, :2]
print(k, "target series at reward, unique values:", np.unique(tsr.round(4), axis=0).tolist(),
"; obstacle series at go, unique:", np.unique(d["obs_details"][T.go.values].round(3), axis=0).shape[0])
print("\nMonkey G end positions: distance of the cursor at reward from the center (cm), 5/50/95%:",
np.percentile(np.concatenate([np.linalg.norm(RAW[k]["cursor"][BB[k]["trials"].reward.values], axis=1) for k in KEYS if k[0] == "G"]), [5, 50, 95]).round(2))
print("Monkey J center-out end positions: distance from the center (cm), 5/50/95%:",
np.percentile(np.concatenate([100 * np.linalg.norm(RAW[k]["cursor"][BB[k]["trials"].query("task == 0").reward.values] - bmi.J_CENTER, axis=1) for k in KEYS if k[0] == "J"]), [5, 50, 95]).round(2))
tj = pd.concat([BB[k]["trials"].assign(session=k[1]) for k in KEYS if k[0] == "J"])
print("\nMonkey J obstacle trials by the code's trial type (from start and end targets; -1 = pair not in the code's list):")
print(tj[tj.task == 1].groupby(["session", "ttype"]).size().unstack(fill_value=0).to_string())
('G', 0) target series at reward, unique values: 1 ; median distance from the cursor at reward 9.38 cm
obstacle trials: obstacle direction minus target direction from the cursor at reward (deg), unique: [ 0. 360.] ; obstacle distance from the center (cm): [5.]
('G', 1) target series at reward, unique values: 1 ; median distance from the cursor at reward 9.4 cm
obstacle trials: obstacle direction minus target direction from the cursor at reward (deg), unique: [ 0. 360.] ; obstacle distance from the center (cm): [5.]
('J', 0) target series at reward, unique values: [[-3.25, 8.1292], [-0.0273, 0.1384], [0.1027, 0.1384]] ; obstacle series at go, unique: 201
('J', 1) target series at reward, unique values: [[0.0, 0.0], [4.25, 2.5]] ; obstacle series at go, unique: 17
('J', 2) target series at reward, unique values: [[-4.25, 2.5], [-0.0273, 0.1384]] ; obstacle series at go, unique: 91
('J', 3) target series at reward, unique values: [[-3.25, -3.1292], [0.1027, 0.1384]] ; obstacle series at go, unique: 81
Monkey G end positions: distance of the cursor at reward from the center (cm), 5/50/95%: [ 8.95 9.83 10.89]
Monkey J center-out end positions: distance from the center (cm), 5/50/95%: [5.58 6.3 7.05]
Monkey J obstacle trials by the code's trial type (from start and end targets; -1 = pair not in the code's list):
ttype 1 2 3 4 5 6 8 9
session
0 11 23 24 21 24 23 23 22
1 10 20 20 20 21 20 21 20
2 5 13 12 13 13 24 25 28
3 5 9 10 8 10 20 20 19
The target series do not give the trial's target. In Monkey G's files the series at reward is not where the cursor ends (the cursor ends 10 cm from the center, the series sits elsewhere), and in Monkey J's files every trial of a session carries one and the same target and obstacle: the code that wrote the files indexes frames with the position of the event in the event list (nwb_utils.py, targets[(go_ - start_index_overall):] with go_ an index into strobed), so the last trial's values fill the file. Targets are therefore read from the cursor at reward, which must lie inside the target: Monkey G's cursors end 8 to 11 cm from the center (targets at 10 cm, radius 1.7 cm) and Monkey J's center-out cursors near 6.5 cm. Monkey G's obstacles lie halfway to the target so read, as the Methods describe. Monkey J's obstacle trial types come from the start and end targets matched to the code's target list; types 0 and 1 are left out as in the code.
In [7]:
# Recordings deposited twice? Compare per-session summaries of the spike counts.
sig = {k: (RAW[k]["spike_counts"].shape, int(RAW[k]["spike_counts"].sum()), RAW[k]["spike_counts"][:2000].sum(0).tolist()[:5]) for k in KEYS}
dups = [(a, b) for i, a in enumerate(KEYS) for b in KEYS[i + 1:] if sig[a][:2] == sig[b][:2]]
print("Pairs of files with identical spike-count shape and total:", dups)
print("Task entries per session:", {f"{k[0]}{k[1]}": list(np.unique(RAW[k]["te_num"])) for k in KEYS})
te_all = [t for k in KEYS for t in np.unique(RAW[k]["te_num"])]
print("Task entries used by more than one session:", sorted({t for t in te_all if te_all.count(t) > 1}))
print("Trial duration (s), median per session:", {f"{k[0]}{k[1]}": round(float(np.median(RAW[k]['stop'] - RAW[k]['start'])), 2) for k in KEYS})
Pairs of files with identical spike-count shape and total: []
Task entries per session: {'G0': [np.str_('4377'), np.str_('4378'), np.str_('4382')], 'G1': [np.str_('4394'), np.str_('4395')], 'G2': [np.str_('4411'), np.str_('4412'), np.str_('4415')], 'G3': [np.str_('4497'), np.str_('4498'), np.str_('4499'), np.str_('4504')], 'G4': [np.str_('4509'), np.str_('4510'), np.str_('4514')], 'G5': [np.str_('4520'), np.str_('4522'), np.str_('4523'), np.str_('4525')], 'G6': [np.str_('4532'), np.str_('4533'), np.str_('4536')], 'G7': [np.str_('4549'), np.str_('4553')], 'G8': [np.str_('4558'), np.str_('4560')], 'J0': [np.str_('jeev082413d'), np.str_('jeev082413e')], 'J1': [np.str_('jeev082413g'), np.str_('jeev082413h')], 'J2': [np.str_('jeev082913c'), np.str_('jeev082913d')], 'J3': [np.str_('jeev083113e'), np.str_('jeev083113f')]}
Task entries used by more than one session: []
Trial duration (s), median per session: {'G0': 3.9, 'G1': 3.65, 'G2': 2.95, 'G3': 3.13, 'G4': 3.35, 'G5': 2.95, 'G6': 2.9, 'G7': 3.35, 'G8': 3.85, 'J0': 2.06, 'J1': 2.97, 'J2': 2.4, 'J3': 2.65}
No file repeats another, and no task entry appears in two sessions. Units are not followed across sessions; every analysis is per session.
Per-session computation¶
Every per-session quantity (tuples, distances and their nulls, dynamics models, left-out models, 1000 shuffle models, behavior-encoding models and pairs of conditions) is computed in one pass per session, in parallel. Parameters and their sources are in the docstrings of analyses.py and bmi.py.
In [8]:
def job(k):
t0 = time.time()
R = an.run_session(BB[k], k[0], NSHUFF, 1000 * (k[0] == "J") + k[1])
return k, R, time.time() - t0
order = sorted(KEYS, key=lambda k: -BB[k]["spk"].shape[1] ** 2 * len(BB[k]["spk"]))
RR = {}
t0 = time.time()
with ProcessPoolExecutor(4) as ex:
for k, R, dt in ex.map(job, order):
RR[k] = R
print(k, f"{dt:.0f} s", flush=True)
print(f"all sessions in {time.time() - t0:.0f} s")
G = [k for k in KEYS if k[0] == "G"]
J = [k for k in KEYS if k[0] == "J"]
def pval_mean(obs_by_session, null_by_session, greater=True):
"""Pooled test (Methods, Statistics Summary): mean over sessions against the per-shuffle mean over sessions."""
o = np.mean(obs_by_session)
nl = np.mean(np.vstack(null_by_session), axis=0)
return o, (nl >= o).mean() if greater else (nl <= o).mean(), nl
('G', 2) 546 s
('G', 4) 456 s
('G', 3) 204 s
('G', 0) 120 s
('G', 5) 116 s
('G', 7) 54 s
('G', 8) 52 s
('G', 6) 43 s
('G', 1) 57 s
('J', 0) 51 s
('J', 3) 32 s
('J', 2) 30 s
('J', 1) 16 s
all sessions in 546 s
conditions-per-session¶
The two tasks elicited up to 24 movement conditions, with an average of 16 to 17 conditions per session.
In [9]:
rows = []
for k in KEYS:
T = BB[k]["trials"]
T = T[T.cond != "dropped"]
n_b = len({x["cond"] for x in RR[k]["tups"]})
allt = bmi.tuples(BB[k], BB[k]["mags"])
rows.append(dict(monkey=k[0], session=k[1], conds_a_any_trial=T.cond.nunique(),
conds_b_with_tuple=len({x["cond"] for x in allt}), co=T[T.task == 0].cond.nunique(),
cw=T[T.cond.str.startswith("cw")].cond.nunique(), ccw=T[T.cond.str.startswith("ccw")].cond.nunique()))
cps = pd.DataFrame(rows)
print(cps.to_string(index=False))
for col in ["conds_a_any_trial", "conds_b_with_tuple"]:
print(f"{col}: mean over 13 sessions {cps[col].mean():.2f} (G {cps[cps.monkey == 'G'][col].mean():.2f}, J {cps[cps.monkey == 'J'][col].mean():.2f}), max {cps[col].max()}; "
f"expected 14.4 to 18.7")
print("Monkey J obstacle trials of types 0 and 1 (left out):", int((tj.ttype.isin([0, 1]) & (tj.task == 1)).sum()), "; unmatched pairs:", int((tj.ttype == -1).sum()))
monkey session conds_a_any_trial conds_b_with_tuple co cw ccw
G 0 17 17 8 2 7
G 1 18 17 8 3 7
G 2 19 18 8 3 8
G 3 19 19 8 3 8
G 4 20 20 8 4 8
G 5 18 17 8 2 8
G 6 18 14 8 2 8
G 7 16 16 8 0 8
G 8 17 16 8 1 8
J 0 16 16 8 4 4
J 1 16 9 8 4 4
J 2 16 16 8 4 4
J 3 17 16 8 5 4
conds_a_any_trial: mean over 13 sessions 17.46 (G 18.00, J 16.25), max 20; expected 14.4 to 18.7
conds_b_with_tuple: mean over 13 sessions 16.23 (G 17.11, J 14.25), max 20; expected 14.4 to 18.7
Monkey J obstacle trials of types 0 and 1 (left out): 31 ; unmatched pairs: 0
Result: reproduced. Under reading (a), conditions with at least one rewarded trial, the mean over the 13 sessions is 17.46 conditions (Monkey G 18.00, Monkey J 16.25; at most 20), inside 14.4 to 18.7. Under reading (b), conditions with a (command, condition) tuple of 15 or more occurrences, it is 16.23, also inside.
figs1b-conditions-per-command¶
Each analysed command was used regularly (15 or more times) in several conditions, about 7 on average, typically 5 to 10.
In [10]:
rows = []
for k in KEYS:
B = BB[k]
mb, ab = bmi.command_bins(B["cmd"], B["mags"])
for m in range(4):
for a in range(8):
ix = (mb == m) & (ab == a)
c, n = np.unique(B["cond"][ix], return_counts=True)
rows.append(dict(monkey=k[0], session=k[1], mag=m, ang=a, n_conds=int((n >= 15).sum())))
cc = pd.DataFrame(rows)
used = cc[cc.n_conds >= 2]
print("Commands with at least two conditions of 15+ occurrences, per session:")
print(used.groupby(["monkey", "session"]).n_conds.agg(["count", "mean", "median", "min", "max"]).round(2).to_string())
print(f"\n(a) mean over all commands pooled: {used.n_conds.mean():.2f}; (b) mean of session means: {used.groupby(['monkey','session']).n_conds.mean().mean():.2f}; "
f"median over commands {used.n_conds.median():.1f}; 25th-75th percentile {used.n_conds.quantile(0.25):.1f}-{used.n_conds.quantile(0.75):.1f}")
print("By monkey:", used.groupby("monkey").n_conds.agg(["mean", "median"]).round(2).to_dict())
used1 = cc[cc.n_conds >= 1]
print(f"Variant, commands with at least one such condition: mean {used1.n_conds.mean():.2f}, median {used1.n_conds.median():.1f}")
Commands with at least two conditions of 15+ occurrences, per session:
count mean median min max
monkey session
G 0 32 7.44 8.0 3 11
1 24 6.17 5.5 2 12
2 30 7.20 8.0 2 11
3 32 9.03 9.0 3 15
4 32 6.69 6.0 2 14
5 32 8.84 9.0 2 11
6 15 4.73 5.0 2 9
7 26 5.00 5.0 2 9
8 22 5.05 5.0 2 9
J 0 16 9.62 9.5 2 16
1 15 5.53 6.0 3 8
2 12 9.08 8.5 5 15
3 20 6.35 6.5 2 11
(a) mean over all commands pooled: 7.06; (b) mean of session means: 6.98; median over commands 7.0; 25th-75th percentile 5.0-9.0
By monkey: {'mean': {'G': 6.94, 'J': 7.51}, 'median': {'G': 7.0, 'J': 7.0}}
Variant, commands with at least one such condition: mean 6.69, median 7.0
Result: reproduced. Over commands used 15 or more times in at least two conditions, a command is used regularly in 7.06 conditions on average over all commands (reading a) and 6.98 over session means (reading b), both inside 6.3 to 7.7; the median is 7 (25th to 75th percentile 5 to 9), inside 4.5 to 11.
figs1d-command-subtrajectories-differ and figs1e-next-commands-differ¶
The same command is used within command subtrajectories that differ between conditions more than expected from subsampling the condition-pooled subtrajectories (normalized distance 1.920 for Monkey G and 2.342 for Monkey J).
The same command is followed by next commands that differ between conditions more than expected from subsampling (normalized distance 1.676 for Monkey G and 2.178 for Monkey J).
Result: figs1d-command-subtrajectories-differ reproduced; figs1e-next-commands-differ partially reproduced. Command subtrajectories: mean normalized distance 1.953 (Monkey G) and 2.330 (Monkey J), both within 10% of 1.920 and 2.342, with p < 0.001 in 9 of 9 and 4 of 4 sessions and pooled (no bootstrap average reached the data). Next commands: p < 0.001 in every session and pooled, and Monkey J's 2.218 is within 10% of 2.178, but Monkey G's 1.873 is above the upper bound of 1.844 (1.676 + 10%). Counting a tuple only when 15 or more of its occurrences remain after the occurrences near a trial's edges are removed gives 2.026 (G) and 2.481 (J) for subtrajectories, within tolerance, and 1.870 and 2.236 for next commands, Monkey G still above 1.844; every session and pooled test stays at p < 0.001.
fig3d-population-distance, fig3b-single-neuron-distance, fig3e-fractions-significant¶
Neural population activity issuing a given command differs between conditions more than expected from shuffled data, by 10% to 200% on session averages (normalized distance 1.222 for Monkey G and 1.724 for Monkey J), in every session.
Single-neuron activity issuing a given command differs between conditions more than expected from shuffled data, aggregating over (command, condition, neuron) tuples, in every session.
Distances exceed the shuffle for about 30% (Monkey G) and 70% (Monkey J) of individual (command, condition) tuples, about 65% and 90% of commands aggregated over conditions, and about 40% and 80% of neurons aggregated over tuples.
In [12]:
print("Tuples per session: 15+ occurrences, dropped by the matching, analysed:")
print(pd.DataFrame([dict(monkey=k[0], session=k[1], tuples_15=RR[k]["n_tuples_min15"], dropped_match=RR[k]["n_tuples_dropped_match"],
analysed=len(RR[k]["tups"]), with_2plus_conds=sum(x["n_conds"] >= 2 for x in RR[k]["tups"])) for k in KEYS]).to_string(index=False))
f3_rows = []
for null in ["shuffle", "subsample"]:
for mc in [2, 1]:
for mk, ks in [("G", G), ("J", J)]:
sess_obs, sess_null = [], []
for k in ks:
f = RR[k]["f3_" + null]
m = np.array([x["n_conds"] >= mc for x in RR[k]["tups"]])
mn = f["null_pop"][:, m].mean(0)
obs = (f["pop"][m] / mn).mean()
nl = (f["null_pop"][:, m] / mn).mean(1)
sess_obs.append(obs)
sess_null.append(nl)
f3_rows.append(dict(null=null, min_conds=mc, monkey=mk, session=k[1], tuples=int(m.sum()), mean_norm=obs,
null_95=np.percentile(nl, 95), p=(nl >= obs).mean()))
o, p, _ = pval_mean(sess_obs, sess_null)
f3_rows.append(dict(null=null, min_conds=mc, monkey=mk, session="pooled", tuples=np.nan, mean_norm=o, null_95=np.nan, p=p))
F3 = pd.DataFrame(f3_rows)
print(F3.round(4).to_string(index=False))
summ = F3[F3.session != "pooled"].groupby(["null", "min_conds", "monkey"]).agg(sess_min=("mean_norm", "min"), sess_max=("mean_norm", "max"),
n_p001=("p", lambda x: (x < 0.001).sum()), n=("p", "size"))
summ["mean"] = F3[F3.session == "pooled"].set_index(["null", "min_conds", "monkey"]).mean_norm
summ["pooled_p"] = F3[F3.session == "pooled"].set_index(["null", "min_conds", "monkey"]).p
print("\nBy reading (expected means: G 1.222 +- 0.122, J 1.724 +- 0.172; every session 0.99 to 3.3):")
print(summ.round(4).to_string())
Tuples per session: 15+ occurrences, dropped by the matching, analysed:
monkey session tuples_15 dropped_match analysed with_2plus_conds
G 0 238 0 238 238
G 1 148 0 148 148
G 2 216 0 216 216
G 3 289 0 289 289
G 4 214 0 214 214
G 5 283 0 283 283
G 6 78 0 78 71
G 7 135 0 135 130
G 8 116 0 116 111
J 0 154 0 154 154
J 1 84 0 84 83
J 2 110 0 110 109
J 3 128 0 128 127
null min_conds monkey session tuples mean_norm null_95 p
shuffle 2 G 0 238.0 1.3646 1.0179 0.0
shuffle 2 G 1 148.0 1.4036 1.0285 0.0
shuffle 2 G 2 216.0 1.3375 1.0147 0.0
shuffle 2 G 3 289.0 1.2178 1.0134 0.0
shuffle 2 G 4 214.0 1.2700 1.0158 0.0
shuffle 2 G 5 283.0 1.1003 1.0142 0.0
shuffle 2 G 6 71.0 1.1520 1.0318 0.0
shuffle 2 G 7 130.0 1.1323 1.0204 0.0
shuffle 2 G 8 111.0 1.1589 1.0230 0.0
shuffle 2 G pooled NaN 1.2374 NaN 0.0
shuffle 2 J 0 154.0 1.8369 1.0302 0.0
shuffle 2 J 1 83.0 1.9757 1.0423 0.0
shuffle 2 J 2 109.0 1.8366 1.0349 0.0
shuffle 2 J 3 127.0 1.4238 1.0330 0.0
shuffle 2 J pooled NaN 1.7683 NaN 0.0
shuffle 1 G 0 238.0 1.3646 1.0179 0.0
shuffle 1 G 1 148.0 1.4036 1.0285 0.0
shuffle 1 G 2 216.0 1.3375 1.0147 0.0
shuffle 1 G 3 289.0 1.2178 1.0134 0.0
shuffle 1 G 4 214.0 1.2700 1.0158 0.0
shuffle 1 G 5 283.0 1.1003 1.0142 0.0
shuffle 1 G 6 78.0 1.1510 1.0289 0.0
shuffle 1 G 7 135.0 1.1282 1.0201 0.0
shuffle 1 G 8 116.0 1.1563 1.0225 0.0
shuffle 1 G pooled NaN 1.2366 NaN 0.0
shuffle 1 J 0 154.0 1.8369 1.0302 0.0
shuffle 1 J 1 84.0 1.9747 1.0425 0.0
shuffle 1 J 2 110.0 1.8282 1.0354 0.0
shuffle 1 J 3 128.0 1.4204 1.0330 0.0
shuffle 1 J pooled NaN 1.7651 NaN 0.0
subsample 2 G 0 238.0 1.3662 1.0172 0.0
subsample 2 G 1 148.0 1.4066 1.0293 0.0
subsample 2 G 2 216.0 1.3377 1.0144 0.0
subsample 2 G 3 289.0 1.2176 1.0128 0.0
subsample 2 G 4 214.0 1.2706 1.0147 0.0
subsample 2 G 5 283.0 1.1005 1.0141 0.0
subsample 2 G 6 71.0 1.1532 1.0301 0.0
subsample 2 G 7 130.0 1.1316 1.0200 0.0
subsample 2 G 8 111.0 1.1580 1.0214 0.0
subsample 2 G pooled NaN 1.2380 NaN 0.0
subsample 2 J 0 154.0 1.8393 1.0301 0.0
subsample 2 J 1 83.0 1.9796 1.0395 0.0
subsample 2 J 2 109.0 1.8424 1.0341 0.0
subsample 2 J 3 127.0 1.4240 1.0320 0.0
subsample 2 J pooled NaN 1.7713 NaN 0.0
subsample 1 G 0 238.0 1.3662 1.0172 0.0
subsample 1 G 1 148.0 1.4066 1.0293 0.0
subsample 1 G 2 216.0 1.3377 1.0144 0.0
subsample 1 G 3 289.0 1.2176 1.0128 0.0
subsample 1 G 4 214.0 1.2706 1.0147 0.0
subsample 1 G 5 283.0 1.1005 1.0141 0.0
subsample 1 G 6 78.0 1.1520 1.0300 0.0
subsample 1 G 7 135.0 1.1276 1.0199 0.0
subsample 1 G 8 116.0 1.1557 1.0213 0.0
subsample 1 G pooled NaN 1.2372 NaN 0.0
subsample 1 J 0 154.0 1.8393 1.0301 0.0
subsample 1 J 1 84.0 1.9782 1.0382 0.0
subsample 1 J 2 110.0 1.8339 1.0350 0.0
subsample 1 J 3 128.0 1.4205 1.0320 0.0
subsample 1 J pooled NaN 1.7680 NaN 0.0
By reading (expected means: G 1.222 +- 0.122, J 1.724 +- 0.172; every session 0.99 to 3.3):
sess_min sess_max n_p001 n mean pooled_p
null min_conds monkey
shuffle 1 G 1.1003 1.4036 9 9 1.2366 0.0
J 1.4204 1.9747 4 4 1.7651 0.0
2 G 1.1003 1.4036 9 9 1.2374 0.0
J 1.4238 1.9757 4 4 1.7683 0.0
subsample 1 G 1.1005 1.4066 9 9 1.2372 0.0
J 1.4205 1.9782 4 4 1.7680 0.0
2 G 1.1005 1.4066 9 9 1.2380 0.0
J 1.4240 1.9796 4 4 1.7713 0.0
Result: reproduced. Under all four combinations of readings (null by shuffle or subsampling; commands with one or with two qualifying conditions) the mean normalized population distance is 1.237 to 1.238 for Monkey G and 1.765 to 1.771 for Monkey J, within 10% of 1.222 and 1.724; session means range from 1.10 to 1.98, inside 0.99 to 3.3; and every session and both pooled tests give p < 0.001 (no shuffle reached the data). No tuple was lost to the matching.
In [13]:
rows = []
for null in ["shuffle", "subsample"]:
for mc in [2, 1]:
for mk, ks in [("G", G), ("J", J)]:
so, sn = [], []
for k in ks:
f = RR[k]["f3_" + null]
m = np.array([x["n_conds"] >= mc for x in RR[k]["tups"]])
obs = f["neur"][m].mean()
nl = f["null_neur_tuple"][:, m].mean(1)
so.append(obs)
sn.append(nl)
rows.append(dict(null=null, min_conds=mc, monkey=mk, session=k[1], mean_hz=obs, null_mean=nl.mean(), p=(nl >= obs).mean()))
o, p, nl = pval_mean(so, sn)
rows.append(dict(null=null, min_conds=mc, monkey=mk, session="pooled", mean_hz=o, null_mean=nl.mean(), p=p))
F3B = pd.DataFrame(rows)
print(F3B.round(4).to_string(index=False))
print("\nExpected: G 1.167 +- 0.117 Hz, J 1.235 +- 0.124 Hz; shuffle means printed by the code 1.004 (G), 0.745 (J)")
null min_conds monkey session mean_hz null_mean p shuffle 2 G 0 1.5786 1.2287 0.0 shuffle 2 G 1 0.9662 0.7733 0.0 shuffle 2 G 2 0.7661 0.6311 0.0 shuffle 2 G 3 1.1501 0.9764 0.0 shuffle 2 G 4 0.8750 0.7543 0.0 shuffle 2 G 5 1.3214 1.2185 0.0 shuffle 2 G 6 1.1237 1.0082 0.0 shuffle 2 G 7 1.3929 1.2670 0.0 shuffle 2 G 8 1.3194 1.1621 0.0 shuffle 2 G pooled 1.1659 1.0022 0.0 shuffle 2 J 0 1.3009 0.7848 0.0 shuffle 2 J 1 1.4388 0.8217 0.0 shuffle 2 J 2 1.3355 0.8041 0.0 shuffle 2 J 3 0.8143 0.5957 0.0 shuffle 2 J pooled 1.2224 0.7516 0.0 shuffle 1 G 0 1.5786 1.2287 0.0 shuffle 1 G 1 0.9662 0.7733 0.0 shuffle 1 G 2 0.7661 0.6311 0.0 shuffle 1 G 3 1.1501 0.9764 0.0 shuffle 1 G 4 0.8750 0.7543 0.0 shuffle 1 G 5 1.3214 1.2185 0.0 shuffle 1 G 6 1.1266 1.0095 0.0 shuffle 1 G 7 1.3968 1.2754 0.0 shuffle 1 G 8 1.3252 1.1723 0.0 shuffle 1 G pooled 1.1673 1.0044 0.0 shuffle 1 J 0 1.3009 0.7848 0.0 shuffle 1 J 1 1.4415 0.8245 0.0 shuffle 1 J 2 1.3296 0.8033 0.0 shuffle 1 J 3 0.8123 0.5958 0.0 shuffle 1 J pooled 1.2211 0.7521 0.0 subsample 2 G 0 1.5786 1.2276 0.0 subsample 2 G 1 0.9662 0.7715 0.0 subsample 2 G 2 0.7661 0.6311 0.0 subsample 2 G 3 1.1501 0.9763 0.0 subsample 2 G 4 0.8750 0.7540 0.0 subsample 2 G 5 1.3214 1.2184 0.0 subsample 2 G 6 1.1237 1.0077 0.0 subsample 2 G 7 1.3929 1.2675 0.0 subsample 2 G 8 1.3194 1.1621 0.0 subsample 2 G pooled 1.1659 1.0018 0.0 subsample 2 J 0 1.3009 0.7840 0.0 subsample 2 J 1 1.4388 0.8213 0.0 subsample 2 J 2 1.3355 0.8016 0.0 subsample 2 J 3 0.8143 0.5956 0.0 subsample 2 J pooled 1.2224 0.7506 0.0 subsample 1 G 0 1.5786 1.2276 0.0 subsample 1 G 1 0.9662 0.7715 0.0 subsample 1 G 2 0.7661 0.6311 0.0 subsample 1 G 3 1.1501 0.9763 0.0 subsample 1 G 4 0.8750 0.7540 0.0 subsample 1 G 5 1.3214 1.2184 0.0 subsample 1 G 6 1.1266 1.0092 0.0 subsample 1 G 7 1.3968 1.2759 0.0 subsample 1 G 8 1.3252 1.1721 0.0 subsample 1 G pooled 1.1673 1.0040 0.0 subsample 1 J 0 1.3009 0.7840 0.0 subsample 1 J 1 1.4415 0.8243 0.0 subsample 1 J 2 1.3296 0.8007 0.0 subsample 1 J 3 0.8123 0.5957 0.0 subsample 1 J pooled 1.2211 0.7512 0.0 Expected: G 1.167 +- 0.117 Hz, J 1.235 +- 0.124 Hz; shuffle means printed by the code 1.004 (G), 0.745 (J)
Result: reproduced. The mean single-neuron distance is 1.166 to 1.167 Hz for Monkey G and 1.221 to 1.222 Hz for Monkey J under the readings, within 10% of 1.167 and 1.235 Hz, and above the shuffle in every session and pooled with p < 0.001. The shuffle means, 1.002 to 1.004 Hz (G) and 0.751 to 0.752 Hz (J), match the 1.004 and 0.745 the paper's code prints.
In [14]:
rows = []
for null in ["shuffle", "subsample"]:
for mc in [2, 1]:
for k in KEYS:
f = RR[k]["f3_" + null]
tu = RR[k]["tups"]
m = np.array([x["n_conds"] >= mc for x in tu])
pop, nullp = f["pop"][m], f["null_pop"][:, m]
frac_t = (pop > np.percentile(nullp, 95, axis=0)).mean()
cmds = np.array([(x["mag"], x["ang"]) for x in tu])[m]
keys_c = sorted({tuple(c) for c in cmds})
sig_c = []
for c in keys_c:
cm = (cmds[:, 0] == c[0]) & (cmds[:, 1] == c[1])
o = pop[cm].mean()
nl = nullp[:, cm].mean(1)
sig_c.append((nl >= o).mean() < 0.05)
res = dict(null=null, min_conds=mc, monkey=k[0], session=k[1], tuples=frac_t, commands=np.mean(sig_c),
n_tuples=int(m.sum()), n_commands=len(keys_c), n_neurons=f["neur"].shape[1])
for nm, mkey in [("neurons_all_tuples", "all" if mc == 1 else "min2"), ("neurons_sig_tuples", "sig_all" if mc == 1 else "sig_min2")]:
if nm == "neurons_all_tuples":
mm = m
else:
mm = m & f["sig"] if (m & f["sig"]).any() else m
o = f["neur"][mm].mean(0)
nl = f["null_neur_mean"][mkey]
res[nm] = ((nl >= o[None, :]).mean(0) < 0.05).mean()
rows.append(res)
F3E = pd.DataFrame(rows)
print(F3E.round(3).to_string(index=False))
tab = F3E.groupby(["null", "min_conds", "monkey"])[["tuples", "commands", "neurons_all_tuples", "neurons_sig_tuples"]].mean()
print("\nMean of session fractions (expected tuples G 0.30, J 0.70; commands 0.65, 0.90; neurons 0.40, 0.80; 10% of each value):")
print(tab.round(3).to_string())
expv = {"G": dict(tuples=0.30, commands=0.65, neurons=0.40), "J": dict(tuples=0.70, commands=0.90, neurons=0.80)}
def within(v, e):
return abs(v - e) <= 0.1 * e + 1e-12
best = []
for (null, mc), sub in tab.groupby(level=[0, 1]):
for nkey in ["neurons_all_tuples", "neurons_sig_tuples"]:
n_ok = 0
for mk in ["G", "J"]:
r = sub.xs(mk, level=2).iloc[0]
n_ok += int(within(r.tuples, expv[mk]["tuples"])) + int(within(r.commands, expv[mk]["commands"])) + int(within(r[nkey], expv[mk]["neurons"]))
best.append(dict(null=null, min_conds=mc, neurons=nkey, values_within_tolerance=n_ok))
print(pd.DataFrame(best).to_string(index=False))
# Second reading: the fraction pooled over all items of a monkey's sessions (session fractions
# weighted by the session's number of tuples, commands or neurons).
pooled_rows = []
for (null, mc, mk), sub in F3E.groupby(["null", "min_conds", "monkey"]):
pooled_rows.append(dict(null=null, min_conds=mc, monkey=mk,
tuples=np.average(sub.tuples, weights=sub.n_tuples),
commands=np.average(sub.commands, weights=sub.n_commands),
neurons_all_tuples=np.average(sub.neurons_all_tuples, weights=sub.n_neurons),
neurons_sig_tuples=np.average(sub.neurons_sig_tuples, weights=sub.n_neurons)))
pooled_frac = pd.DataFrame(pooled_rows).set_index(["null", "min_conds", "monkey"])
print("\nSecond reading, fraction pooled over all items of the monkey's sessions:")
print(pooled_frac.round(3).to_string())
best = []
for (null, mc), sub in pooled_frac.groupby(level=[0, 1]):
for nkey in ["neurons_all_tuples", "neurons_sig_tuples"]:
n_ok = 0
for mk in ["G", "J"]:
r = sub.xs(mk, level=2).iloc[0]
n_ok += int(within(r.tuples, expv[mk]["tuples"])) + int(within(r.commands, expv[mk]["commands"])) + int(within(r[nkey], expv[mk]["neurons"]))
best.append(dict(null=null, min_conds=mc, neurons=nkey, values_within_tolerance=n_ok))
print(pd.DataFrame(best).to_string(index=False))
null min_conds monkey session tuples commands n_tuples n_commands n_neurons neurons_all_tuples neurons_sig_tuples
shuffle 2 G 0 0.576 0.938 238 32 44 0.636 0.705
shuffle 2 G 1 0.446 0.708 148 24 35 0.400 0.400
shuffle 2 G 2 0.593 0.833 216 30 148 0.392 0.432
shuffle 2 G 3 0.394 0.969 289 32 53 0.623 0.585
shuffle 2 G 4 0.397 0.656 214 32 151 0.272 0.325
shuffle 2 G 5 0.201 0.531 283 32 42 0.357 0.571
shuffle 2 G 6 0.239 0.600 71 15 50 0.220 0.240
shuffle 2 G 7 0.208 0.423 130 26 48 0.312 0.354
shuffle 2 G 8 0.306 0.636 111 22 48 0.417 0.479
shuffle 2 J 0 0.734 1.000 154 16 20 0.950 0.950
shuffle 2 J 1 0.819 1.000 83 15 20 0.900 0.900
shuffle 2 J 2 0.780 1.000 109 12 20 0.950 0.900
shuffle 2 J 3 0.504 0.800 127 20 21 0.762 0.810
shuffle 1 G 0 0.576 0.938 238 32 44 0.636 0.705
shuffle 1 G 1 0.446 0.708 148 24 35 0.400 0.400
shuffle 1 G 2 0.593 0.833 216 30 148 0.392 0.432
shuffle 1 G 3 0.394 0.969 289 32 53 0.623 0.585
shuffle 1 G 4 0.397 0.656 214 32 151 0.272 0.325
shuffle 1 G 5 0.201 0.531 283 32 42 0.357 0.571
shuffle 1 G 6 0.244 0.500 78 22 50 0.200 0.260
shuffle 1 G 7 0.200 0.355 135 31 48 0.333 0.354
shuffle 1 G 8 0.302 0.556 116 27 48 0.396 0.521
shuffle 1 J 0 0.734 1.000 154 16 20 0.950 0.950
shuffle 1 J 1 0.821 1.000 84 16 20 0.900 0.900
shuffle 1 J 2 0.773 0.923 110 13 20 1.000 0.900
shuffle 1 J 3 0.500 0.762 128 21 21 0.762 0.810
subsample 2 G 0 0.576 0.938 238 32 44 0.636 0.705
subsample 2 G 1 0.446 0.708 148 24 35 0.429 0.457
subsample 2 G 2 0.593 0.833 216 30 148 0.385 0.453
subsample 2 G 3 0.391 0.969 289 32 53 0.604 0.585
subsample 2 G 4 0.397 0.625 214 32 151 0.252 0.351
subsample 2 G 5 0.201 0.531 283 32 42 0.405 0.619
subsample 2 G 6 0.254 0.600 71 15 50 0.220 0.260
subsample 2 G 7 0.215 0.462 130 26 48 0.333 0.333
subsample 2 G 8 0.306 0.636 111 22 48 0.417 0.479
subsample 2 J 0 0.734 1.000 154 16 20 0.950 0.950
subsample 2 J 1 0.819 1.000 83 15 20 0.900 0.900
subsample 2 J 2 0.789 1.000 109 12 20 0.950 0.900
subsample 2 J 3 0.496 0.800 127 20 21 0.762 0.762
subsample 1 G 0 0.576 0.938 238 32 44 0.636 0.705
subsample 1 G 1 0.446 0.708 148 24 35 0.429 0.457
subsample 1 G 2 0.593 0.833 216 30 148 0.385 0.453
subsample 1 G 3 0.391 0.969 289 32 53 0.604 0.585
subsample 1 G 4 0.397 0.625 214 32 151 0.252 0.351
subsample 1 G 5 0.201 0.531 283 32 42 0.405 0.619
subsample 1 G 6 0.256 0.500 78 22 50 0.200 0.240
subsample 1 G 7 0.207 0.387 135 31 48 0.333 0.333
subsample 1 G 8 0.302 0.556 116 27 48 0.396 0.500
subsample 1 J 0 0.734 1.000 154 16 20 0.950 0.950
subsample 1 J 1 0.821 1.000 84 16 20 0.900 0.900
subsample 1 J 2 0.782 0.923 110 13 20 1.000 0.900
subsample 1 J 3 0.492 0.762 128 21 21 0.762 0.762
Mean of session fractions (expected tuples G 0.30, J 0.70; commands 0.65, 0.90; neurons 0.40, 0.80; 10% of each value):
tuples commands neurons_all_tuples neurons_sig_tuples
null min_conds monkey
shuffle 1 G 0.373 0.672 0.401 0.461
J 0.707 0.921 0.903 0.890
2 G 0.373 0.699 0.403 0.455
J 0.709 0.950 0.890 0.890
subsample 1 G 0.374 0.672 0.404 0.471
J 0.707 0.921 0.903 0.878
2 G 0.375 0.700 0.409 0.471
J 0.710 0.950 0.890 0.878
null min_conds neurons values_within_tolerance
shuffle 1 neurons_all_tuples 4
shuffle 1 neurons_sig_tuples 3
shuffle 2 neurons_all_tuples 4
shuffle 2 neurons_sig_tuples 3
subsample 1 neurons_all_tuples 4
subsample 1 neurons_sig_tuples 4
subsample 2 neurons_all_tuples 4
subsample 2 neurons_sig_tuples 4
Second reading, fraction pooled over all items of the monkey's sessions:
tuples commands neurons_all_tuples neurons_sig_tuples
null min_conds monkey
shuffle 1 G 0.389 0.679 0.378 0.433
J 0.695 0.909 0.901 0.889
2 G 0.391 0.714 0.380 0.428
J 0.698 0.937 0.889 0.889
subsample 1 G 0.390 0.679 0.375 0.446
J 0.695 0.909 0.901 0.877
2 G 0.392 0.714 0.378 0.446
J 0.698 0.937 0.889 0.877
null min_conds neurons values_within_tolerance
shuffle 1 neurons_all_tuples 4
shuffle 1 neurons_sig_tuples 4
shuffle 2 neurons_all_tuples 4
shuffle 2 neurons_sig_tuples 4
subsample 1 neurons_all_tuples 4
subsample 1 neurons_sig_tuples 4
subsample 2 neurons_all_tuples 4
subsample 2 neurons_sig_tuples 4
Result: partially reproduced. Under the best combinations of readings, four of the six fractions are within tolerance: commands 0.67 to 0.70 (G) and 0.92 to 0.95 (J), neurons 0.40 to 0.41 (G), tuples 0.71 (J). The fraction of significant tuples for Monkey G is 0.37, above the 0.27 to 0.33 of the claim, and the fraction of neurons for Monkey J is 0.89 or 0.90 when neurons are aggregated over all tuples (0.878 over population-significant tuples with the subsampling null, inside 0.72 to 0.88, where Monkey G's neurons are then 0.47, outside). Every fraction is far above the 0.05 chance level. The second reading, the fraction pooled over all items of a monkey's sessions, also puts four of the six values within tolerance under every combination of readings; Monkey G's fraction of tuples is 0.389 to 0.392 under it.
Dynamics models: fig4c, figs4b, figs4g, fig5c, fig5e¶
A linear dynamics model fit on all conditions predicts, on held-out data, the activity issuing a given command from the previous activity better than shuffle dynamics, and generalizes to left-out commands and left-out conditions.
Dynamics models trained without all commands of one angle bin, or of one magnitude bin, still predict the left-out activity (given the command) better than shuffle dynamics.
Dynamics models trained without whole classes of conditions (vertical, horizontal or diagonal targets, both tasks) still predict the left-out activity better than shuffle dynamics.
The full dynamics model, its command and condition left-out versions, and dynamics restricted to the decoder-null space all predict the next neural activity from the current activity better than shuffle dynamics.
The full dynamics model and its command and condition left-out versions predict the next command from the current activity better than shuffle dynamics, while decoder-null dynamics give no prediction of it.
Which pairs an R2 is scored on. The main reading is the code's: each model's test pairs whose commands at t and t-1 are both below the 95th percentile of magnitude (generalization_plots.check_ix_lo, "Remove mag command 4"), with the full model and every shuffle model scored on exactly the same pairs (plot_loo_r2_overall, which scores the full model and the shuffles on the left-out model's test bins). Two other readings are shown: "analysed", pairs whose command at t (or at t-1 for next predictions) is below the 95th percentile, and "all", every test pair. Under every reading the shuffle is scored on the model's own pairs.
In [15]:
print(pd.DataFrame([dict(monkey=k[0], session=k[1], alpha=RR[k]["alpha"], alpha_one_shuffle=RR[k]["alpha_shuffle"], alpha_null=RR[k]["null_alpha"],
pairs=len(RR[k]["resid"])) for k in KEYS]).to_string(index=False))
READINGS = ["code", "analysed", "all"]
def r2_table(model, metric, full_model=False):
"""R2 of a model per session and pooled, against the shuffle models' R2 on exactly the pairs the
model is scored on, under the three readings of the scored pairs ("code": test pairs whose
commands at t and t-1 are both analysed, generalization_plots.check_ix_lo; "analysed": the
command at t, or at t-1 for next predictions, analysed; "all": every test pair). With
full_model=True, the full model's R2 on the same pairs as the named model."""
name = model + (" [full model]" if full_model else "")
rows = []
for reading in READINGS:
for mk, ks in [("G", G), ("J", J)]:
so, sn = [], []
for k in ks:
v = RR[k]["r2tab"].get((name, metric, reading), np.nan)
if not np.isfinite(v):
continue
nl = RR[k]["shuffle"]["masked"][(model, metric, reading)]
so.append(v)
sn.append(nl)
rows.append(dict(model=name, reading=reading, monkey=mk, session=k[1], r2=v,
pairs=int(RR[k]["score_masks"][(model, metric, reading)].sum()),
shuffle_95=np.percentile(nl, 95), shuffle_mean=nl.mean(), p=(nl >= v).mean()))
if not so:
continue
o, p, nl = pval_mean(so, sn)
rows.append(dict(model=name, reading=reading, monkey=mk, session="pooled", r2=o,
shuffle_95=np.percentile(nl, 95), shuffle_mean=nl.mean(), p=p))
return rows
def summarize(rows):
df = pd.DataFrame(rows)
s = df[df.session != "pooled"].groupby(["model", "reading", "monkey"], sort=False).agg(n_p001=("p", lambda x: (x < 0.001).sum()), n_p05=("p", lambda x: (x < 0.05).sum()), n=("p", "size"))
pooled = df[df.session == "pooled"].set_index(["model", "reading", "monkey"])
s["mean_r2"] = pooled.r2
s["pooled_shuffle_95"] = pooled.shuffle_95
s["pooled_p"] = pooled.p
return df, s
monkey session alpha alpha_one_shuffle alpha_null pairs
G 0 2000.0 1000.0 1000.0 11684
G 1 750.0 1000.0 500.0 7059
G 2 5000.0 7000.0 5000.0 15978
G 3 2000.0 3000.0 1000.0 15396
G 4 3000.0 10000.0 3000.0 12039
G 5 1000.0 3000.0 1000.0 11298
G 6 750.0 2000.0 500.0 4380
G 7 1000.0 3000.0 1000.0 5528
G 8 1000.0 3000.0 1000.0 5273
J 0 250.0 500.0 100.0 8489
J 1 250.0 750.0 100.0 4952
J 2 100.0 750.0 50.0 5307
J 3 100.0 500.0 75.0 5839
In [16]:
# Shuffle-model ridge parameter, reading (a): selected on one shuffled dataset and applied to 200 shuffles
# (R2 only); reading (b), the real data's parameter, is the 1000-shuffle distribution used in the tables.
rows = []
for k in KEYS:
R = RR[k]
sa = R["shuffle_a"]
for metric, real in [("cond", R["r2"]["cond"]), ("next", R["r2"]["next"]), ("cmd", R["r2"]["cmd"])]:
nb_ = R["shuffle"]["r2_" + metric]
na_ = nb_ if sa is None else sa["r2_" + metric]
rows.append(dict(monkey=k[0], session=k[1], metric=metric, r2=real, alpha_b=R["alpha"], shuffle95_b=np.percentile(nb_, 95),
alpha_a=R["alpha_shuffle"], shuffle95_a=np.percentile(na_, 95), p_a=(na_ >= real).mean()))
SA = pd.DataFrame(rows)
print(SA.round(4).to_string(index=False))
print("Sessions where the full model does not beat the 95th percentile under reading (a):", SA[SA.r2 <= SA.shuffle95_a][["monkey", "session", "metric"]].values.tolist())
monkey session metric r2 alpha_b shuffle95_b alpha_a shuffle95_a p_a
G 0 cond 0.3657 2000.0 0.3195 1000.0 0.3192 0.0
G 0 next 0.2406 2000.0 0.1791 1000.0 0.1779 0.0
G 0 cmd 0.4594 2000.0 0.4204 1000.0 0.4193 0.0
G 1 cond 0.4258 750.0 0.3757 1000.0 0.3760 0.0
G 1 next 0.2837 750.0 0.2137 1000.0 0.2141 0.0
G 1 cmd 0.4181 750.0 0.3820 1000.0 0.3817 0.0
G 2 cond 0.1463 5000.0 0.1075 7000.0 0.1079 0.0
G 2 next 0.0916 5000.0 0.0412 7000.0 0.0412 0.0
G 2 cmd 0.3189 5000.0 0.2608 7000.0 0.2566 0.0
G 3 cond 0.1411 2000.0 0.1072 3000.0 0.1075 0.0
G 3 next 0.0842 2000.0 0.0421 3000.0 0.0419 0.0
G 3 cmd 0.3460 2000.0 0.2952 3000.0 0.2923 0.0
G 4 cond 0.1182 3000.0 0.0753 10000.0 0.0770 0.0
G 4 next 0.0860 3000.0 0.0280 10000.0 0.0279 0.0
G 4 cmd 0.3071 3000.0 0.2266 10000.0 0.2112 0.0
G 5 cond 0.1363 1000.0 0.1153 3000.0 0.1159 0.0
G 5 next 0.0607 1000.0 0.0367 3000.0 0.0365 0.0
G 5 cmd 0.3334 1000.0 0.3096 3000.0 0.3023 0.0
G 6 cond 0.1881 750.0 0.1481 2000.0 0.1494 0.0
G 6 next 0.1042 750.0 0.0533 2000.0 0.0525 0.0
G 6 cmd 0.3234 750.0 0.2849 2000.0 0.2690 0.0
G 7 cond 0.1157 1000.0 0.0871 3000.0 0.0889 0.0
G 7 next 0.0642 1000.0 0.0297 3000.0 0.0295 0.0
G 7 cmd 0.3178 1000.0 0.2786 3000.0 0.2566 0.0
G 8 cond 0.1104 1000.0 0.0855 3000.0 0.0870 0.0
G 8 next 0.0596 1000.0 0.0292 3000.0 0.0280 0.0
G 8 cmd 0.3029 1000.0 0.2660 3000.0 0.2404 0.0
J 0 cond 0.3297 250.0 0.2704 500.0 0.2706 0.0
J 0 next 0.1640 250.0 0.0938 500.0 0.0935 0.0
J 0 cmd 0.3109 250.0 0.2811 500.0 0.2787 0.0
J 1 cond 0.3104 250.0 0.2451 750.0 0.2455 0.0
J 1 next 0.1610 250.0 0.0847 750.0 0.0835 0.0
J 1 cmd 0.3090 250.0 0.2790 750.0 0.2726 0.0
J 2 cond 0.3179 100.0 0.2616 750.0 0.2625 0.0
J 2 next 0.1672 100.0 0.0972 750.0 0.0951 0.0
J 2 cmd 0.3310 100.0 0.3007 750.0 0.2904 0.0
J 3 cond 0.2074 100.0 0.1814 500.0 0.1824 0.0
J 3 next 0.0674 100.0 0.0366 500.0 0.0360 0.0
J 3 cmd 0.1823 100.0 0.1698 500.0 0.1609 0.0
Sessions where the full model does not beat the 95th percentile under reading (a): []
In [17]:
rows = sum([r2_table(m, "cond") for m in ["full", "command left out", "condition left out"]], [])
df4c, s4c = summarize(rows)
print(df4c.round(4).to_string(index=False))
print("\nExpected R2: full G 0.167, J 0.252; command left out G 0.163, J 0.243; condition left out G 0.163, J 0.240 (10%)")
print(s4c.round(4).to_string())
print("\nReference, prediction from the command alone (mean over sessions):", {mk: round(np.mean([RR[k]["r2"]["conly"] for k in ks]), 4) for mk, ks in [("G", G), ("J", J)]})
_, sfull = summarize(sum([r2_table(m, "cond", full_model=True) for m in ["command left out", "condition left out"]], []))
print("\nThe full model scored on each left-out model's pairs, against the shuffle on the same pairs:")
print(sfull.round(4).to_string())
model reading monkey session r2 pairs shuffle_95 shuffle_mean p
full code G 0 0.2998 10734.0 0.2505 0.2485 0.0
full code G 1 0.3571 6433.0 0.3036 0.3009 0.0
full code G 2 0.1244 14510.0 0.0851 0.0839 0.0
full code G 3 0.1241 13963.0 0.0916 0.0905 0.0
full code G 4 0.1055 10911.0 0.0648 0.0633 0.0
full code G 5 0.1148 10265.0 0.0936 0.0928 0.0
full code G 6 0.1629 3965.0 0.1249 0.1224 0.0
full code G 7 0.1034 5006.0 0.0754 0.0737 0.0
full code G 8 0.0974 4778.0 0.0737 0.0721 0.0
full code G pooled 0.1655 NaN 0.1282 0.1276 0.0
full code J 0 0.2906 7754.0 0.2293 0.2264 0.0
full code J 1 0.2682 4529.0 0.2044 0.2006 0.0
full code J 2 0.2706 4842.0 0.2149 0.2115 0.0
full code J 3 0.1675 5333.0 0.1424 0.1403 0.0
full code J pooled 0.2492 NaN 0.1963 0.1947 0.0
full analysed G 0 0.3101 11085.0 0.2606 0.2585 0.0
full analysed G 1 0.3730 6697.0 0.3181 0.3154 0.0
full analysed G 2 0.1273 15152.0 0.0876 0.0864 0.0
full analysed G 3 0.1268 14604.0 0.0934 0.0923 0.0
full analysed G 4 0.1079 11425.0 0.0663 0.0649 0.0
full analysed G 5 0.1172 10720.0 0.0959 0.0950 0.0
full analysed G 6 0.1673 4153.0 0.1282 0.1256 0.0
full analysed G 7 0.1052 5243.0 0.0766 0.0748 0.0
full analysed G 8 0.0997 5004.0 0.0753 0.0737 0.0
full analysed G pooled 0.1705 NaN 0.1324 0.1318 0.0
full analysed J 0 0.2984 8088.0 0.2356 0.2326 0.0
full analysed J 1 0.2772 4717.0 0.2119 0.2080 0.0
full analysed J 2 0.2805 5046.0 0.2225 0.2190 0.0
full analysed J 3 0.1735 5564.0 0.1469 0.1447 0.0
full analysed J pooled 0.2574 NaN 0.2027 0.2011 0.0
full all G 0 0.3657 11684.0 0.3195 0.3176 0.0
full all G 1 0.4258 7059.0 0.3757 0.3732 0.0
full all G 2 0.1463 15978.0 0.1075 0.1063 0.0
full all G 3 0.1411 15396.0 0.1072 0.1061 0.0
full all G 4 0.1182 12039.0 0.0753 0.0738 0.0
full all G 5 0.1363 11298.0 0.1153 0.1144 0.0
full all G 6 0.1881 4380.0 0.1481 0.1454 0.0
full all G 7 0.1157 5528.0 0.0871 0.0854 0.0
full all G 8 0.1104 5273.0 0.0855 0.0838 0.0
full all G pooled 0.1942 NaN 0.1568 0.1562 0.0
full all J 0 0.3297 8489.0 0.2704 0.2676 0.0
full all J 1 0.3104 4952.0 0.2451 0.2413 0.0
full all J 2 0.3179 5307.0 0.2616 0.2583 0.0
full all J 3 0.2074 5839.0 0.1814 0.1793 0.0
full all J pooled 0.2913 NaN 0.2382 0.2366 0.0
command left out code G 0 0.2962 10734.0 0.2505 0.2485 0.0
command left out code G 1 0.3500 6433.0 0.3036 0.3009 0.0
command left out code G 2 0.1227 14510.0 0.0851 0.0839 0.0
command left out code G 3 0.1217 13963.0 0.0916 0.0905 0.0
command left out code G 4 0.1040 10911.0 0.0648 0.0633 0.0
command left out code G 5 0.1127 10265.0 0.0936 0.0928 0.0
command left out code G 6 0.1588 3965.0 0.1249 0.1224 0.0
command left out code G 7 0.1009 5006.0 0.0754 0.0737 0.0
command left out code G 8 0.0950 4778.0 0.0737 0.0721 0.0
command left out code G pooled 0.1624 NaN 0.1282 0.1276 0.0
command left out code J 0 0.2838 7754.0 0.2293 0.2264 0.0
command left out code J 1 0.2604 4529.0 0.2044 0.2006 0.0
command left out code J 2 0.2641 4842.0 0.2149 0.2115 0.0
command left out code J 3 0.1588 5333.0 0.1424 0.1403 0.0
command left out code J pooled 0.2418 NaN 0.1963 0.1947 0.0
command left out analysed G 0 0.3064 11085.0 0.2606 0.2585 0.0
command left out analysed G 1 0.3657 6697.0 0.3181 0.3154 0.0
command left out analysed G 2 0.1257 15152.0 0.0876 0.0864 0.0
command left out analysed G 3 0.1244 14604.0 0.0934 0.0923 0.0
command left out analysed G 4 0.1064 11425.0 0.0663 0.0649 0.0
command left out analysed G 5 0.1150 10720.0 0.0959 0.0950 0.0
command left out analysed G 6 0.1631 4153.0 0.1282 0.1256 0.0
command left out analysed G 7 0.1025 5243.0 0.0766 0.0748 0.0
command left out analysed G 8 0.0973 5004.0 0.0753 0.0737 0.0
command left out analysed G pooled 0.1674 NaN 0.1324 0.1318 0.0
command left out analysed J 0 0.2916 8088.0 0.2356 0.2326 0.0
command left out analysed J 1 0.2690 4717.0 0.2119 0.2080 0.0
command left out analysed J 2 0.2739 5046.0 0.2225 0.2190 0.0
command left out analysed J 3 0.1645 5564.0 0.1469 0.1447 0.0
command left out analysed J pooled 0.2497 NaN 0.2027 0.2011 0.0
command left out all G 0 0.3064 11085.0 0.2606 0.2585 0.0
command left out all G 1 0.3657 6697.0 0.3181 0.3154 0.0
command left out all G 2 0.1257 15152.0 0.0876 0.0864 0.0
command left out all G 3 0.1244 14604.0 0.0934 0.0923 0.0
command left out all G 4 0.1064 11425.0 0.0663 0.0649 0.0
command left out all G 5 0.1150 10720.0 0.0959 0.0950 0.0
command left out all G 6 0.1631 4153.0 0.1282 0.1256 0.0
command left out all G 7 0.1025 5243.0 0.0766 0.0748 0.0
command left out all G 8 0.0973 5004.0 0.0753 0.0737 0.0
command left out all G pooled 0.1674 NaN 0.1324 0.1318 0.0
command left out all J 0 0.2916 8088.0 0.2356 0.2326 0.0
command left out all J 1 0.2690 4717.0 0.2119 0.2080 0.0
command left out all J 2 0.2739 5046.0 0.2225 0.2190 0.0
command left out all J 3 0.1645 5564.0 0.1469 0.1447 0.0
command left out all J pooled 0.2497 NaN 0.2027 0.2011 0.0
condition left out code G 0 0.2972 10734.0 0.2505 0.2485 0.0
condition left out code G 1 0.3523 6433.0 0.3036 0.3009 0.0
condition left out code G 2 0.1223 14510.0 0.0851 0.0839 0.0
condition left out code G 3 0.1222 13963.0 0.0916 0.0905 0.0
condition left out code G 4 0.1036 10911.0 0.0648 0.0633 0.0
condition left out code G 5 0.1135 10265.0 0.0936 0.0928 0.0
condition left out code G 6 0.1599 3965.0 0.1249 0.1224 0.0
condition left out code G 7 0.1013 5006.0 0.0754 0.0737 0.0
condition left out code G 8 0.0952 4778.0 0.0737 0.0721 0.0
condition left out code G pooled 0.1631 NaN 0.1282 0.1276 0.0
condition left out code J 0 0.2817 7754.0 0.2293 0.2264 0.0
condition left out code J 1 0.2521 4529.0 0.2044 0.2006 0.0
condition left out code J 2 0.2612 4842.0 0.2149 0.2115 0.0
condition left out code J 3 0.1613 5333.0 0.1424 0.1403 0.0
condition left out code J pooled 0.2391 NaN 0.1963 0.1947 0.0
condition left out analysed G 0 0.3074 11085.0 0.2606 0.2585 0.0
condition left out analysed G 1 0.3679 6697.0 0.3181 0.3154 0.0
condition left out analysed G 2 0.1251 15152.0 0.0876 0.0864 0.0
condition left out analysed G 3 0.1249 14604.0 0.0934 0.0923 0.0
condition left out analysed G 4 0.1060 11425.0 0.0663 0.0649 0.0
condition left out analysed G 5 0.1159 10720.0 0.0959 0.0950 0.0
condition left out analysed G 6 0.1642 4153.0 0.1282 0.1256 0.0
condition left out analysed G 7 0.1030 5243.0 0.0766 0.0748 0.0
condition left out analysed G 8 0.0974 5004.0 0.0753 0.0737 0.0
condition left out analysed G pooled 0.1680 NaN 0.1324 0.1318 0.0
condition left out analysed J 0 0.2894 8088.0 0.2356 0.2326 0.0
condition left out analysed J 1 0.2610 4717.0 0.2119 0.2080 0.0
condition left out analysed J 2 0.2710 5046.0 0.2225 0.2190 0.0
condition left out analysed J 3 0.1670 5564.0 0.1469 0.1447 0.0
condition left out analysed J pooled 0.2471 NaN 0.2027 0.2011 0.0
condition left out all G 0 0.3628 11684.0 0.3195 0.3176 0.0
condition left out all G 1 0.4206 7059.0 0.3757 0.3732 0.0
condition left out all G 2 0.1438 15978.0 0.1075 0.1063 0.0
condition left out all G 3 0.1391 15396.0 0.1072 0.1061 0.0
condition left out all G 4 0.1162 12039.0 0.0753 0.0738 0.0
condition left out all G 5 0.1350 11298.0 0.1153 0.1144 0.0
condition left out all G 6 0.1852 4380.0 0.1481 0.1454 0.0
condition left out all G 7 0.1135 5528.0 0.0871 0.0854 0.0
condition left out all G 8 0.1079 5273.0 0.0855 0.0838 0.0
condition left out all G pooled 0.1916 NaN 0.1568 0.1562 0.0
condition left out all J 0 0.3206 8489.0 0.2704 0.2676 0.0
condition left out all J 1 0.2938 4952.0 0.2451 0.2413 0.0
condition left out all J 2 0.3081 5307.0 0.2616 0.2583 0.0
condition left out all J 3 0.2007 5839.0 0.1814 0.1793 0.0
condition left out all J pooled 0.2808 NaN 0.2382 0.2366 0.0
Expected R2: full G 0.167, J 0.252; command left out G 0.163, J 0.243; condition left out G 0.163, J 0.240 (10%)
n_p001 n_p05 n mean_r2 pooled_shuffle_95 pooled_p
model reading monkey
full code G 9 9 9 0.1655 0.1282 0.0
J 4 4 4 0.2492 0.1963 0.0
analysed G 9 9 9 0.1705 0.1324 0.0
J 4 4 4 0.2574 0.2027 0.0
all G 9 9 9 0.1942 0.1568 0.0
J 4 4 4 0.2913 0.2382 0.0
command left out code G 9 9 9 0.1624 0.1282 0.0
J 4 4 4 0.2418 0.1963 0.0
analysed G 9 9 9 0.1674 0.1324 0.0
J 4 4 4 0.2497 0.2027 0.0
all G 9 9 9 0.1674 0.1324 0.0
J 4 4 4 0.2497 0.2027 0.0
condition left out code G 9 9 9 0.1631 0.1282 0.0
J 4 4 4 0.2391 0.1963 0.0
analysed G 9 9 9 0.1680 0.1324 0.0
J 4 4 4 0.2471 0.2027 0.0
all G 9 9 9 0.1916 0.1568 0.0
J 4 4 4 0.2808 0.2382 0.0
Reference, prediction from the command alone (mean over sessions): {'G': np.float64(0.1557), 'J': np.float64(0.2375)}
The full model scored on each left-out model's pairs, against the shuffle on the same pairs:
n_p001 n_p05 n mean_r2 pooled_shuffle_95 pooled_p
model reading monkey
command left out [full model] code G 9 9 9 0.1655 0.1282 0.0
J 4 4 4 0.2492 0.1963 0.0
analysed G 9 9 9 0.1705 0.1324 0.0
J 4 4 4 0.2574 0.2027 0.0
all G 9 9 9 0.1705 0.1324 0.0
J 4 4 4 0.2574 0.2027 0.0
condition left out [full model] code G 9 9 9 0.1655 0.1282 0.0
J 4 4 4 0.2492 0.1963 0.0
analysed G 9 9 9 0.1705 0.1324 0.0
J 4 4 4 0.2574 0.2027 0.0
all G 9 9 9 0.1942 0.1568 0.0
J 4 4 4 0.2913 0.2382 0.0
Result: reproduced. Scored as the code scores it, R2 is 0.1655 (G) and 0.2492 (J) for the full model, 0.1624 and 0.2418 with the command left out, and 0.1631 and 0.2391 with the condition left out, all within 10% of the paper's 0.167, 0.252, 0.163, 0.243, 0.163 and 0.240. Each model beats the shuffle, scored on the same pairs, with p < 0.001 in 9 of 9 and 4 of 4 sessions and pooled. The shuffle's pooled 95th percentiles under this reading, 0.1282 and 0.1963, match Table S1's 0.130 and 0.196. On "analysed" pairs the tests are the same and the R2 values are within tolerance (0.1705 and 0.2574 for the full model). On all pairs the tests are the same (the command left-out model beats the shuffle in 9 of 9 and 4 of 4 sessions), but the full and condition left-out R2, 0.1942 and 0.2913, 0.1916 and 0.2808, are outside tolerance.
In [18]:
rows = sum([r2_table(m, "cond") for m in ["angle bin left out", "magnitude bin left out"]], [])
dfs4b, ss4b = summarize(rows)
print(dfs4b.round(4).to_string(index=False))
print("\nExpected R2: angle G 0.155, J 0.216; magnitude G 0.159, J 0.238 (10%)")
print(ss4b.round(4).to_string())
_, sfull = summarize(sum([r2_table(m, "cond", full_model=True) for m in ["angle bin left out", "magnitude bin left out"]], []))
print("\nThe full model scored on each left-out model's pairs, against the shuffle on the same pairs:")
print(sfull.round(4).to_string())
model reading monkey session r2 pairs shuffle_95 shuffle_mean p
angle bin left out code G 0 0.2928 10734.0 0.2505 0.2485 0.000
angle bin left out code G 1 0.3398 6433.0 0.3036 0.3009 0.000
angle bin left out code G 2 0.1185 14510.0 0.0851 0.0839 0.000
angle bin left out code G 3 0.1168 13963.0 0.0916 0.0905 0.000
angle bin left out code G 4 0.0995 10911.0 0.0648 0.0633 0.000
angle bin left out code G 5 0.1093 10265.0 0.0936 0.0928 0.000
angle bin left out code G 6 0.1519 3965.0 0.1249 0.1224 0.000
angle bin left out code G 7 0.0958 5006.0 0.0754 0.0737 0.000
angle bin left out code G 8 0.0894 4778.0 0.0737 0.0721 0.000
angle bin left out code G pooled 0.1571 NaN 0.1282 0.1276 0.000
angle bin left out code J 0 0.2658 7754.0 0.2293 0.2264 0.000
angle bin left out code J 1 0.2397 4529.0 0.2044 0.2006 0.000
angle bin left out code J 2 0.2498 4842.0 0.2149 0.2115 0.000
angle bin left out code J 3 0.1409 5333.0 0.1424 0.1403 0.321
angle bin left out code J pooled 0.2241 NaN 0.1963 0.1947 0.000
angle bin left out analysed G 0 0.3030 11085.0 0.2606 0.2585 0.000
angle bin left out analysed G 1 0.3559 6697.0 0.3181 0.3154 0.000
angle bin left out analysed G 2 0.1214 15152.0 0.0876 0.0864 0.000
angle bin left out analysed G 3 0.1195 14604.0 0.0934 0.0923 0.000
angle bin left out analysed G 4 0.1019 11425.0 0.0663 0.0649 0.000
angle bin left out analysed G 5 0.1116 10720.0 0.0959 0.0950 0.000
angle bin left out analysed G 6 0.1561 4153.0 0.1282 0.1256 0.000
angle bin left out analysed G 7 0.0973 5243.0 0.0766 0.0748 0.000
angle bin left out analysed G 8 0.0917 5004.0 0.0753 0.0737 0.000
angle bin left out analysed G pooled 0.1621 NaN 0.1324 0.1318 0.000
angle bin left out analysed J 0 0.2739 8088.0 0.2356 0.2326 0.000
angle bin left out analysed J 1 0.2488 4717.0 0.2119 0.2080 0.000
angle bin left out analysed J 2 0.2600 5046.0 0.2225 0.2190 0.000
angle bin left out analysed J 3 0.1470 5564.0 0.1469 0.1447 0.044
angle bin left out analysed J pooled 0.2324 NaN 0.2027 0.2011 0.000
angle bin left out all G 0 0.3030 11085.0 0.2606 0.2585 0.000
angle bin left out all G 1 0.3559 6697.0 0.3181 0.3154 0.000
angle bin left out all G 2 0.1214 15152.0 0.0876 0.0864 0.000
angle bin left out all G 3 0.1195 14604.0 0.0934 0.0923 0.000
angle bin left out all G 4 0.1019 11425.0 0.0663 0.0649 0.000
angle bin left out all G 5 0.1116 10720.0 0.0959 0.0950 0.000
angle bin left out all G 6 0.1561 4153.0 0.1282 0.1256 0.000
angle bin left out all G 7 0.0973 5243.0 0.0766 0.0748 0.000
angle bin left out all G 8 0.0917 5004.0 0.0753 0.0737 0.000
angle bin left out all G pooled 0.1621 NaN 0.1324 0.1318 0.000
angle bin left out all J 0 0.2739 8088.0 0.2356 0.2326 0.000
angle bin left out all J 1 0.2488 4717.0 0.2119 0.2080 0.000
angle bin left out all J 2 0.2600 5046.0 0.2225 0.2190 0.000
angle bin left out all J 3 0.1470 5564.0 0.1469 0.1447 0.044
angle bin left out all J pooled 0.2324 NaN 0.2027 0.2011 0.000
magnitude bin left out code G 0 0.2912 10734.0 0.2505 0.2485 0.000
magnitude bin left out code G 1 0.3438 6433.0 0.3036 0.3009 0.000
magnitude bin left out code G 2 0.1198 14510.0 0.0851 0.0839 0.000
magnitude bin left out code G 3 0.1194 13963.0 0.0916 0.0905 0.000
magnitude bin left out code G 4 0.1023 10911.0 0.0648 0.0633 0.000
magnitude bin left out code G 5 0.1100 10265.0 0.0936 0.0928 0.000
magnitude bin left out code G 6 0.1546 3965.0 0.1249 0.1224 0.000
magnitude bin left out code G 7 0.0971 5006.0 0.0754 0.0737 0.000
magnitude bin left out code G 8 0.0924 4778.0 0.0737 0.0721 0.000
magnitude bin left out code G pooled 0.1589 NaN 0.1282 0.1276 0.000
magnitude bin left out code J 0 0.2777 7754.0 0.2293 0.2264 0.000
magnitude bin left out code J 1 0.2532 4529.0 0.2044 0.2006 0.000
magnitude bin left out code J 2 0.2592 4842.0 0.2149 0.2115 0.000
magnitude bin left out code J 3 0.1566 5333.0 0.1424 0.1403 0.000
magnitude bin left out code J pooled 0.2367 NaN 0.1963 0.1947 0.000
magnitude bin left out analysed G 0 0.3016 11085.0 0.2606 0.2585 0.000
magnitude bin left out analysed G 1 0.3596 6697.0 0.3181 0.3154 0.000
magnitude bin left out analysed G 2 0.1227 15152.0 0.0876 0.0864 0.000
magnitude bin left out analysed G 3 0.1221 14604.0 0.0934 0.0923 0.000
magnitude bin left out analysed G 4 0.1046 11425.0 0.0663 0.0649 0.000
magnitude bin left out analysed G 5 0.1124 10720.0 0.0959 0.0950 0.000
magnitude bin left out analysed G 6 0.1592 4153.0 0.1282 0.1256 0.000
magnitude bin left out analysed G 7 0.0987 5243.0 0.0766 0.0748 0.000
magnitude bin left out analysed G 8 0.0947 5004.0 0.0753 0.0737 0.000
magnitude bin left out analysed G pooled 0.1640 NaN 0.1324 0.1318 0.000
magnitude bin left out analysed J 0 0.2859 8088.0 0.2356 0.2326 0.000
magnitude bin left out analysed J 1 0.2619 4717.0 0.2119 0.2080 0.000
magnitude bin left out analysed J 2 0.2689 5046.0 0.2225 0.2190 0.000
magnitude bin left out analysed J 3 0.1626 5564.0 0.1469 0.1447 0.000
magnitude bin left out analysed J pooled 0.2448 NaN 0.2027 0.2011 0.000
magnitude bin left out all G 0 0.3016 11085.0 0.2606 0.2585 0.000
magnitude bin left out all G 1 0.3596 6697.0 0.3181 0.3154 0.000
magnitude bin left out all G 2 0.1227 15152.0 0.0876 0.0864 0.000
magnitude bin left out all G 3 0.1221 14604.0 0.0934 0.0923 0.000
magnitude bin left out all G 4 0.1046 11425.0 0.0663 0.0649 0.000
magnitude bin left out all G 5 0.1124 10720.0 0.0959 0.0950 0.000
magnitude bin left out all G 6 0.1592 4153.0 0.1282 0.1256 0.000
magnitude bin left out all G 7 0.0987 5243.0 0.0766 0.0748 0.000
magnitude bin left out all G 8 0.0947 5004.0 0.0753 0.0737 0.000
magnitude bin left out all G pooled 0.1640 NaN 0.1324 0.1318 0.000
magnitude bin left out all J 0 0.2859 8088.0 0.2356 0.2326 0.000
magnitude bin left out all J 1 0.2619 4717.0 0.2119 0.2080 0.000
magnitude bin left out all J 2 0.2689 5046.0 0.2225 0.2190 0.000
magnitude bin left out all J 3 0.1626 5564.0 0.1469 0.1447 0.000
magnitude bin left out all J pooled 0.2448 NaN 0.2027 0.2011 0.000
Expected R2: angle G 0.155, J 0.216; magnitude G 0.159, J 0.238 (10%)
n_p001 n_p05 n mean_r2 pooled_shuffle_95 pooled_p
model reading monkey
angle bin left out code G 9 9 9 0.1571 0.1282 0.0
J 3 3 4 0.2241 0.1963 0.0
analysed G 9 9 9 0.1621 0.1324 0.0
J 3 4 4 0.2324 0.2027 0.0
all G 9 9 9 0.1621 0.1324 0.0
J 3 4 4 0.2324 0.2027 0.0
magnitude bin left out code G 9 9 9 0.1589 0.1282 0.0
J 4 4 4 0.2367 0.1963 0.0
analysed G 9 9 9 0.1640 0.1324 0.0
J 4 4 4 0.2448 0.2027 0.0
all G 9 9 9 0.1640 0.1324 0.0
J 4 4 4 0.2448 0.2027 0.0
The full model scored on each left-out model's pairs, against the shuffle on the same pairs:
n_p001 n_p05 n mean_r2 pooled_shuffle_95 pooled_p
model reading monkey
angle bin left out [full model] code G 9 9 9 0.1655 0.1282 0.0
J 4 4 4 0.2492 0.1963 0.0
analysed G 9 9 9 0.1705 0.1324 0.0
J 4 4 4 0.2574 0.2027 0.0
all G 9 9 9 0.1705 0.1324 0.0
J 4 4 4 0.2574 0.2027 0.0
magnitude bin left out [full model] code G 9 9 9 0.1655 0.1282 0.0
J 4 4 4 0.2492 0.1963 0.0
analysed G 9 9 9 0.1705 0.1324 0.0
J 4 4 4 0.2574 0.2027 0.0
all G 9 9 9 0.1705 0.1324 0.0
J 4 4 4 0.2574 0.2027 0.0
Result: reproduced. Scored as the code scores it, models with an angle bin left out give R2 0.1571 (G) and 0.2241 (J) and with a magnitude bin left out 0.1589 and 0.2367, all within 10% of 0.155, 0.216, 0.159 and 0.238. Against the shuffle scored on the same pairs, the angle models beat it with p < 0.001 in 9 of 9 Monkey G and 3 of 4 Monkey J sessions (session 3, p = 0.321), the magnitude models in 9 of 9 and 4 of 4, and both pooled with p < 0.001. The other two readings give the same session counts and R2 values within tolerance (0.1621 to 0.2448).
In [19]:
rows = sum([r2_table(m, "cond") for m in ["vertical", "horizontal", "diagonal 1", "diagonal 2"]], [])
print("Pairs per class:", {f"{k[0]}{k[1]}": RR[k]["class_pairs"] for k in KEYS})
dfs4g, ss4g = summarize(rows)
print(dfs4g.round(4).to_string(index=False))
print("\nExpected R2: vertical G 0.169, J 0.251; horizontal G 0.161, J 0.251; diagonal 1 G 0.172, J 0.212; diagonal 2 G 0.168 (10%)")
print(ss4g.round(4).to_string())
_, sfull = summarize(sum([r2_table(m, "cond", full_model=True) for m in ["vertical", "horizontal", "diagonal 1", "diagonal 2"]], []))
print("\nThe full model scored on each left-out model's pairs, against the shuffle on the same pairs:")
print(sfull.round(4).to_string())
Pairs per class: {'G0': {'vertical': 2738, 'horizontal': 2923, 'diagonal 1': 3502, 'diagonal 2': 2521}, 'G1': {'vertical': 2056, 'horizontal': 1591, 'diagonal 1': 1917, 'diagonal 2': 1495}, 'G2': {'vertical': 4082, 'horizontal': 4424, 'diagonal 1': 4066, 'diagonal 2': 3406}, 'G3': {'vertical': 3941, 'horizontal': 3912, 'diagonal 1': 3995, 'diagonal 2': 3548}, 'G4': {'vertical': 3189, 'horizontal': 2710, 'diagonal 1': 3256, 'diagonal 2': 2884}, 'G5': {'vertical': 2750, 'horizontal': 2959, 'diagonal 1': 2850, 'diagonal 2': 2739}, 'G6': {'vertical': 1065, 'horizontal': 1101, 'diagonal 1': 1143, 'diagonal 2': 1071}, 'G7': {'vertical': 1405, 'horizontal': 1293, 'diagonal 1': 1551, 'diagonal 2': 1279}, 'G8': {'vertical': 1407, 'horizontal': 1087, 'diagonal 1': 1453, 'diagonal 2': 1326}, 'J0': {'vertical': 3836, 'horizontal': 2019, 'diagonal 1': 1789}, 'J1': {'vertical': 2541, 'horizontal': 1379, 'diagonal 1': 927}, 'J2': {'vertical': 1889, 'horizontal': 1787, 'diagonal 1': 1234}, 'J3': {'vertical': 1741, 'horizontal': 1952, 'diagonal 1': 1743}}
model reading monkey session r2 pairs shuffle_95 shuffle_mean p
vertical code G 0 0.2983 2512.0 0.2429 0.2405 0.0
vertical code G 1 0.3614 1864.0 0.3075 0.3040 0.0
vertical code G 2 0.1259 3652.0 0.0853 0.0840 0.0
vertical code G 3 0.1162 3621.0 0.0858 0.0846 0.0
vertical code G 4 0.1057 2898.0 0.0656 0.0641 0.0
vertical code G 5 0.1147 2518.0 0.0905 0.0894 0.0
vertical code G 6 0.1582 973.0 0.1194 0.1166 0.0
vertical code G 7 0.0910 1286.0 0.0700 0.0682 0.0
vertical code G 8 0.0875 1277.0 0.0673 0.0655 0.0
vertical code G pooled 0.1621 NaN 0.1248 0.1241 0.0
vertical code J 0 0.2812 3538.0 0.2285 0.2256 0.0
vertical code J 1 0.2660 2279.0 0.2203 0.2161 0.0
vertical code J 2 0.2633 1758.0 0.2180 0.2145 0.0
vertical code J 3 0.1725 1533.0 0.1482 0.1454 0.0
vertical code J pooled 0.2458 NaN 0.2022 0.2004 0.0
vertical analysed G 0 0.3133 2593.0 0.2573 0.2549 0.0
vertical analysed G 1 0.3800 1936.0 0.3253 0.3217 0.0
vertical analysed G 2 0.1303 3831.0 0.0893 0.0880 0.0
vertical analysed G 3 0.1186 3766.0 0.0872 0.0860 0.0
vertical analysed G 4 0.1088 3028.0 0.0668 0.0652 0.0
vertical analysed G 5 0.1172 2622.0 0.0929 0.0918 0.0
vertical analysed G 6 0.1619 1018.0 0.1227 0.1197 0.0
vertical analysed G 7 0.0918 1341.0 0.0713 0.0694 0.0
vertical analysed G 8 0.0897 1335.0 0.0688 0.0668 0.0
vertical analysed G pooled 0.1680 NaN 0.1300 0.1293 0.0
vertical analysed J 0 0.2842 3676.0 0.2311 0.2280 0.0
vertical analysed J 1 0.2714 2400.0 0.2259 0.2217 0.0
vertical analysed J 2 0.2683 1821.0 0.2219 0.2183 0.0
vertical analysed J 3 0.1776 1629.0 0.1524 0.1495 0.0
vertical analysed J pooled 0.2503 NaN 0.2062 0.2044 0.0
vertical all G 0 0.3814 2738.0 0.3295 0.3271 0.0
vertical all G 1 0.4554 2056.0 0.4085 0.4055 0.0
vertical all G 2 0.1590 4082.0 0.1198 0.1185 0.0
vertical all G 3 0.1321 3941.0 0.0995 0.0983 0.0
vertical all G 4 0.1208 3189.0 0.0780 0.0764 0.0
vertical all G 5 0.1358 2750.0 0.1119 0.1108 0.0
vertical all G 6 0.1790 1065.0 0.1402 0.1373 0.0
vertical all G 7 0.1015 1405.0 0.0814 0.0795 0.0
vertical all G 8 0.0997 1407.0 0.0788 0.0768 0.0
vertical all G pooled 0.1961 NaN 0.1596 0.1589 0.0
vertical all J 0 0.3058 3836.0 0.2559 0.2530 0.0
vertical all J 1 0.3028 2541.0 0.2600 0.2561 0.0
vertical all J 2 0.2895 1889.0 0.2437 0.2402 0.0
vertical all J 3 0.2207 1741.0 0.1977 0.1949 0.0
vertical all J pooled 0.2797 NaN 0.2377 0.2360 0.0
horizontal code G 0 0.2907 2672.0 0.2539 0.2520 0.0
horizontal code G 1 0.3171 1433.0 0.2719 0.2693 0.0
horizontal code G 2 0.1071 4123.0 0.0711 0.0697 0.0
horizontal code G 3 0.1229 3526.0 0.0914 0.0903 0.0
horizontal code G 4 0.0972 2470.0 0.0609 0.0593 0.0
horizontal code G 5 0.1061 2684.0 0.0943 0.0934 0.0
horizontal code G 6 0.1586 986.0 0.1223 0.1193 0.0
horizontal code G 7 0.1023 1137.0 0.0734 0.0712 0.0
horizontal code G 8 0.0985 972.0 0.0777 0.0755 0.0
horizontal code G pooled 0.1556 NaN 0.1229 0.1222 0.0
horizontal code J 0 0.2820 1738.0 0.2394 0.2364 0.0
horizontal code J 1 0.2364 1239.0 0.1878 0.1830 0.0
horizontal code J 2 0.2602 1546.0 0.2141 0.2101 0.0
horizontal code J 3 0.1576 1759.0 0.1439 0.1417 0.0
horizontal code J pooled 0.2341 NaN 0.1946 0.1928 0.0
horizontal analysed G 0 0.2989 2767.0 0.2631 0.2612 0.0
horizontal analysed G 1 0.3298 1504.0 0.2817 0.2790 0.0
horizontal analysed G 2 0.1092 4258.0 0.0730 0.0717 0.0
horizontal analysed G 3 0.1252 3700.0 0.0935 0.0923 0.0
horizontal analysed G 4 0.0986 2582.0 0.0623 0.0607 0.0
horizontal analysed G 5 0.1082 2807.0 0.0962 0.0953 0.0
horizontal analysed G 6 0.1638 1038.0 0.1264 0.1231 0.0
horizontal analysed G 7 0.1054 1208.0 0.0745 0.0723 0.0
horizontal analysed G 8 0.1031 1023.0 0.0804 0.0783 0.0
horizontal analysed G pooled 0.1602 NaN 0.1267 0.1260 0.0
horizontal analysed J 0 0.2975 1861.0 0.2523 0.2491 0.0
horizontal analysed J 1 0.2474 1295.0 0.1958 0.1909 0.0
horizontal analysed J 2 0.2765 1645.0 0.2271 0.2229 0.0
horizontal analysed J 3 0.1641 1844.0 0.1480 0.1456 0.0
horizontal analysed J pooled 0.2464 NaN 0.2041 0.2021 0.0
horizontal all G 0 0.3508 2923.0 0.3182 0.3164 0.0
horizontal all G 1 0.3659 1591.0 0.3208 0.3183 0.0
horizontal all G 2 0.1209 4424.0 0.0847 0.0834 0.0
horizontal all G 3 0.1389 3912.0 0.1076 0.1065 0.0
horizontal all G 4 0.1057 2710.0 0.0679 0.0663 0.0
horizontal all G 5 0.1269 2959.0 0.1151 0.1141 0.0
horizontal all G 6 0.1899 1101.0 0.1498 0.1465 0.0
horizontal all G 7 0.1187 1293.0 0.0871 0.0850 0.0
horizontal all G 8 0.1164 1087.0 0.0926 0.0904 0.0
horizontal all G pooled 0.1816 NaN 0.1481 0.1474 0.0
horizontal all J 0 0.3490 2019.0 0.3083 0.3053 0.0
horizontal all J 1 0.2849 1379.0 0.2336 0.2285 0.0
horizontal all J 2 0.3374 1787.0 0.2940 0.2900 0.0
horizontal all J 3 0.1983 1952.0 0.1847 0.1825 0.0
horizontal all J pooled 0.2924 NaN 0.2535 0.2516 0.0
diagonal 1 code G 0 0.3071 3319.0 0.2596 0.2575 0.0
diagonal 1 code G 1 0.3548 1793.0 0.3172 0.3147 0.0
diagonal 1 code G 2 0.1214 3682.0 0.0863 0.0851 0.0
diagonal 1 code G 3 0.1291 3639.0 0.0939 0.0926 0.0
diagonal 1 code G 4 0.1064 2875.0 0.0622 0.0605 0.0
diagonal 1 code G 5 0.1195 2581.0 0.0923 0.0910 0.0
diagonal 1 code G 6 0.1631 1039.0 0.1250 0.1222 0.0
diagonal 1 code G 7 0.1090 1416.0 0.0779 0.0757 0.0
diagonal 1 code G 8 0.0974 1326.0 0.0746 0.0726 0.0
diagonal 1 code G pooled 0.1675 NaN 0.1308 0.1302 0.0
diagonal 1 code J 0 0.2734 1689.0 0.1998 0.1962 0.0
diagonal 1 code J 1 0.1746 912.0 0.1237 0.1186 0.0
diagonal 1 code J 2 0.2363 1180.0 0.1786 0.1744 0.0
diagonal 1 code J 3 0.1370 1676.0 0.1188 0.1164 0.0
diagonal 1 code J pooled 0.2053 NaN 0.1534 0.1514 0.0
diagonal 1 analysed G 0 0.3143 3401.0 0.2652 0.2629 0.0
diagonal 1 analysed G 1 0.3595 1850.0 0.3224 0.3200 0.0
diagonal 1 analysed G 2 0.1237 3855.0 0.0880 0.0868 0.0
diagonal 1 analysed G 3 0.1311 3797.0 0.0952 0.0939 0.0
diagonal 1 analysed G 4 0.1103 3047.0 0.0652 0.0634 0.0
diagonal 1 analysed G 5 0.1216 2700.0 0.0945 0.0932 0.0
diagonal 1 analysed G 6 0.1697 1084.0 0.1291 0.1261 0.0
diagonal 1 analysed G 7 0.1101 1479.0 0.0783 0.0761 0.0
diagonal 1 analysed G 8 0.0991 1384.0 0.0760 0.0739 0.0
diagonal 1 analysed G pooled 0.1710 NaN 0.1336 0.1329 0.0
diagonal 1 analysed J 0 0.2855 1736.0 0.2097 0.2061 0.0
diagonal 1 analysed J 1 0.1800 920.0 0.1289 0.1238 0.0
diagonal 1 analysed J 2 0.2445 1205.0 0.1853 0.1811 0.0
diagonal 1 analysed J 3 0.1411 1707.0 0.1232 0.1209 0.0
diagonal 1 analysed J pooled 0.2128 NaN 0.1600 0.1580 0.0
diagonal 1 all G 0 0.3381 3502.0 0.2894 0.2872 0.0
diagonal 1 all G 1 0.3875 1917.0 0.3519 0.3497 0.0
diagonal 1 all G 2 0.1392 4066.0 0.1045 0.1034 0.0
diagonal 1 all G 3 0.1455 3995.0 0.1097 0.1083 0.0
diagonal 1 all G 4 0.1222 3256.0 0.0751 0.0734 0.0
diagonal 1 all G 5 0.1407 2850.0 0.1140 0.1126 0.0
diagonal 1 all G 6 0.1895 1143.0 0.1487 0.1457 0.0
diagonal 1 all G 7 0.1194 1551.0 0.0881 0.0860 0.0
diagonal 1 all G 8 0.1079 1453.0 0.0845 0.0825 0.0
diagonal 1 all G pooled 0.1878 NaN 0.1505 0.1499 0.0
diagonal 1 all J 0 0.3160 1789.0 0.2440 0.2405 0.0
diagonal 1 all J 1 0.1893 927.0 0.1382 0.1332 0.0
diagonal 1 all J 2 0.2656 1234.0 0.2072 0.2030 0.0
diagonal 1 all J 3 0.1593 1743.0 0.1413 0.1390 0.0
diagonal 1 all J pooled 0.2325 NaN 0.1809 0.1789 0.0
diagonal 2 code G 0 0.2696 2231.0 0.2192 0.2171 0.0
diagonal 2 code G 1 0.3490 1343.0 0.2880 0.2845 0.0
diagonal 2 code G 2 0.1370 3053.0 0.0973 0.0960 0.0
diagonal 2 code G 3 0.1220 3177.0 0.0936 0.0924 0.0
diagonal 2 code G 4 0.1040 2668.0 0.0686 0.0671 0.0
diagonal 2 code G 5 0.1096 2482.0 0.0948 0.0938 0.0
diagonal 2 code G 6 0.1573 967.0 0.1289 0.1263 0.0
diagonal 2 code G 7 0.0941 1167.0 0.0749 0.0732 0.0
diagonal 2 code G 8 0.0912 1203.0 0.0688 0.0669 0.0
diagonal 2 code G pooled 0.1593 NaN 0.1248 0.1241 0.0
diagonal 2 analysed G 0 0.2856 2324.0 0.2377 0.2355 0.0
diagonal 2 analysed G 1 0.3793 1407.0 0.3188 0.3153 0.0
diagonal 2 analysed G 2 0.1397 3208.0 0.0997 0.0984 0.0
diagonal 2 analysed G 3 0.1266 3341.0 0.0964 0.0952 0.0
diagonal 2 analysed G 4 0.1049 2768.0 0.0694 0.0679 0.0
diagonal 2 analysed G 5 0.1122 2591.0 0.0973 0.0963 0.0
diagonal 2 analysed G 6 0.1593 1013.0 0.1307 0.1281 0.0
diagonal 2 analysed G 7 0.0958 1215.0 0.0767 0.0749 0.0
diagonal 2 analysed G 8 0.0926 1262.0 0.0701 0.0683 0.0
diagonal 2 analysed G pooled 0.1662 NaN 0.1318 0.1311 0.0
diagonal 2 all G 0 0.3852 2521.0 0.3447 0.3428 0.0
diagonal 2 all G 1 0.4452 1495.0 0.3929 0.3899 0.0
diagonal 2 all G 2 0.1596 3406.0 0.1218 0.1205 0.0
diagonal 2 all G 3 0.1421 3548.0 0.1109 0.1097 0.0
diagonal 2 all G 4 0.1132 2884.0 0.0775 0.0760 0.0
diagonal 2 all G 5 0.1324 2739.0 0.1177 0.1167 0.0
diagonal 2 all G 6 0.1794 1071.0 0.1491 0.1466 0.0
diagonal 2 all G 7 0.1062 1279.0 0.0865 0.0846 0.0
diagonal 2 all G 8 0.1035 1326.0 0.0808 0.0790 0.0
diagonal 2 all G pooled 0.1963 NaN 0.1635 0.1629 0.0
Expected R2: vertical G 0.169, J 0.251; horizontal G 0.161, J 0.251; diagonal 1 G 0.172, J 0.212; diagonal 2 G 0.168 (10%)
n_p001 n_p05 n mean_r2 pooled_shuffle_95 pooled_p
model reading monkey
vertical code G 9 9 9 0.1621 0.1248 0.0
J 4 4 4 0.2458 0.2022 0.0
analysed G 9 9 9 0.1680 0.1300 0.0
J 4 4 4 0.2503 0.2062 0.0
all G 9 9 9 0.1961 0.1596 0.0
J 4 4 4 0.2797 0.2377 0.0
horizontal code G 9 9 9 0.1556 0.1229 0.0
J 4 4 4 0.2341 0.1946 0.0
analysed G 9 9 9 0.1602 0.1267 0.0
J 4 4 4 0.2464 0.2041 0.0
all G 9 9 9 0.1816 0.1481 0.0
J 4 4 4 0.2924 0.2535 0.0
diagonal 1 code G 9 9 9 0.1675 0.1308 0.0
J 4 4 4 0.2053 0.1534 0.0
analysed G 9 9 9 0.1710 0.1336 0.0
J 4 4 4 0.2128 0.1600 0.0
all G 9 9 9 0.1878 0.1505 0.0
J 4 4 4 0.2325 0.1809 0.0
diagonal 2 code G 9 9 9 0.1593 0.1248 0.0
analysed G 9 9 9 0.1662 0.1318 0.0
all G 9 9 9 0.1963 0.1635 0.0
The full model scored on each left-out model's pairs, against the shuffle on the same pairs:
n_p001 n_p05 n mean_r2 pooled_shuffle_95 pooled_p
model reading monkey
vertical [full model] code G 9 9 9 0.1645 0.1248 0.0
J 4 4 4 0.2577 0.2022 0.0
analysed G 9 9 9 0.1705 0.1300 0.0
J 4 4 4 0.2627 0.2062 0.0
all G 9 9 9 0.1988 0.1596 0.0
J 4 4 4 0.2930 0.2377 0.0
horizontal [full model] code G 9 9 9 0.1573 0.1229 0.0
J 4 4 4 0.2433 0.1946 0.0
analysed G 9 9 9 0.1620 0.1267 0.0
J 4 4 4 0.2566 0.2041 0.0
all G 9 9 9 0.1833 0.1481 0.0
J 4 4 4 0.3040 0.2535 0.0
diagonal 1 [full model] code G 9 9 9 0.1695 0.1308 0.0
J 4 4 4 0.2147 0.1534 0.0
analysed G 9 9 9 0.1730 0.1336 0.0
J 4 4 4 0.2221 0.1600 0.0
all G 9 9 9 0.1900 0.1505 0.0
J 4 4 4 0.2415 0.1809 0.0
diagonal 2 [full model] code G 9 9 9 0.1615 0.1248 0.0
analysed G 9 9 9 0.1685 0.1318 0.0
all G 9 9 9 0.1986 0.1635 0.0
Result: reproduced. Scored as the code scores it, R2 is 0.1621 (G) and 0.2458 (J) for the vertical class, 0.1556 and 0.2341 for horizontal, 0.1675 and 0.2053 for diagonal 1, and 0.1593 for diagonal 2 (Monkey G only), all within 10% of the paper's values. Against the shuffle scored on the same pairs, every class beats it with p < 0.001 in 9 of 9 Monkey G and 4 of 4 Monkey J sessions and pooled. The "analysed" and "all" readings give the same tests; on all pairs five of the seven R2 values are outside tolerance (for example vertical 0.1961 and 0.2797).
Result: reproduced. Scored as the code scores it, R2 of the next activity is 0.0995 (G) and 0.1115 (J) for the full model, 0.0984 and 0.1083 with the command left out, 0.0966 and 0.0992 with the condition left out, and 0.0811 and 0.0793 for decoder-null dynamics, all within 10% of the paper's values. Every model beats the shuffle, scored on the same pairs, pooled with p < 0.001, and in every session with p < 0.001 except decoder-null dynamics in Monkey G session 1 (p = 0.007): 8 of 9, as the criterion needs. Under the "analysed" and "all" readings decoder-null dynamics beat the shuffle with p < 0.001 in only 7 of 9 Monkey G sessions, and on all pairs the full model's R2, 0.1194 and 0.1399, is outside tolerance.
Result: partially reproduced. Scored as the code scores it, every R2 is within tolerance: 0.3095 (G) and 0.2297 (J) for the full model, 0.3064 and 0.2231 with the command left out, 0.3008 and 0.2096 with the condition left out, and -0.0083 and -0.0228 for decoder-null dynamics (within 0.05 of 0). The full and command left-out models beat the shuffle with p < 0.001 in 9 of 9 and 4 of 4 sessions and pooled, and decoder-null dynamics in no session. The condition left-out model beats it with p < 0.001 in 9 of 9 Monkey G sessions and pooled, but for Monkey J in 2 of 4 sessions (p = 0.560 and 1.000 in sessions 1 and 3) and pooled with p = 0.003, where the paper gives p < 0.001 pooled. On "analysed" pairs that pooled test gives p < 0.001 but Monkey J's R2 (0.2535, 0.2372, 0.2326) are above the upper bounds 0.233, 0.232 and 0.212; on all pairs they are higher still (0.2833 for the full model). Under every reading one part of the claim falls short.
fig4de-condition-specific-activity-predicted, fig4f-condition-specific-component, fig4g-majority-predicted¶
The dynamics model predicts the condition-specific average activity for a command, of single neurons and of the population, with lower error than shuffle dynamics.
Invariant dynamics predict 20% to 40% of the condition-specific component of neural activity (the part not explained by the command alone), significantly more than shuffle dynamics.
The dynamics model predicts the condition-specific activity better than shuffle dynamics for the vast majority of (command, condition) tuples, commands and neurons.
In [22]:
rows = []
for reading in ["all", "sig"]:
for mk, ks in [("G", G), ("J", J)]:
so_n, sn_n, so_p, sn_p = [], [], [], []
for k in ks:
R = RR[k]
m = R["ts_masks"][reading]
sh = R["shuffle"]
o_n = R["ts_real"]["neur"][m].mean()
n_n = sh["neur_tuple"][:, m].mean(1)
mn = sh["pop"][:, m].mean(0)
o_p = (R["ts_real"]["pop"][m] / mn).mean()
n_p = (sh["pop"][:, m] / mn).mean(1)
so_n.append(o_n); sn_n.append(n_n); so_p.append(o_p); sn_p.append(n_p)
rows.append(dict(tuples=reading, monkey=mk, session=k[1], n_tuples=int(m.sum()), neuron_err_hz=o_n, neuron_shuffle_5=np.percentile(n_n, 5),
p_neuron=(n_n <= o_n).mean(), pop_err_norm=o_p, pop_shuffle_5=np.percentile(n_p, 5), p_pop=(n_p <= o_p).mean()))
on = np.mean(so_n); nn = np.mean(np.vstack(sn_n), 0); op = np.mean(so_p); npp = np.mean(np.vstack(sn_p), 0)
rows.append(dict(tuples=reading, monkey=mk, session="pooled", neuron_err_hz=on, neuron_shuffle_5=np.percentile(nn, 5), p_neuron=(nn <= on).mean(),
pop_err_norm=op, pop_shuffle_5=np.percentile(npp, 5), p_pop=(npp <= op).mean()))
D4 = pd.DataFrame(rows)
print(D4.round(4).to_string(index=False))
print("\nExpected: single-neuron error G 1.232, J 1.182 Hz; population error over shuffle mean G 0.883, J 0.809 (10%)")
print("Tuples with a prediction (of the analysed tuples):", {f"{k[0]}{k[1]}": f"{len(RR[k]['ts_keep'])}/{len(RR[k]['tups'])}" for k in KEYS})
tuples monkey session n_tuples neuron_err_hz neuron_shuffle_5 p_neuron pop_err_norm pop_shuffle_5 p_pop
all G 0 238.0 1.5965 1.8532 0.0 0.8246 0.9954 0.0
all G 1 148.0 1.0327 1.1767 0.0 0.8371 0.9944 0.0
all G 2 216.0 0.7945 0.8662 0.0 0.8833 0.9974 0.0
all G 3 289.0 1.2273 1.3372 0.0 0.9016 0.9976 0.0
all G 4 214.0 0.8969 0.9968 0.0 0.8617 0.9962 0.0
all G 5 283.0 1.4497 1.5255 0.0 0.9446 0.9979 0.0
all G 6 78.0 1.2033 1.2834 0.0 0.9225 0.9955 0.0
all G 7 135.0 1.5447 1.6307 0.0 0.9308 0.9964 0.0
all G 8 116.0 1.4191 1.5118 0.0 0.9254 0.9959 0.0
all G pooled NaN 1.2405 1.3561 0.0 0.8924 0.9987 0.0
all J 0 154.0 1.2734 1.5405 0.0 0.8166 0.9940 0.0
all J 1 84.0 1.3458 1.6632 0.0 0.7936 0.9923 0.0
all J 2 110.0 1.2028 1.4717 0.0 0.8045 0.9925 0.0
all J 3 128.0 0.9682 1.0467 0.0 0.9183 0.9948 0.0
all J pooled NaN 1.1975 1.4347 0.0 0.8332 0.9965 0.0
sig G 0 137.0 1.6482 1.9452 0.0 0.8024 0.9950 0.0
sig G 1 66.0 0.9817 1.1571 0.0 0.8081 0.9930 0.0
sig G 2 128.0 0.7831 0.8632 0.0 0.8660 0.9970 0.0
sig G 3 114.0 1.2788 1.4040 0.0 0.8915 0.9973 0.0
sig G 4 85.0 0.8434 0.9838 0.0 0.7993 0.9936 0.0
sig G 5 57.0 1.6124 1.7465 0.0 0.9115 0.9963 0.0
sig G 6 19.0 1.3404 1.4534 0.0 0.8957 0.9933 0.0
sig G 7 27.0 1.6587 1.7543 0.0 0.9210 0.9950 0.0
sig G 8 35.0 1.5039 1.6035 0.0 0.9147 0.9949 0.0
sig G pooled NaN 1.2945 1.4385 0.0 0.8678 0.9982 0.0
sig J 0 113.0 1.2572 1.5539 0.0 0.7938 0.9936 0.0
sig J 1 69.0 1.3162 1.6506 0.0 0.7748 0.9914 0.0
sig J 2 85.0 1.2091 1.5141 0.0 0.7811 0.9920 0.0
sig J 3 64.0 1.0167 1.1242 0.0 0.8918 0.9938 0.0
sig J pooled NaN 1.1998 1.4652 0.0 0.8104 0.9962 0.0
Expected: single-neuron error G 1.232, J 1.182 Hz; population error over shuffle mean G 0.883, J 0.809 (10%)
Tuples with a prediction (of the analysed tuples): {'G0': '238/238', 'G1': '148/148', 'G2': '216/216', 'G3': '289/289', 'G4': '214/214', 'G5': '283/283', 'G6': '78/78', 'G7': '135/135', 'G8': '116/116', 'J0': '154/154', 'J1': '84/84', 'J2': '110/110', 'J3': '128/128'}
Result: reproduced. Over all tuples (reading a), the single-neuron error is 1.241 Hz (G) and 1.198 Hz (J), and the population error over the shuffle mean 0.892 and 0.833, all within 10% of the paper's values; both errors are below the shuffle with p < 0.001 in every session and pooled. Reading (b), tuples significant in fig3d, gives 1.295 and 1.200 Hz and 0.868 and 0.810, also within tolerance.
In [23]:
rows = []
for reading in ["all", "sig"]:
for mk, ks in [("G", G), ("J", J)]:
so, sn = [], []
for k in ks:
R = RR[k]
m = R["ts_masks"][reading]
v = an.component_r2(R["true_mean"][m], R["ts_real"]["pred_mean"][m], R["conly_mean"][m], BB[k]["Nb"])
nl = R["shuffle"]["comp_r2"][reading]
so.append(v); sn.append(nl)
rows.append(dict(tuples=reading, monkey=mk, session=k[1], n_tuples=int(m.sum()), r2=v, shuffle_mean=np.nanmean(nl), shuffle_95=np.nanpercentile(nl, 95), p=(nl >= v).mean()))
o, p, nl = pval_mean(so, sn)
rows.append(dict(tuples=reading, monkey=mk, session="pooled", r2=o, shuffle_mean=np.nanmean(nl), shuffle_95=np.nanpercentile(nl, 95), p=p))
D4F = pd.DataFrame(rows)
print(D4F.round(4).to_string(index=False))
print("\nExpected: G 0.226, J 0.330 (10%); every session 0.18 to 0.44; shuffle about 0 (Table S1: -0.006 G, -0.016 J)")
tuples monkey session n_tuples r2 shuffle_mean shuffle_95 p all G 0 238.0 0.3554 0.0301 0.0399 0.0 all G 1 148.0 0.3161 0.0144 0.0253 0.0 all G 2 216.0 0.2203 0.0025 0.0071 0.0 all G 3 289.0 0.1881 0.0020 0.0067 0.0 all G 4 214.0 0.2747 0.0073 0.0151 0.0 all G 5 283.0 0.1206 0.0026 0.0072 0.0 all G 6 78.0 0.1493 -0.0151 -0.0052 0.0 all G 7 135.0 0.1496 -0.0022 0.0058 0.0 all G 8 116.0 0.1600 0.0072 0.0156 0.0 all G pooled NaN 0.2149 0.0054 0.0081 0.0 all J 0 154.0 0.3261 -0.0002 0.0116 0.0 all J 1 84.0 0.3652 -0.0036 0.0132 0.0 all J 2 110.0 0.3661 -0.0094 0.0055 0.0 all J 3 128.0 0.1546 -0.0033 0.0072 0.0 all J pooled NaN 0.3030 -0.0041 0.0031 0.0 sig G 0 137.0 0.3914 0.0339 0.0442 0.0 sig G 1 66.0 0.3782 0.0215 0.0353 0.0 sig G 2 128.0 0.2459 0.0031 0.0086 0.0 sig G 3 114.0 0.2013 -0.0101 -0.0045 0.0 sig G 4 85.0 0.3519 0.0005 0.0129 0.0 sig G 5 57.0 0.1670 -0.0124 -0.0044 0.0 sig G 6 19.0 0.1678 -0.0651 -0.0492 0.0 sig G 7 27.0 0.1562 -0.0241 -0.0126 0.0 sig G 8 35.0 0.1299 -0.0331 -0.0229 0.0 sig G pooled NaN 0.2433 -0.0095 -0.0059 0.0 sig J 0 113.0 0.3516 -0.0032 0.0094 0.0 sig J 1 69.0 0.3881 0.0000 0.0175 0.0 sig J 2 85.0 0.3958 -0.0092 0.0075 0.0 sig J 3 64.0 0.1958 -0.0085 0.0045 0.0 sig J pooled NaN 0.3328 -0.0052 0.0024 0.0 Expected: G 0.226, J 0.330 (10%); every session 0.18 to 0.44; shuffle about 0 (Table S1: -0.006 G, -0.016 J)
Result: partially reproduced. The dynamics model predicts 0.215 (G) and 0.303 (J) of the condition-specific component over all tuples (reading a) and 0.243 and 0.333 over significant tuples (reading b), within 10% of 0.226 and 0.330, against shuffle means of about 0 (-0.010 to 0.005), with p < 0.001 in every session and pooled. Session values do not all lie in 0.18 to 0.44: under reading (a) Monkey G sessions 5 to 8 and Monkey J session 3 are 0.12 to 0.16, and under reading (b) Monkey G sessions 5, 6, 7 and 8 are 0.13 to 0.17.
In [24]:
rows = []
for k in KEYS:
R = RR[k]
sh = R["shuffle"]
pop, spop = R["ts_real"]["pop"], sh["pop"]
frac_t = (pop < np.percentile(spop, 5, axis=0)).mean()
tk = [R["tups"][i] for i in R["ts_keep"]]
cmds = np.array([(x["mag"], x["ang"]) for x in tk])
sig_c = []
for c in sorted({tuple(c) for c in cmds}):
cm = (cmds[:, 0] == c[0]) & (cmds[:, 1] == c[1])
sig_c.append((spop[:, cm].mean(1) <= pop[cm].mean()).mean() < 0.05)
o = R["ts_real"]["neur"].mean(0)
nl = sh["neur_by_neuron"]["all"]
frac_n = ((nl <= o[None, :]).mean(0) < 0.05).mean()
rows.append(dict(monkey=k[0], session=k[1], tuples=frac_t, commands=np.mean(sig_c), neurons=frac_n))
D4G = pd.DataFrame(rows)
print(D4G.round(3).to_string(index=False))
print("\nMean of session fractions (criterion: above 0.5 for all three in both monkeys; chance 0.05):")
print(D4G.groupby("monkey")[["tuples", "commands", "neurons"]].mean().round(3).to_string())
monkey session tuples commands neurons
G 0 0.929 1.000 0.909
G 1 0.919 1.000 0.686
G 2 0.940 1.000 0.709
G 3 0.948 1.000 0.943
G 4 0.902 0.969 0.649
G 5 0.781 1.000 0.810
G 6 0.833 1.000 0.580
G 7 0.852 0.903 0.729
G 8 0.888 0.963 0.708
J 0 0.929 1.000 0.950
J 1 0.964 1.000 0.900
J 2 0.936 1.000 0.900
J 3 0.805 0.857 0.905
Mean of session fractions (criterion: above 0.5 for all three in both monkeys; chance 0.05):
tuples commands neurons
monkey
G 0.888 0.982 0.747
J 0.908 0.964 0.914
Result: reproduced. The mean of session fractions predicted better than shuffle dynamics is 0.888 (G) and 0.908 (J) for tuples, 0.982 and 0.964 for commands, and 0.747 and 0.914 for neurons, all above 0.5 (chance 0.05).
figs3-low-dimensional-decaying-dynamics¶
The fitted dynamics are low-dimensional, with 2 to 4 eigenvalues per session whose decay time exceeds the 0.1 s bin (2 to 5 for decoder-null dynamics), and decay to a fixed point with little rotation (mean frequency about 0.1 Hz for Monkey G and 0 Hz for Monkey J).
Result: not reproduced. The dynamics are low-dimensional and decaying: 2 to 5 eigenvalues per session have decay times above 0.1 s (longest 0.38 s), and Monkey J's long-timescale eigenvalues are real in 3 of 4 sessions. But the counts fall outside the claim's ranges in some sessions (5 for the full model in Monkey G session 4, where the claim gives 2 to 4; 1 for decoder-null dynamics in Monkey J session 3, where it gives 2 to 5), Monkey G's mean frequencies, 0.083 Hz (full) and 0.068 Hz (decoder-null), are outside 0.09 to 0.11 and 0.054 to 0.066, and one Monkey J eigenvalue pair rotates at 0.036 Hz. The 99% step changes no count. Fitting the dynamics on rates in Hz (counts times 10) with the same ridge grid leaves every count of the full model's long-timescale eigenvalues unchanged and Monkey G's mean frequency at 0.083 Hz (the variant refits only the full model, so the decoder-null counts are not tested under it); the selected parameter is the counts' parameter times 100 except in Monkey G sessions 0, 3 and 4, where the grid holds a neighbouring value.
fig5d-four-dimensions¶
The full dynamics model predicts about 4 dimensions of neural activity (dynamics eigenvectors) more accurately than the rest.
In [26]:
rows = []
for k in KEYS:
r2d = RR[k]["dims_r2"]
rows.append(dict(monkey=k[0], session=k[1], n_dims=len(r2d), above_pooled_r2=int((r2d > RR[k]["r2"]["next"]).sum()),
above_shuffle_95=int((r2d > np.percentile(RR[k]["shuffle"]["r2_next"], 95)).sum()),
top8_r2=np.round(r2d[:8], 3).tolist()))
D5D = pd.DataFrame(rows)
print(D5D.to_string(index=False))
print("\nMedian over sessions:", D5D.groupby("monkey")[["above_pooled_r2", "above_shuffle_95"]].median().to_dict())
monkey session n_dims above_pooled_r2 above_shuffle_95 top8_r2
G 0 44 3 4 [0.583, 0.542, 0.282, 0.18, 0.045, 0.04, 0.029, 0.014]
G 1 35 2 3 [0.618, 0.414, 0.266, 0.076, 0.034, 0.029, 0.027, 0.005]
G 2 148 5 6 [0.504, 0.417, 0.193, 0.133, 0.1, 0.086, 0.034, 0.031]
G 3 53 4 6 [0.464, 0.437, 0.169, 0.138, 0.073, 0.053, 0.03, 0.024]
G 4 151 5 9 [0.503, 0.443, 0.256, 0.211, 0.178, 0.055, 0.05, 0.042]
G 5 42 5 7 [0.391, 0.317, 0.252, 0.115, 0.079, 0.043, 0.041, 0.026]
G 6 50 4 6 [0.382, 0.371, 0.211, 0.142, 0.099, 0.064, 0.052, 0.03]
G 7 48 4 6 [0.451, 0.391, 0.136, 0.13, 0.058, 0.04, 0.028, 0.017]
G 8 48 5 7 [0.428, 0.373, 0.134, 0.101, 0.06, 0.05, 0.037, 0.025]
J 0 20 3 4 [0.379, 0.277, 0.26, 0.118, 0.092, 0.064, 0.051, 0.025]
J 1 20 3 5 [0.416, 0.278, 0.275, 0.147, 0.104, 0.08, 0.067, 0.041]
J 2 20 3 4 [0.459, 0.332, 0.18, 0.152, 0.088, 0.062, 0.053, 0.039]
J 3 21 3 4 [0.261, 0.171, 0.083, 0.046, 0.034, 0.031, 0.031, 0.017]
Median over sessions: {'above_pooled_r2': {'G': 4.0, 'J': 3.0}, 'above_shuffle_95': {'G': 6.0, 'J': 4.0}}
Result: reproduced. Under reading (a), dimensions predicted better than the session's pooled R2, the median over sessions is 4 for Monkey G and 3 for Monkey J, inside 3 to 5. Under reading (b), dimensions above the shuffle's 95th percentile, it is 6 and 4.
fig5f-condition-specific-next-command and fig5h-turn-predicted¶
The dynamics model predicts the condition-specific average next command better than shuffle dynamics, aggregated over tuples, and for the majority of (command, condition) tuples.
The dynamics model predicts whether the next command turns clockwise or counterclockwise relative to the condition-pooled next command, and the angle of the turn, better than shuffle dynamics.
Result: reproduced. The model predicts the condition-specific next command better than shuffle dynamics with p < 0.001 in every session and pooled, and for 0.805 (G) and 0.617 (J) of tuples on average, a majority. The errors, taken as the code prints them (ten times the command's own units, as the amended criterion says), are 3.973 (G) and 7.052 (J), within tolerance of 3.956 and 7.324; in the command's own units they are 0.397 and 0.705. A reader's answer (questions/000404.yaml) settled the units; before it the claim was recorded as partially reproduced.
In [28]:
rows = []
for mk, ks in [("G", G), ("J", J)]:
s_o, s_n, a_o, a_n = [], [], [], []
for k in ks:
R = RR[k]
so_, sn_ = R["ts_real"]["sign_ok"].mean(), R["shuffle"]["sign_ok"].mean(1)
ao_, an_ = R["ts_real"]["ang_err"].mean(), R["shuffle"]["ang_err"].mean(1)
s_o.append(so_); s_n.append(sn_); a_o.append(ao_); a_n.append(an_)
rows.append(dict(monkey=mk, session=k[1], sign_frac=so_, sign_shuffle_95=np.percentile(sn_, 95), p_sign=(sn_ >= so_).mean(),
angle_err=ao_, angle_shuffle_5=np.percentile(an_, 5), p_angle=(an_ <= ao_).mean()))
o1, p1, n1 = pval_mean(s_o, s_n)
o2, p2, n2 = pval_mean(a_o, a_n, greater=False)
rows.append(dict(monkey=mk, session="pooled", sign_frac=o1, sign_shuffle_95=np.percentile(n1, 95), p_sign=p1, angle_err=o2, angle_shuffle_5=np.percentile(n2, 5), p_angle=p2))
D5H = pd.DataFrame(rows)
print(D5H.round(4).to_string(index=False))
print("\nExpected: sign G 0.708, J 0.617; angle error G 19.503, J 9.608 degrees (10%)")
monkey session sign_frac sign_shuffle_95 p_sign angle_err angle_shuffle_5 p_angle
G 0 0.7437 0.6639 0.000 22.6074 26.9089 0.0
G 1 0.6419 0.4324 0.000 20.3634 28.3347 0.0
G 2 0.7130 0.5139 0.000 21.7422 30.1117 0.0
G 3 0.7785 0.6021 0.000 22.9939 28.9508 0.0
G 4 0.7383 0.5093 0.000 21.1596 32.3617 0.0
G 5 0.6184 0.5795 0.000 24.4155 25.5577 0.0
G 6 0.7308 0.6923 0.005 13.9898 14.9543 0.0
G 7 0.7333 0.6889 0.001 15.3963 16.7989 0.0
G 8 0.6293 0.4655 0.000 18.4349 22.7073 0.0
G pooled 0.7030 0.5587 0.000 20.1226 25.4560 0.0
J 0 0.6039 0.5260 0.000 8.4344 9.7722 0.0
J 1 0.5952 0.5476 0.000 9.0653 10.6203 0.0
J 2 0.7091 0.5818 0.000 7.6921 9.9536 0.0
J 3 0.5156 0.4297 0.000 11.7122 12.3426 0.0
J pooled 0.6060 0.5167 0.000 9.2260 10.7885 0.0
Expected: sign G 0.708, J 0.617; angle error G 19.503, J 9.608 degrees (10%)
Result: partially reproduced. The four values are within 10% of the paper's: the sign of the turn is predicted for 0.703 (G) and 0.606 (J) of tuples, and the angle error is 20.1 and 9.2 degrees. Both beat the shuffle pooled with p < 0.001 and the angle in every session; the sign beats it with p < 0.001 in 7 of 9 Monkey G sessions (p = 0.005 and 0.001 in sessions 6 and 7), where the criterion needs 8, and in 4 of 4 Monkey J sessions.
figs5b-dynamics-beyond-behavior¶
Adding previous neural activity (invariant dynamics) to a model that encodes the command, cursor state, target and task significantly improves the held-out prediction of neural activity.
In [29]:
rows = [dict(monkey=k[0], session=k[1], behavior=RR[k]["s5b"]["behavior"], behavior_plus_dynamics=RR[k]["s5b"]["behavior_plus_dynamics"]) for k in KEYS]
D5B = pd.DataFrame(rows)
print(D5B.round(4).to_string(index=False))
for mk in ["G", "J"]:
x = D5B[D5B.monkey == mk]
t, p = scipy.stats.ttest_rel(x.behavior, x.behavior_plus_dynamics)
print(f"Monkey {mk}: paired t-test over {len(x)} sessions, T = {t:.3f}, p = {p:.3g}; mean gain {np.mean(x.behavior_plus_dynamics - x.behavior):.4f}; "
f"sessions improved {(x.behavior_plus_dynamics > x.behavior).sum()} of {len(x)}")
monkey session behavior behavior_plus_dynamics
G 0 0.3399 0.3711
G 1 0.4086 0.4346
G 2 0.1271 0.1450
G 3 0.1238 0.1470
G 4 0.1009 0.1186
G 5 0.1238 0.1393
G 6 0.1686 0.1907
G 7 0.1052 0.1177
G 8 0.1023 0.1128
J 0 0.3168 0.3485
J 1 0.3017 0.3324
J 2 0.3148 0.3384
J 3 0.2097 0.2225
Monkey G: paired t-test over 9 sessions, T = -8.852, p = 2.09e-05; mean gain 0.0196; sessions improved 9 of 9
Monkey J: paired t-test over 4 sessions, T = -5.634, p = 0.0111; mean gain 0.0247; sessions improved 4 of 4
Result: reproduced. Adding the previous activity raises the held-out R2 in 9 of 9 Monkey G sessions (paired t-test T = -8.852, p = 2.09e-05) and 4 of 4 Monkey J sessions (T = -5.634, p = 0.0111).
figs3f-residual-early-late¶
The residual of the dynamics model is larger early in the trial than late, consistent with large input at movement start.
In [30]:
rows = []
for k in KEYS:
R = RR[k]
T = BB[k]["trials"]
T = T[T.cond != "dropped"].set_index("trl")
med = int(np.median(T.nbins))
lo, hi = med / 3.0, 2 * med / 3.0
for tr in np.unique(R["resid_trl"]):
if T.loc[tr, "nbins"] < med:
continue
m = R["resid_trl"] == tr
b, r = R["resid_bin"][m], R["resid"][m]
early, late = r[b < lo], r[(b >= hi) & (b < med)]
if len(early) and len(late):
rows.append(dict(monkey=k[0], session=f"{k[0]}{k[1]}", epoch=0, resid=early.mean()))
rows.append(dict(monkey=k[0], session=f"{k[0]}{k[1]}", epoch=1, resid=late.mean()))
DS3F = pd.DataFrame(rows)
print(DS3F.groupby(["session", "epoch"]).resid.agg(["mean", "size"]).unstack().round(4).to_string())
for mk in ["G", "J"]:
x = DS3F[DS3F.monkey == mk]
fit = smf.mixedlm("resid ~ epoch", x, groups=x["session"]).fit()
print(f"Monkey {mk}: LME late minus early slope {fit.params['epoch']:.4f}, z = {fit.tvalues['epoch']:.2f}, p = {fit.pvalues['epoch']:.3g}, N = {len(x)}")
mean size epoch 0 1 0 1 session G0 1.1435 0.9855 157 157 G1 0.7224 0.6169 98 98 G2 0.7054 0.6449 277 277 G3 0.8852 0.7936 240 240 G4 0.7596 0.7184 176 176 G5 0.9607 0.9086 204 204 G6 0.7135 0.6569 81 81 G7 0.8980 0.8354 84 84 G8 0.8523 0.8144 70 70 J0 0.7440 0.7078 183 183 J1 0.7623 0.7656 88 88 J2 0.7110 0.6788 114 114 J3 0.5135 0.4869 89 89 Monkey G: LME late minus early slope -0.0752, z = -24.55, p = 3.97e-133, N = 2774 Monkey J: LME late minus early slope -0.0261, z = -4.63, p = 3.65e-06, N = 948
Result: reproduced. The residual is smaller late in the trial than early: LME slope -0.0752 (z = -24.55, p = 3.97e-133, N = 2774 trial epochs) for Monkey G and -0.0261 (z = -4.63, p = 3.65e-06, N = 948) for Monkey J.
figs6cd-pairwise-distances-predicted, figs6g-distance-tracks-behavior, figs6k-next-command-pairs¶
Across pairs of conditions issuing the same command, the dynamics-predicted distance between average activity correlates with the true distance, for single neurons and for the population.
Across pairs of conditions, the distance between activity patterns issuing the same command correlates with the distance between the command subtrajectories, and the dynamics model's predicted distances keep this relationship.
Across pairs of conditions, the dynamics-predicted distance between the next commands that follow the same command correlates with the true distance.
In [31]:
# Main reading: pairwise command matching as in the code (plot_pred_fr_diffs.py 2016-2050), the
# predicted distance regressed on the true one as the code does (linregress(true, pred), lines
# 1632-1751). Variant: each condition's mean over its own bins, without the matching.
rows = []
for reading, sfx in [("matched", "_m"), ("unmatched", "")]:
for mk, ks in [("G", G), ("J", J)]:
acc = {"neur": ([], []), "pop": ([], [])}
for k in ks:
P = [p for p in RR[k]["pairs"] if sfx == "" or np.isfinite(p["pop_true_m"])]
nt = np.concatenate([p["neur_true" + sfx] for p in P]); npd = np.concatenate([p["neur_pred" + sfx] for p in P])
pt = np.array([p["pop_true" + sfx] for p in P]); pp = np.array([p["pop_pred" + sfx] for p in P])
for nm, (a, b) in [("neur", (nt, npd)), ("pop", (pt, pp))]:
lr = scipy.stats.linregress(a, b)
acc[nm][0].append(a); acc[nm][1].append(b)
rows.append(dict(reading=reading, monkey=mk, session=k[1], measure=nm, n=len(a), slope=lr.slope, r=lr.rvalue, p=lr.pvalue))
for nm in ["neur", "pop"]:
a, b = np.concatenate(acc[nm][0]), np.concatenate(acc[nm][1])
lr = scipy.stats.linregress(a, b)
rows.append(dict(reading=reading, monkey=mk, session="pooled", measure=nm, n=len(a), slope=lr.slope, r=lr.rvalue, p=lr.pvalue))
npairs = sum(len(RR[k]["pairs"]) for k in KEYS)
nkept = sum(np.isfinite(p["pop_true_m"]) for k in KEYS for p in RR[k]["pairs"])
print(f"Pairs of conditions: {npairs}; kept by the pairwise matching (15 or more of each after it): {nkept}")
D6 = pd.DataFrame(rows)
print(D6.to_string(index=False))
print("\nExpected pooled r: neurons G 0.584, J 0.693; population G 0.638, J 0.800 (10%); paper slopes 0.221, 0.317, 0.350, 0.449; "
"paper N: neurons 377489 (G), 34643 (J); population 5112, 1714")
print("Regression of predicted on true distance (scipy linregress), pooled:")
print(D6[D6.session == "pooled"].to_string(index=False))
print("Sessions with p < 0.001 and positive slope:")
print(D6[D6.session != "pooled"].assign(ok=lambda x: (x.p < 0.001) & (x.slope > 0)).groupby(["reading", "monkey", "measure"]).ok.agg(["sum", "size"]).to_string())
Pairs of conditions: 7925; kept by the pairwise matching (15 or more of each after it): 7310 reading monkey session measure n slope r p matched G 0 neur 32648 0.316888 0.633623 0.000000e+00 matched G 0 pop 742 0.425443 0.659142 1.139262e-93 matched G 1 neur 13930 0.334220 0.665359 0.000000e+00 matched G 1 pop 398 0.344100 0.633555 4.717626e-46 matched G 2 neur 102860 0.232353 0.624660 0.000000e+00 matched G 2 pop 695 0.308714 0.729578 1.848121e-116 matched G 3 neur 61374 0.136778 0.464073 0.000000e+00 matched G 3 pop 1158 0.218651 0.586793 4.433900e-108 matched G 4 neur 112646 0.273871 0.669884 0.000000e+00 matched G 4 pop 746 0.444241 0.757251 8.774358e-140 matched G 5 neur 43764 0.088636 0.361801 0.000000e+00 matched G 5 pop 1042 0.198235 0.577567 9.516257e-94 matched G 6 neur 7150 0.151696 0.505818 0.000000e+00 matched G 6 pop 143 0.222311 0.614252 3.386521e-16 matched G 7 neur 13632 0.110025 0.381474 0.000000e+00 matched G 7 pop 284 0.229565 0.522394 2.751692e-21 matched G 8 neur 11520 0.098168 0.371507 0.000000e+00 matched G 8 pop 240 0.198327 0.539559 1.565791e-19 matched G pooled neur 399524 0.224955 0.586003 0.000000e+00 matched G pooled pop 5448 0.356734 0.650386 0.000000e+00 matched J 0 neur 16520 0.306962 0.632661 0.000000e+00 matched J 0 pop 826 0.418152 0.678515 1.702018e-112 matched J 1 neur 3980 0.319764 0.694648 0.000000e+00 matched J 1 pop 199 0.397605 0.759622 1.158575e-38 matched J 2 neur 9000 0.307060 0.678261 0.000000e+00 matched J 2 pop 450 0.383364 0.751844 4.577614e-83 matched J 3 neur 8127 0.152116 0.523749 0.000000e+00 matched J 3 pop 387 0.227583 0.621094 1.173098e-42 matched J pooled neur 37627 0.300947 0.655128 0.000000e+00 matched J pooled pop 1862 0.433786 0.766134 0.000000e+00 unmatched G 0 neur 35596 0.338835 0.644628 0.000000e+00 unmatched G 0 pop 809 0.455537 0.669789 1.987355e-106 unmatched G 1 neur 16345 0.549590 0.777644 0.000000e+00 unmatched G 1 pop 467 0.737349 0.771876 1.641574e-93 unmatched G 2 neur 110556 0.256033 0.627395 0.000000e+00 unmatched G 2 pop 747 0.373585 0.737765 2.857490e-129 unmatched G 3 neur 65773 0.143023 0.467263 0.000000e+00 unmatched G 3 pop 1241 0.225935 0.592035 2.969851e-118 unmatched G 4 neur 116119 0.286653 0.677825 0.000000e+00 unmatched G 4 pop 769 0.466960 0.759024 4.240768e-145 unmatched G 5 neur 49392 0.096509 0.373079 0.000000e+00 unmatched G 5 pop 1176 0.218069 0.608100 7.623373e-120 unmatched G 6 neur 7900 0.160716 0.518077 0.000000e+00 unmatched G 6 pop 158 0.251256 0.687242 2.044918e-23 unmatched G 7 neur 15456 0.105273 0.361083 0.000000e+00 unmatched G 7 pop 322 0.233458 0.528714 1.380930e-24 unmatched G 8 neur 13008 0.104333 0.376254 0.000000e+00 unmatched G 8 pop 271 0.189338 0.507707 3.666721e-19 unmatched G pooled neur 430145 0.247824 0.594209 0.000000e+00 unmatched G pooled pop 5960 0.406979 0.645806 0.000000e+00 unmatched J 0 neur 16940 0.342209 0.655101 0.000000e+00 unmatched J 0 pop 847 0.464020 0.719407 6.930580e-136 unmatched J 1 neur 4280 0.342233 0.698647 0.000000e+00 unmatched J 1 pop 214 0.388444 0.740952 1.616260e-38 unmatched J 2 neur 9840 0.369205 0.715041 0.000000e+00 unmatched J 2 pop 492 0.459164 0.774951 1.059515e-99 unmatched J 3 neur 8652 0.173207 0.527723 0.000000e+00 unmatched J 3 pop 412 0.272251 0.645053 7.668150e-50 unmatched J pooled neur 39712 0.342816 0.680156 0.000000e+00 unmatched J pooled pop 1965 0.478859 0.789433 0.000000e+00 Expected pooled r: neurons G 0.584, J 0.693; population G 0.638, J 0.800 (10%); paper slopes 0.221, 0.317, 0.350, 0.449; paper N: neurons 377489 (G), 34643 (J); population 5112, 1714 Regression of predicted on true distance (scipy linregress), pooled: reading monkey session measure n slope r p matched G pooled neur 399524 0.224955 0.586003 0.0 matched G pooled pop 5448 0.356734 0.650386 0.0 matched J pooled neur 37627 0.300947 0.655128 0.0 matched J pooled pop 1862 0.433786 0.766134 0.0 unmatched G pooled neur 430145 0.247824 0.594209 0.0 unmatched G pooled pop 5960 0.406979 0.645806 0.0 unmatched J pooled neur 39712 0.342816 0.680156 0.0 unmatched J pooled pop 1965 0.478859 0.789433 0.0 Sessions with p < 0.001 and positive slope:
sum size
reading monkey measure
matched G neur 9 9
pop 9 9
J neur 4 4
pop 4 4
unmatched G neur 9 9
pop 9 9
J neur 4 4
pop 4 4
Result: reproduced. With the code's pairwise command matching (7310 of 7925 pairs of conditions keep 15 or more occurrences of each), the pooled r of predicted on true distance is 0.586 (G) and 0.655 (J) for single neurons and 0.650 and 0.766 for the population, within 10% of 0.584, 0.693, 0.638 and 0.800; the slopes, 0.225, 0.301, 0.357 and 0.434, are close to the paper's 0.221, 0.317, 0.350 and 0.449. Every session and both pooled regressions give a positive slope with p < 0.001. Without the matching r is 0.594, 0.680, 0.646 and 0.789, with every session at p < 0.001.
In [32]:
rows = []
for k in KEYS:
for p in RR[k]["pairs"]:
if p["n_conds"] > 5 and np.isfinite(p["beh"]):
rows.append(dict(monkey=k[0], session=f"{k[0]}{k[1]}", command=f"{k[0]}{k[1]}_{p['command']}", beh=p["beh"], true=p["pop_true"], pred=p["pop_pred"]))
D6G = pd.DataFrame(rows)
for mk in ["G", "J"]:
x = D6G[D6G.monkey == mk].reset_index(drop=True)
for y in ["true", "pred"]:
md = smf.mixedlm(f"{y} ~ beh", x, groups=x["session"], re_formula="1", vc_formula={"command": "0 + C(command)"}).fit()
print(f"Monkey {mk}, {y} distance: N = {len(x)}, slope {md.params['beh']:.4f}, z = {md.tvalues['beh']:.2f}, p = {md.pvalues['beh']:.3g}")
per = []
for s, xs in x.groupby("session"):
if xs.command.nunique() > 1:
f1 = smf.mixedlm("true ~ beh", xs, groups=xs["command"]).fit()
per.append((s, len(xs), round(f1.params["beh"], 4), f1.pvalues["beh"]))
print(" per session (true distances, command as random effect):", per)
print("Paper: G N = 4674, z = 31.04 (true) and 33.1 (predicted); J N = 1629, z = 14.18 and 22.5")
Monkey G, true distance: N = 5509, slope 0.2512, z = 38.69, p = 0 Monkey G, pred distance: N = 5509, slope 0.1573, z = 42.27, p = 0
per session (true distances, command as random effect): [('G0', 784, np.float64(0.3923), np.float64(1.9069637496534024e-75)), ('G1', 383, np.float64(0.369), np.float64(2.3076200243627976e-55)), ('G2', 705, np.float64(0.149), np.float64(3.773426384055983e-112)), ('G3', 1238, np.float64(0.4469), np.float64(1.8452352454633536e-68)), ('G4', 698, np.float64(0.5895), np.float64(4.319897982584192e-57)), ('G5', 1169, np.float64(0.2133), np.float64(5.1543728691741025e-53)), ('G6', 102, np.float64(0.3135), np.float64(1.5280220790209052e-05)), ('G7', 238, np.float64(0.4711), np.float64(1.4932269013892397e-10)), ('G8', 192, np.float64(0.2462), np.float64(6.758843684361057e-20))]
Monkey J, true distance: N = 1844, slope 0.1292, z = 16.66, p = 2.69e-62
Monkey J, pred distance: N = 1844, slope 0.1092, z = 24.37, p = 3.72e-131
per session (true distances, command as random effect): [('J0', 814, np.float64(0.1089), np.float64(2.510989614523188e-22)), ('J1', 177, np.float64(0.0542), np.float64(0.020740971078605095)), ('J2', 482, np.float64(0.1803), np.float64(1.267364469572404e-29)), ('J3', 371, np.float64(0.1834), np.float64(5.814333749616757e-26))]
Paper: G N = 4674, z = 31.04 (true) and 33.1 (predicted); J N = 1629, z = 14.18 and 22.5
Result: reproduced. Over commands used in more than five conditions, population distance rises with command subtrajectory distance pooled over sessions, for true distances (Monkey G slope 0.251, z = 38.69; Monkey J 0.129, z = 16.66) and for dynamics-predicted distances (z = 42.27 and 24.37), all with p < 0.001.
Result: reproduced. Pooled r is 0.774 (G) and 0.654 (J), within 10% of 0.76 and 0.63, with positive slopes and p < 0.001; all 9 Monkey G sessions (session 7 at p = 0.046) and 4 of 4 Monkey J sessions are significant at p < 0.05.
monkey-j-linear-command¶
For Monkey J, a time-invariant linear map of neural activity accounts for 99.6, 99.6, 99.5 and 99.8% of the variance of the point process filter's command in the four sessions.
In [34]:
for k in J:
dec = DEC[k]
print(f"Monkey J session {k[1]}: variance of the velocity update accounted for {100 * dec['vaf_command']:.2f}% "
f"(velocity itself {100 * dec['vaf_velocity']:.3f}%); per task entry {[round(100 * v['vaf_command'], 2) for v in dec['per_te'].values()]}")
print("Paper: 99.6, 99.6, 99.5, 99.8%; criterion: every session at least 99% and within 0.5 points of the reported value")
# Readings of the fit (the criterion: rewarded trials, no intercept, the PPF's own update, which the files do not hold)
rows = []
for k in J:
for r in bmi.j_command_fit_variants(RAW[k]):
rows.append(dict(session=k[1], **r))
JV = pd.DataFrame(rows)
JV["vaf_percent"] = 100 * JV.vaf
print(JV.pivot_table(index=["frames", "dynamics"], columns="session", values="vaf_percent").round(2).to_string())
print("frames per session:", JV.groupby(["session", "frames"]).n.first().to_dict())
Monkey J session 0: variance of the velocity update accounted for 98.92% (velocity itself 99.985%); per task entry [np.float64(99.17), np.float64(98.8)] Monkey J session 1: variance of the velocity update accounted for 98.89% (velocity itself 99.985%); per task entry [np.float64(98.94), np.float64(98.41)] Monkey J session 2: variance of the velocity update accounted for 99.57% (velocity itself 99.995%); per task entry [np.float64(99.5), np.float64(99.59)] Monkey J session 3: variance of the velocity update accounted for 99.25% (velocity itself 99.992%); per task entry [np.float64(99.63), np.float64(99.16)] Paper: 99.6, 99.6, 99.5, 99.8%; criterion: every session at least 99% and within 0.5 points of the reported value
session 0 1 2 3
frames dynamics
all fitted, no offset 98.87 98.84 98.90 99.15
fitted, with offset 98.92 98.89 99.57 99.25
identity on velocity 87.05 87.10 87.45 92.43
rewarded trials fitted, no offset 99.10 99.09 99.09 99.63
fitted, with offset 99.63 99.63 99.57 99.80
identity on velocity 87.33 87.35 87.66 93.89
frames per session: {(0, 'all'): 851999, (0, 'rewarded trials'): 181137, (1, 'all'): 544996, (1, 'rewarded trials'): 105940, (2, 'all'): 311999, (2, 'rewarded trials'): 111354, (3, 'all'): 431999, (3, 'rewarded trials'): 121203}
Result: blocked, data missing. The criterion needs the point process filter's own velocity update or its parameters and says the claim is blocked when the dataset holds neither, as here. The readings shown estimate the update with cursor dynamics fitted jointly with the map, which favours the fit, so none is the criterion's test. Over all frames, with the offset column in the fitted dynamics, the map accounts for 98.92, 98.89, 99.57 and 99.25% of the variance of the update; over the frames of rewarded trials, as the criterion says, 99.63, 99.63, 99.57 and 99.80% with the offset column and 99.10, 99.09, 99.09 and 99.63% without it (the paper gives 99.6, 99.6, 99.5 and 99.8%). With the update taken as the change in velocity it is 87.05 to 93.89%.
Summary¶
| Claim | Outcome |
|---|---|
| conditions-per-session | reproduced |
| figs1b-conditions-per-command | reproduced |
| figs1d-command-subtrajectories-differ | reproduced |
| figs1e-next-commands-differ | partially reproduced |
| fig3d-population-distance | reproduced |
| fig3b-single-neuron-distance | reproduced |
| fig3e-fractions-significant | partially reproduced |
| fig4c-dynamics-predict-activity-given-command | reproduced |
| fig4de-condition-specific-activity-predicted | reproduced |
| fig4f-condition-specific-component | partially reproduced |
| fig4g-majority-predicted | reproduced |
| figs3-low-dimensional-decaying-dynamics | not reproduced |
| fig5c-next-activity | reproduced |
| fig5d-four-dimensions | reproduced |
| fig5e-next-command | partially reproduced |
| fig5f-condition-specific-next-command | reproduced |
| fig5h-turn-predicted | partially reproduced |
| figs5b-dynamics-beyond-behavior | reproduced |
| figs5f-nonlinear-no-better | blocked (not runnable) |
| figs3f-residual-early-late | reproduced |
| figs4b-generalize-command-sets | reproduced |
| figs4g-generalize-condition-classes | reproduced |
| figs6cd-pairwise-distances-predicted | reproduced |
| figs6g-distance-tracks-behavior | reproduced |
| figs6k-next-command-pairs | reproduced |
| fig6c-dynamics-reduce-input | reproduced, in ofc_simulations.ipynb (assumption) |
| fig6d-decoder-null-dynamics-no-saving | reproduced, in ofc_simulations.ipynb (assumption) |
| fig6g-ofc-uses-distinct-patterns | partially reproduced, in ofc_simulations.ipynb (assumption) |
| monkey-j-linear-command | 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 decoder matrices are not in the files; the gain is recovered by least squares from the decoder state and spike counts (exact for Monkey G, the paper's own linear approximation for Monkey J). Targets come from the cursor at reward, not the target series. The R2 of the dynamics models is scored as the paper's code scores it, on each model's test pairs whose commands at t and t-1 are both below the 95th percentile of magnitude, with the full model and the shuffle models on the same pairs; two other readings are shown. Folds are random over bins with a fixed seed; the behavior-preserving shuffle uses 1000 datasets with the real data's ridge parameter (the parameter chosen on one shuffled dataset is shown on 200 shuffles), and, as in the code, each shuffle permutes the whole session before the folds are split. The rsLDS comparison was not run: the package cannot be installed in the entry's locked environment. Monkey J's linear command is blocked as its criterion directs, the files holding neither the filter's update nor its parameters; the estimated readings are shown. The optimal feedback control simulations of Figure 6 are in ofc_simulations.ipynb.
What was found about the dataset. The files hold no decoder parameters, though the paper's Data availability lists them; the gain can be recovered from the decoder state exactly for Monkey G. Monkey J's target and obstacle series carry one value for every trial of a session (a bug in the code that wrote them), and Monkey G's target series does not hold the trial's target; targets have to be read from the cursor at reward. Units differ between monkeys (cm against m) and Monkey J's target and obstacle series are in cm while its cursor is in m. The trials table has no task, target or outcome columns.
In [35]:
print(f"notebook ran in {(time.time() - started) / 60:.1f} min")
notebook ran in 10.4 min