Optimal feedback control simulations of Athalye et al. (2023), Figure 6, on 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. Claims file written from the accepted author manuscript (PMC10527529) with its supplement.
Dandiset. DANDI:000404, version 0.230605.2024, the version the paper cites and the only published one. Loading, the decoder and the dynamics fits are those of reproduce_athalye_2023.ipynb (helper bmi.py); the simulation is in ofc.py.
The simulations need, per session, the fitted dynamics (A, b), the decoder's gain and its cursor dynamics at the 100 ms step. The files hold no decoder matrices; Monkey G's cursor update is recovered exactly from the decoder state, and Monkey J's 100 ms update is estimated by least squares from the decoder state over each 100 ms bin (an assumption: the paper builds it from the PPF's parameters, which are not in the files).
| Claim | Tested here |
|---|---|
| fig6c-dynamics-reduce-input | tested |
| fig6d-decoder-null-dynamics-no-saving | tested |
| fig6g-ofc-uses-distinct-patterns | tested |
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 scipy.stats
import statsmodels.formula.api as smf
import bmi
import ofc
warnings.filterwarnings("ignore")
pd.set_option("display.width", 200)
pd.set_option("display.max_rows", 300)
NSHUFF = 1000
started = time.time()
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)
return (m, s), d, bmi.prepare(B, dec, m)
with ThreadPoolExecutor(4) as ex:
loaded = list(ex.map(load, sessions))
RAW = {k: d for k, d, _ 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 50 s
Data check¶
rows = []
for k in KEYS:
F, K, vaf = ofc.cursor_update(RAW[k], BB[k], k[0])
rows.append(dict(monkey=k[0], session=k[1], units=BB[k]["spk"].shape[1], bins=len(BB[k]["spk"]),
vaf_pos_x=vaf[0], vaf_pos_y=vaf[1], vaf_vel_x=vaf[2], vaf_vel_y=vaf[3],
velocity_decay=round(F[2, 2], 4), pos_from_vel=round(F[0, 2], 4),
max_rel_diff_Kvel_vs_command_matrix=np.abs(K[2:4] - BB[k]["Kc"]).max() / np.abs(BB[k]["Kc"]).max()))
print(pd.DataFrame(rows).round(5).to_string(index=False))
monkey session units bins vaf_pos_x vaf_pos_y vaf_vel_x vaf_vel_y velocity_decay pos_from_vel max_rel_diff_Kvel_vs_command_matrix
G 0 44 11971 1.00000 0.99978 0.99983 0.99973 0.5000 0.0700 0.00000
G 1 35 7252 0.99920 0.99982 0.99997 1.00000 0.5000 0.0700 0.00001
G 2 148 16505 0.99961 0.99955 0.98196 0.96672 0.5000 0.0700 0.00000
G 3 53 15870 0.99983 0.99969 0.99557 0.99923 0.7764 0.0918 0.00001
G 4 151 12383 0.99916 0.99984 0.98409 0.99634 0.7750 0.0914 0.00001
G 5 42 11696 0.99986 0.99949 0.99937 0.99942 0.7074 0.0776 0.00001
G 6 50 4536 0.99925 0.99910 0.99888 0.99853 0.7375 0.0825 0.00001
G 7 48 5692 0.99999 1.00000 0.99986 0.99999 0.7854 0.0947 0.00002
G 8 48 5412 0.99978 0.99974 0.99833 0.99932 0.5000 0.0700 0.00000
J 0 20 8848 0.99985 0.99990 0.99511 0.99685 0.6065 0.0803 0.26361
J 1 20 5110 0.99987 0.99991 0.99531 0.99690 0.6130 0.0807 0.23435
J 2 20 5525 0.99987 0.99990 0.99558 0.99745 0.6190 0.0811 0.19161
J 3 21 6015 0.99991 0.99991 0.99781 0.99776 0.7019 0.0855 0.14412
Monkey G's 100 ms cursor update is the Kalman filter's (variance accounted for 0.967 to 1 on the updates; its velocity rows equal the command matrix). Monkey J's 100 ms update, estimated over 100 ms bins, accounts for over 99% of the position and velocity; its velocity gain differs from the sum of the 5 ms gains because the cursor dynamics act within the bin.
def job(k):
t0 = time.time()
return k, ofc.run_session(RAW[k], BB[k], k[0], 2000 * (k[0] == "J") + k[1], NSHUFF), time.time() - t0
order = sorted(KEYS, key=lambda k: -BB[k]["spk"].shape[1])
SIM = {}
with ProcessPoolExecutor(4) as ex:
for k, o, dt in ex.map(job, order):
SIM[k] = o
print(k, f"{dt:.0f} s", flush=True)
rows = []
for k in KEYS:
for name, r in SIM[k]["models"].items():
rows.append(dict(monkey=k[0], session=k[1], model=name, noise_sd=SIM[k]["sigmas"][name], tuples=r["n_tuples"],
final_error_cm_median=np.median(r["final_err"]), final_error_cm_max=r["final_err"].max()))
print(pd.DataFrame(rows).round(4).to_string(index=False))
print("Largest final distance from the target over all trials (cm):", round(max(r["final_err"].max() for k in KEYS for r in SIM[k]["models"].values()), 3),
"(target radius 1.7 cm for Monkey G, 1.2 cm for Monkey J)")
('G', 4) 135 s
('G', 2) 122 s
('G', 3) 70 s
('G', 6) 65 s
('G', 7) 65 s
('G', 8) 55 s
('G', 0) 56 s
('G', 5) 58 s
('G', 1) 52 s
('J', 3) 50 s
('J', 0) 42 s
('J', 1) 31 s
('J', 2) 41 s
monkey session model noise_sd tuples final_error_cm_median final_error_cm_max
G 0 Full Dynamics 0.2613 301 0.0625 0.1963
G 0 No Dynamics 0.2613 237 0.0866 0.2027
G 0 Decoder-null Dynamics 0.2613 239 0.0835 0.2566
G 0 No Dynamics (decoder-null offset) 0.2613 237 0.0804 0.2354
G 1 Full Dynamics 0.3377 299 0.0773 0.3420
G 1 No Dynamics 0.3377 212 0.1293 0.3649
G 1 Decoder-null Dynamics 0.3377 219 0.1204 0.3136
G 1 No Dynamics (decoder-null offset) 0.3377 211 0.1217 0.3316
G 2 Full Dynamics 0.0210 141 0.0263 0.0987
G 2 No Dynamics 0.0210 138 0.0273 0.0848
G 2 Decoder-null Dynamics 0.0210 139 0.0287 0.0791
G 2 No Dynamics (decoder-null offset) 0.0210 137 0.0260 0.0818
G 3 Full Dynamics 0.1730 334 0.0640 0.1926
G 3 No Dynamics 0.1730 295 0.0666 0.1965
G 3 Decoder-null Dynamics 0.1730 294 0.0617 0.1947
G 3 No Dynamics (decoder-null offset) 0.1730 294 0.0647 0.1860
G 4 Full Dynamics 0.0226 221 0.0346 0.1038
G 4 No Dynamics 0.0226 222 0.0353 0.1023
G 4 Decoder-null Dynamics 0.0226 220 0.0369 0.1064
G 4 No Dynamics (decoder-null offset) 0.0226 224 0.0362 0.1139
G 5 Full Dynamics 0.1682 324 0.0632 0.1840
G 5 No Dynamics 0.1682 296 0.0646 0.1868
G 5 Decoder-null Dynamics 0.1682 283 0.0648 0.2214
G 5 No Dynamics (decoder-null offset) 0.1682 296 0.0610 0.2032
G 6 Full Dynamics 0.1183 305 0.0638 0.1906
G 6 No Dynamics 0.1183 263 0.0589 0.2074
G 6 Decoder-null Dynamics 0.1183 269 0.0598 0.1778
G 6 No Dynamics (decoder-null offset) 0.1183 261 0.0593 0.1763
G 7 Full Dynamics 0.2326 346 0.0746 0.2488
G 7 No Dynamics 0.2326 326 0.0707 0.2407
G 7 Decoder-null Dynamics 0.2326 324 0.0784 0.2396
G 7 No Dynamics (decoder-null offset) 0.2326 332 0.0759 0.2212
G 8 Full Dynamics 0.1602 247 0.0532 0.1872
G 8 No Dynamics 0.1602 189 0.0517 0.1448
G 8 Decoder-null Dynamics 0.1602 188 0.0562 0.1649
G 8 No Dynamics (decoder-null offset) 0.1602 187 0.0551 0.1651
J 0 Full Dynamics 0.2103 184 0.0619 0.1802
J 0 No Dynamics 0.2103 145 0.0581 0.1890
J 0 Decoder-null Dynamics 0.2103 146 0.0667 0.1832
J 0 No Dynamics (decoder-null offset) 0.2103 143 0.0648 0.1936
J 1 Full Dynamics 0.2202 181 0.0568 0.1735
J 1 No Dynamics 0.2202 150 0.0639 0.2198
J 1 Decoder-null Dynamics 0.2203 153 0.0724 0.1857
J 1 No Dynamics (decoder-null offset) 0.2203 149 0.0715 0.2027
J 2 Full Dynamics 0.2028 182 0.0550 0.1902
J 2 No Dynamics 0.2028 158 0.0544 0.1925
J 2 Decoder-null Dynamics 0.2028 156 0.0592 0.1824
J 2 No Dynamics (decoder-null offset) 0.2028 157 0.0608 0.1650
J 3 Full Dynamics 0.1342 224 0.0434 0.1282
J 3 No Dynamics 0.1342 201 0.0401 0.1542
J 3 Decoder-null Dynamics 0.1342 200 0.0416 0.1149
J 3 No Dynamics (decoder-null offset) 0.1342 202 0.0442 0.1404
Largest final distance from the target over all trials (cm): 0.365 (target radius 1.7 cm for Monkey G, 1.2 cm for Monkey J)
Every simulated trial ends within the printed distance of its target, less than a target radius, so every simulated trial completes its task.
fig6c-dynamics-reduce-input and fig6d-decoder-null-dynamics-no-saving¶
In an optimal feedback control model that drives the recorded decoders, a neural population with the fitted invariant dynamics needs significantly less input to perform the tasks than one with no dynamics.
When the model's invariant dynamics are restricted to the decoder-null space, the input needed is not different from a population with no dynamics.
def input_table(model, base):
rows, lme = [], []
for k in KEYS:
a = pd.DataFrame(SIM[k]["models"][model]["mags"]).groupby(["task", "target"]).mag.mean()
c = pd.DataFrame(SIM[k]["models"][base]["mags"]).groupby(["task", "target"]).mag.mean()
w = scipy.stats.wilcoxon(a.values, c.values)
rows.append(dict(monkey=k[0], session=k[1], input_model=a.mean(), input_no_dynamics=c.mean(), ratio=a.mean() / c.mean(),
conditions_lower=int((a.values < c.values).sum()), wilcoxon_p=w.pvalue))
for v in a.values:
lme.append(dict(monkey=k[0], session=f"{k[0]}{k[1]}", model=1, mag=v))
for v in c.values:
lme.append(dict(monkey=k[0], session=f"{k[0]}{k[1]}", model=0, mag=v))
df = pd.DataFrame(rows)
L = pd.DataFrame(lme)
print(f"{model} against {base}:")
print(df.round(4).to_string(index=False))
for mk in ["G", "J"]:
x = L[L.monkey == mk]
fit = smf.mixedlm("mag ~ model", x, groups=x["session"]).fit()
print(f" Monkey {mk}: LME effect of the dynamics model on input {fit.params['model']:.4f}, z = {fit.tvalues['model']:.2f}, p = {fit.pvalues['model']:.3g}, N = {len(x)}; "
f"sessions with Wilcoxon p < 0.05: {(df[df.monkey == mk].wilcoxon_p < 0.05).sum()} of {(df.monkey == mk).sum()}")
return df
C6 = input_table("Full Dynamics", "No Dynamics")
print("Paper: Monkey G z = 10.49, N = 432; Monkey J z = 5.20, N = 192; Wilcoxon p < 0.05 in 9 of 9 and 4 of 4 sessions\n")
D6 = input_table("Decoder-null Dynamics", "No Dynamics (decoder-null offset)")
print("Paper: Monkey G z = 0.002, p = 0.998; Monkey J z = -0.003, p = 0.990; no session significant")
Full Dynamics against No Dynamics:
monkey session input_model input_no_dynamics ratio conditions_lower wilcoxon_p
G 0 6.9567 11.0494 0.6296 23 0.0000
G 1 6.2980 13.4545 0.4681 23 0.0000
G 2 1.4886 1.9522 0.7625 18 0.0002
G 3 4.8973 7.5439 0.6492 24 0.0000
G 4 1.5970 1.7647 0.9050 18 0.0025
G 5 3.9824 7.0781 0.5626 23 0.0000
G 6 3.5832 5.1222 0.6995 22 0.0000
G 7 6.1644 9.4908 0.6495 24 0.0000
G 8 4.4683 6.9326 0.6445 24 0.0000
J 0 3.4846 7.4997 0.4646 24 0.0000
J 1 3.7585 8.5066 0.4418 24 0.0000
J 2 3.2724 7.2974 0.4484 24 0.0000
J 3 2.6735 4.3449 0.6153 23 0.0000
Monkey G: LME effect of the dynamics model on input -2.7725, z = -11.27, p = 1.87e-29, N = 432; sessions with Wilcoxon p < 0.05: 9 of 9
Monkey J: LME effect of the dynamics model on input -3.6149, z = -14.94, p = 1.87e-50, N = 192; sessions with Wilcoxon p < 0.05: 4 of 4
Paper: Monkey G z = 10.49, N = 432; Monkey J z = 5.20, N = 192; Wilcoxon p < 0.05 in 9 of 9 and 4 of 4 sessions
Decoder-null Dynamics against No Dynamics (decoder-null offset):
monkey session input_model input_no_dynamics ratio conditions_lower wilcoxon_p
G 0 10.7000 10.6894 1.0010 12 0.7469
G 1 13.5668 13.5746 0.9994 13 0.6033
G 2 1.9418 1.9410 1.0004 11 0.7257
G 3 7.6219 7.6325 0.9986 12 0.3902
G 4 1.9007 1.9000 1.0004 11 0.6840
G 5 6.6492 6.6518 0.9996 16 0.5838
G 6 5.2373 5.2370 1.0001 10 0.8774
G 7 9.4891 9.4742 1.0016 10 0.3596
G 8 7.0262 7.0371 0.9985 15 0.1434
J 0 12.6770 12.6572 1.0016 11 0.2292
J 1 13.4485 13.4750 0.9980 14 0.1140
J 2 11.5580 11.5491 1.0008 11 0.4732
J 3 5.8232 5.8058 1.0030 9 0.0691
Monkey G: LME effect of the dynamics model on input -0.0005, z = -0.00, p = 0.998, N = 432; sessions with Wilcoxon p < 0.05: 0 of 9
Monkey J: LME effect of the dynamics model on input 0.0049, z = 0.02, p = 0.988, N = 192; sessions with Wilcoxon p < 0.05: 0 of 4
Paper: Monkey G z = 0.002, p = 0.998; Monkey J z = -0.003, p = 0.990; no session significant
Result: fig6c-dynamics-reduce-input reproduced; fig6d-decoder-null-dynamics-no-saving reproduced; both resting on an assumption for Monkey J. With the fitted dynamics the population needs less input than without dynamics in every session (ratios 0.44 to 0.91; Wilcoxon p < 0.05 in 9 of 9 Monkey G and 4 of 4 Monkey J sessions), and the LME finds the reduction in each monkey (z = -11.27, p = 1.87e-29, N = 432; z = -14.94, p = 1.87e-50, N = 192). Decoder-null dynamics need the same input as their No Dynamics model (ratios 0.998 to 1.003): the LME finds no difference (p = 0.998 and 0.988) and no session's Wilcoxon test does (smallest p 0.069). Monkey J's simulations use a 100 ms decoder update estimated from the files. With Monkey J's velocity gain fixed to the summed 5 ms gain (sensitivity reading, below), the Full Dynamics effect holds (LME z = -16.61, Wilcoxon p < 0.05 in 4 of 4 sessions) and the decoder-null LME finds no difference (p = 0.995), but Monkey J session 1's Wilcoxon test does (p = 0.0016), so under that reading fig6d-decoder-null-dynamics-no-saving's criterion, which fails on any session's Wilcoxon p < 0.05, would give not reproduced.
fig6g-ofc-uses-distinct-patterns¶
In the OFC simulations, activity issuing the same command differs between conditions more than shuffle with the fitted full dynamics, but not without dynamics or with decoder-null dynamics.
for null in ["shuffle", "subsample"]:
rows, lme = [], []
for k in KEYS:
for name, r in SIM[k]["models"].items():
if r["n_tuples"] == 0:
continue
obs, nm = r["dist_" + null], r["null_mean_" + null]
u = scipy.stats.mannwhitneyu(obs, nm)
rows.append(dict(model=name, monkey=k[0], session=k[1], tuples=len(obs), mean_norm=(obs / nm).mean(), mw_p=u.pvalue))
for v in obs:
lme.append(dict(model=name, monkey=k[0], session=f"{k[0]}{k[1]}", source=1, dist=v))
for v in nm:
lme.append(dict(model=name, monkey=k[0], session=f"{k[0]}{k[1]}", source=0, dist=v))
df = pd.DataFrame(rows)
L = pd.DataFrame(lme)
print(f"Null: {null} ({NSHUFF} draws)")
print(df.round(4).to_string(index=False))
out = []
for name in df.model.unique():
for mk in ["G", "J"]:
x = L[(L.model == name) & (L.monkey == mk)]
fit = smf.mixedlm("dist ~ source", x, groups=x["session"]).fit()
sub = df[(df.model == name) & (df.monkey == mk)]
out.append(dict(model=name, monkey=mk, N=len(x), lme_effect=fit.params["source"], z=fit.tvalues["source"], p=fit.pvalues["source"],
sessions_mw_p05=int((sub.mw_p < 0.05).sum()), sessions=len(sub), mean_norm=sub.mean_norm.mean()))
print(pd.DataFrame(out).to_string(index=False))
print()
print("Paper (LME, simulation against shuffle): Full G z = -23.09, J z = -16.68; No Dynamics G z = 0.168, J 0.462; Decoder-null G 0.932, J -1.490; "
"No Dynamics (decoder-null) G 0.611, J 0.449; Mann-Whitney p < 0.05 in 9 of 9 and 4 of 4 sessions for Full Dynamics only")
Null: shuffle (1000 draws)
model monkey session tuples mean_norm mw_p
Full Dynamics G 0 301 2.2141 0.0000
No Dynamics G 0 237 1.0106 0.9626
Decoder-null Dynamics G 0 239 1.0491 0.3860
No Dynamics (decoder-null offset) G 0 237 1.0017 0.8795
Full Dynamics G 1 299 1.5894 0.0000
No Dynamics G 1 212 0.9959 0.8346
Decoder-null Dynamics G 1 219 1.0439 0.4954
No Dynamics (decoder-null offset) G 1 211 1.0008 0.9231
Full Dynamics G 2 141 1.5243 0.0001
No Dynamics G 2 138 1.0069 0.9441
Decoder-null Dynamics G 2 139 1.0136 0.8767
No Dynamics (decoder-null offset) G 2 137 0.9978 0.9295
Full Dynamics G 3 334 1.4170 0.0000
No Dynamics G 3 295 0.9888 0.6474
Decoder-null Dynamics G 3 294 1.0218 0.7192
No Dynamics (decoder-null offset) G 3 294 0.9940 0.7874
Full Dynamics G 4 221 1.4028 0.0000
No Dynamics G 4 222 1.0047 0.9354
Decoder-null Dynamics G 4 220 1.0122 0.8416
No Dynamics (decoder-null offset) G 4 224 1.0009 0.9834
Full Dynamics G 5 324 1.1751 0.0000
No Dynamics G 5 296 0.9998 0.9152
Decoder-null Dynamics G 5 283 1.0181 0.8387
No Dynamics (decoder-null offset) G 5 296 1.0066 0.8912
Full Dynamics G 6 305 1.4045 0.0000
No Dynamics G 6 263 1.0072 0.8769
Decoder-null Dynamics G 6 269 1.0153 0.8050
No Dynamics (decoder-null offset) G 6 261 1.0035 0.9362
Full Dynamics G 7 346 1.3255 0.0000
No Dynamics G 7 326 0.9994 0.8550
Decoder-null Dynamics G 7 324 1.0083 0.9948
No Dynamics (decoder-null offset) G 7 332 1.0031 0.9708
Full Dynamics G 8 247 1.3948 0.0000
No Dynamics G 8 189 1.0058 0.9250
Decoder-null Dynamics G 8 188 0.9997 0.8782
No Dynamics (decoder-null offset) G 8 187 0.9968 0.8889
Full Dynamics J 0 184 1.6181 0.0000
No Dynamics J 0 145 1.0070 0.8654
Decoder-null Dynamics J 0 146 1.1123 0.1621
No Dynamics (decoder-null offset) J 0 143 0.9662 0.5227
Full Dynamics J 1 181 1.5803 0.0000
No Dynamics J 1 150 0.9941 0.7600
Decoder-null Dynamics J 1 153 1.1331 0.0930
No Dynamics (decoder-null offset) J 1 149 1.0226 0.8015
Full Dynamics J 2 182 1.7121 0.0000
No Dynamics J 2 158 0.9827 0.6904
Decoder-null Dynamics J 2 156 1.0818 0.3745
No Dynamics (decoder-null offset) J 2 157 1.0195 0.7969
Full Dynamics J 3 224 1.2970 0.0000
No Dynamics J 3 201 1.0116 0.8502
Decoder-null Dynamics J 3 200 1.0320 0.7816
No Dynamics (decoder-null offset) J 3 202 1.0021 0.8184
model monkey N lme_effect z p sessions_mw_p05 sessions mean_norm
Full Dynamics G 5036 0.154909 25.445118 7.994965e-143 9 9 1.494180
Full Dynamics J 1542 0.133990 15.511577 2.896782e-54 4 4 1.551894
No Dynamics G 4356 0.000003 0.000923 9.992639e-01 0 9 1.002103
No Dynamics J 1308 0.000414 0.064128 9.488680e-01 0 4 0.998852
Decoder-null Dynamics G 4350 0.004852 1.340602 1.800496e-01 0 9 1.020218
Decoder-null Dynamics J 1310 0.015775 2.328034 1.991029e-02 0 4 1.089798
No Dynamics (decoder-null offset) G 4358 0.000190 0.055496 9.557430e-01 0 9 1.000579
No Dynamics (decoder-null offset) J 1302 -0.001300 -0.210459 8.333091e-01 0 4 1.002596
Null: subsample (1000 draws)
model monkey session tuples mean_norm mw_p
Full Dynamics G 0 301 2.2605 0.0000
No Dynamics G 0 237 1.0104 0.9487
Decoder-null Dynamics G 0 239 1.0488 0.3846
No Dynamics (decoder-null offset) G 0 237 1.0016 0.8742
Full Dynamics G 1 299 1.6123 0.0000
No Dynamics G 1 212 0.9963 0.8488
Decoder-null Dynamics G 1 219 1.0440 0.4926
No Dynamics (decoder-null offset) G 1 211 1.0006 0.9174
Full Dynamics G 2 141 1.5402 0.0001
No Dynamics G 2 138 1.0069 0.9513
Decoder-null Dynamics G 2 139 1.0139 0.8556
No Dynamics (decoder-null offset) G 2 137 0.9973 0.9223
Full Dynamics G 3 334 1.4352 0.0000
No Dynamics G 3 295 0.9889 0.6387
Decoder-null Dynamics G 3 294 1.0216 0.7221
No Dynamics (decoder-null offset) G 3 294 0.9941 0.7885
Full Dynamics G 4 221 1.4070 0.0000
No Dynamics G 4 222 1.0046 0.9566
Decoder-null Dynamics G 4 220 1.0127 0.8381
No Dynamics (decoder-null offset) G 4 224 1.0008 0.9811
Full Dynamics G 5 324 1.1811 0.0000
No Dynamics G 5 296 0.9991 0.9004
Decoder-null Dynamics G 5 283 1.0180 0.8451
No Dynamics (decoder-null offset) G 5 296 1.0065 0.8874
Full Dynamics G 6 305 1.4164 0.0000
No Dynamics G 6 263 1.0069 0.8864
Decoder-null Dynamics G 6 269 1.0154 0.7990
No Dynamics (decoder-null offset) G 6 261 1.0035 0.9394
Full Dynamics G 7 346 1.3377 0.0000
No Dynamics G 7 326 0.9993 0.8527
Decoder-null Dynamics G 7 324 1.0084 0.9918
No Dynamics (decoder-null offset) G 7 332 1.0033 0.9798
Full Dynamics G 8 247 1.4036 0.0000
No Dynamics G 8 189 1.0057 0.9280
Decoder-null Dynamics G 8 188 1.0001 0.8789
No Dynamics (decoder-null offset) G 8 187 0.9966 0.8844
Full Dynamics J 0 184 1.6360 0.0000
No Dynamics J 0 145 1.0081 0.8853
Decoder-null Dynamics J 0 146 1.1123 0.1588
No Dynamics (decoder-null offset) J 0 143 0.9671 0.5330
Full Dynamics J 1 181 1.5993 0.0000
No Dynamics J 1 150 0.9937 0.7509
Decoder-null Dynamics J 1 153 1.1319 0.0966
No Dynamics (decoder-null offset) J 1 149 1.0228 0.7994
Full Dynamics J 2 182 1.7355 0.0000
No Dynamics J 2 158 0.9821 0.6858
Decoder-null Dynamics J 2 156 1.0806 0.3923
No Dynamics (decoder-null offset) J 2 157 1.0180 0.8375
Full Dynamics J 3 224 1.3005 0.0000
No Dynamics J 3 201 1.0106 0.8441
Decoder-null Dynamics J 3 200 1.0320 0.7723
No Dynamics (decoder-null offset) J 3 202 1.0020 0.8137
model monkey N lme_effect z p sessions_mw_p05 sessions mean_norm
Full Dynamics G 5036 0.158369 25.973970 9.749962e-149 9 9 1.510435
Full Dynamics J 1542 0.136019 15.747259 7.172291e-56 4 4 1.567814
No Dynamics G 4356 -0.000040 -0.011576 9.907637e-01 0 9 1.002006
No Dynamics J 1308 0.000443 0.068658 9.452616e-01 0 4 0.998640
Decoder-null Dynamics G 4350 0.004810 1.328311 1.840755e-01 0 9 1.020303
Decoder-null Dynamics J 1310 0.015606 2.301731 2.135037e-02 0 4 1.089189
No Dynamics (decoder-null offset) G 4358 0.000187 0.054506 9.565317e-01 0 9 1.000463
No Dynamics (decoder-null offset) J 1302 -0.001283 -0.207832 8.353599e-01 0 4 1.002505
Paper (LME, simulation against shuffle): Full G z = -23.09, J z = -16.68; No Dynamics G z = 0.168, J 0.462; Decoder-null G 0.932, J -1.490; No Dynamics (decoder-null) G 0.611, J 0.449; Mann-Whitney p < 0.05 in 9 of 9 and 4 of 4 sessions for Full Dynamics only
Result: partially reproduced, resting on an assumption for Monkey J. Under both nulls, Full Dynamics activity for a command differs between conditions more than the null (mean normalized distance 1.49 to 1.51 for Monkey G and 1.55 to 1.57 for Monkey J; LME z = 25.4 to 26.0 and 15.5 to 15.7; Mann-Whitney p < 0.05 in 9 of 9 and 4 of 4 sessions). Both No Dynamics models show no difference (LME p 0.83 to 0.999, no session significant). Decoder-null dynamics show none for Monkey G (LME p = 0.18) and in no session, but the Monkey J LME finds a small difference (effect 0.016, z = 2.30 to 2.33, p = 0.020 to 0.021; mean normalized distance 1.09), so that absence is not reproduced. With the summed 5 ms gain for Monkey J the result is the same: Full Dynamics distances exceed the null (z = 15.53 and 15.79), the No Dynamics models show no difference, and the decoder-null LME again finds one (p = 0.0059 and 0.0051, Mann-Whitney p < 0.05 in 1 of 4 sessions).
Monkey J: sensitivity to the 100 ms decoder¶
Monkey J's 100 ms velocity gain, fitted over 100 ms bins, differs from the sum of the 5 ms gains (cell above). Here the simulations of Monkey J are repeated with the velocity gain fixed to the summed 5 ms gain (the command matrix of the other analyses) and the velocity rows of the cursor dynamics refitted with it, and the tests of fig6c, fig6d and fig6g are repeated on them.
JK = [k for k in KEYS if k[0] == "J"]
def job_summed(k):
return k, ofc.run_session(RAW[k], BB[k], k[0], 2000 + k[1], NSHUFF, j_gain="summed")
with ProcessPoolExecutor(4) as ex:
SIM_S = dict(ex.map(job_summed, JK))
for k in JK:
F, K, vaf = ofc.cursor_update(RAW[k], BB[k], "J", j_gain="summed")
print(k, "summed-gain decoder: variance accounted for (pos x, pos y, vel x, vel y)", np.round(vaf, 4), "velocity decay", round(F[2, 2], 4),
"max final error (cm)", round(max(r["final_err"].max() for r in SIM_S[k]["models"].values()), 3))
for model, base in [("Full Dynamics", "No Dynamics"), ("Decoder-null Dynamics", "No Dynamics (decoder-null offset)")]:
rows, lme = [], []
for k in JK:
a = pd.DataFrame(SIM_S[k]["models"][model]["mags"]).groupby(["task", "target"]).mag.mean()
c = pd.DataFrame(SIM_S[k]["models"][base]["mags"]).groupby(["task", "target"]).mag.mean()
rows.append(dict(session=k[1], ratio=a.mean() / c.mean(), wilcoxon_p=scipy.stats.wilcoxon(a.values, c.values).pvalue))
lme += [dict(session=k[1], model=1, mag=v) for v in a.values] + [dict(session=k[1], model=0, mag=v) for v in c.values]
L = pd.DataFrame(lme)
fit = smf.mixedlm("mag ~ model", L, groups=L["session"]).fit()
print(f"{model} against {base}, Monkey J, summed gain:")
print(pd.DataFrame(rows).round(4).to_string(index=False))
print(f" LME effect {fit.params['model']:.4f}, z = {fit.tvalues['model']:.2f}, p = {fit.pvalues['model']:.3g}, N = {len(L)}")
out = []
for null in ["shuffle", "subsample"]:
for name in SIM_S[JK[0]]["models"]:
rows, lme = [], []
for k in JK:
r = SIM_S[k]["models"][name]
if r["n_tuples"] == 0:
continue
obs, nm = r["dist_" + null], r["null_mean_" + null]
rows.append(dict(p=scipy.stats.mannwhitneyu(obs, nm).pvalue, norm=(obs / nm).mean()))
lme += [dict(session=k[1], source=1, dist=v) for v in obs] + [dict(session=k[1], source=0, dist=v) for v in nm]
L = pd.DataFrame(lme)
fit = smf.mixedlm("dist ~ source", L, groups=L["session"]).fit()
out.append(dict(null=null, model=name, N=len(L), lme_effect=fit.params["source"], z=fit.tvalues["source"], p=fit.pvalues["source"],
sessions_mw_p05=sum(x["p"] < 0.05 for x in rows), sessions=len(rows), mean_norm=np.mean([x["norm"] for x in rows])))
print("\nFigure 6G tests, Monkey J, summed gain:")
print(pd.DataFrame(out).round(4).to_string(index=False))
('J', 0) summed-gain decoder: variance accounted for (pos x, pos y, vel x, vel y) [0.9998 0.9999 0.9839 0.99 ] velocity decay 0.543 max final error (cm) 0.194
('J', 1) summed-gain decoder: variance accounted for (pos x, pos y, vel x, vel y) [0.9999 0.9999 0.9844 0.9913] velocity decay 0.5583 max final error (cm) 0.218
('J', 2) summed-gain decoder: variance accounted for (pos x, pos y, vel x, vel y) [0.9999 0.9999 0.9865 0.9923] velocity decay 0.5519 max final error (cm) 0.197
('J', 3) summed-gain decoder: variance accounted for (pos x, pos y, vel x, vel y) [0.9999 0.9999 0.994 0.9935] velocity decay 0.6755 max final error (cm) 0.159
Full Dynamics against No Dynamics, Monkey J, summed gain:
session ratio wilcoxon_p
0 0.4389 0.0
1 0.4160 0.0
2 0.4271 0.0
3 0.5972 0.0
LME effect -3.7533, z = -16.61, p = 5.9e-62, N = 192
Decoder-null Dynamics against No Dynamics (decoder-null offset), Monkey J, summed gain:
session ratio wilcoxon_p
0 1.0009 0.4223
1 0.9964 0.0016
2 1.0012 0.2522
3 1.0028 0.0738
LME effect -0.0018, z = -0.01, p = 0.995, N = 192
Figure 6G tests, Monkey J, summed gain:
null model N lme_effect z p sessions_mw_p05 sessions mean_norm
shuffle Full Dynamics 1460 0.1286 15.5297 0.0000 4 4 1.5887
shuffle No Dynamics 1234 -0.0017 -0.2752 0.7831 0 4 0.9969
shuffle Decoder-null Dynamics 1222 0.0182 2.7532 0.0059 1 4 1.1047
shuffle No Dynamics (decoder-null offset) 1218 -0.0014 -0.2319 0.8166 0 4 0.9983
subsample Full Dynamics 1460 0.1307 15.7895 0.0000 4 4 1.6069
subsample No Dynamics 1234 -0.0017 -0.2881 0.7733 0 4 0.9965
subsample Decoder-null Dynamics 1222 0.0185 2.7981 0.0051 1 4 1.1059
subsample No Dynamics (decoder-null offset) 1218 -0.0013 -0.2233 0.8233 0 4 0.9983
Sensitivity. The summed 5 ms gain accounts for 0.984 to 0.994 of the variance of the 100 ms velocity, against 0.995 to 0.998 for the fitted 100 ms update (data check above), and every simulated trial still ends within 0.218 cm of its target. With it, fig6c-dynamics-reduce-input holds as before; the absence of fig6d-decoder-null-dynamics-no-saving fails in one Monkey J session (Wilcoxon p = 0.0016, LME p = 0.995); fig6g-ofc-uses-distinct-patterns is partially reproduced as before. The fitted 100 ms update, the closer approximation of the decoder at the simulation's step, stays the adopted reading.
Summary¶
| Claim | Outcome |
|---|---|
| fig6c-dynamics-reduce-input | reproduced (assumption for Monkey J) |
| fig6d-decoder-null-dynamics-no-saving | reproduced (assumption for Monkey J) |
| fig6g-ofc-uses-distinct-patterns | partially reproduced (assumption for Monkey J) |
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 decoders' 100 ms cursor update is recovered from the files (exactly for Monkey G; by least squares over 100 ms bins for Monkey J, where the paper derived it from the PPF's parameters). The dynamics are the in-sample ridge fits on counts per bin, with the session's cross-validated parameter. The terminal weight (100) is the code's. For Monkey J a second reading, the velocity gain fixed to the summed 5 ms gain, is shown; under it the absence of fig6d-decoder-null-dynamics-no-saving fails in one session. Commands in the simulations are discretized with the simulation's own magnitude percentiles.
What was found about the dataset. The decoder parameters the simulations need are not in the files; Monkey G's can be recovered exactly from the decoder state.
print(f"notebook ran in {(time.time() - started) / 60:.1f} min")
notebook ran in 5.5 min