Reproducing Neupane, Fiete and Jazayeri (2024) from dandiset 000897¶
Paper. Neupane S, Fiete I, Jazayeri M (2024). Mental navigation in the primate entorhinal
cortex. Nature 630, 704-711. doi:10.1038/s41586-024-07557-z.
The claims file (claims/000897/10.1038-s41586-024-07557-z.yaml) was written from the version of
record in Europe PMC (PMC11224022): main text, Methods and Extended Data legends, with the
correction of 4 July 2024, and the paper's code (github.com/jazlab/Neupane_Fiete_Jazayeri_mental_navigation,
commit 4bee109).
Dandiset. DANDI:000897, version
0.240605.1710 (doi:10.48324/dandi.000897/0.240605.1710). The paper links the dandiset without a
version; this is its only published version, and the draft lists the same 30 files with the same
total size. Only the *_behavior+ecephys.nwb file of each session is read (trials, sorted units,
eye and joystick traces); the *_ecephys.nwb files hold the raw voltage.
What the dandiset holds. 15 sessions (monkey A 6, monkey M 9) of the paper's 32 (A 17, M 15), with 364 and 729 units against the paper's 614 and 864. Proportions, not counts, can therefore be compared with the paper's.
| Claim | Tested here |
|---|---|
| fig1bc-produced-matches-actual | tested (recorded sessions only) |
| fig1e-immediate-generalization | blocked: the generalization sessions are not in the dandiset |
| fig2-periodic-proportion | tested (proportions) |
| fig2d-periodicity-at-065 | tested |
| ed3-no-eye-hand-periodicity | tested; in 10 of 15 sessions the eye and joystick clock offset is estimated from the joystick trace |
| ed2b-pi-navigation-above-iti | blocked: the ITI cannot be located in the files |
| ed2c-pi-error-below-correct | tested |
| fig2e-first-peak-near-065 | tested |
| fig2f-earlier-peaks-multiples | tested |
| ed4d-overshoot-longer-period | tested |
| fig3b-onset-and-offset-coding | tested |
| fig3d-offset-coding-stronger | tested |
| ed4e-rate-codes-distance-over-slope | tested |
| ed4fg-distance-over-landmark | tested |
| fig3e-xcorr-conserved-iti | blocked: the ITI cannot be located in the files |
| ed5-xcorr-conserved-other-epochs | tested (its comparisons do not use the ITI) |
| fig4h-periodic-fano-factor | tested (proportions) |
| ed8-fano-periodic-without-rate-periodicity | tested |
| ed10-motor-preparation-signal | tested (proportions) |
| fig4g-reset-model-fits-better | tested (the 15 recorded sessions) |
| fig4b-learned-landmark-input | blocked: a model simulation, not attempted in the runs' time |
| fig4d-model-subscalar-variability | tested (a simulation of the paper's model, ported from its code) |
| ed6d-model-fano-periodic | blocked: a model simulation, not attempted in the runs' time |
| ed5efg-tighter-phase-clustering | blocked: the model's phase distribution was not simulated |
Methods come from the claims file, and details it lacks from the paper, cited in mnav.py and
in each section. The notebook takes about 30 minutes, most of it the PI of every unit, the
surrogate nulls, the ANOVAs, streaming the eye and joystick traces and the CAN model.
import time
import warnings
import numpy as np
import pandas as pd
import statsmodels.api as sm
import statsmodels.formula.api as smf
from scipy import stats
import mnav
pd.set_option("display.width", 220)
pd.set_option("display.max_columns", 40)
pd.set_option("display.max_rows", 300)
warnings.filterwarnings("ignore")
started = time.time()
sessions = mnav.load_all()
print(f"loaded {len(sessions)} sessions from dandiset {mnav.DANDISET_ID} version {mnav.DANDISET_VERSION} "
f"in {time.time() - started:.0f} s")
loaded 15 sessions from dandiset 000897 version 0.240605.1710 in 41 s
Data check¶
Each check of the metadata against the values is stated under its output.
rows = []
for s in sessions:
tr = s["trials"]
q = s["units"].get("quality", s["units"].get("KSLabel"))
rows.append({"monkey": s["monkey"], "session": s["session"], "units": len(s["spikes"]),
"label_column": "quality" if "quality" in s["units"] else "KSLabel",
**{f"label_{k}": int(v) for k, v in q.value_counts().items()},
"trials": len(tr), "mnav": int((tr.trial_type == 3).sum()),
"nts_type1": int((tr.trial_type == 1).sum()), "type2": int((tr.trial_type == 2).sum()),
"seqq": "/".join(f"{int(k)}:{v}" for k, v in tr.seqq.value_counts().sort_index().items()),
"attempt": "/".join(f"{int(k)}:{v}" for k, v in tr.attempt.value_counts().sort_index().items())})
info = pd.DataFrame(rows).fillna(0)
print(info.to_string(index=False))
print("\nSessions per monkey:", info.groupby("monkey").size().to_dict(), "(paper: A 17, M 15)")
print("Units per monkey, all labels:", info.groupby("monkey").units.sum().to_dict(), "(paper: A 614, M 864)")
noise = info.get("label_noise", pd.Series(0, index=info.index))
print("Units per monkey without those labelled noise:", (info.units - noise).groupby(info.monkey).sum().to_dict())
monkey session units label_column label_noise label_mua label_good label_unsorted trials mnav nts_type1 type2 seqq attempt label_nan
A 08152019 59 quality 23.0 16.0 12.0 8.0 2078 1710 368 0 1:1314/2:764 1:1018/2:881/99:179 0.0
A 08212019 57 quality 27.0 20.0 8.0 1.0 2629 1915 714 0 1:1685/2:944 1:1481/2:1013/99:135 1.0
A 08232019 43 quality 12.0 3.0 19.0 9.0 2741 2195 546 0 1:1713/2:1028 1:1368/2:1107/99:266 0.0
A 08272019 46 quality 6.0 2.0 29.0 9.0 2988 2105 883 0 1:1600/2:1388 1:1547/2:1157/99:284 0.0
A 08292019 95 quality 35.0 27.0 30.0 0.0 2434 1800 634 0 1:801/2:1633 1:1390/2:813/99:231 3.0
A 09172019 64 quality 0.0 0.0 0.0 64.0 2257 1261 996 0 1:2257 1:1322/2:832/99:103 0.0
M 04122021 98 KSLabel 0.0 87.0 11.0 0.0 3254 2372 0 882 1:1607/2:2/3:1645 1:2450/2:749/999:55 0.0
M 04132021 93 KSLabel 0.0 69.0 24.0 0.0 2223 1877 0 346 1:1102/2:3/3:1118 1:1534/2:612/999:77 0.0
M 04142021 112 KSLabel 0.0 56.0 56.0 0.0 2055 1716 0 339 1:1012/2:1/3:1042 1:1392/2:590/999:73 0.0
M 04212021 81 KSLabel 0.0 48.0 33.0 0.0 2641 1970 19 652 1:1425/2:2/3:1214 1:2057/2:514/999:70 0.0
M 04222021 60 KSLabel 0.0 52.0 8.0 0.0 1445 1314 0 131 1:861/2:3/3:581 1:823/2:538/999:84 0.0
M 04252021 64 KSLabel 0.0 50.0 14.0 0.0 2888 2146 0 742 1:1344/2:7/3:1537 1:2115/2:720/999:53 0.0
M 05042021 64 KSLabel 0.0 41.0 23.0 0.0 2979 2397 0 582 1:1491/2:1/3:1487 1:2181/2:653/999:145 0.0
M 05062021 103 KSLabel 0.0 71.0 32.0 0.0 3015 2372 0 643 1:1247/2:1/3:1767 1:2341/2:521/999:153 0.0
M 05072021 54 KSLabel 0.0 47.0 7.0 0.0 2589 1887 617 85 1:1142/2:4/3:1443 1:1850/2:614/999:125 0.0
Sessions per monkey: {'A': 6, 'M': 9} (paper: A 17, M 15)
Units per monkey, all labels: {'A': 364, 'M': 729} (paper: A 614, M 864)
Units per monkey without those labelled noise: {'A': 261.0, 'M': 729.0}
The dandiset holds 15 sessions (A 6, M 9) of the paper's 32 (A 17, M 15), with 364 units in
monkey A and 729 in M, against 614 and 864 in the paper. The monkey A files label units with
phy's quality (good, mua, noise, unsorted); 103 of the 364 are labelled noise, which the
paper's "single units and multi-units" would not include. The monkey M files carry only
Kilosort's KSLabel (good, mua). Units are analysed here as the files hold them, all labels
included; the count without noise units is given as a variant where it matters.
rows = []
for s in sessions:
tr = s["trials"]
m = (tr.trial_type == 3) & (tr.attempt <= 2)
d = tr[m]
joy = d.joy1offttl - d.joy1onttl
by_seq = {}
for q, g in d.groupby("seqq"):
by_seq[f"ta/0.65==target-curr, seqq {int(q)}"] = round(float(np.mean(np.isclose(g.ta / 0.65, g.target - g.curr))), 3)
sp_all = np.concatenate([x for x in s["spikes"] if len(x)])
rows.append({"monkey": s["monkey"], "session": s["session"],
"ta_values": len(np.unique(d.ta)), "ta_range": f"{d.ta.min():.2f}..{d.ta.max():.2f}",
"|joy off-on| - |tp| p1,p99 (ms)": f"{1000 * np.percentile(joy - d.tp.abs(), 1):.0f},{1000 * np.percentile(joy - d.tp.abs(), 99):.0f}",
"sign(tp)==sign(ta)": round(float(np.mean(np.sign(d.tp) == np.sign(d.ta))), 3),
**by_seq,
"trials_s": f"{tr.start_time.min():.0f}..{tr.stop_time.max():.0f}",
"spikes_s": f"{sp_all.min():.0f}..{sp_all.max():.0f}",
"eye_ts_s": f"{s['info']['eye_ts_first']:.0f}..{s['info']['eye_ts_last']:.0f}",
"hand_ts_s": f"{s['info']['hand_ts_first']:.0f}..{s['info']['hand_ts_last']:.0f}"})
print(pd.DataFrame(rows).to_string(index=False))
print("\nDistances (|ta|/0.65) of MNAV trials, all sessions:",
pd.concat([s["trials"][s["trials"].trial_type == 3].distance for s in sessions]).value_counts().sort_index().to_dict())
for m in ("A", "M"):
t3 = pd.concat([s["trials"] for s in sessions if s["monkey"] == m])
t3 = t3[t3.trial_type == 3]
print(f"monkey {m}, |ta| values (s) per seqq:", {int(q): sorted(np.unique(np.round(g.ta.abs(), 3)).tolist()) for q, g in t3.groupby("seqq")})
print("tp == 99 or 999 (no production) rows:", sum(int((s["trials"].tp.abs() >= 99).sum()) for s in sessions),
"; attempt codes 99/999:", sum(int((s["trials"].attempt >= 99).sum()) for s in sessions))
monkey session ta_values ta_range |joy off-on| - |tp| p1,p99 (ms) sign(tp)==sign(ta) ta/0.65==target-curr, seqq 1 ta/0.65==target-curr, seqq 2 trials_s spikes_s eye_ts_s hand_ts_s ta/0.65==target-curr, seqq 3
A 08152019 10 -3.25..3.25 -16,20 0.949 1.0 1.0 5..18974 0..19081 775418..795185 775418..795184 NaN
A 08212019 10 -3.25..3.25 -18,17 0.974 1.0 1.0 5..16780 0..16817 -99..16805 -99..16804 NaN
A 08232019 10 -3.25..3.25 -15,14 0.939 1.0 1.0 5..19035 0..19040 -99..19214 -99..19214 NaN
A 08272019 10 -3.25..3.25 -15,14 0.944 1.0 1.0 4..20261 0..20275 -58..20265 -58..20265 NaN
A 08292019 10 -3.25..3.25 -16,16 0.977 1.0 1.0 5..15686 0..15729 -37..15716 -37..15715 NaN
A 09172019 10 -3.25..3.25 -15,15 0.919 1.0 NaN 4..13778 0..13867 -37..13855 -34..13855 NaN
M 04122021 18 -4.88..4.88 -16,17 0.977 1.0 NaN 4..20655 0..20674 430976..452273 431522..452273 0.0
M 04132021 18 -4.88..4.88 -16,15 0.965 1.0 NaN 4..18597 0..9306 515768..538443 515768..538443 0.0
M 04142021 18 -4.88..4.88 -15,16 0.950 1.0 1.0 4..19761 0..19792 603985..623867 603985..623867 0.0
M 04212021 18 -4.88..4.88 -16,15 0.985 1.0 NaN 6..14041 0..14060 1207146..1221274 1207146..1221274 0.0
M 04222021 18 -4.88..4.88 -16,16 0.945 1.0 NaN 4..16679 0..16700 1298303..1315006 1298303..1315006 0.0
M 04252021 18 -4.88..4.88 -16,16 0.972 1.0 NaN 5..15468 0..15490 1559866..1576931 1559866..1576931 0.0
M 05042021 18 -4.88..4.88 -17,17 0.976 1.0 NaN 5..16585 0..16603 3209..19917 3297..19917 0.0
M 05062021 18 -4.88..4.88 -16,15 0.990 1.0 NaN 4..16677 0..16696 172389..189122 172389..189121 0.0
M 05072021 18 -4.88..4.88 -16,15 0.983 1.0 NaN 5..18282 0..18293 258417..276725 258417..276725 0.0
Distances (|ta|/0.65) of MNAV trials, all sessions: {1.0: 8570, 2.0: 7931, 3.0: 6082, 4.0: 4001, 5.0: 2453}
monkey A, |ta| values (s) per seqq: {1: [0.65, 1.3, 1.95, 2.6, 3.25], 2: [0.65, 1.3, 1.95, 2.6, 3.25]}
monkey M, |ta| values (s) per seqq: {1: [0.65, 1.3, 1.95, 2.6, 3.25], 2: [1.3], 3: [0.975, 1.95, 2.925, 3.9, 4.875]}
tp == 99 or 999 (no production) rows: 2033 ; attempt codes 99/999: 2033
The trial table's ta takes the ten values +-0.65 to +-3.25 s, as the six-landmark line with
0.65 s spacing has, and the joystick offset minus onset equals |tp| within about 20 ms, so
tp is the produced vector in seconds, signed as ta. In monkey A ta/0.65 equals target -
curr on every trial. In monkey M, about half of the MNAV trials have seqq 3, whose vectors
come in steps of 0.975 s (0.975 to 4.875 s, printed above), not the 0.65 s the paper describes for
both monkeys; the paper and the files do not describe this sequence. Analyses that use distance
keep seqq < 3, as the paper's code does for the PI (claims, opening comment).
Rows with attempt 99 or 999 have tp 99 or 999 (no production) and are not used.
Spike times and trial times share one range in every session. The eye and joystick timestamps do not in several sessions: they start hundreds of thousands of seconds away from the trial times (A 08152019; M 04122021, 04132021, 04142021, 04212021, 04222021, 04252021, 05062021, 05072021), so in those sessions the behavioural traces are on another clock, which the files do not describe. The eye and joystick traces are used only for claim ed3, in the sessions whose timestamps span the trial times.
# The trial table's columns, what start_time marks, and which attempt the joystick times of a
# trial with two attempts (attempt == 2) belong to.
tr0 = sessions[0]["trials"]
print("intervals/trials columns:", [c for c in tr0.columns if c not in ("distance", "direction")])
print("start_time == stim1onttl in every session:",
all(np.allclose(s["trials"].start_time, s["trials"].stim1onttl) for s in sessions))
rows = []
for s in sessions:
tr = s["trials"]
for a in (1, 2):
g = tr[(tr.trial_type == 3) & (tr.attempt == a)]
rows.append({"monkey": s["monkey"], "session": s["session"], "attempt": a, "rows": len(g),
"go cue to joy1on, median s": round((g.joy1onttl - g.gocuettl).median(), 2),
"joy1off to row end, median s": round((g.stop_time - g.joy1offttl).median(), 2),
"| |tp| - |ta| | / |ta|, median": round(((g.tp.abs() - g.ta.abs()).abs() / g.ta.abs()).median(), 2),
"succ mean": round(g.succ.mean(), 2)})
print(pd.DataFrame(rows).to_string(index=False))
intervals/trials columns: ['attempt', 'curr', 'delay', 'gocuettl', 'joy1offttl', 'joy1onttl', 'numrepeat', 'rt', 'seqq', 'start_time', 'stim1onttl', 'stop_time', 'succ', 'ta', 'target', 'tp', 'trial_type', 'validtrials_mm'] start_time == stim1onttl in every session: True
monkey session attempt rows go cue to joy1on, median s joy1off to row end, median s | |tp| - |ta| | / |ta|, median succ mean
A 08152019 1 754 0.44 0.98 0.08 1.00
A 08152019 2 804 0.53 2.76 0.38 0.31
A 08212019 1 936 0.55 0.91 0.09 1.00
A 08212019 2 882 0.56 2.75 0.33 0.27
A 08232019 1 982 0.41 0.87 0.10 1.00
A 08232019 2 999 0.51 2.73 0.41 0.26
A 08272019 1 953 0.48 0.93 0.10 1.00
A 08272019 2 950 0.53 2.73 0.34 0.22
A 08292019 1 978 0.47 0.87 0.09 1.00
A 08292019 2 666 0.46 2.68 0.31 0.29
A 09172019 1 610 0.57 0.87 0.08 1.00
A 09172019 2 590 0.66 2.68 0.35 0.36
M 04122021 1 1676 0.39 1.24 0.10 1.00
M 04122021 2 658 0.50 3.00 0.37 0.27
M 04132021 1 1239 0.39 1.26 0.12 1.00
M 04132021 2 576 0.51 2.68 0.38 0.21
M 04142021 1 1085 0.40 1.27 0.12 1.00
M 04142021 2 567 0.52 4.23 0.45 0.21
M 04212021 1 1473 0.39 1.22 0.10 1.00
M 04212021 2 457 0.43 2.48 0.31 0.32
M 04222021 1 728 0.41 1.30 0.10 1.00
M 04222021 2 505 0.53 5.42 0.50 0.22
M 04252021 1 1435 0.40 1.25 0.10 1.00
M 04252021 2 669 0.43 2.39 0.29 0.35
M 05042021 1 1731 0.43 1.23 0.10 1.00
M 05042021 2 550 0.49 2.63 0.30 0.45
M 05062021 1 1815 0.41 1.23 0.10 1.00
M 05062021 2 445 0.44 2.53 0.29 0.39
M 05072021 1 1299 0.43 1.22 0.10 1.00
M 05072021 2 518 0.49 2.93 0.37 0.25
The trial table has no fixation or ITI column: start_time equals stim1onttl, the onset of the
start landmark. Each row is one trial. On rows with attempt 2 the joystick is pressed as soon
after the go cue as on attempt 1 rows (medians printed above), and the row then runs longer after
the joystick offset (medians 2.4 to 5.4 s against 0.9 to 1.3 s), with relative errors of about 0.3 to 0.5 against 0.1.
This is consistent with the row's tp and joystick times being the first attempt, followed by a
second attempt the table does not time; the files do not say so.
# Clock check: the paper's motor response before joystick offset (Extended Data Fig. 10). For each
# session, the fraction of units whose mean rate differs (Welch t-test, p < 0.01) between 0.4-0.1 s
# and 1.0-0.6 s before joystick offset on first-attempt MNAV trials, with the true offset times and
# with the offset times of each trial replaced by those of a trial 5 rows later (a shifted label).
rows = []
for s in sessions:
tr = s["trials"]
d = tr[(tr.trial_type == 3) & (tr.attempt == 1) & (tr.tp.abs() > 1.5)]
off = d.joy1offttl.values
out = {"monkey": s["monkey"], "session": s["session"], "trials": len(d)}
for name, t in (("true", off), ("shifted", np.roll(off, 5) + 7.0)):
sig = 0
for sp in s["spikes"]:
late = mnav.mean_rate(sp, t - 0.4, t - 0.1)
early = mnav.mean_rate(sp, t - 1.0, t - 0.6)
if np.std(late - early) > 0 and stats.ttest_rel(late, early).pvalue < 0.01:
sig += 1
out[f"frac_units_{name}"] = round(sig / len(s["spikes"]), 3)
rows.append(out)
clock = pd.DataFrame(rows)
print(clock.to_string(index=False))
monkey session trials frac_units_true frac_units_shifted
A 08152019 319 0.102 0.000
A 08212019 393 0.053 0.000
A 08232019 418 0.023 0.000
A 08272019 401 0.065 0.000
A 08292019 416 0.432 0.011
A 09172019 258 0.000 0.000
M 04122021 710 0.490 0.020
M 04132021 456 0.000 0.000
M 04142021 419 0.000 0.000
M 04212021 698 0.444 0.062
M 04222021 307 0.350 0.000
M 04252021 807 0.266 0.047
M 05042021 739 0.141 0.000
M 05062021 897 0.505 0.184
M 05072021 594 0.889 0.019
In 12 of the 15 sessions more units change their rate before the true joystick offset than before the shifted times (printed above), which is consistent with spikes and trials sharing a clock there. In three sessions (A 09172019, M 04132021, M 04142021) no unit does so with either timing, and while preparing this notebook the population rate of these sessions was flat (within a few percent) around the stimulus onset, go cue, joystick onset and joystick offset (not run here). In M 04132021 the spikes also stop at 9,306 s while trials run to 18,599 s. For these three sessions a shared clock could not be established; they are kept, as the paper analysed all units, and results are broken down by session where it matters.
fig1bc-produced-matches-actual¶
In MNAV the produced vectors match the actual vectors in direction and magnitude; performance is quantified by the slope of the regression of v_p on v_a.
Criterion (claims file): per recorded session, OLS of signed v_p on signed v_a over all MNAV
trials with a recorded production (code: 0.01 < |tp| and tp < 99 s); reproduced when, in each
monkey, the median session slope is at least 0.8 and every slope is positive. Variant: the
validtrials_mm trials only.
rows = []
for s in sessions:
tr = s["trials"]
d = tr[(tr.trial_type == 3) & (tr.tp.abs() > 0.01) & (tr.tp < 99)]
v = d[d.validtrials_mm == 1]
fit = stats.linregress(d.ta, d.tp)
fitv = stats.linregress(v.ta, v.tp)
rows.append({"monkey": s["monkey"], "session": s["session"], "trials": len(d), "slope_all": fit.slope,
"r_all": fit.rvalue, "valid_trials": len(v), "slope_valid": fitv.slope})
fig1 = pd.DataFrame(rows)
print(fig1.round(3).to_string(index=False))
print("\nPer monkey:")
print(fig1.groupby("monkey").agg(sessions=("slope_all", "size"), median_slope_all=("slope_all", "median"),
min_slope_all=("slope_all", "min"), max_slope_all=("slope_all", "max"),
median_slope_valid=("slope_valid", "median")).round(3).to_string())
for m, g in fig1.groupby("monkey"):
print(f"monkey {m}: median slope >= 0.8: {g.slope_all.median() >= 0.8}; all slopes positive: {(g.slope_all > 0).all()}")
monkey session trials slope_all r_all valid_trials slope_valid
A 08152019 1550 0.889 0.811 1252 0.942
A 08212019 1803 0.855 0.933 1321 0.867
A 08232019 1965 0.843 0.808 876 1.015
A 08272019 1892 0.825 0.857 1558 0.897
A 08292019 1639 0.907 0.944 1407 0.930
A 09172019 1132 0.946 0.929 986 0.997
M 04122021 2325 0.867 0.832 2035 0.869
M 04132021 1799 0.785 0.797 1592 0.815
M 04142021 1638 0.687 0.569 1422 0.753
M 04212021 1926 0.879 0.936 1716 0.882
M 04222021 1210 0.893 0.617 958 0.904
M 04252021 2071 0.777 0.924 1904 0.813
M 05042021 2269 0.786 0.907 2049 0.838
M 05062021 2258 0.866 0.949 2063 0.886
M 05072021 1814 0.860 0.900 1628 0.889
Per monkey:
sessions median_slope_all min_slope_all max_slope_all median_slope_valid
monkey
A 6 0.872 0.825 0.946 0.936
M 9 0.860 0.687 0.893 0.869
monkey A: median slope >= 0.8: True; all slopes positive: True
monkey M: median slope >= 0.8: True; all slopes positive: True
Result: reproduced. The median session slope is 0.872 in monkey A (6 sessions, 0.825 to
0.946) and 0.860 in monkey M (9 sessions, 0.687 to 0.893), so at least 0.8 in both, and every
slope is positive. One monkey M session (04142021) is below 0.8. On the validtrials_mm trials
the medians are 0.936 and 0.869. These are the recorded sessions only, a subset of the
behavioural sessions of Fig. 1c.
The periodicity index (PI) of every unit¶
Computed once, as the claims file's opening comment and mnav.py describe, for both lag ranges
(Methods 0 to 1,300 ms, code 100 to 1,600 ms), per unit and direction, on first-attempt MNAV
trials over at least three landmarks with |tp| > 1 s and seqq < 3.
t = time.time()
pit = mnav.periodicity_table(sessions)
print(f"PI curves for {len(pit)} unit-direction pairs in {time.time() - t:.0f} s")
print("Units with more than 15 kept trials, per monkey and direction:")
print(pit[pit.k > mnav.MIN_TRIALS].groupby(["monkey", "direction"]).size().to_string())
print("Trials in the PI per session and direction (median over sessions):",
pit.drop_duplicates(["session", "direction"]).groupby(["monkey", "direction"]).n_trials.median().to_dict())
print("Monkey M trials with seqq 3 left out by seqq < 3:",
sum(int(((s["trials"].trial_type == 3) & (s["trials"].seqq == 3)).sum()) for s in sessions if s["monkey"] == "M"))
PI curves for 2186 unit-direction pairs in 112 s
Units with more than 15 kept trials, per monkey and direction:
monkey direction
A -1 291
1 292
M -1 439
1 397
Trials in the PI per session and direction (median over sessions): {('A', -1): 189.5, ('A', 1): 191.0, ('M', -1): 122.0, ('M', 1): 117.0}
Monkey M trials with seqq 3 left out by seqq < 3: 9332
# Null distributions (Methods: GP surrogate, Poisson, smoothed as the data; 1,000 repeats) and the
# code's variant without the Poisson step ("gp"). The own-duration null of a unit averages as many
# surrogate trials as the unit keeps (k, matched to the nearest level of a grid in steps of 25%),
# with durations drawn from its session and direction's trials. The global null averages the
# monkey's median k trials, all at its median duration, from a bank generated at that duration.
rng_null = np.random.default_rng(897)
LAGSETS = {"methods": mnav.LAGS_METHODS, "code": mnav.LAGS_CODE}
NULL_TYPES = {"poisson": True, "gp": False}
def thresholds(null):
m, sd = null.mean(0), null.std(0)
scores = mnav.pi_scores(null)
return {"per_lag": m + 2 * sd, "a1": float((m + 2 * sd).max()), "a2": float(scores.mean() + 2 * scores.std()),
"score_mean": float(scores.mean()), "score_sd": float(scores.std())}
kept = pit.k > mnav.MIN_TRIALS
pit["k_level"] = 0
pit.loc[kept, "k_level"] = mnav.k_levels(pit.loc[kept, "k"].values)[1]
thr, banks, surr_rate = {}, {}, {}
t = time.time()
for monkey in ("A", "M"):
sub = pit[(pit.monkey == monkey)]
groups = sub.drop_duplicates(["session", "direction"])
all_durations = np.concatenate(groups.durations.values)
rates = []
for s in sessions:
if s["monkey"] != monkey:
continue
trs = mnav.pi_trials(s["trials"])
t0, t1 = mnav.nav_windows(trs)
rates += [float(np.mean(mnav.mean_rate(sp, t0, t1))) for sp in s["spikes"]]
rate = surr_rate[monkey] = float(np.median(rates))
k_med = int(np.median(sub[sub.k > mnav.MIN_TRIALS].k))
d_med = float(np.median(all_durations))
for null_type, poisson in NULL_TYPES.items():
bank, bank_d = mnav.surrogate_bank(rng_null, all_durations, rate, 4000, poisson=poisson)
gbank, gbank_d = mnav.surrogate_bank(rng_null, [d_med], rate, 3000, poisson=poisson)
banks[(monkey, null_type)] = (bank, bank_d)
for lag, curves in mnav.null_from_bank(rng_null, gbank, gbank_d, [d_med], k_med, LAGSETS).items():
thr[(monkey, "global", null_type, lag)] = thresholds(curves)
for _, g in groups.iterrows():
levels = sub[(sub.session == g.session) & (sub.direction == g.direction) & (sub.k > mnav.MIN_TRIALS)].k_level.unique()
for lev in levels:
for lag, curves in mnav.null_from_bank(rng_null, bank, bank_d, g.durations, int(lev), LAGSETS).items():
thr[(monkey, g.session, g.direction, int(lev), null_type, lag)] = thresholds(curves)
print(f"monkey {monkey}: surrogate rate {rate:.2f} spikes/s (median unit rate); global null with {k_med} trials "
f"of {d_med:.0f} ms; banks of {len(bank)} surrogate trials (durations of the PI trials) and {len(gbank)} (median duration)")
print(f"{len(thr)} nulls in {time.time() - t:.0f} s")
print("Kept trials per unit (k) as a fraction of its session and direction's PI trials, units with k > 15:",
pit[kept].assign(f=pit.k / pit.n_trials).groupby("monkey").f.quantile([0.1, 0.25, 0.5, 0.75]).round(2).to_dict())
monkey A: surrogate rate 1.65 spikes/s (median unit rate); global null with 106 trials of 1748 ms; banks of 4000 surrogate trials (durations of the PI trials) and 3000 (median duration)
monkey M: surrogate rate 0.51 spikes/s (median unit rate); global null with 72 trials of 1672 ms; banks of 3998 surrogate trials (durations of the PI trials) and 2999 (median duration)
1060 nulls in 274 s
Kept trials per unit (k) as a fraction of its session and direction's PI trials, units with k > 15: {('A', 0.1): 0.2, ('A', 0.25): 0.39, ('A', 0.5): 0.69, ('A', 0.75): 0.92, ('M', 0.1): 0.21, ('M', 0.25): 0.37, ('M', 0.5): 0.68, ('M', 0.75): 0.93}
# Classification under every reading. Rules: (a1) score = max PI - mean PI from the first
# positive lag on, above the maximum over lags of the surrogate PI mean + 2 SD; (a2) the same score
# above the surrogate scores' mean + 2 SD; (b) PI at any lag above the surrogate mean + 2 SD at that
# lag. Null durations: own (the session and direction's trials, k-matched) or global. Lags: methods
# or code. Null: poisson (Methods) or gp (code's variant, no Poisson step).
READINGS = [(rule, dur, lag, nt) for nt in NULL_TYPES for rule in ("a1", "a2", "b") for dur in ("own", "global")
for lag in ("methods", "code")]
def thr_key(row, dur, lag, nt):
if dur == "global":
return (row.monkey, "global", nt, lag)
return (row.monkey, row.session, row.direction, int(row.k_level), nt, lag)
def classify(row, rule, dur, lag, nt):
pi = getattr(row, "pi_" + lag)
if pi is None or row.k <= mnav.MIN_TRIALS:
return False
th = thr[thr_key(row, dur, lag, nt)]
if rule == "b":
return bool(np.any(pi > th["per_lag"]))
return mnav.pi_score(pi) > th[rule]
def rname(r):
return "/".join(r)
for r in READINGS:
pit["per_" + "_".join(r)] = [classify(row, *r) for row in pit.itertuples()]
for lag, lags in LAGSETS.items():
pit[f"period_{lag}"] = [np.nan if p is None else lags[int(np.argmax(p))] / 1000 for p in pit["pi_" + lag]]
pit[f"score_{lag}"] = [np.nan if p is None else mnav.pi_score(p) for p in pit["pi_" + lag]]
MAIN = ("a2", "own", "methods", "poisson")
main_col = "per_" + "_".join(MAIN)
# periodic sets carried to the claims that use periodic units: the main reading and the readings
# that differ from it in the rule or the null
SETS = {"a2 (main)": main_col, "a1": "per_a1_own_methods_poisson", "b": "per_b_own_methods_poisson",
"a2, gp null": "per_a2_own_methods_gp"}
i650 = int(np.nonzero(mnav.LAGS_METHODS == 650)[0][0])
base = pit[["monkey", "session", "unit", "direction", "k", "k_level", "pi_methods", "period_methods", "score_methods"]
+ list(dict.fromkeys(SETS.values()))].copy()
base["above650"] = [False if (p is None or r.k <= mnav.MIN_TRIALS) else
p[i650] > thr[thr_key(r, "own", "methods", "poisson")]["per_lag"][i650]
for p, r in zip(base.pi_methods, base.itertuples())]
print("Thresholds, global null (rule a1, rule a2; score mean and SD):")
for key in [k for k in thr if k[1] == "global"]:
print(" ", key, {k: round(v, 3) for k, v in thr[key].items() if k != "per_lag"})
own = pd.DataFrame([{"monkey": k[0], "session": k[1], "direction": k[2], "k_level": k[3], "null": k[4], "lags": k[5],
"a1": v["a1"], "a2": v["a2"]} for k, v in thr.items() if k[1] != "global"])
print("\nOwn-duration null thresholds (rule a2), by k level, median over sessions and directions:")
print(own[own.lags == "methods"].pivot_table(index="k_level", columns=["monkey", "null"], values="a2", aggfunc="median").round(3).to_string())
print("\nOwn-duration null thresholds (rule a2, Methods lags, Poisson) per session and direction, range over k levels:")
print(own[(own.lags == "methods") & (own.null == "poisson")].groupby(["monkey", "session", "direction"]).a2
.agg(["min", "max", "size"]).round(3).to_string())
print("\nNeuron scores (methods lags), median and 90th percentile per monkey:",
pit[pit.k > mnav.MIN_TRIALS].groupby("monkey").score_methods.quantile([0.5, 0.9]).round(3).to_dict())
Thresholds, global null (rule a1, rule a2; score mean and SD):
('A', 'global', 'poisson', 'methods') {'a1': 0.355, 'a2': 0.133, 'score_mean': 0.08, 'score_sd': 0.027}
('A', 'global', 'poisson', 'code') {'a1': 0.355, 'a2': 0.162, 'score_mean': 0.095, 'score_sd': 0.033}
('A', 'global', 'gp', 'methods') {'a1': 0.548, 'a2': 0.181, 'score_mean': 0.112, 'score_sd': 0.034}
('A', 'global', 'gp', 'code') {'a1': 0.548, 'a2': 0.217, 'score_mean': 0.136, 'score_sd': 0.04}
('M', 'global', 'poisson', 'methods') {'a1': 0.409, 'a2': 0.191, 'score_mean': 0.097, 'score_sd': 0.047}
('M', 'global', 'poisson', 'code') {'a1': 0.409, 'a2': 0.234, 'score_mean': 0.123, 'score_sd': 0.055}
('M', 'global', 'gp', 'methods') {'a1': 0.585, 'a2': 0.192, 'score_mean': 0.12, 'score_sd': 0.036}
('M', 'global', 'gp', 'code') {'a1': 0.585, 'a2': 0.253, 'score_mean': 0.156, 'score_sd': 0.049}
Own-duration null thresholds (rule a2), by k level, median over sessions and directions:
monkey A M
null gp poisson gp poisson
k_level
16 0.407 0.357 0.412 0.291
20 0.367 0.328 0.386 0.268
25 0.329 0.293 0.340 0.256
31 0.310 0.271 0.314 0.222
39 0.283 0.241 0.292 0.209
49 0.264 0.219 0.267 0.193
61 0.249 0.200 0.248 0.187
76 0.233 0.188 0.245 0.173
95 0.222 0.172 0.231 0.164
119 0.210 0.161 0.226 0.159
149 0.198 0.150 0.223 0.153
186 0.190 0.140 NaN NaN
Own-duration null thresholds (rule a2, Methods lags, Poisson) per session and direction, range over k levels:
min max size
monkey session direction
A 08152019 -1 0.152 0.362 11
1 0.150 0.358 10
08212019 -1 0.140 0.324 10
1 0.140 0.354 11
08232019 -1 0.144 0.356 9
1 0.139 0.364 10
08272019 -1 0.136 0.349 8
1 0.141 0.363 8
08292019 -1 0.132 0.349 12
1 0.151 0.324 10
09172019 -1 0.154 0.343 10
1 0.161 0.364 10
M 04122021 -1 0.132 0.283 10
1 0.138 0.291 10
04132021 -1 0.170 0.283 7
1 0.232 0.299 5
04142021 -1 0.159 0.290 8
1 0.232 0.323 6
04212021 -1 0.145 0.301 11
1 0.144 0.292 10
04222021 -1 0.139 0.282 8
1 0.197 0.295 4
04252021 -1 0.166 0.296 8
1 0.166 0.304 9
05042021 -1 0.162 0.286 7
1 0.159 0.258 6
05062021 -1 0.164 0.315 8
1 0.159 0.278 9
05072021 -1 0.151 0.275 9
1 0.131 0.268 9
Neuron scores (methods lags), median and 90th percentile per monkey: {('A', 0.5): 0.089, ('A', 0.9): 0.21, ('M', 0.5): 0.094, ('M', 0.9): 0.167}
# Check of the bank: one group's null built from the bank against a null generated afresh (each
# surrogate trial a new GP and Poisson draw), at the group's largest k level, 300 surrogates.
rng_chk = np.random.default_rng(8970)
g0 = pit[(pit.monkey == "A") & (pit.k > mnav.MIN_TRIALS)].sort_values(["session", "direction"]).iloc[0]
lev = int(pit[(pit.session == g0.session) & (pit.direction == g0.direction) & (pit.k > mnav.MIN_TRIALS)].k_level.max())
t = time.time()
fresh = mnav.surrogate_pis(rng_chk, g0.durations, lev, surr_rate["A"], mnav.LAGS_METHODS, n=300)
bank, bank_d = banks[("A", "poisson")]
from_bank = mnav.null_from_bank(rng_chk, bank, bank_d, g0.durations, lev, {"methods": mnav.LAGS_METHODS})["methods"]
for name, null in (("fresh", fresh), ("bank", from_bank), ("bank, as used", None)):
th = thresholds(null) if null is not None else thr[("A", g0.session, g0.direction, lev, "poisson", "methods")]
print(f"monkey A {g0.session} direction {int(g0.direction)}, k = {lev}, {name}: a2 threshold {th['a2']:.3f} "
f"(score mean {th['score_mean']:.3f}, SD {th['score_sd']:.3f}); a1 threshold {th['a1']:.3f}; "
f"per-lag threshold at 650 ms {th['per_lag'][i650]:.3f}")
print(f"in {time.time() - t:.0f} s")
monkey A 08152019 direction -1, k = 149, fresh: a2 threshold 0.141 (score mean 0.081, SD 0.030); a1 threshold 0.360; per-lag threshold at 650 ms 0.349 monkey A 08152019 direction -1, k = 149, bank: a2 threshold 0.155 (score mean 0.089, SD 0.033); a1 threshold 0.377; per-lag threshold at 650 ms 0.365 monkey A 08152019 direction -1, k = 149, bank, as used: a2 threshold 0.152 (score mean 0.088, SD 0.032); a1 threshold 0.377; per-lag threshold at 650 ms 0.358 in 55 s
The nulls are now drawn at random from the whole neighbourhood of 200 bank entries nearest in duration (all entries for the global null, whose bank is generated at the median duration), and each unit's own-duration null averages as many surrogate trials as the unit keeps (k). Units keep a median of 68 to 69% of their session and direction's PI trials (10th percentile 20%), and the a2 threshold falls from about 0.36 (A) and 0.29 (M) at k = 16 to 0.14 to 0.15 at the largest k (printed above), so matching k matters. For one group the bank null gives an a2 threshold of 0.155 (0.152 as used) against 0.141 for a null generated afresh, a difference of about 0.4 SD of the null scores (0.03), so the bank thresholds are slightly conservative.
fig2-periodic-proportion¶
A substantial fraction of EC neurons have significantly periodic activity during mental navigation, 37% (231/614) in monkey A and 36% (311/864) in monkey M, and the proportion per session ranged from 0% to 45% depending on recording site.
A unit counts as periodic when it is periodic in either direction; the denominator is all units of the monkey. The dandiset holds 6 of the 17 monkey A sessions and 9 of the 15 monkey M sessions, so the counts cannot reach 231 and 311; they are compared as proportions of the units the dandiset holds (231/614 and 311/864, each with the claims file's 10% relative tolerance), and the raw counts are printed.
units_total = info.groupby("monkey").units.sum()
noise_by_monkey = (noise).groupby(info.monkey).sum()
paper = {"A": (231, 614), "M": (311, 864)}
labels = pit.drop_duplicates(["monkey", "session", "unit"])[["monkey", "session", "unit", "label"]]
rows, sess_rows = [], []
for r in READINGS:
col = "per_" + "_".join(r)
per_unit = pit.groupby(["monkey", "session", "unit"])[col].any().reset_index()
out = {"reading": rname(r)}
for m in ("A", "M"):
n = int(per_unit[per_unit.monkey == m][col].sum())
p_ref = paper[m][0] / paper[m][1]
out[f"{m}_count"] = n
out[f"{m}_pct"] = round(100 * n / units_total[m], 1)
nn = per_unit[(per_unit.monkey == m)].merge(labels, on=["monkey", "session", "unit"])
nn = nn[nn.label != "noise"]
out[f"{m}_pct_no_noise"] = round(100 * nn[col].mean(), 1)
out[f"{m}_in_tol"] = abs(n / units_total[m] - p_ref) <= 0.1 * p_ref
ses = per_unit.groupby("session")[col].mean() * 100
out["session_pct_min"] = round(ses.min(), 1)
out["session_pct_max"] = round(ses.max(), 1)
out["sessions_in_0_49.5"] = bool((ses <= 49.5).all())
out["conditions_met"] = int(out["A_in_tol"]) + int(out["M_in_tol"]) + int(out["sessions_in_0_49.5"])
rows.append(out)
sess_rows.append(ses.rename(rname(r)))
fig2 = pd.DataFrame(rows)
print("Paper: A 231/614 = 37.6%, M 311/864 = 36.0%; tolerance 10% of each proportion (A 33.9-41.4%, M 32.4-39.6%).")
print(fig2.to_string(index=False))
print("\nPer-session percentage of periodic units under each reading:")
print(pd.concat(sess_rows, axis=1).round(1).T.to_string())
no_clock = ["09172019", "04132021", "04142021"]
pu = pit[~pit.session.isin(no_clock)].groupby(["monkey", "session", "unit"])[main_col].any().reset_index()
print("\nMain reading without the three sessions whose clock could not be checked (data check):",
{m: f"{int(g[main_col].sum())}/{len(g)} = {100 * g[main_col].mean():.1f}%" for m, g in pu.groupby("monkey")})
Paper: A 231/614 = 37.6%, M 311/864 = 36.0%; tolerance 10% of each proportion (A 33.9-41.4%, M 32.4-39.6%).
reading A_count A_pct A_pct_no_noise A_in_tol M_count M_pct M_pct_no_noise M_in_tol session_pct_min session_pct_max sessions_in_0_49.5 conditions_met
a1/own/methods/poisson 10 2.7 3.8 False 0 0.0 0.0 False 0.0 9.5 True 1
a1/own/code/poisson 15 4.1 5.7 False 0 0.0 0.0 False 0.0 12.6 True 1
a1/global/methods/poisson 14 3.8 5.4 False 1 0.1 0.1 False 0.0 9.5 True 1
a1/global/code/poisson 19 5.2 7.3 False 3 0.4 0.4 False 0.0 11.6 True 1
a2/own/methods/poisson 56 15.4 20.7 False 57 7.8 7.8 False 1.0 30.6 True 1
a2/own/code/poisson 44 12.1 16.5 False 43 5.9 5.9 False 0.0 24.5 True 1
a2/global/methods/poisson 111 30.5 35.6 False 43 5.9 5.9 False 0.0 43.8 True 1
a2/global/code/poisson 114 31.3 35.6 False 51 7.0 7.0 False 0.0 43.8 True 1
b/own/methods/poisson 239 65.7 82.4 False 379 52.0 52.0 False 21.5 73.9 False 0
b/own/code/poisson 241 66.2 83.1 False 385 52.8 52.8 False 22.6 76.1 False 0
b/global/methods/poisson 244 67.0 79.7 False 421 57.8 57.8 False 40.9 74.1 False 0
b/global/code/poisson 246 67.6 80.5 False 420 57.6 57.6 False 40.9 74.1 False 0
a1/own/methods/gp 4 1.1 1.5 False 0 0.0 0.0 False 0.0 4.2 True 1
a1/own/code/gp 5 1.4 1.9 False 0 0.0 0.0 False 0.0 5.3 True 1
a1/global/methods/gp 4 1.1 1.5 False 0 0.0 0.0 False 0.0 4.2 True 1
a1/global/code/gp 4 1.1 1.5 False 0 0.0 0.0 False 0.0 4.2 True 1
a2/own/methods/gp 32 8.8 12.3 False 8 1.1 1.1 False 0.0 21.1 True 1
a2/own/code/gp 26 7.1 10.0 False 1 0.1 0.1 False 0.0 20.0 True 1
a2/global/methods/gp 57 15.7 19.2 False 42 5.8 5.8 False 0.0 26.3 True 1
a2/global/code/gp 59 16.2 19.9 False 35 4.8 4.8 False 0.0 26.3 True 1
b/own/methods/gp 296 81.3 98.5 False 436 59.8 59.8 False 41.9 95.3 False 0
b/own/code/gp 296 81.3 98.5 False 436 59.8 59.8 False 41.9 95.3 False 0
b/global/methods/gp 296 81.3 98.5 False 440 60.4 60.4 False 41.9 95.3 False 0
b/global/code/gp 296 81.3 98.5 False 440 60.4 60.4 False 41.9 95.3 False 0
Per-session percentage of periodic units under each reading:
session 04122021 04132021 04142021 04212021 04222021 04252021 05042021 05062021 05072021 08152019 08212019 08232019 08272019 08292019 09172019
a1/own/methods/poisson 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.8 0.0 0.0 9.5 0.0
a1/own/code/poisson 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.7 3.5 0.0 0.0 12.6 0.0
a1/global/methods/poisson 0.0 0.0 0.9 0.0 0.0 0.0 0.0 0.0 0.0 1.7 3.5 0.0 0.0 9.5 3.1
a1/global/code/poisson 0.0 0.0 1.8 0.0 1.7 0.0 0.0 0.0 0.0 1.7 5.3 0.0 0.0 11.6 6.2
a2/own/methods/poisson 30.6 2.2 7.1 6.2 3.3 6.2 6.2 1.0 1.9 6.8 17.5 2.3 8.7 28.4 15.6
a2/own/code/poisson 24.5 0.0 7.1 6.2 3.3 3.1 1.6 1.0 0.0 6.8 14.0 0.0 6.5 22.1 12.5
a2/global/methods/poisson 14.3 4.3 16.1 1.2 1.7 3.1 3.1 1.0 0.0 22.0 31.6 14.0 23.9 36.8 43.8
a2/global/code/poisson 13.3 5.4 18.8 3.7 5.0 4.7 1.6 1.9 0.0 32.2 33.3 14.0 17.4 35.8 43.8
b/own/methods/poisson 62.2 21.5 37.5 66.7 40.0 67.2 70.3 49.5 72.2 64.4 54.4 65.1 73.9 71.6 62.5
b/own/code/poisson 62.2 22.6 37.5 66.7 41.7 70.3 73.4 49.5 72.2 64.4 54.4 65.1 76.1 71.6 64.1
b/global/methods/poisson 57.1 40.9 48.2 71.6 55.0 71.9 68.8 50.5 74.1 66.1 57.9 60.5 71.7 72.6 68.8
b/global/code/poisson 56.1 40.9 48.2 71.6 55.0 71.9 68.8 50.5 74.1 66.1 61.4 60.5 71.7 72.6 68.8
a1/own/methods/gp 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 4.2 0.0
a1/own/code/gp 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 5.3 0.0
a1/global/methods/gp 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 4.2 0.0
a1/global/code/gp 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 4.2 0.0
a2/own/methods/gp 5.1 0.0 1.8 1.2 0.0 0.0 0.0 0.0 0.0 1.7 10.5 0.0 2.2 21.1 6.2
a2/own/code/gp 0.0 0.0 0.0 1.2 0.0 0.0 0.0 0.0 0.0 1.7 7.0 0.0 0.0 20.0 3.1
a2/global/methods/gp 13.3 4.3 16.1 1.2 1.7 3.1 3.1 1.0 0.0 6.8 17.5 0.0 4.3 26.3 25.0
a2/global/code/gp 9.2 2.2 13.4 1.2 5.0 3.1 1.6 1.9 0.0 8.5 19.3 0.0 4.3 26.3 25.0
b/own/methods/gp 66.3 41.9 50.9 70.4 55.0 75.0 73.4 50.5 70.4 76.3 71.9 76.7 91.3 77.9 95.3
b/own/code/gp 66.3 41.9 50.9 70.4 55.0 75.0 73.4 50.5 70.4 76.3 71.9 76.7 91.3 77.9 95.3
b/global/methods/gp 67.3 41.9 51.8 70.4 55.0 75.0 73.4 50.5 74.1 76.3 71.9 76.7 91.3 77.9 95.3
b/global/code/gp 67.3 41.9 51.8 70.4 55.0 75.0 73.4 50.5 74.1 76.3 71.9 76.7 91.3 77.9 95.3
Main reading without the three sessions whose clock could not be checked (data check): {'A': '46/300 = 15.3%', 'M': '47/524 = 9.0%'}
Result: not reproduced, resting on an assumption. Under no reading is either monkey's proportion
within 10% of the paper's (A 37.6%, M 36.0%). Under the main reading (rule a2, own-duration null
matched to each unit's kept trials, Methods lags, Poisson null) 56 of 364 units in monkey A (15.4%)
and 57 of 729 in monkey M (7.8%) are periodic, with session proportions of 1.0% to 30.6%. The a2
readings give 7.1% to 31.3% (A) and 0.1% to 7.8% (M); a1 at most 5.2% and 0.4%; b 65.7% to 81.3% and
52.0% to 60.4%, with session maxima of 73.9% to 95.3%, above 49.5%. Only the per-session range is
met, under the a1 and a2 readings, so at most one of the three conditions holds under any reading.
The code's GP null without the Poisson step gives fewer periodic units than the Methods' null under
a1 and a2 (A 1.1% to 16.2%, M 0% to 5.8%). Without units labelled noise, the main reading gives
20.7% in monkey A; without the three sessions whose clock could not be checked, 15.3% (A, 46/300)
and 9.0% (M, 47/524). The conversion of the GP into a firing rate is an assumption (mnav.py), and
the counts are judged as proportions because the dandiset holds 15 of the 32 sessions.
fig2d-periodicity-at-065¶
Across neurons with significant PI, the lag of the peak PI was at or near 0.65 s, the temporal distance between landmarks.
Criterion: median periodicity (lag of the maximum PI) of the periodic units within 0.60 to 0.70 s in every monkey and direction; the mode in 50 ms bins is reported as well.
rows = []
for r in READINGS:
col, lag = "per_" + "_".join(r), r[2]
for (m, d), g in pit[pit[col]].groupby(["monkey", "direction"]):
p = g[f"period_{lag}"]
hist, edges = np.histogram(p, bins=np.arange(0, 1.651, 0.05))
rows.append({"reading": rname(r), "monkey": m, "direction": int(d), "n": len(p),
"median_s": round(float(p.median()), 3), "iqr_s": f"{p.quantile(.25):.2f}-{p.quantile(.75):.2f}",
"mode_bin_s": f"{edges[hist.argmax()]:.2f}-{edges[hist.argmax() + 1]:.2f}",
"frac_0.6_0.7": round(float(((p >= 0.6) & (p <= 0.7)).mean()), 3),
"in_0.60_0.70": 0.60 <= p.median() <= 0.70})
fig2d = pd.DataFrame(rows)
print(fig2d.to_string(index=False))
print("\nGroups (monkey x direction) with median in 0.60-0.70 s, per reading:")
print(fig2d.groupby("reading")["in_0.60_0.70"].agg(["sum", "size"]).to_string())
main_sel = pit[pit[main_col]]
print("\nMain reading (a2/own/methods/poisson), all periodic unit-directions pooled: median periodicity",
round(float(main_sel.period_methods.median()), 3), "s; per monkey:", main_sel.groupby("monkey").period_methods.median().round(3).to_dict())
allp = pit[pit.k > mnav.MIN_TRIALS]
print("\nReference: median periodicity of all units with >15 trials (methods lags):",
allp.groupby(["monkey", "direction"]).period_methods.median().round(3).to_dict())
reading monkey direction n median_s iqr_s mode_bin_s frac_0.6_0.7 in_0.60_0.70
a1/own/methods/poisson A -1 10 0.565 0.49-0.58 0.55-0.60 0.200 False
a1/own/methods/poisson A 1 6 0.525 0.49-0.54 0.50-0.55 0.000 False
a1/own/code/poisson A -1 14 0.570 0.50-0.62 0.55-0.60 0.286 False
a1/own/code/poisson A 1 7 0.530 0.50-0.54 0.50-0.55 0.000 False
a1/global/methods/poisson A -1 11 0.560 0.48-0.57 0.55-0.60 0.182 False
a1/global/methods/poisson A 1 9 0.540 0.52-0.55 0.50-0.55 0.111 False
a1/global/methods/poisson M -1 1 0.520 0.52-0.52 0.50-0.55 0.000 False
a1/global/code/poisson A -1 14 0.565 0.49-0.61 0.55-0.60 0.214 False
a1/global/code/poisson A 1 12 0.545 0.53-0.61 0.50-0.55 0.333 False
a1/global/code/poisson M -1 1 0.520 0.52-0.52 0.50-0.55 0.000 False
a1/global/code/poisson M 1 2 0.685 0.62-0.75 0.55-0.60 0.000 True
a2/own/methods/poisson A -1 37 0.580 0.53-0.66 0.55-0.60 0.324 False
a2/own/methods/poisson A 1 31 0.590 0.54-0.66 0.55-0.60 0.387 False
a2/own/methods/poisson M -1 30 0.650 0.59-0.70 0.65-0.70 0.500 True
a2/own/methods/poisson M 1 34 0.640 0.62-0.67 0.60-0.65 0.765 True
a2/own/code/poisson A -1 26 0.595 0.54-0.66 0.55-0.60 0.423 False
a2/own/code/poisson A 1 27 0.600 0.54-0.67 0.50-0.55 0.407 True
a2/own/code/poisson M -1 18 0.665 0.63-0.74 0.65-0.70 0.500 True
a2/own/code/poisson M 1 29 0.640 0.63-0.69 0.60-0.65 0.655 True
a2/global/methods/poisson A -1 67 0.620 0.57-0.70 0.55-0.60 0.313 True
a2/global/methods/poisson A 1 65 0.650 0.58-0.73 0.55-0.60 0.400 True
a2/global/methods/poisson M -1 22 0.635 0.57-0.67 0.55-0.60 0.455 True
a2/global/methods/poisson M 1 27 0.620 0.60-0.67 0.60-0.65 0.667 True
a2/global/code/poisson A -1 72 0.630 0.57-0.72 0.55-0.60 0.319 True
a2/global/code/poisson A 1 68 0.665 0.57-0.74 0.65-0.70 0.353 True
a2/global/code/poisson M -1 22 0.640 0.58-0.69 0.65-0.70 0.455 True
a2/global/code/poisson M 1 35 0.670 0.62-0.71 0.65-0.70 0.600 True
b/own/methods/poisson A -1 186 0.720 0.59-0.98 0.55-0.60 0.210 False
b/own/methods/poisson A 1 184 0.810 0.61-1.04 0.55-0.60 0.174 False
b/own/methods/poisson M -1 333 0.960 0.77-1.09 1.00-1.05 0.126 False
b/own/methods/poisson M 1 285 0.920 0.74-1.08 0.90-0.95 0.147 False
b/own/code/poisson A -1 190 0.720 0.59-0.97 0.55-0.60 0.216 False
b/own/code/poisson A 1 189 0.800 0.60-1.04 0.55-0.60 0.180 False
b/own/code/poisson M -1 350 0.950 0.76-1.09 1.00-1.05 0.134 False
b/own/code/poisson M 1 298 0.920 0.71-1.09 0.90-0.95 0.164 False
b/global/methods/poisson A -1 202 0.695 0.59-0.94 0.55-0.60 0.228 True
b/global/methods/poisson A 1 193 0.740 0.60-0.99 0.55-0.60 0.223 False
b/global/methods/poisson M -1 390 0.945 0.76-1.09 1.00-1.05 0.121 False
b/global/methods/poisson M 1 332 0.900 0.74-1.07 0.90-0.95 0.136 False
b/global/code/poisson A -1 203 0.700 0.59-0.94 0.55-0.60 0.227 True
b/global/code/poisson A 1 193 0.740 0.60-1.00 0.55-0.60 0.218 False
b/global/code/poisson M -1 391 0.950 0.77-1.10 1.00-1.05 0.118 False
b/global/code/poisson M 1 336 0.900 0.74-1.07 0.90-0.95 0.134 False
a1/own/methods/gp A -1 4 0.565 0.45-0.57 0.55-0.60 0.000 False
a1/own/methods/gp A 1 4 0.535 0.43-0.54 0.50-0.55 0.000 False
a1/own/code/gp A -1 5 0.560 0.48-0.57 0.55-0.60 0.000 False
a1/own/code/gp A 1 4 0.535 0.43-0.54 0.50-0.55 0.000 False
a1/global/methods/gp A -1 4 0.565 0.45-0.57 0.55-0.60 0.000 False
a1/global/methods/gp A 1 4 0.535 0.43-0.54 0.50-0.55 0.000 False
a1/global/code/gp A -1 4 0.565 0.45-0.57 0.55-0.60 0.000 False
a1/global/code/gp A 1 4 0.535 0.43-0.54 0.50-0.55 0.000 False
a2/own/methods/gp A -1 24 0.575 0.53-0.64 0.55-0.60 0.417 False
a2/own/methods/gp A 1 17 0.540 0.52-0.60 0.50-0.55 0.294 False
a2/own/methods/gp M -1 5 0.610 0.52-0.64 0.50-0.55 0.400 True
a2/own/methods/gp M 1 3 0.620 0.61-0.63 0.60-0.65 1.000 True
a2/own/code/gp A -1 23 0.580 0.53-0.64 0.55-0.60 0.435 False
a2/own/code/gp A 1 11 0.540 0.52-0.57 0.50-0.55 0.182 False
a2/own/code/gp M -1 1 0.770 0.77-0.77 0.75-0.80 0.000 False
a2/global/methods/gp A -1 38 0.580 0.55-0.63 0.55-0.60 0.395 False
a2/global/methods/gp A 1 30 0.595 0.54-0.65 0.55-0.60 0.400 False
a2/global/methods/gp M -1 22 0.635 0.57-0.67 0.55-0.60 0.455 True
a2/global/methods/gp M 1 25 0.630 0.60-0.67 0.60-0.65 0.640 True
a2/global/code/gp A -1 38 0.590 0.55-0.65 0.55-0.60 0.395 False
a2/global/code/gp A 1 35 0.600 0.55-0.68 0.55-0.60 0.400 True
a2/global/code/gp M -1 15 0.630 0.58-0.69 0.50-0.55 0.467 True
a2/global/code/gp M 1 23 0.640 0.61-0.71 0.60-0.65 0.565 True
b/own/methods/gp A -1 290 0.730 0.62-0.94 0.55-0.60 0.241 False
b/own/methods/gp A 1 288 0.760 0.63-0.98 0.55-0.60 0.247 False
b/own/methods/gp M -1 432 0.910 0.72-1.07 0.65-0.70 0.160 False
b/own/methods/gp M 1 387 0.860 0.68-1.05 0.60-0.65 0.222 False
b/own/code/gp A -1 290 0.730 0.62-0.94 0.55-0.60 0.241 False
b/own/code/gp A 1 288 0.760 0.63-0.98 0.55-0.60 0.247 False
b/own/code/gp M -1 432 0.910 0.72-1.08 0.65-0.70 0.157 False
b/own/code/gp M 1 387 0.860 0.69-1.05 0.60-0.65 0.217 False
b/global/methods/gp A -1 291 0.730 0.62-0.94 0.55-0.60 0.241 False
b/global/methods/gp A 1 291 0.760 0.63-0.98 0.55-0.60 0.244 False
b/global/methods/gp M -1 436 0.910 0.72-1.07 0.65-0.70 0.158 False
b/global/methods/gp M 1 394 0.860 0.68-1.05 0.60-0.65 0.221 False
b/global/code/gp A -1 291 0.730 0.62-0.94 0.55-0.60 0.241 False
b/global/code/gp A 1 291 0.760 0.63-0.98 0.55-0.60 0.244 False
b/global/code/gp M -1 436 0.910 0.72-1.08 0.65-0.70 0.156 False
b/global/code/gp M 1 394 0.860 0.69-1.05 0.65-0.70 0.216 False
Groups (monkey x direction) with median in 0.60-0.70 s, per reading:
sum size
reading
a1/global/code/gp 0 2
a1/global/code/poisson 1 4
a1/global/methods/gp 0 2
a1/global/methods/poisson 0 3
a1/own/code/gp 0 2
a1/own/code/poisson 0 2
a1/own/methods/gp 0 2
a1/own/methods/poisson 0 2
a2/global/code/gp 3 4
a2/global/code/poisson 4 4
a2/global/methods/gp 2 4
a2/global/methods/poisson 4 4
a2/own/code/gp 0 3
a2/own/code/poisson 3 4
a2/own/methods/gp 2 4
a2/own/methods/poisson 2 4
b/global/code/gp 0 4
b/global/code/poisson 1 4
b/global/methods/gp 0 4
b/global/methods/poisson 1 4
b/own/code/gp 0 4
b/own/code/poisson 0 4
b/own/methods/gp 0 4
b/own/methods/poisson 0 4
Main reading (a2/own/methods/poisson), all periodic unit-directions pooled: median periodicity 0.63 s; per monkey: {'A': 0.585, 'M': 0.65}
Reference: median periodicity of all units with >15 trials (methods lags): {('A', -1): 0.73, ('A', 1): 0.76, ('M', -1): 0.91, ('M', 1): 0.86}
Result: partially reproduced, resting on an assumption. Under the main reading the median periodicity of the periodic units is 0.580 s (A, direction -1), 0.590 s (A, 1), 0.650 s (M, -1) and 0.640 s (M, 1): within 0.60 to 0.70 s in the two monkey M groups and not in monkey A. The a2 readings put 0 to 4 of the groups within the range (4 of 4 with the global null and the Poisson null, 0 of 3 with the own-duration GP null and the code's lags); under a1, which selects 0 to 14 units per group, 0 or 1; under b, which selects most units, 0 or 1, with medians of 0.695 to 0.96 s, close to those of all units (0.73 to 0.91 s). The readings disagree; reading a2 is adopted, as it compares the code's score with the same score on the surrogates. Pooling all periodic unit-directions of the main reading, the median is 0.63 s.
ed2c-pi-error-below-correct¶
Neurons significantly periodic at 0.65 s had a lower PI on error trials than on correct trials, suggesting weaker periodicity contributed to larger errors.
PI measure (claims, ed2b and ed2c): PI at the 650 ms lag minus the mean PI from the first lag at
which PI is positive. Correct trials are the PI trials; error trials are the first attempt of
trials with more than one attempt (attempt == 2; Methods, MNAV; the data check shows the row's
joystick times are the first attempt's), with the same distance, |tp| and seqq rules; the
code's variant takes succ == 0 trials with |tp| > 1.8 s. Neurons, two readings: (i) code,
periodic units whose PI at 650 ms exceeds the surrogate PI mean + 2 SD at 650 ms; (ii) Results,
periodic units with periodicity in 0.60 to 0.70 s. Periodic units are those of each periodic set
(the main reading and the a1, b and GP-null readings). Both conditions need more than 15 trials.
One-tailed paired t-test, correct > error. As a control for the number of trials (correct trials
outnumber error trials, and a mean ACG over fewer trials is noisier), the main reading is also
run with the correct trials of each session and direction subsampled at random to the number of
error trials.
def pi650_measure(pi, lags=mnav.LAGS_METHODS):
i650 = int(np.nonzero(lags == 650)[0][0])
pos = np.nonzero(pi > 0)[0]
base_ = pi[pos[0]:].mean() if len(pos) else pi.mean()
return float(pi[i650] - base_)
rng2c = np.random.default_rng(2002)
def matched_correct(tr):
"""Control for trial count: per direction, a random subset of the correct PI trials as large as
the error PI trials (Methods' error trials)."""
out = pd.Series(False, index=tr.index)
for d in (1, -1):
c = mnav.pi_trials(tr, d)
e = mnav.pi_trials(tr, d, attempt=2)
out[rng2c.choice(c.index, min(len(c), len(e)), replace=False)] = True
return out
t = time.time()
pit_err = mnav.periodicity_table(sessions, attempt=2, lag_sets={"methods": mnav.LAGS_METHODS})
pit_err_code = mnav.periodicity_table(
sessions, attempt=2, lag_sets={"methods": mnav.LAGS_METHODS},
extra=lambda tr: (tr.succ == 0) & (tr.tp.abs() > 1.8))
pit_cmatch = mnav.periodicity_table(sessions, extra=matched_correct, lag_sets={"methods": mnav.LAGS_METHODS})
print(f"error-trial and matched correct-trial PI in {time.time() - t:.0f} s")
rows = []
for set_name, set_col in SETS.items():
for err_name, err in (("attempt2", pit_err), ("code_succ0_tp1.8", pit_err_code)):
for corr_name in ("all correct", "correct matched to error count"):
if corr_name != "all correct" and (err_name != "attempt2" or set_name != "a2 (main)"):
continue
mrg = base.merge(err[["monkey", "session", "unit", "direction", "k", "pi_methods"]],
on=["monkey", "session", "unit", "direction"], suffixes=("", "_err"))
if corr_name != "all correct":
mrg = mrg.drop(columns=["k", "pi_methods"]).merge(
pit_cmatch[["monkey", "session", "unit", "direction", "k", "pi_methods"]], on=["monkey", "session", "unit", "direction"])
for sel_name, sel in (("i_code_pi650", mrg[set_col] & mrg.above650),
("ii_results_0.6-0.7", mrg[set_col] & mrg.period_methods.between(0.6, 0.7))):
g = mrg[sel & (mrg.k_err > mnav.MIN_TRIALS) & (mrg.k > mnav.MIN_TRIALS)]
for (m, d), gg in g.groupby(["monkey", "direction"]):
c = np.array([pi650_measure(p) for p in gg.pi_methods])
e = np.array([pi650_measure(p) for p in gg.pi_methods_err])
tt = stats.ttest_rel(c, e, alternative="greater") if len(gg) > 1 else None
rows.append({"periodic_set": set_name, "error_trials": err_name, "correct_trials": corr_name, "neurons": sel_name,
"monkey": m, "direction": int(d), "n": len(gg),
"dropped_<=15_trials": int((sel & (mrg.monkey == m) & (mrg.direction == d)).sum() - len(gg)),
"k_correct_median": int(np.median(gg.k)), "k_error_median": int(np.median(gg.k_err)),
"correct_mean": round(c.mean(), 3), "error_mean": round(e.mean(), 3),
"t": None if tt is None else round(tt.statistic, 2), "p_one_tailed": None if tt is None else tt.pvalue,
"meets_p<1e-4": bool(tt is not None and tt.pvalue < 1e-4 and c.mean() > e.mean())})
ed2c = pd.DataFrame(rows)
print(ed2c.to_string(index=False))
print("\nGroups meeting the criterion per reading (of 4; groups with fewer than 2 neurons are not tested):")
print(ed2c.groupby(["periodic_set", "error_trials", "correct_trials", "neurons"], sort=False)["meets_p<1e-4"].agg(["sum", "size"]).to_string())
error-trial and matched correct-trial PI in 237 s
periodic_set error_trials correct_trials neurons monkey direction n dropped_<=15_trials k_correct_median k_error_median correct_mean error_mean t p_one_tailed meets_p<1e-4
a2 (main) attempt2 all correct i_code_pi650 A -1 18 0 180 109 0.319 0.225 4.80 8.399531e-05 True
a2 (main) attempt2 all correct i_code_pi650 A 1 18 0 102 92 0.266 0.155 3.38 1.776904e-03 False
a2 (main) attempt2 all correct i_code_pi650 M -1 10 0 59 26 0.203 0.121 2.71 1.193537e-02 False
a2 (main) attempt2 all correct i_code_pi650 M 1 24 0 93 62 0.209 0.139 4.40 1.040497e-04 False
a2 (main) attempt2 all correct ii_results_0.6-0.7 A -1 12 0 184 125 0.276 0.166 4.12 8.463947e-04 False
a2 (main) attempt2 all correct ii_results_0.6-0.7 A 1 12 0 78 85 0.248 0.149 2.10 2.963580e-02 False
a2 (main) attempt2 all correct ii_results_0.6-0.7 M -1 15 0 74 28 0.194 0.129 2.37 1.638907e-02 False
a2 (main) attempt2 all correct ii_results_0.6-0.7 M 1 26 0 93 66 0.206 0.130 5.05 1.646197e-05 True
a2 (main) attempt2 correct matched to error count i_code_pi650 A -1 18 0 120 109 0.327 0.225 4.38 2.036114e-04 False
a2 (main) attempt2 correct matched to error count i_code_pi650 A 1 18 0 102 92 0.267 0.155 3.42 1.635246e-03 False
a2 (main) attempt2 correct matched to error count i_code_pi650 M -1 10 0 28 26 0.156 0.121 1.15 1.395894e-01 False
a2 (main) attempt2 correct matched to error count i_code_pi650 M 1 24 0 63 62 0.195 0.139 3.08 2.638833e-03 False
a2 (main) attempt2 correct matched to error count ii_results_0.6-0.7 A -1 12 0 127 125 0.295 0.166 4.39 5.408360e-04 False
a2 (main) attempt2 correct matched to error count ii_results_0.6-0.7 A 1 12 0 78 85 0.250 0.149 2.14 2.789731e-02 False
a2 (main) attempt2 correct matched to error count ii_results_0.6-0.7 M -1 15 0 29 28 0.173 0.129 1.54 7.285219e-02 False
a2 (main) attempt2 correct matched to error count ii_results_0.6-0.7 M 1 26 0 63 66 0.190 0.130 3.57 7.435453e-04 False
a2 (main) code_succ0_tp1.8 all correct i_code_pi650 A -1 18 0 180 61 0.319 0.191 4.72 9.815036e-05 True
a2 (main) code_succ0_tp1.8 all correct i_code_pi650 A 1 16 2 104 46 0.256 0.081 4.86 1.039185e-04 False
a2 (main) code_succ0_tp1.8 all correct i_code_pi650 M -1 3 7 151 16 0.150 -0.137 15.15 2.164339e-03 False
a2 (main) code_succ0_tp1.8 all correct i_code_pi650 M 1 23 1 94 32 0.208 -0.023 10.02 5.824205e-10 True
a2 (main) code_succ0_tp1.8 all correct ii_results_0.6-0.7 A -1 12 0 184 68 0.276 0.114 4.20 7.382230e-04 False
a2 (main) code_succ0_tp1.8 all correct ii_results_0.6-0.7 A 1 9 3 130 69 0.220 0.053 3.49 4.083698e-03 False
a2 (main) code_succ0_tp1.8 all correct ii_results_0.6-0.7 M -1 7 8 140 16 0.159 -0.036 2.62 1.981284e-02 False
a2 (main) code_succ0_tp1.8 all correct ii_results_0.6-0.7 M 1 25 1 94 32 0.205 -0.037 11.33 2.035494e-11 True
a1 attempt2 all correct i_code_pi650 A -1 7 0 184 106 0.416 0.310 6.52 3.106248e-04 False
a1 attempt2 all correct i_code_pi650 A 1 4 0 153 159 0.326 0.292 0.80 2.417293e-01 False
a1 attempt2 all correct ii_results_0.6-0.7 A -1 2 0 196 126 0.388 0.227 8.25 3.840639e-02 False
a1 code_succ0_tp1.8 all correct i_code_pi650 A -1 7 0 184 58 0.416 0.284 2.73 1.705288e-02 False
a1 code_succ0_tp1.8 all correct i_code_pi650 A 1 4 0 153 63 0.326 0.293 0.93 2.097695e-01 False
a1 code_succ0_tp1.8 all correct ii_results_0.6-0.7 A -1 2 0 196 70 0.388 0.098 2.80 1.090133e-01 False
b attempt2 all correct i_code_pi650 A -1 18 0 180 109 0.319 0.225 4.80 8.399531e-05 True
b attempt2 all correct i_code_pi650 A 1 18 0 102 92 0.266 0.155 3.38 1.776904e-03 False
b attempt2 all correct i_code_pi650 M -1 12 0 57 29 0.202 0.118 3.35 3.248783e-03 False
b attempt2 all correct i_code_pi650 M 1 24 0 93 62 0.209 0.139 4.40 1.040497e-04 False
b attempt2 all correct ii_results_0.6-0.7 A -1 39 0 171 128 0.145 0.118 1.72 4.688857e-02 False
b attempt2 all correct ii_results_0.6-0.7 A 1 32 0 132 138 0.150 0.130 0.83 2.067168e-01 False
b attempt2 all correct ii_results_0.6-0.7 M -1 39 3 111 48 0.122 0.124 -0.10 5.399914e-01 False
b attempt2 all correct ii_results_0.6-0.7 M 1 41 1 114 65 0.157 0.099 4.05 1.132892e-04 False
b code_succ0_tp1.8 all correct i_code_pi650 A -1 18 0 180 61 0.319 0.191 4.72 9.815036e-05 True
b code_succ0_tp1.8 all correct i_code_pi650 A 1 16 2 104 46 0.256 0.081 4.86 1.039185e-04 False
b code_succ0_tp1.8 all correct i_code_pi650 M -1 5 7 43 18 0.170 -0.075 5.22 3.220971e-03 False
b code_succ0_tp1.8 all correct i_code_pi650 M 1 23 1 94 32 0.208 -0.023 10.02 5.824205e-10 True
b code_succ0_tp1.8 all correct ii_results_0.6-0.7 A -1 39 0 171 70 0.145 0.051 5.12 4.502578e-06 True
b code_succ0_tp1.8 all correct ii_results_0.6-0.7 A 1 30 2 136 69 0.137 0.026 4.61 3.735079e-05 True
b code_succ0_tp1.8 all correct ii_results_0.6-0.7 M -1 24 18 140 20 0.107 -0.075 6.73 3.666758e-07 True
b code_succ0_tp1.8 all correct ii_results_0.6-0.7 M 1 36 6 118 35 0.166 -0.041 9.75 8.283834e-12 True
a2, gp null attempt2 all correct i_code_pi650 A -1 17 0 181 112 0.321 0.241 5.37 3.123191e-05 True
a2, gp null attempt2 all correct i_code_pi650 A 1 11 0 103 93 0.302 0.202 2.29 2.239872e-02 False
a2, gp null attempt2 all correct i_code_pi650 M -1 3 0 60 19 0.218 0.084 2.06 8.800553e-02 False
a2, gp null attempt2 all correct i_code_pi650 M 1 3 0 47 30 0.297 0.089 53.08 1.774018e-04 False
a2, gp null attempt2 all correct ii_results_0.6-0.7 A -1 10 0 191 125 0.290 0.203 4.01 1.531706e-03 False
a2, gp null attempt2 all correct ii_results_0.6-0.7 A 1 5 0 44 41 0.311 0.172 1.61 9.118127e-02 False
a2, gp null attempt2 all correct ii_results_0.6-0.7 M -1 2 0 54 19 0.277 0.095 2.43 1.243834e-01 False
a2, gp null attempt2 all correct ii_results_0.6-0.7 M 1 3 0 47 30 0.297 0.089 53.08 1.774018e-04 False
a2, gp null code_succ0_tp1.8 all correct i_code_pi650 A -1 17 0 181 63 0.321 0.208 4.70 1.210041e-04 False
a2, gp null code_succ0_tp1.8 all correct i_code_pi650 A 1 9 2 106 46 0.294 0.164 3.08 7.522004e-03 False
a2, gp null code_succ0_tp1.8 all correct i_code_pi650 M -1 1 2 152 16 0.098 -0.189 NaN NaN False
a2, gp null code_succ0_tp1.8 all correct i_code_pi650 M 1 3 0 47 20 0.297 -0.041 8.54 6.713604e-03 False
a2, gp null code_succ0_tp1.8 all correct ii_results_0.6-0.7 A -1 10 0 191 68 0.290 0.159 3.43 3.730183e-03 False
a2, gp null code_succ0_tp1.8 all correct ii_results_0.6-0.7 A 1 2 3 116 57 0.276 0.159 3.91 7.962109e-02 False
a2, gp null code_succ0_tp1.8 all correct ii_results_0.6-0.7 M 1 3 0 47 20 0.297 -0.041 8.54 6.713604e-03 False
Groups meeting the criterion per reading (of 4; groups with fewer than 2 neurons are not tested):
sum size
periodic_set error_trials correct_trials neurons
a2 (main) attempt2 all correct i_code_pi650 1 4
ii_results_0.6-0.7 1 4
correct matched to error count i_code_pi650 0 4
ii_results_0.6-0.7 0 4
code_succ0_tp1.8 all correct i_code_pi650 2 4
ii_results_0.6-0.7 1 4
a1 attempt2 all correct i_code_pi650 0 2
ii_results_0.6-0.7 0 1
code_succ0_tp1.8 all correct i_code_pi650 0 2
ii_results_0.6-0.7 0 1
b attempt2 all correct i_code_pi650 1 4
ii_results_0.6-0.7 0 4
code_succ0_tp1.8 all correct i_code_pi650 2 4
ii_results_0.6-0.7 4 4
a2, gp null attempt2 all correct i_code_pi650 1 4
ii_results_0.6-0.7 0 4
code_succ0_tp1.8 all correct i_code_pi650 0 4
ii_results_0.6-0.7 0 3
Result: not reproduced, resting on an assumption. Under the main reading, with the Methods' error trials, correct-trial PI exceeds error-trial PI in all four groups (0.203 to 0.319 against 0.121 to 0.225 under the code's neuron selection), but p < 0.0001 in only 1 of 4 groups under either neuron selection (code: A -1, t = 4.80, n = 18, while A 1, M -1 and M 1 have p = 0.0018, 0.0119 and 0.00010; Results: M 1, t = 5.05, n = 26). With the correct trials subsampled to the error-trial count, no group reaches it (p 2.0e-04 to 0.14). The other readings differ: with the code's error trials, 2 of 4 groups under the code's selection, and 4 of 4 under rule b with the Results' selection; under rule a1 too few units are periodic (2 to 7 per group in monkey A, none in M). The groups hold 10 to 26 neurons under the main reading against the paper's 66 to 104.
ed4d-overshoot-longer-period¶
For trials of the same desired distance, the periodicity of EC neurons is significantly longer (lower frequency) on overshoot trials (top tertile of produced vector) than on undershoot trials (bottom tertile).
Periodic units (each periodic set: the main reading and the a1, b and GP-null readings), per direction in which they are periodic; correct trials of distance 3 (with the PI rules); per session and direction, trials sorted by |tp| - |ta|, bottom and top tertiles; periodicity (Methods lags) on each; units with fewer than 15 kept trials in a tertile dropped, and periodicities of 0 or above 1 s dropped. Rank-sum test of overshoot against undershoot across units, per monkey; the code's variant lags (100 to 1,090 ms) are also shown.
def tertile_mask(which):
def f(tr):
m = (tr.trial_type == 3) & (tr.attempt == 1) & (tr.distance == 3) & (tr.tp.abs() > 1) & (tr.seqq < 3)
out = pd.Series(False, index=tr.index)
for d in (1, -1):
g = tr[m & (tr.direction == d)]
err = g.tp.abs() - g.ta.abs()
lo, hi = err.quantile(1 / 3), err.quantile(2 / 3)
out[g.index[err <= lo] if which == "under" else g.index[err >= hi]] = True
return out
return f
lagsets = {"methods": mnav.LAGS_METHODS, "code": np.arange(100, 1091, 10)}
t = time.time()
under = mnav.periodicity_table(sessions, extra=tertile_mask("under"), lag_sets=lagsets)
over = mnav.periodicity_table(sessions, extra=tertile_mask("over"), lag_sets=lagsets)
print(f"tertile PI in {time.time() - t:.0f} s")
rows = []
for (set_name, set_col), (lag, lags) in [(sc, ll) for sc in SETS.items() for ll in lagsets.items()]:
mrg = under.merge(over, on=["monkey", "session", "unit", "direction"], suffixes=("_u", "_o"))
mrg = mrg.merge(base[["monkey", "session", "unit", "direction", set_col]], on=["monkey", "session", "unit", "direction"])
mrg = mrg[mrg[set_col]]
n_periodic = mrg.groupby("monkey").size()
ok = (mrg.k_u >= 15) & (mrg.k_o >= 15) & mrg[f"pi_{lag}_u"].notna() & mrg[f"pi_{lag}_o"].notna()
mrg = mrg[ok].copy()
mrg["per_u"] = [lags[int(np.argmax(p))] / 1000 for p in mrg[f"pi_{lag}_u"]]
mrg["per_o"] = [lags[int(np.argmax(p))] / 1000 for p in mrg[f"pi_{lag}_o"]]
keep = (mrg.per_u > 0) & (mrg.per_u <= 1) & (mrg.per_o > 0) & (mrg.per_o <= 1)
for m in ("A", "M"):
g = mrg[keep & (mrg.monkey == m)]
if len(g) < 2:
rows.append({"periodic_set": set_name, "lags": lag, "monkey": m, "n": len(g)})
continue
rs = stats.ranksums(g.per_o, g.per_u)
rows.append({"periodic_set": set_name, "lags": lag, "monkey": m, "periodic_unit_dirs": int(n_periodic.get(m, 0)),
"kept_>=15_trials": int((mrg.monkey == m).sum()), "n": len(g),
"median_under_s": g.per_u.median(), "median_over_s": g.per_o.median(),
"mean_diff_over-under_s": round((g.per_o - g.per_u).mean(), 4),
"ranksum_Z(over-under)": round(rs.statistic, 2), "p": rs.pvalue,
"meets(over>under,p<0.05)": bool(rs.pvalue < 0.05 and g.per_o.median() >= g.per_u.median() and rs.statistic > 0)})
ed4d = pd.DataFrame(rows)
print(ed4d.to_string(index=False))
print("Trials per tertile per session and direction (median):",
int(np.median([len(d) for d in under.drop_duplicates(['session', 'direction']).durations])))
tertile PI in 31 s
periodic_set lags monkey periodic_unit_dirs kept_>=15_trials n median_under_s median_over_s mean_diff_over-under_s ranksum_Z(over-under) p meets(over>under,p<0.05)
a2 (main) methods A 68.0 45.0 41 0.570 0.590 0.0207 0.78 4.332431e-01 False
a2 (main) methods M 64.0 34.0 29 0.580 0.710 0.1131 3.83 1.263869e-04 True
a2 (main) code A 68.0 45.0 41 0.570 0.590 0.0207 0.78 4.332431e-01 False
a2 (main) code M 64.0 34.0 30 0.580 0.710 0.1147 3.87 1.072748e-04 True
a1 methods A 16.0 14.0 14 0.535 0.545 0.0107 0.44 6.624719e-01 False
a1 methods M NaN NaN 0 NaN NaN NaN NaN NaN NaN
a1 code A 16.0 14.0 14 0.535 0.545 0.0107 0.44 6.624719e-01 False
a1 code M NaN NaN 0 NaN NaN NaN NaN NaN NaN
b methods A 370.0 253.0 179 0.630 0.650 0.0260 1.45 1.482408e-01 False
b methods M 618.0 274.0 143 0.680 0.770 0.0948 5.10 3.373963e-07 True
b code A 370.0 253.0 188 0.640 0.650 0.0282 1.53 1.254864e-01 False
b code M 618.0 274.0 157 0.680 0.770 0.0906 5.12 3.125705e-07 True
a2, gp null methods A 41.0 27.0 26 0.540 0.560 0.0315 1.16 2.451845e-01 False
a2, gp null methods M NaN NaN 1 NaN NaN NaN NaN NaN NaN
a2, gp null code A 41.0 27.0 26 0.540 0.560 0.0315 1.16 2.451845e-01 False
a2, gp null code M NaN NaN 1 NaN NaN NaN NaN NaN NaN
Trials per tertile per session and direction (median): 26
Result: partially reproduced, resting on an assumption. Under the main reading, in monkey M the overshoot tertile has a longer periodicity (median 0.71 against 0.58 s, rank-sum Z = 3.83, p = 1.3e-04, n = 29 unit-directions); in monkey A the difference is in the same direction but not significant (0.590 against 0.570 s, Z = 0.78, p = 0.43, n = 41). The code's lag range gives the same (M p = 1.1e-04, n = 30). Rule b gives the same split (M p = 3.4e-07, A p = 0.15); rule a1 has no periodic units in monkey M, and the GP-null reading leaves 1 in M; neither reaches p < 0.05 in A (p = 0.66 and 0.25). Of the periodic unit-directions of the main reading, 45 of 68 (A) and 34 of 64 (M) had at least 15 trials in both tertiles.
ed10-motor-preparation-signal¶
Most EC neurons show a motor-related response before joystick offset, 77% (477/614) in monkey A and 70% (607/864) in monkey M, including most non-periodic neurons (81%, 368/452, in A and 89%, 432/481, in M).
Criterion (claims, from the code): direction 1, first-attempt MNAV trials of distance at least 3 with |tp| > 1.5 s, at least 6 trials; spike-count rate in tp - 1 to tp - 0.6 s and in tp - 0.4 to tp - 0.1 s after joystick onset; 200 bootstrap means of each; two-sample t-test at p < 0.0001 (two-sided; increase only as a variant). Non-periodic: not periodic (each periodic set) with more than 15 PI trials. As for fig2, the dandiset holds a subset of the sessions, so the counts are compared as proportions of the units it holds; the raw counts are printed. Variants: either direction, and a paired t-test across the trials themselves (without the bootstrap, whose t-test on 200 resampled means has a p-value that shrinks with the number of resamples).
rng10 = np.random.default_rng(10)
rows = []
for s in sessions:
tr = s["trials"]
for d in (1, -1):
g = tr[(tr.trial_type == 3) & (tr.attempt == 1) & (tr.distance >= 3) & (tr.tp.abs() > 1.5) & (tr.seqq < 3) & (tr.direction == d)]
on, tp = g.joy1onttl.values, g.tp.abs().values
for u, sp in enumerate(s["spikes"]):
row = {"monkey": s["monkey"], "session": s["session"], "unit": u, "direction": d, "n": len(g)}
if len(g) >= 6:
early = mnav.mean_rate(sp, on + tp - 1.0, on + tp - 0.6)
late = mnav.mean_rate(sp, on + tp - 0.4, on + tp - 0.1)
idx = rng10.integers(0, len(g), size=(200, len(g)))
be, bl = early[idx].mean(1), late[idx].mean(1)
p = stats.ttest_ind(bl, be).pvalue if (be.std() + bl.std()) > 0 else 1.0
# variant without the bootstrap: paired t-test across the trials themselves
p_tr = stats.ttest_rel(late, early).pvalue if np.std(late - early) > 0 else 1.0
row.update(p=p, p_trials=p_tr, increase=bl.mean() > be.mean())
rows.append(row)
mot = pd.DataFrame(rows)
mot["motor"] = mot.p.fillna(1) < 1e-4
mot["motor_increase"] = mot.motor & mot.increase.fillna(False).astype(bool)
mot["motor_trial_ttest"] = mot.p_trials.fillna(1) < 1e-4
per_dir = mot.merge(base[["monkey", "session", "unit", "direction", "k"] + list(dict.fromkeys(SETS.values()))],
on=["monkey", "session", "unit", "direction"])
unit_k = base.groupby(["monkey", "session", "unit"]).k.max()
paper10 = {"A": dict(motor=(477, 614), nonper_motor=(368, 452), nonper=(452, 614)),
"M": dict(motor=(607, 864), nonper_motor=(432, 481), nonper=(481, 864))}
rows = []
for set_name, set_col, variant, sel, col in (
[(sn, sc, "direction 1 (code)", per_dir.direction == 1, "motor") for sn, sc in SETS.items()]
+ [("a2 (main)", main_col, v, sl, c) for v, sl in (("direction 1 (code)", per_dir.direction == 1), ("either direction", per_dir.direction != 0))
for c in ("motor", "motor_increase", "motor_trial_ttest") if (v, c) != ("direction 1 (code)", "motor")]):
unit_periodic = base.groupby(["monkey", "session", "unit"])[set_col].any()
u = per_dir[sel].groupby(["monkey", "session", "unit"])[col].any()
df = pd.DataFrame({"motor": u, "periodic": unit_periodic, "k": unit_k}).reset_index()
df["nonper"] = (~df.periodic) & (df.k > mnav.MIN_TRIALS)
for m in ("A", "M"):
g = df[df.monkey == m]
n_all, n_mot = len(g), int(g.motor.sum())
n_np, n_np_mot = int(g.nonper.sum()), int((g.nonper & g.motor).sum())
ref = paper10[m]
rows.append({"periodic_set": set_name, "variant": variant, "test": col, "monkey": m, "units": n_all, "motor": n_mot,
"motor_pct": round(100 * n_mot / n_all, 1), "paper_pct": round(100 * ref["motor"][0] / ref["motor"][1], 1),
"nonperiodic": n_np, "nonper_pct_of_units": round(100 * n_np / n_all, 1),
"paper_nonper_pct": round(100 * ref["nonper"][0] / ref["nonper"][1], 1),
"nonper_motor": n_np_mot, "nonper_motor_pct": round(100 * n_np_mot / max(n_np, 1), 1),
"paper_nonper_motor_pct": round(100 * ref["nonper_motor"][0] / ref["nonper_motor"][1], 1)})
ed10 = pd.DataFrame(rows)
print(ed10.to_string(index=False))
main10 = ed10[(ed10.periodic_set == "a2 (main)") & (ed10.variant == "direction 1 (code)") & (ed10.test == "motor")]
for r in main10.itertuples():
within = [abs(r.motor_pct - r.paper_pct) <= 0.1 * r.paper_pct,
abs(r.nonper_pct_of_units - r.paper_nonper_pct) <= 0.1 * r.paper_nonper_pct,
abs(r.nonper_motor_pct - r.paper_nonper_motor_pct) <= 0.1 * r.paper_nonper_motor_pct]
print(f"monkey {r.monkey}: within 10% of the paper's proportion (motor, non-periodic, motor among non-periodic): {within}; "
f"motor and motor-among-non-periodic above 50%: {r.motor_pct > 50}, {r.nonper_motor_pct > 50}")
print("\nMotor units (direction 1, two-sided) per session (%):")
print((per_dir[per_dir.direction == 1].groupby(["monkey", "session"]).motor.mean() * 100).round(1).to_string())
periodic_set variant test monkey units motor motor_pct paper_pct nonperiodic nonper_pct_of_units paper_nonper_pct nonper_motor nonper_motor_pct paper_nonper_motor_pct
a2 (main) direction 1 (code) motor A 364 275 75.5 77.7 240 65.9 73.6 203 84.6 81.4
a2 (main) direction 1 (code) motor M 729 495 67.9 70.3 385 52.8 55.7 322 83.6 89.8
a1 direction 1 (code) motor A 364 275 75.5 77.7 286 78.6 73.6 243 85.0 81.4
a1 direction 1 (code) motor M 729 495 67.9 70.3 442 60.6 55.7 373 84.4 89.8
b direction 1 (code) motor A 364 275 75.5 77.7 57 15.7 73.6 46 80.7 81.4
b direction 1 (code) motor M 729 495 67.9 70.3 63 8.6 55.7 44 69.8 89.8
a2, gp null direction 1 (code) motor A 364 275 75.5 77.7 264 72.5 73.6 222 84.1 81.4
a2, gp null direction 1 (code) motor M 729 495 67.9 70.3 434 59.5 55.7 366 84.3 89.8
a2 (main) direction 1 (code) motor_increase A 364 120 33.0 77.7 240 65.9 73.6 86 35.8 81.4
a2 (main) direction 1 (code) motor_increase M 729 241 33.1 70.3 385 52.8 55.7 171 44.4 89.8
a2 (main) direction 1 (code) motor_trial_ttest A 364 21 5.8 77.7 240 65.9 73.6 12 5.0 81.4
a2 (main) direction 1 (code) motor_trial_ttest M 729 50 6.9 70.3 385 52.8 55.7 43 11.2 89.8
a2 (main) either direction motor A 364 314 86.3 77.7 240 65.9 73.6 229 95.4 81.4
a2 (main) either direction motor M 729 600 82.3 70.3 385 52.8 55.7 374 97.1 89.8
a2 (main) either direction motor_increase A 364 188 51.6 77.7 240 65.9 73.6 142 59.2 81.4
a2 (main) either direction motor_increase M 729 340 46.6 70.3 385 52.8 55.7 227 59.0 89.8
a2 (main) either direction motor_trial_ttest A 364 31 8.5 77.7 240 65.9 73.6 16 6.7 81.4
a2 (main) either direction motor_trial_ttest M 729 85 11.7 70.3 385 52.8 55.7 76 19.7 89.8
monkey A: within 10% of the paper's proportion (motor, non-periodic, motor among non-periodic): [True, False, True]; motor and motor-among-non-periodic above 50%: True, True
monkey M: within 10% of the paper's proportion (motor, non-periodic, motor among non-periodic): [True, True, True]; motor and motor-among-non-periodic above 50%: True, True
Motor units (direction 1, two-sided) per session (%):
monkey session
A 08152019 76.3
08212019 68.4
08232019 69.8
08272019 78.3
08292019 80.0
09172019 76.6
M 04122021 71.4
04132021 64.5
04142021 47.3
04212021 84.0
04222021 68.3
04252021 71.9
05042021 67.2
05062021 61.2
05072021 94.4
Result: partially reproduced, resting on an assumption. Under the code's test (two-sided, direction 1), 75.5% of the monkey A units (275/364) and 67.9% of the monkey M units (495/729) are motor-related, within 10% of the paper's 77.7% and 70.3%. Under the main periodic reading, 65.9% (A, 240/364) and 52.8% (M, 385/729) of the units are non-periodic, against 73.6% and 55.7%: M is within 10% and A just outside (the bound is 66.2%). Among them 84.6% (203/240) and 83.6% (322/385) are motor-related, within 10% of 81.4% and 89.8%. So five of the six values are within tolerance, and the four proportions are each above 50%. Under the a1 and GP-null readings all six are within tolerance (non-periodic 78.6% and 72.5% in A, 60.6% and 59.5% in M); under b too few units are non-periodic (15.7% and 8.6%). The test says little about a motor signal: counting increases only gives 33.0% and 33.1%, and a paired t-test across the trials themselves at p < 0.0001 finds 5.8% (A) and 6.9% (M) of the units; in M 04142021 and A 09172019, where no clock check succeeded, the bootstrap test still finds 47.3% and 76.6%.
Ramping and distance coding: shared selection¶
Correct (first-attempt) MNAV trials with seqq < 3 (distance from ta in 0.65 s steps).
Units with at least 10 such trials per distance (1 to 5) in both directions (claims, fig3b and
fig3d, from the code); direction 1 (ta > 0) is analysed (code line 93).
def coding_trials(tr):
return tr[(tr.trial_type == 3) & (tr.attempt == 1) & (tr.seqq < 3) & (tr.tp.abs() > 0.01) & (tr.tp < 99)]
sel_rows = []
for s in sessions:
g = coding_trials(s["trials"])
cnt = g.groupby(["direction", "distance"]).size()
ok = all(cnt.get((d, k), 0) >= 10 for d in (1, -1) for k in range(1, 6))
sel_rows.append({"monkey": s["monkey"], "session": s["session"], "min_trials_per_dir_dist": int(cnt.min()), "eligible": ok,
"units": len(s["spikes"])})
elig = pd.DataFrame(sel_rows)
print(elig.to_string(index=False))
print("Eligible units per monkey:", elig[elig.eligible].groupby("monkey").units.sum().to_dict())
monkey session min_trials_per_dir_dist eligible units
A 08152019 25 True 59
A 08212019 30 True 57
A 08232019 32 True 43
A 08272019 32 True 46
A 08292019 32 True 95
A 09172019 19 True 64
M 04122021 18 True 98
M 04132021 5 False 93
M 04142021 4 False 112
M 04212021 17 True 81
M 04222021 2 False 60
M 04252021 13 True 64
M 05042021 17 True 64
M 05062021 14 True 103
M 05072021 14 True 54
Eligible units per monkey: {'A': 364, 'M': 464}
fig3b-onset-and-offset-coding¶
Single EC neurons encode temporal distance in their firing rate before joystick onset, before joystick offset, or, in a small number, at both.
Cross-validated R2 (Methods, Cross-validated R2; claims): mean rate in the window, each distance subsampled to the smallest count, random halves, linear fit of rate on distance 1 to 5 on the training half scored on the held-out half, 50 splits averaged; positive coding is mean R2 > 0; units with R2 below -1 dropped. Windows: (a) 100 ms before onset and offset, (b) 200 ms before. Null: the same with the distances shuffled.
rng3b = np.random.default_rng(3)
def cv_r2(y, x, rng, n_split=50, shuffle=False):
out = []
for _ in range(n_split):
xx = rng.permutation(x) if shuffle else x
perm = rng.permutation(len(y))
tr_i, te_i = perm[: len(y) // 2], perm[len(y) // 2:]
b = np.polyfit(xx[tr_i], y[tr_i], 1)
pred = np.polyval(b, xx[te_i])
yt = y[te_i]
ss = np.sum((yt - yt.mean()) ** 2)
out.append(1 - np.sum((yt - pred) ** 2) / ss if ss > 0 else np.nan)
return np.nanmean(out) if np.any(np.isfinite(out)) else np.nan
rows = []
for s in sessions:
if not elig.set_index("session").eligible[s["session"]]:
continue
g = coding_trials(s["trials"])
g = g[g.direction == 1]
nmin = g.groupby("distance").size().loc[1:5].min()
g = pd.concat([gg.sample(nmin, random_state=int(rng3b.integers(1e9))) for k, gg in g.groupby("distance") if 1 <= k <= 5])
x = g.distance.values.astype(float)
on, off = g.joy1onttl.values, g.joy1offttl.values
for u, sp in enumerate(s["spikes"]):
row = {"monkey": s["monkey"], "session": s["session"], "unit": u}
for w in (0.1, 0.2):
for ev, t in (("onset", on), ("offset", off)):
y = mnav.mean_rate(sp, t - w, t)
row[f"r2_{ev}_{int(w * 1000)}"] = cv_r2(y, x, rng3b)
row[f"r2shuf_{ev}_{int(w * 1000)}"] = cv_r2(y, x, rng3b, shuffle=True)
rows.append(row)
r2 = pd.DataFrame(rows)
res = []
for w in (100, 200):
for m, g in r2.groupby("monkey"):
on, off = g[f"r2_onset_{w}"], g[f"r2_offset_{w}"]
valid = (on >= -1) & (off >= -1)
on, off = on[valid], off[valid]
son, soff = g[f"r2shuf_onset_{w}"][valid], g[f"r2shuf_offset_{w}"][valid]
p_on_only = float(((on > 0) & (off <= 0)).mean())
p_off_only = float(((off > 0) & (on <= 0)).mean())
p_both = float(((on > 0) & (off > 0)).mean())
p_shuf = float(((son > 0) | (soff > 0)).mean())
p_shuf_on, p_shuf_off = float((son > 0).mean()), float((soff > 0).mean())
res.append({"window_ms": w, "monkey": m, "units": int(valid.sum()), "dropped_r2<-1": int((~valid).sum()),
"onset_only": round(p_on_only, 3), "offset_only": round(p_off_only, 3), "both": round(p_both, 3),
"shuffled_onset>0": round(p_shuf_on, 3), "shuffled_offset>0": round(p_shuf_off, 3),
"meets": bool(p_on_only > p_shuf_on and p_off_only > p_shuf_off and p_both < p_on_only and p_both < p_off_only)})
fig3b = pd.DataFrame(res)
print(fig3b.to_string(index=False))
print("\nShuffled-regressor proportions are compared with the coding proportions at the same epoch (onset only against shuffled onset, offset only against shuffled offset).")
window_ms monkey units dropped_r2<-1 onset_only offset_only both shuffled_onset>0 shuffled_offset>0 meets
100 A 288 76 0.007 0.083 0.000 0.0 0.0 True
100 M 263 201 0.080 0.023 0.011 0.0 0.0 True
200 A 301 63 0.010 0.060 0.000 0.0 0.0 True
200 M 279 185 0.082 0.036 0.029 0.0 0.0 True
Shuffled-regressor proportions are compared with the coding proportions at the same epoch (onset only against shuffled onset, offset only against shuffled offset).
Result: reproduced. In both monkeys and both windows, the proportions of units coding distance at onset only (A 0.7% and 1.0%, M 8.0% and 8.2% for 100 and 200 ms) and at offset only (A 8.3% and 6.0%, M 2.3% and 3.6%) exceed the proportion with positive cross-validated R2 under the shuffled regressor (0 in every case), and fewer units code at both (A 0%, M 1.1% and 2.9%) than at either alone. The shuffled proportions are 0 because R2 averaged over 50 held-out halves is always negative for a shuffled regressor, so this comparison is weak; the coding proportions themselves are small. Units with R2 below -1 were dropped (A 76 and 63, M 201 and 185). Monkey M sessions with fewer than 10 trials per distance in a direction (3 of 9) are not used.
fig3d-offset-coding-stronger¶
At the population level distance is encoded weakly before joystick onset and robustly and linearly before joystick offset; offset R2 exceeds onset R2 in both monkeys.
Targeted dimensionality reduction (Methods; claims): pseudo-population of each monkey's eligible units, direction 1. Per resample, each unit's trials of each distance are split at random into halves. Rates in 20 ms bins (a standard choice for the time resolution of the regression, which the Methods do not give) in the windows -200 to 0 ms before onset and -500 to -300 ms before offset (Methods); per bin, the regression r = beta * dist across training trials; the axis is the coefficient vector at the bin of largest norm; held-out trials of each distance are averaged, averaged over the window and projected on the axis; R2 of a linear fit of the projections on distance; 50 resamples. Two-sided rank-sum test of offset against onset R2.
Centring, two readings. The Methods centre each unit "by subtracting its mean response across the two windows of interest" and regress without an intercept; a difference between the unit's mean rates in the two windows then loads onto beta. The method the Methods cite for the regression (ref. 63, Mante et al. 2013) has a constant term at each time point, which is equivalent to centring each unit in each window and time bin. Both are run ("windows", the main reading, as the Methods state it, and "bins", a variant). The bin of largest norm is found in each window separately, which the Methods leave open. The criterion's variant with the code's windows (-200 to -100 ms and -400 to -300 ms on time-warped bins) is not run. Other variants: the coefficients projected onto the principal components holding 80% of the training variance (the code's denoising), and distances 2 to 5 only (the Methods' control for distance 1).
rng3d = np.random.default_rng(33)
BIN = 0.02
win = {"onset": (-0.2, 0.0), "offset": (-0.5, -0.3)}
pop = {"A": [], "M": []}
for s in sessions:
if not elig.set_index("session").eligible[s["session"]]:
continue
g = coding_trials(s["trials"])
g = g[(g.direction == 1) & g.distance.between(1, 5)]
for u, sp in enumerate(s["spikes"]):
unit = {}
for ev, col in (("onset", "joy1onttl"), ("offset", "joy1offttl")):
a, b = win[ev]
edges = np.arange(a, b + 1e-9, BIN)
t = g[col].values
unit[ev] = np.stack([mnav.mean_rate(sp, t + e0, t + e0 + BIN) for e0 in edges[:-1]], axis=1) # trials x bins
unit["dist"] = g.distance.values.astype(int)
pop[s["monkey"]].append(unit)
def tdr_once(units, rng, denoise=False, centring="windows", dists=(1, 2, 3, 4, 5)):
"""One resample. centring "windows": each unit's training mean across both windows (Methods,
TDR: "subtracting its mean response across the two windows"); "bins": its training mean in
each window and time bin (variant). dists: the distances used (the Methods' control leaves out 1)."""
dists = list(dists)
train = {"onset": [], "offset": []}
test = {"onset": [], "offset": []}
for unit in units:
tr_idx, te_idx = [], []
for k in dists:
idx = rng.permutation(np.nonzero(unit["dist"] == k)[0])
tr_idx.append(idx[: len(idx) // 2])
te_idx.append(idx[len(idx) // 2:])
tr_all = np.concatenate(tr_idx)
mu = np.mean([unit[ev][tr_all].mean() for ev in ("onset", "offset")])
for ev in ("onset", "offset"):
X = unit[ev] - (mu if centring == "windows" else unit[ev][tr_all].mean(0))
train[ev].append([X[i] for i in tr_idx])
test[ev].append(np.array([X[i].mean() for i in te_idx])) # per distance, averaged over bins and trials
out = {}
for ev in ("onset", "offset"):
nb = train[ev][0][0].shape[1]
beta = np.zeros((len(units), nb))
cond = np.zeros((len(units), len(dists) * nb))
for i, per_d in enumerate(train[ev]):
num = sum(k * per_d[j].sum(0) for j, k in enumerate(dists))
den = sum(k ** 2 * len(per_d[j]) for j, k in enumerate(dists))
beta[i] = num / den
cond[i] = np.concatenate([per_d[j].mean(0) for j in range(len(dists))])
if denoise:
c = cond - cond.mean(1, keepdims=True)
U, S, _ = np.linalg.svd(c, full_matrices=False)
npc = int(np.searchsorted(np.cumsum(S ** 2) / np.sum(S ** 2), 0.8) + 1)
beta = U[:, :npc] @ (U[:, :npc].T @ beta)
tstar = int(np.argmax(np.linalg.norm(beta, axis=0)))
axis = beta[:, tstar] / np.linalg.norm(beta[:, tstar])
proj = np.array(test[ev]).T @ axis # one value per distance
fit = stats.linregress(dists, proj)
out[ev] = fit.rvalue ** 2
return out
rows = []
for centring, dists, denoise in [("windows", (1, 2, 3, 4, 5), False), ("windows", (1, 2, 3, 4, 5), True),
("bins", (1, 2, 3, 4, 5), False), ("bins", (1, 2, 3, 4, 5), True),
("windows", (2, 3, 4, 5), False), ("bins", (2, 3, 4, 5), False)]:
for m in ("A", "M"):
if not pop[m]:
continue
rs = [tdr_once(pop[m], rng3d, denoise, centring, dists) for _ in range(50)]
on = np.array([r["onset"] for r in rs])
off = np.array([r["offset"] for r in rs])
test = stats.ranksums(on, off)
rows.append({"centring": centring, "distances": "1-5" if 1 in dists else "2-5", "denoise_80pct_PCs": denoise,
"monkey": m, "units": len(pop[m]),
"onset_R2_median": round(float(np.median(on)), 3), "offset_R2_median": round(float(np.median(off)), 3),
"onset_R2_iqr": f"{np.percentile(on, 25):.2f}-{np.percentile(on, 75):.2f}",
"offset_R2_iqr": f"{np.percentile(off, 25):.2f}-{np.percentile(off, 75):.2f}",
"ranksum_Z(onset-offset)": round(test.statistic, 2), "p": test.pvalue,
"meets(offset>onset,p<1e-4)": bool(np.median(off) > np.median(on) and test.pvalue < 1e-4)})
fig3d = pd.DataFrame(rows)
print(fig3d.to_string(index=False))
diff = {m: np.mean([np.mean([u["offset"].mean() - u["onset"].mean()]) for u in pop[m]]) for m in pop if pop[m]}
print("\nMean over units of the offset-window rate minus the onset-window rate (spikes/s):", {m: round(v, 2) for m, v in diff.items()})
centring distances denoise_80pct_PCs monkey units onset_R2_median offset_R2_median onset_R2_iqr offset_R2_iqr ranksum_Z(onset-offset) p meets(offset>onset,p<1e-4)
windows 1-5 False A 364 0.709 0.854 0.55-0.82 0.82-0.88 -5.14 2.808032e-07 True
windows 1-5 False M 464 0.787 0.684 0.75-0.83 0.64-0.73 6.85 7.259223e-12 False
windows 1-5 True A 364 0.568 0.853 0.39-0.73 0.83-0.87 -7.53 5.150316e-14 True
windows 1-5 True M 464 0.790 0.676 0.75-0.82 0.62-0.73 6.53 6.645493e-11 False
bins 1-5 False A 364 0.502 0.853 0.24-0.79 0.83-0.87 -6.03 1.618502e-09 True
bins 1-5 False M 464 0.808 0.935 0.77-0.86 0.91-0.97 -6.90 5.174192e-12 True
bins 1-5 True A 364 0.530 0.861 0.23-0.81 0.81-0.88 -5.48 4.239175e-08 True
bins 1-5 True M 464 0.798 0.944 0.74-0.84 0.90-0.96 -7.56 3.953331e-14 True
windows 2-5 False A 364 0.839 0.833 0.58-0.94 0.79-0.88 -0.12 9.067056e-01 False
windows 2-5 False M 464 0.711 0.928 0.52-0.83 0.89-0.94 -6.91 4.928954e-12 True
bins 2-5 False A 364 0.377 0.871 0.12-0.68 0.82-0.94 -6.68 2.387606e-11 True
bins 2-5 False M 464 0.562 0.832 0.48-0.69 0.78-0.90 -7.56 4.168482e-14 True
Mean over units of the offset-window rate minus the onset-window rate (spikes/s): {'A': np.float64(-0.81), 'M': np.float64(1.36)}
Result: partially reproduced. With each unit centred by its mean across the two windows, as the Methods state ("We first centred the responses of each neuron by subtracting its mean response across the two windows of interest before joystick onset and offset"), offset R2 exceeds onset R2 at p < 0.0001 in monkey A (p = 2.8e-07, denoised 5.2e-14) and monkey M reverses (onset 0.787 against offset 0.684, Z = 6.85, p = 7.3e-12; denoised 6.6e-11). With the distances 1 to 5 of the criterion this holds in one monkey, which the criterion calls partially reproduced. Without distance 1, the Methods' control and outside the criterion, A no longer differs (p = 0.91) and M gives offset above onset (p = 4.9e-12).
With the variant centring per window and time bin (the cited method's constant term), offset R2 exceeds onset R2 in both monkeys: A medians 0.853 against 0.502 (Z = -6.03, p = 1.6e-09), M 0.935 against 0.808 (Z = -6.90, p = 5.2e-12); with denoising p = 4.2e-08 and 4.0e-14, and without distance 1 p = 2.4e-11 and 4.2e-14. The mean rate in the offset window is 1.36 spikes/s above the onset window in M and 0.81 below it in A, the difference that the Methods' centring loads onto the coefficients. The bin width of the regression (20 ms) is a standard choice.
ed4fg-distance-over-landmark and ed4e-rate-codes-distance-over-slope: the ANOVAs¶
Four-way ANOVA (distance, start landmark, target landmark, direction; Methods), per unit with at least 10 correct trials in each of the 30 start and target conditions (claims, from the code's landmark ANOVA), on the mean rate (spike count over the window) and on the slope of a line fitted to the smoothed rate in the window. Type II sums of squares (each factor adjusted for the others; the Methods do not say which, a standard choice). Units with zero variance in a window are skipped.
def anova_F(df, y):
model = smf.ols(f"{y} ~ C(distance) + C(curr) + C(target) + C(direction)", data=df).fit()
a = sm.stats.anova_lm(model, typ=2)
return {f: (a.loc[f"C({f})", "F"], a.loc[f"C({f})", "PR(>F)"]) for f in ("distance", "curr", "target", "direction")}
def slope_of(sp, t0, t1):
rates = mnav.smoothed_rates(sp, t0, t1)
return np.array([np.polyfit(np.arange(len(r)) / 1000, r, 1)[0] if len(r) > 2 else np.nan for r in rates])
WINDOWS = {"rate_-500_-300": ("rate", -0.5, -0.3), "rate_-200_0": ("rate", -0.2, 0.0), "rate_-400_-300": ("rate", -0.4, -0.3),
"rate_-550_-250": ("rate", -0.55, -0.25), "slope_-550_-250": ("slope", -0.55, -0.25), "slope_-300_0": ("slope", -0.3, 0.0)}
t = time.time()
rows = []
anova_sessions = []
for s in sessions:
g = coding_trials(s["trials"])
cnt = g.groupby(["curr", "target"]).size()
if len(cnt) < 30 or cnt.min() < 10:
anova_sessions.append({"session": s["session"], "monkey": s["monkey"], "conditions": len(cnt), "min_trials": int(cnt.min()), "used": False})
continue
anova_sessions.append({"session": s["session"], "monkey": s["monkey"], "conditions": len(cnt), "min_trials": int(cnt.min()), "used": True})
off = g.joy1offttl.values
for u, sp in enumerate(s["spikes"]):
row = {"monkey": s["monkey"], "session": s["session"], "unit": u}
df = g[["distance", "curr", "target", "direction"]].copy()
for name, (kind, a, b) in WINDOWS.items():
y = mnav.mean_rate(sp, off + a, off + b) if kind == "rate" else slope_of(sp, off + a, off + b)
if np.nanstd(y) == 0:
continue
df["y"] = y
res_ = anova_F(df.dropna(), "y")
for f, (F, p) in res_.items():
row[f"{name}_{f}_F"] = F
row[f"{name}_{f}_p"] = p
rows.append(row)
anova = pd.DataFrame(rows)
print(pd.DataFrame(anova_sessions).to_string(index=False))
print(f"ANOVAs for {len(anova)} units in {time.time() - t:.0f} s; units per monkey:", anova.groupby("monkey").size().to_dict())
session monkey conditions min_trials used
08152019 A 30 24 True
08212019 A 30 30 True
08232019 A 30 31 True
08272019 A 30 31 True
08292019 A 30 32 True
09172019 A 30 19 True
04122021 M 30 13 True
04132021 M 30 5 False
04142021 M 30 4 False
04212021 M 30 11 True
04222021 M 30 1 False
04252021 M 30 13 True
05042021 M 30 16 True
05062021 M 30 14 True
05072021 M 30 4 False
ANOVAs for 774 units in 274 s; units per monkey: {'A': 364, 'M': 410}
ed4fg-distance-over-landmark¶
Significantly less variance of firing rate near joystick offset is explained by the start or the target landmark than by distance.
Units with a significant effect (p < 0.05) of distance or of the landmark; distance F against landmark F, by (a) a t-test on the differences (Methods; paired) and (b) a two-sample t-test (legend). Windows (a) -500 to -300 ms (Methods), (b) the last 200 ms (legend); the code's -400 to -300 ms is a variant.
rows = []
for wname in ("rate_-500_-300", "rate_-200_0", "rate_-400_-300"):
for m, g in anova.groupby("monkey"):
for lm in ("curr", "target"):
cols = [f"{wname}_distance_F", f"{wname}_{lm}_F", f"{wname}_distance_p", f"{wname}_{lm}_p"]
gg = g.dropna(subset=cols)
sig = gg[(gg[cols[2]] < 0.05) | (gg[cols[3]] < 0.05)]
fd, fl = sig[cols[0]], sig[cols[1]]
tp_ = stats.ttest_1samp(fd - fl, 0) if len(sig) > 1 else None
t2 = stats.ttest_ind(fd, fl) if len(sig) > 1 else None
rows.append({"window": wname, "monkey": m, "landmark": "start" if lm == "curr" else "target", "units_tested": len(gg),
"units_sig_either": len(sig), "mean_F_distance": round(fd.mean(), 2), "mean_F_landmark": round(fl.mean(), 2),
"paired_t": None if tp_ is None else round(tp_.statistic, 2), "paired_p": None if tp_ is None else tp_.pvalue,
"two_sample_t": None if t2 is None else round(t2.statistic, 2), "two_sample_p": None if t2 is None else t2.pvalue,
"meets_paired": bool(tp_ is not None and tp_.pvalue < 1e-4 and fd.mean() > fl.mean()),
"meets_two_sample": bool(t2 is not None and t2.pvalue < 1e-4 and fd.mean() > fl.mean())})
ed4fg = pd.DataFrame(rows)
print(ed4fg.to_string(index=False))
main_fg = ed4fg[ed4fg.window != "rate_-400_-300"]
print("\nMain readings (windows a, b; tests a, b): comparisons meeting p < 1e-4 with distance F larger:",
int(main_fg.meets_paired.sum() + main_fg.meets_two_sample.sum()), "of", 2 * len(main_fg))
window monkey landmark units_tested units_sig_either mean_F_distance mean_F_landmark paired_t paired_p two_sample_t two_sample_p meets_paired meets_two_sample rate_-500_-300 A start 324 70 6.12 2.55 4.04 1.391301e-04 3.68 3.384088e-04 False False rate_-500_-300 A target 324 69 6.29 3.40 3.35 1.311916e-03 2.89 4.531491e-03 False False rate_-500_-300 M start 356 73 4.88 1.34 5.66 2.858284e-07 5.77 4.637736e-08 True True rate_-500_-300 M target 356 90 4.22 2.21 3.16 2.144890e-03 3.30 1.150169e-03 False False rate_-200_0 A start 329 67 7.47 3.79 5.01 4.321812e-06 3.04 2.843883e-03 True False rate_-200_0 A target 329 73 7.06 4.17 4.23 6.758202e-05 2.53 1.262763e-02 True False rate_-200_0 M start 330 105 8.17 1.57 10.11 3.761316e-17 9.65 1.911632e-18 True True rate_-200_0 M target 330 110 7.85 1.47 9.08 5.398348e-15 9.46 5.264819e-18 True True rate_-400_-300 A start 319 61 5.39 2.21 4.09 1.286246e-04 3.74 2.865748e-04 False False rate_-400_-300 A target 319 68 5.01 2.94 2.82 6.295716e-03 2.60 1.025075e-02 False False rate_-400_-300 M start 337 58 4.24 1.13 6.19 7.007675e-08 6.34 4.698853e-09 True True rate_-400_-300 M target 337 68 3.82 2.01 3.27 1.725284e-03 3.40 8.948051e-04 False False Main readings (windows a, b; tests a, b): comparisons meeting p < 1e-4 with distance F larger: 8 of 16
Result: partially reproduced. Distance F exceeds landmark F in every comparison (mean F 4.2 to 8.2 against 1.3 to 4.2), but p < 0.0001 is reached in 8 of the 16 main comparisons (2 windows, 2 tests, 2 landmarks, 2 monkeys). In the last 200 ms (legend window) with the t-test on differences (Methods), all four monkey and landmark comparisons meet it (p 6.8e-05 to 3.8e-17); in -500 to -300 ms (Methods window) only monkey M's start landmark does (p 1.4e-04 to 4.5e-03 for the others). Monkey M sessions with fewer than 10 trials in a condition (4 of 9) are not used.
ed4e-rate-codes-distance-over-slope¶
Near joystick offset the distance effect on the mean firing rate is significantly larger than that on the slope of the firing rate.
Distance F of the rate ANOVA against distance F of the slope ANOVA across units, two-sample t-test (legend; the code's paired test is a variant). Window readings: (a) Methods, rate -500 to -300 ms and slope -550 to -250 ms; (b) legend, rate in the last 200 ms and slope in the last 300 ms; (c) Results, both in -550 to -250 ms.
rows = []
for reading, rw, sw in (("a_methods", "rate_-500_-300", "slope_-550_-250"), ("b_legend", "rate_-200_0", "slope_-300_0"),
("c_results", "rate_-550_-250", "slope_-550_-250")):
for m, g in anova.groupby("monkey"):
gg = g.dropna(subset=[f"{rw}_distance_F", f"{sw}_distance_F"])
fr, fs = gg[f"{rw}_distance_F"], gg[f"{sw}_distance_F"]
t2 = stats.ttest_ind(fr, fs)
tp_ = stats.ttest_rel(fr, fs)
rows.append({"reading": reading, "monkey": m, "units": len(gg), "mean_F_rate": round(fr.mean(), 2), "mean_F_slope": round(fs.mean(), 2),
"median_F_rate": round(fr.median(), 2), "median_F_slope": round(fs.median(), 2),
"two_sample_t": round(t2.statistic, 2), "two_sample_p": t2.pvalue, "paired_t": round(tp_.statistic, 2), "paired_p": tp_.pvalue,
"meets(two-sample)": bool(t2.pvalue < 1e-4 and fr.mean() > fs.mean())})
ed4e = pd.DataFrame(rows)
print(ed4e.to_string(index=False))
reading monkey units mean_F_rate mean_F_slope median_F_rate median_F_slope two_sample_t two_sample_p paired_t paired_p meets(two-sample) a_methods A 324 2.01 1.37 0.90 1.02 2.56 0.010602 2.74 6.452344e-03 False a_methods M 356 1.83 1.69 1.17 1.03 0.71 0.474941 1.40 1.620364e-01 False b_legend A 329 2.29 2.77 1.02 1.13 -1.26 0.207268 -2.04 4.254132e-02 False b_legend M 330 3.27 1.86 1.23 1.05 4.31 0.000019 5.43 1.109513e-07 True c_results A 331 2.25 1.37 0.95 1.02 2.97 0.003054 3.17 1.641541e-03 False c_results M 363 1.99 1.67 1.17 1.02 1.48 0.138149 2.82 5.023089e-03 False
Result: not reproduced. The two-sample t-test finds rate F larger than slope F at p < 0.0001 only in monkey M with the legend's windows (mean F 3.27 against 1.86, p = 1.9e-05). Under the Methods' and the Results' windows neither monkey reaches it (p 0.003 to 0.47), and in monkey A with the legend's windows slope F is the larger (2.77 against 2.29, p = 0.21). The median F is about 1 in every case. No monkey meets the criterion under all readings.
The inter-trial interval (claims ed2b, fig3e)¶
The paper compares the navigation epoch with ITI segments longer than 4 s (Methods, NTS). The trial table has no ITI or fixation event: each row starts at the stimulus onset and stops at the next row's start. The only interval that can hold the ITI is from the joystick offset to the next trial's start, printed here.
rows = []
for s in sessions:
tr = s["trials"]
nxt = tr.start_time.shift(-1)
a1 = tr[(tr.trial_type == 3) & (tr.attempt == 1)]
gap = (a1.stop_time - a1.joy1offttl).dropna()
rows.append({"monkey": s["monkey"], "session": s["session"],
"stop==next_start": round(float(np.mean(np.isclose(tr.stop_time, nxt)[:-1])), 4),
"offset_to_next_start_median_s": round(gap.median(), 2), "p99_s": round(gap.quantile(0.99), 2),
"n_longer_than_4s": int((gap > 4).sum()), "n_trials": len(gap)})
iti = pd.DataFrame(rows)
print(iti.to_string(index=False))
ab = []
for s in sessions:
tr = s["trials"]
a = tr[tr.attempt >= 99]
ab.append({"monkey": s["monkey"], "session": s["session"], "aborts": len(a),
"abort_to_next_start_median_s": round(float((a.stop_time - a[["gocuettl", "joy1onttl"]].max(axis=1)).median()), 2)})
print("\nAborted trials (attempt 99 or 999): time from the joystick press to the next trial's start")
print(pd.DataFrame(ab).to_string(index=False))
monkey session stop==next_start offset_to_next_start_median_s p99_s n_longer_than_4s n_trials
A 08152019 1.0 0.98 2.20 0 754
A 08212019 1.0 0.91 1.80 0 936
A 08232019 1.0 0.87 1.67 2 982
A 08272019 1.0 0.93 1.73 1 953
A 08292019 1.0 0.87 1.98 0 978
A 09172019 1.0 0.87 3.51 2 610
M 04122021 1.0 1.24 1.99 2 1676
M 04132021 1.0 1.26 1.86 2 1239
M 04142021 1.0 1.27 3.58 6 1085
M 04212021 1.0 1.22 2.01 3 1473
M 04222021 1.0 1.30 1.95 2 728
M 04252021 1.0 1.25 2.56 3 1435
M 05042021 1.0 1.23 1.97 2 1731
M 05062021 1.0 1.23 1.84 4 1815
M 05072021 1.0 1.22 2.14 3 1299
Aborted trials (attempt 99 or 999): time from the joystick press to the next trial's start
monkey session aborts abort_to_next_start_median_s
A 08152019 179 2.39
A 08212019 135 2.29
A 08232019 266 2.31
A 08272019 284 2.28
A 08292019 231 2.28
A 09172019 103 2.28
M 04122021 55 6.11
M 04132021 77 6.19
M 04142021 73 6.25
M 04212021 70 6.07
M 04222021 84 6.16
M 04252021 53 6.15
M 05042021 145 6.11
M 05062021 153 6.10
M 05072021 125 6.27
Every row ends where the next begins, and the time from a completed trial's joystick offset to the next trial's start is about 0.9 s (A) or 1.2 s (M), longer than 4 s in 0 to 6 trials per session. After aborted trials (the Methods add a 5 s time-out) the next trial starts a median 6.1 to 6.3 s after the abort in monkey M sessions and 2.3 to 2.4 s in monkey A (printed for each session above); the paper does not say that its ITI segments came from there. The ITI segments longer than 4 s that claims ed2b and fig3e use cannot be located in the files, so those claims are not tested. Claim ed5 uses the selection of fig3e but none of its comparisons uses the ITI, so it is tested below.
ed5-xcorr-conserved-other-epochs¶
The conserved correlation structure generalizes to the inference epoch, error trials and both directions of movement.
This claim's comparisons do not use the ITI: the navigation epoch against the inference epoch (first image onset -0.2 to 1.2 s), against error trials (first attempt of two-attempt trials, the same window as navigation) and between the two directions. Selection as fig3e: the session readings are the chronological fifth monkey A and first monkey M session among the sessions the dandiset holds (the paper's list holds more sessions, so these need not be the paper's) and every session with at least 20 periodic units (main reading); the 20 (legend) or 25 (code) units with the highest PI score. Per pair, the trial-averaged normalised cross-correlation of the detrended smoothed rates (as for the PI), averaged over lags -5 to 5 ms; navigation trials are first-attempt MNAV trials of distance above 3, from joystick onset to 500 ms before offset. Pearson r across pairs; the p-value at lag 0 is Bonferroni-corrected over the lags of Extended Data Fig. 5c (-900 to 900 ms), counted as 1,801 lags (1 ms steps, the stricter) and as 181 (10 ms steps).
sess_by = {s["session"]: s for s in sessions}
rng5 = np.random.default_rng(55)
def pair_xcorr(sp_list, t0, t1, max_lag=5):
"""Per pair of units, the trial-averaged normalised cross-correlation of the detrended smoothed
rates, averaged over lags -max_lag..max_lag ms (claims, fig3e). Units whose rate on a trial is
at or below the PI's 0.05 spikes/s do not enter that trial."""
n = len(sp_list)
acc, cnt = np.zeros((n, n)), np.zeros((n, n))
rates = [mnav.smoothed_rates(sp, t0, t1) for sp in sp_list]
for k in range(len(t0)):
X = np.array([mnav._detrend(r[k]) for r in rates])
ok = np.array([r[k].mean() > mnav.MIN_TRIAL_RATE for r in rates]) & (np.linalg.norm(X, axis=1) > 0)
Xn = np.where(ok[:, None], X / np.where(ok, np.linalg.norm(X, axis=1), 1)[:, None], 0)
T = X.shape[1]
c = sum(Xn[:, :T - l] @ Xn[:, l:].T + (Xn[:, l:] @ Xn[:, :T - l].T if l else 0) for l in range(max_lag + 1))
c = c / (2 * max_lag + 1)
both = np.outer(ok, ok)
acc += np.where(both, c, 0)
cnt += both
iu = np.triu_indices(n, 1)
return np.where(cnt[iu] > 0, acc[iu] / np.maximum(cnt[iu], 1), np.nan)
per_unit_score = base[base[main_col]].groupby(["monkey", "session", "unit"]).score_methods.max().reset_index()
n_per_session = per_unit_score.groupby(["monkey", "session"]).size()
print("Periodic units per session (main reading):", n_per_session.to_dict())
chrono = {"A": sorted(s["session"] for s in sessions if s["monkey"] == "A")[4],
"M": sorted(s["session"] for s in sessions if s["monkey"] == "M")[0]}
chosen = sorted(set([(m, ss) for m, ss in chrono.items()] + [k for k, v in n_per_session.items() if v >= 20]))
print("Sessions: chronological fifth (A) and first (M) in the dandiset:", chrono, "; all with at least 20 periodic units:",
[k for k, v in n_per_session.items() if v >= 20])
N_LAGS_BONF = {"10 ms steps": 181, "1 ms steps": 1801}
rows = []
t = time.time()
for m, ses in chosen:
s = sess_by[ses]
tr = s["trials"]
mn = tr[(tr.trial_type == 3) & (tr.attempt == 1) & (tr.distance > 3) & (tr.seqq < 3) & (tr.tp.abs() < 99)]
er = tr[(tr.trial_type == 3) & (tr.attempt == 2) & (tr.distance > 3) & (tr.seqq < 3) & (tr.tp.abs() > 1) & (tr.tp.abs() < mnav.MAX_TP_S)]
epochs = {"MNAV": (mn.joy1onttl.values, mn.joy1offttl.values - 0.5),
"MNAV direction 1": (mn[mn.direction == 1].joy1onttl.values, mn[mn.direction == 1].joy1offttl.values - 0.5),
"MNAV direction -1": (mn[mn.direction == -1].joy1onttl.values, mn[mn.direction == -1].joy1offttl.values - 0.5),
"inference": (mn.stim1onttl.values - 0.2, mn.stim1onttl.values + 1.2),
"error": (er.joy1onttl.values, er.joy1offttl.values - 0.5)}
cand = per_unit_score[(per_unit_score.monkey == m) & (per_unit_score.session == ses)].sort_values("score_methods", ascending=False)
for n_units in (20, 25):
if len(cand) < n_units:
rows.append({"monkey": m, "session": ses, "neurons": n_units, "note": f"only {len(cand)} periodic units"})
continue
units = cand.unit.values[:n_units]
xc = {name: pair_xcorr([s["spikes"][u] for u in units], a, b) for name, (a, b) in epochs.items()}
for ref, other in (("MNAV", "inference"), ("MNAV", "error"), ("MNAV direction 1", "MNAV direction -1")):
ok = np.isfinite(xc[ref]) & np.isfinite(xc[other])
r_, p_ = stats.pearsonr(xc[ref][ok], xc[other][ok])
rows.append({"monkey": m, "session": ses, "neurons": n_units, "comparison": f"{ref} vs {other}",
"trials": f"{len(epochs[ref][0])} / {len(epochs[other][0])}", "pairs": int(ok.sum()),
"r": round(r_, 3), "p": p_, **{f"p_bonferroni_{k}": min(1.0, p_ * v) for k, v in N_LAGS_BONF.items()},
"meets(r>0, p_bonf<0.05)": bool(r_ > 0 and p_ * N_LAGS_BONF["1 ms steps"] < 0.05)})
ed5 = pd.DataFrame(rows)
print(f"cross-correlations in {time.time() - t:.0f} s")
print(ed5.to_string(index=False))
Periodic units per session (main reading): {('A', '08152019'): 4, ('A', '08212019'): 10, ('A', '08232019'): 1, ('A', '08272019'): 4, ('A', '08292019'): 27, ('A', '09172019'): 10, ('M', '04122021'): 30, ('M', '04132021'): 2, ('M', '04142021'): 8, ('M', '04212021'): 5, ('M', '04222021'): 2, ('M', '04252021'): 4, ('M', '05042021'): 4, ('M', '05062021'): 1, ('M', '05072021'): 1}
Sessions: chronological fifth (A) and first (M) in the dandiset: {'A': '08292019', 'M': '04122021'} ; all with at least 20 periodic units: [('A', '08292019'), ('M', '04122021')]
cross-correlations in 29 s
monkey session neurons comparison trials pairs r p p_bonferroni_10 ms steps p_bonferroni_1 ms steps meets(r>0, p_bonf<0.05)
A 08292019 20 MNAV vs inference 196 / 196 186 0.934 3.695438e-84 6.688742e-82 6.655483e-81 True
A 08292019 20 MNAV vs error 196 / 180 186 0.954 3.381101e-98 6.119792e-96 6.089362e-95 True
A 08292019 20 MNAV direction 1 vs MNAV direction -1 98 / 98 186 0.797 3.705772e-42 6.707448e-40 6.674096e-39 True
A 08292019 25 MNAV vs inference 196 / 196 295 0.902 1.311507e-108 2.373827e-106 2.362024e-105 True
A 08292019 25 MNAV vs error 196 / 180 295 0.940 5.162994e-139 9.345020e-137 9.298553e-136 True
A 08292019 25 MNAV direction 1 vs MNAV direction -1 98 / 98 294 0.815 2.645642e-71 4.788611e-69 4.764800e-68 True
M 04122021 20 MNAV vs inference 136 / 136 190 0.904 2.302415e-71 4.167371e-69 4.146649e-68 True
M 04122021 20 MNAV vs error 136 / 87 190 0.883 1.680464e-63 3.041639e-61 3.026515e-60 True
M 04122021 20 MNAV direction 1 vs MNAV direction -1 62 / 74 190 0.812 6.250833e-46 1.131401e-43 1.125775e-42 True
M 04122021 25 MNAV vs inference 136 / 136 300 0.904 4.748689e-112 8.595128e-110 8.552389e-109 True
M 04122021 25 MNAV vs error 136 / 87 300 0.870 9.650196e-94 1.746686e-91 1.738000e-90 True
M 04122021 25 MNAV direction 1 vs MNAV direction -1 62 / 74 300 0.782 2.869264e-63 5.193368e-61 5.167544e-60 True
Result: reproduced. In the one session per monkey that both session readings select (A 08292019 with 27 periodic units, M 04122021 with 30), every comparison gives a positive r far below the corrected threshold: navigation against the inference epoch r = 0.934 and 0.902 (A, 20 and 25 units) and 0.904 and 0.904 (M); against error trials 0.954 and 0.940 (A) and 0.883 and 0.870 (M); direction 1 against direction -1 0.797 and 0.815 (A) and 0.812 and 0.782 (M); the largest p after correction for 1,801 lags is 6.7e-39. The pairs are 186 to 300 (190 or 300 when every pair has trials in both epochs). The chronological readings rest on the sessions the dandiset holds, which need not be the paper's fifth and first sessions.
fig4h-periodic-fano-factor¶
Many EC neurons have a Fano factor that is periodic at about 0.65 s: 55% (341/614) in monkey A and 43% (378/864) in monkey M.
Extended Data Fig. 8 legend and claims (from the code): correct MNAV trials with |tp| >= 2.3 s
(seqq < 3), per direction, at least 12 trials; spike counts in 200 ms windows stepped by 30 ms
from joystick onset. The end of the time course is not given; it is taken here as 500 ms before
the shortest trial's offset, as for the PI (an assumption). Fano factor across trials,
bootstrapped 100 times (trials resampled); Poisson null from the unit's mean count time course,
100 times; each Fano time course linearly detrended, autocorrelated to 1.6 s, and its PI computed
as for the rate, in 30 ms steps (650 ms is read at the nearest step, 660 ms). Selection readings:
(a) legend, maximum PI in the 300 to 900 ms lag window above the null by a two-sample t-test at
p < 0.0001/20; (b) code, PI at 650 ms against the null by a one-tailed Welch t-test at 0.01.
Direction readings: (i) either direction, (ii) direction 1. Periodicity: lag of the maximum mean
bootstrap PI over 100 to 1,600 ms. Counts are compared as proportions of the units the dandiset
holds, as for fig2.
rng4h = np.random.default_rng(48)
STEP, WIN = 0.03, 0.2
max_lag = int(round(1.6 / STEP))
flags = np.arange(0, max_lag + 1)
i660 = int(round(0.65 / STEP))
win_lo, win_hi = int(np.ceil(0.3 / STEP)), int(np.floor(0.9 / STEP))
def series_pi(x):
"""PI curve (lags 0..1.6 s in 30 ms steps) of a time course, detrended, ACG normalised at 0."""
t_ = np.arange(len(x))
d = x - np.polyval(np.polyfit(t_, x, 1), t_)
n = len(d)
nfft = 1 << int(np.ceil(np.log2(2 * n + 2 * max_lag + 2)))
F = np.fft.rfft(d, nfft)
ac = np.fft.irfft(F * np.conj(F), nfft)
if ac[0] <= 0:
return None
ac = ac[:max_lag + 1] / ac[0]
full = np.concatenate([ac[:0:-1], ac])
return mnav.pi_curve(full, flags)
def fano_curves(counts, rng, n=100, poisson=False):
out = []
mu = counts.mean(0)
for _ in range(n):
c = rng.poisson(np.broadcast_to(mu, counts.shape)) if poisson else counts[rng.integers(0, len(counts), len(counts))]
m_, v_ = c.mean(0), c.var(0, ddof=1)
ff = np.where(m_ > 0, v_ / np.where(m_ > 0, m_, 1), 0.0)
p = series_pi(ff)
if p is not None:
out.append(p)
return np.array(out)
t = time.time()
rows = []
for s in sessions:
tr = s["trials"]
for d in (1, -1):
g = tr[(tr.trial_type == 3) & (tr.attempt == 1) & (tr.tp.abs() >= 2.3) & (tr.tp.abs() < 99) & (tr.seqq < 3) & (tr.direction == d)]
if len(g) < 12:
continue
on = g.joy1onttl.values
end = g.tp.abs().min() - mnav.OFFSET_CUT_S
starts = np.arange(0, end - WIN + 1e-9, STEP)
for u, sp in enumerate(s["spikes"]):
counts = np.stack([np.searchsorted(sp, on + a + WIN) - np.searchsorted(sp, on + a) for a in starts], axis=1)
if counts.sum() == 0:
continue
boot = fano_curves(counts, rng4h)
null = fano_curves(counts, rng4h, poisson=True)
if len(boot) < 10 or len(null) < 10:
continue
bmax = boot[:, win_lo:win_hi + 1].max(1)
nmax = null[:, win_lo:win_hi + 1].max(1)
ta_ = stats.ttest_ind(bmax, nmax)
tb_ = stats.ttest_ind(boot[:, i660], null[:, i660], equal_var=False, alternative="greater")
mean_pi = boot.mean(0)
lo = int(np.ceil(0.1 / STEP))
rows.append({"monkey": s["monkey"], "session": s["session"], "unit": u, "direction": d, "trials": len(g),
"sel_a": bool(ta_.pvalue < 1e-4 / 20 and bmax.mean() > nmax.mean()),
"sel_b": bool(tb_.pvalue < 0.01),
"fano_period": (lo + int(np.argmax(mean_pi[lo:]))) * STEP})
ff = pd.DataFrame(rows)
print(f"Fano factor PI for {len(ff)} unit-direction pairs in {time.time() - t:.0f} s")
paper4h = {"A": (341, 614), "M": (378, 864)}
rows = []
for sel in ("sel_a", "sel_b"):
for dname, dsel in (("either", ff.direction != 0), ("direction 1", ff.direction == 1)):
g = ff[dsel]
per_unit = g.groupby(["monkey", "session", "unit"])[sel].any()
periods = g[g[sel]].groupby(["monkey", "session", "unit"]).fano_period.first()
out = {"selection": sel, "direction": dname}
ok = True
for m in ("A", "M"):
n = int(per_unit.loc[m].sum()) if m in per_unit.index.get_level_values(0) else 0
pct = 100 * n / units_total[m]
ref = 100 * paper4h[m][0] / paper4h[m][1]
out[f"{m}_count"] = n
out[f"{m}_pct"] = round(pct, 1)
out[f"{m}_paper_pct"] = round(ref, 1)
out[f"{m}_in_tol"] = abs(pct - ref) <= 0.1 * ref
pm = periods.loc[m] if m in periods.index.get_level_values(0) else pd.Series(dtype=float)
out[f"{m}_median_period_s"] = round(float(pm.median()), 3) if len(pm) else None
ok &= bool(len(pm) and 0.6 <= pm.median() <= 0.7)
out["medians_in_0.60_0.70"] = ok
rows.append(out)
fig4h = pd.DataFrame(rows)
print(fig4h.to_string(index=False))
print("\nUnit-direction pairs tested per monkey:", ff.groupby("monkey").size().to_dict())
Fano factor PI for 1809 unit-direction pairs in 108 s
selection direction A_count A_pct A_paper_pct A_in_tol A_median_period_s M_count M_pct M_paper_pct M_in_tol M_median_period_s medians_in_0.60_0.70
sel_a either 238 65.4 55.5 False 0.69 381 52.3 43.8 False 0.66 True
sel_a direction 1 161 44.2 55.5 False 0.66 238 32.6 43.8 False 0.69 True
sel_b either 218 59.9 55.5 True 0.69 394 54.0 43.8 False 0.72 False
sel_b direction 1 129 35.4 55.5 False 0.69 215 29.5 43.8 False 0.72 False
Unit-direction pairs tested per monkey: {'A': 651, 'M': 1158}
Result: partially reproduced. Under the legend's selection (a) and either direction, the periodic Fano factor units have median periodicities of 0.69 s (A) and 0.66 s (M), within 0.60 to 0.70 s, but they are 65.4% (238/364) and 52.3% (381/729) of the units, above the paper's 55.5% and 43.8% by more than 10%. With direction 1 only, 44.2% and 32.6%, below it; medians 0.66 and 0.69 s. Under the code's selection (b), monkey A's 59.9% is within tolerance but monkey M's 54.0% is not, and M's median is 0.72 s. No combination puts both proportions within tolerance; under (a) the medians hold. The end of the Fano time course is an assumption (500 ms before the shortest trial's offset).
ed8-fano-periodic-without-rate-periodicity¶
The periodic Fano factor is not a by-product of firing-rate periodicity; it is also present in neurons whose firing rate is not periodic.
Units not periodic in rate (each periodic set, either direction) that have a periodic Fano factor (fig4h readings). One-sided binomial test of their count against 5% of all units (main) and, as variants, 5% of the non-periodic units and the false-positive rate of the Fano selection itself, measured on simulated Poisson units with each sampled unit's own mean count time course.
# False-positive rate of the Fano selections: units whose trials are Poisson with the unit's own mean
# count time course (so with no Fano periodicity), passed through the same selection; 150 random
# unit-directions per monkey.
rng8 = np.random.default_rng(808)
sess_by = {s["session"]: s for s in sessions}
fp_rows = []
for m in ("A", "M"):
sample = ff[ff.monkey == m].sample(150, random_state=int(rng8.integers(1e9)))
for r in sample.itertuples():
s = sess_by[r.session]
tr = s["trials"]
g = tr[(tr.trial_type == 3) & (tr.attempt == 1) & (tr.tp.abs() >= 2.3) & (tr.tp.abs() < 99) & (tr.seqq < 3) & (tr.direction == r.direction)]
on = g.joy1onttl.values
starts = np.arange(0, g.tp.abs().min() - mnav.OFFSET_CUT_S - WIN + 1e-9, STEP)
sp = s["spikes"][r.unit]
counts = np.stack([np.searchsorted(sp, on + a + WIN) - np.searchsorted(sp, on + a) for a in starts], axis=1)
fake = rng8.poisson(np.broadcast_to(counts.mean(0), counts.shape))
if fake.sum() == 0:
continue
boot, null = fano_curves(fake, rng8), fano_curves(fake, rng8, poisson=True)
if len(boot) < 10 or len(null) < 10:
continue
bmax, nmax = boot[:, win_lo:win_hi + 1].max(1), null[:, win_lo:win_hi + 1].max(1)
fp_rows.append({"monkey": m, "sel_a": bool(stats.ttest_ind(bmax, nmax).pvalue < 1e-4 / 20 and bmax.mean() > nmax.mean()),
"sel_b": bool(stats.ttest_ind(boot[:, i660], null[:, i660], equal_var=False, alternative="greater").pvalue < 0.01)})
fp = pd.DataFrame(fp_rows)
fp_rate = fp.groupby("monkey")[["sel_a", "sel_b"]].mean()
print("False-positive rate of each Fano selection per unit-direction, Poisson units without Fano periodicity:")
print(fp.groupby("monkey").size().rename("units simulated").to_frame().join(fp_rate.round(3)).to_string())
rows = []
for set_name, set_col in SETS.items():
rate_periodic = base.groupby(["monkey", "session", "unit"])[set_col].any()
for sel in ("sel_a", "sel_b"):
for dname, dsel in (("either", ff.direction != 0), ("direction 1", ff.direction == 1)):
if set_name != "a2 (main)" and (sel, dname) != ("sel_a", "either"):
continue
per_unit = ff[dsel].groupby(["monkey", "session", "unit"])[sel].any()
df = pd.DataFrame({"fano": per_unit}).join(rate_periodic.rename("rate_periodic"), how="right").fillna(False)
for m in ("A", "M"):
g = df.loc[m]
n_np = int((~g.rate_periodic).sum())
k = int((g.fano & ~g.rate_periodic).sum())
p_all = stats.binomtest(k, int(units_total[m]), 0.05, alternative="greater").pvalue
p_np = stats.binomtest(k, n_np, 0.05, alternative="greater").pvalue
# chance per unit from the measured false-positive rate per unit-direction (either
# direction: 1 - (1 - rate)^2)
fpr = fp_rate.loc[m, sel]
fpr_unit = 1 - (1 - fpr) ** 2 if dname == "either" else fpr
p_fp = stats.binomtest(k, n_np, min(fpr_unit, 1.0), alternative="greater").pvalue
rows.append({"periodic_set": set_name, "selection": sel, "direction": dname, "monkey": m, "units": int(units_total[m]),
"nonperiodic_rate": n_np, "fano_periodic_among_nonperiodic": k,
"binom_p_vs_5pct_all": p_all, "binom_p_vs_5pct_nonperiodic": p_np,
"measured_chance_per_unit": round(fpr_unit, 3), "binom_p_vs_measured_chance": p_fp,
"meets_p<0.05": bool(p_all < 0.05)})
ed8 = pd.DataFrame(rows)
print(ed8.to_string(index=False))
print("Code's unlabelled counts, for information: A 99/614, M 74/864.")
False-positive rate of each Fano selection per unit-direction, Poisson units without Fano periodicity:
units simulated sel_a sel_b
monkey
A 150 0.033 0.220
M 150 0.020 0.167
periodic_set selection direction monkey units nonperiodic_rate fano_periodic_among_nonperiodic binom_p_vs_5pct_all binom_p_vs_5pct_nonperiodic measured_chance_per_unit binom_p_vs_measured_chance meets_p<0.05
a2 (main) sel_a either A 364 308 193 1.048235e-147 2.544874e-167 0.066 1.950753e-145 True
a2 (main) sel_a either M 729 672 334 3.146055e-227 5.312973e-242 0.040 3.108544e-274 True
a2 (main) sel_a direction 1 A 364 308 135 1.511747e-78 7.841919e-90 0.033 2.619374e-112 True
a2 (main) sel_a direction 1 M 729 672 205 1.923708e-92 9.470757e-100 0.020 4.674957e-175 True
a2 (main) sel_b either A 364 308 179 1.560362e-129 7.389577e-147 0.392 1.425101e-11 True
a2 (main) sel_b either M 729 672 347 4.126752e-243 8.885945e-259 0.306 5.142059e-30 True
a2 (main) sel_b direction 1 A 364 308 104 1.595764e-48 2.596021e-56 0.220 1.428463e-06 True
a2 (main) sel_b direction 1 M 729 672 188 4.797432e-78 1.612811e-84 0.167 1.630873e-13 True
a1 sel_a either A 364 354 230 2.752326e-200 1.683540e-204 0.066 2.502845e-178 True
a1 sel_a either M 729 729 381 1.492681e-286 1.492681e-286 0.040 1.482197e-323 True
b sel_a either A 364 125 55 3.952729e-13 1.014333e-37 0.066 9.653625e-32 True
b sel_a either M 729 350 119 1.629619e-29 1.480043e-64 0.040 1.586416e-75 True
a2, gp null sel_a either A 364 332 211 1.893976e-172 1.084594e-184 0.066 9.865382e-161 True
a2, gp null sel_a either M 729 721 373 4.395921e-276 2.055249e-278 0.040 1.506875e-314 True
Code's unlabelled counts, for information: A 99/614, M 74/864.
Result: reproduced. Among units not periodic in rate under the main reading, 193 of 308 (A) and 334 of 672 (M) have a periodic Fano factor under the legend's selection and either direction, far above 5% of all units (binomial p 1.0e-147 and 3.1e-227). The measured false-positive rate of that selection is 3.3% (A) and 2.0% (M) per unit-direction (6.6% and 4.0% per unit over two directions), and the counts are far above it too (p 2.0e-145 and 3.1e-274). The code's selection has false-positive rates of 22.0% and 16.7% per unit-direction; its counts are still above them (p 1.4e-06 or below). Every selection, direction and periodic-set reading gives p < 1e-12 against 5% of all units (largest 4.0e-13, rule b, monkey A). The counts are larger than the code's unlabelled 99/614 and 74/864, as the Fano selection here admits more units (fig4h).
fig2e-first-peak-near-065 and fig2f-earlier-peaks-multiples: phase of the local peaks¶
Methods (Quantification of the phase of local bumps) and claims (from the code): per periodic
unit and direction in which it is periodic, for each periodic set; first-attempt trials not flagged as
lapses (validtrials_mm), distance 4 or 5, seqq < 3; trials with no spikes from 1.5 to 0.5 s
before offset dropped (the claims give "low rate" without a threshold; the PI's 0.05 spikes/s is
used, an assumption); at least 25 trials. 100 times, the smoothed rate averaged over a random
third of the trials; the first peak is the highest local maximum (findpeaks, prominence 0.05 on
the min-max normalised window) in the 1 s before offset (Methods), the second in the 1 s before
the first, the third in the 1 s before the second. Null: each trial's spike train circularly
shifted by a random length within the 3.5 s analysed, then the same. Per unit, the KS test of
bootstrap against null first-peak phases, the mode of the bootstrap phases (10 ms bins, a
standard choice) and the median intervals between peaks. Across units whose KS p < 0.001, the
median of the modes; readings for the KS P << 0.0001: (a) per unit, (b) pooled phases. Variant:
the code's -1 to -0.1 s window, on the main reading's units.
from scipy.signal import find_peaks
rng2e = np.random.default_rng(265)
SPAN = 3.5
def peak_phases(rate, first_window=(-1.0, 0.0)):
"""Phases (s, negative) of the first three peaks of a rate on the window [-SPAN, 0) before offset."""
out = []
hi = first_window[1]
lo = first_window[0]
for _ in range(3):
a, b = int(round((SPAN + lo) * 1000)), int(round((SPAN + hi) * 1000))
if a < 0 or b - a < 3:
out.append(np.nan)
continue
seg = rate[a:b]
rngv = seg.max() - seg.min()
if rngv <= 0:
out.append(np.nan)
continue
x = (seg - seg.min()) / rngv
pk, _ = find_peaks(x, prominence=0.05)
if len(pk) == 0:
out.append(np.nan)
break
best = pk[np.argmax(x[pk])]
phase = (a + best) / 1000 - SPAN
out.append(phase)
hi, lo = phase, phase - 1.0
while len(out) < 3:
out.append(np.nan)
return out
def phase_unit(counts, rng, n=100, window=(-1.0, 0.0)):
k = max(1, len(counts) // 3)
boot, null = [], []
for _ in range(n):
idx = rng.choice(len(counts), k, replace=False)
r = gaussian(counts[idx].mean(0))
boot.append(peak_phases(r, window))
sh = np.stack([np.roll(counts[i], int(rng.integers(0, counts.shape[1]))) for i in idx])
null.append(peak_phases(gaussian(sh.mean(0)), window))
return np.array(boot), np.array(null)
from scipy.ndimage import gaussian_filter1d
def gaussian(x):
return gaussian_filter1d(x, mnav.SMOOTH_SD_MS, truncate=mnav.SMOOTH_TRUNCATE, mode="nearest") * 1000
set_cols = list(dict.fromkeys(SETS.values()))
periodic_pairs = base[base[set_cols].any(axis=1)][["monkey", "session", "unit", "direction"] + set_cols]
sess_by = {s["session"]: s for s in sessions}
t = time.time()
rows, pooled = [], {}
for window_name, window in (("methods_-1_0", (-1.0, 0.0)), ("code_-1_-0.1", (-1.0, -0.1))):
# the code's window is run on the main reading's units only
pairs = periodic_pairs if window_name == "methods_-1_0" else periodic_pairs[periodic_pairs[main_col]]
for r in pairs.itertuples():
s = sess_by[r.session]
tr = s["trials"]
g = tr[(tr.trial_type == 3) & (tr.attempt == 1) & (tr.validtrials_mm == 1) & tr.distance.isin([4, 5])
& (tr.seqq < 3) & (tr.direction == r.direction)]
sp = s["spikes"][r.unit]
off = g.joy1offttl.values
keep = mnav.mean_rate(sp, off - 1.5, off - 0.5) > mnav.MIN_TRIAL_RATE
off = off[keep]
if len(off) < 25:
rows.append({"window": window_name, "monkey": r.monkey, "session": r.session, "unit": r.unit, "direction": r.direction, "trials": len(off)})
continue
counts, _n = mnav.trial_counts(sp, off - SPAN, off)
boot, null = phase_unit(counts, rng2e, window=window)
b1, n1 = boot[:, 0][np.isfinite(boot[:, 0])], null[:, 0][np.isfinite(null[:, 0])]
ks = stats.ks_2samp(b1, n1).pvalue if len(b1) > 1 and len(n1) > 1 else np.nan
row = {"window": window_name, "monkey": r.monkey, "session": r.session, "unit": r.unit, "direction": r.direction,
"trials": len(off), "ks_p": ks}
for j, name in enumerate(("first", "second", "third")):
v = boot[:, j][np.isfinite(boot[:, j])]
if len(v):
h, e = np.histogram(v, bins=np.arange(-SPAN, 0.0001, 0.01))
row[f"mode_{name}"] = round(e[h.argmax()] + 0.005, 3)
pooled[(window_name, r.monkey, r.session, r.unit, r.direction)] = (b1, n1)
row["interval_1_2"] = float(np.nanmedian(boot[:, 0] - boot[:, 1])) if np.isfinite(boot[:, 1]).any() else np.nan
row["interval_2_3"] = float(np.nanmedian(boot[:, 1] - boot[:, 2])) if np.isfinite(boot[:, 2]).any() else np.nan
rows.append(row)
phase = pd.DataFrame(rows).merge(periodic_pairs, on=["monkey", "session", "unit", "direction"])
print(f"phase analysis in {time.time() - t:.0f} s; periodic unit-directions: {len(periodic_pairs)}; "
f"with at least 25 trials: {int(phase[phase.window == 'methods_-1_0'].ks_p.notna().sum())}")
phase analysis in 114 s; periodic unit-directions: 995; with at least 25 trials: 660
def peak_set_rows(variant_rows, exclude=()):
rows = []
for set_name, set_col in SETS.items():
for (w, m), g in phase[phase.ks_p.notna() & phase[set_col] & ~phase.session.isin(exclude)].groupby(["window", "monkey"]):
if w != "methods_-1_0" and set_col != main_col: # the code's window is run on the main reading only
continue
sig = g[g.ks_p < 0.001]
keys = [(w, r.monkey, r.session, r.unit, r.direction) for r in g.itertuples()]
b = np.concatenate([pooled[k][0] for k in keys])
nn = np.concatenate([pooled[k][1] for k in keys])
rows.append({"periodic_set": set_name, "window": w, "monkey": m, "units_tested": len(g), "units_ks_p<0.001": len(sig),
"units_ks_p<0.0001": int((g.ks_p < 1e-4).sum()), "frac_p<0.0001": round(float((g.ks_p < 1e-4).mean()), 3),
"pooled_ks_p": stats.ks_2samp(b, nn).pvalue, **variant_rows(sig)})
return pd.DataFrame(rows)
def first_peak(sig):
return {"median_first_mode_s": sig.mode_first.median(),
"first_mode_iqr": f"{sig.mode_first.quantile(.25):.2f}..{sig.mode_first.quantile(.75):.2f}",
"in_-0.70_-0.60": bool(-0.70 <= sig.mode_first.median() <= -0.60)}
fig2e = peak_set_rows(first_peak)
print(fig2e.to_string(index=False))
sig_m = phase[(phase.window == "methods_-1_0") & (phase.ks_p < 0.001) & phase[main_col]]
print("\nMain reading, Methods window, both monkeys pooled: median first-peak mode", sig_m.mode_first.median(), "s over", len(sig_m), "units")
print("\nBy direction (main reading, Methods window, units with KS p < 0.001):")
print(sig_m.groupby(["monkey", "direction"]).mode_first.agg(["size", "median"]).to_string())
print("Code's choice of direction for Fig. 2f (A direction -1, M direction 1):",
sig_m[((sig_m.monkey == "A") & (sig_m.direction == -1)) | ((sig_m.monkey == "M") & (sig_m.direction == 1))]
.groupby("monkey").mode_first.median().to_dict())
print("\nWithout A 09172019, M 04132021 and M 04142021 (clock not established):")
print(peak_set_rows(first_peak, exclude=no_clock).query("window == 'methods_-1_0'").to_string(index=False))
periodic_set window monkey units_tested units_ks_p<0.001 units_ks_p<0.0001 frac_p<0.0001 pooled_ks_p median_first_mode_s first_mode_iqr in_-0.70_-0.60
a2 (main) code_-1_-0.1 A 55 45 38 0.691 5.984811e-94 -0.675 -0.77..-0.60 True
a2 (main) code_-1_-0.1 M 34 27 24 0.706 3.043076e-39 -0.825 -0.93..-0.49 False
a2 (main) methods_-1_0 A 55 45 36 0.655 1.700576e-71 -0.655 -0.77..-0.43 True
a2 (main) methods_-1_0 M 34 30 26 0.765 1.565572e-40 -0.780 -0.90..-0.54 False
a1 methods_-1_0 A 16 12 9 0.562 1.256105e-24 -0.645 -0.70..-0.47 True
b methods_-1_0 A 298 227 184 0.617 6.639578e-105 -0.625 -0.80..-0.20 True
b methods_-1_0 M 358 288 265 0.740 6.759537e-73 -0.595 -0.87..-0.34 False
a2, gp null methods_-1_0 A 36 28 23 0.639 1.925503e-66 -0.665 -0.74..-0.61 True
a2, gp null methods_-1_0 M 1 1 1 1.000 5.562326e-08 -0.555 -0.56..-0.56 False
Main reading, Methods window, both monkeys pooled: median first-peak mode -0.695 s over 75 units
By direction (main reading, Methods window, units with KS p < 0.001):
size median
monkey direction
A -1 29 -0.685
1 16 -0.520
M -1 11 -0.575
1 19 -0.835
Code's choice of direction for Fig. 2f (A direction -1, M direction 1): {'A': -0.685, 'M': -0.835}
Without A 09172019, M 04132021 and M 04142021 (clock not established):
periodic_set window monkey units_tested units_ks_p<0.001 units_ks_p<0.0001 frac_p<0.0001 pooled_ks_p median_first_mode_s first_mode_iqr in_-0.70_-0.60
a2 (main) methods_-1_0 A 54 44 35 0.648 1.971944e-69 -0.655 -0.75..-0.43 True
a2 (main) methods_-1_0 M 34 30 26 0.765 1.565572e-40 -0.780 -0.90..-0.54 False
a1 methods_-1_0 A 16 12 9 0.562 1.256105e-24 -0.645 -0.70..-0.47 True
b methods_-1_0 A 266 210 171 0.643 4.111383e-95 -0.615 -0.78..-0.19 True
b methods_-1_0 M 346 280 260 0.751 9.706447e-81 -0.595 -0.87..-0.35 False
a2, gp null methods_-1_0 A 36 28 23 0.639 1.925503e-66 -0.665 -0.74..-0.61 True
a2, gp null methods_-1_0 M 1 1 1 1.000 5.562326e-08 -0.555 -0.56..-0.56 False
Result: partially reproduced. Under the main reading and the Methods' window, the median first-peak mode of the units whose KS test gives p < 0.001 is -0.655 s in monkey A (45 units), within -0.70 to -0.60 s, and -0.780 s in monkey M (30 units), 0.08 s outside it. The KS test separates bootstrap from null in most units (p < 0.0001 in 36 of 55 in A and 26 of 34 in M) and pooled (p 1.7e-71 and 1.6e-40), so both readings of the test hold. The code's window gives -0.675 (A) and -0.825 s (M). Rule b gives -0.625 (A) and -0.595 s (M), the GP-null reading -0.665 s in A (one unit in M), and rule a1 -0.645 s in A (no unit in M); no reading puts M within the range. By direction the medians are -0.685 and -0.520 s (A) and -0.575 and -0.835 s (M). Pooled over both monkeys the median is -0.695 s (75 units). Without the three sessions whose clock could not be checked the main reading gives the same (-0.655 and -0.780 s).
def later_peaks(g):
return {"median_second_s": g.mode_second.median(), "n_second": int(g.mode_second.notna().sum()),
"median_third_s": g.mode_third.median(), "n_third": int(g.mode_third.notna().sum()),
"median_interval_1_2_s": round(g.interval_1_2.median(), 3), "median_interval_2_3_s": round(g.interval_2_3.median(), 3),
"second_in_-1.35_-1.25": bool(abs(g.mode_second.median() + 1.30) <= 0.05),
"third_in_-2.00_-1.90": bool(abs(g.mode_third.median() + 1.95) <= 0.05)}
fig2f = peak_set_rows(later_peaks)
print(fig2f.drop(columns=["units_ks_p<0.0001", "frac_p<0.0001", "pooled_ks_p"]).to_string(index=False))
print("\nMain reading, Methods window, both monkeys pooled: median second-peak mode", sig_m.mode_second.median(),
"s; median third-peak mode", sig_m.mode_third.median(), "s")
print("\nPer-unit intervals between peaks (median over bootstraps), main reading, Methods window, units with KS p < 0.001:")
print(sig_m.groupby("monkey")[["interval_1_2", "interval_2_3"]].describe(percentiles=[0.25, 0.5, 0.75]).round(3).T.to_string())
print("\nWithout A 09172019, M 04132021 and M 04142021:")
print(peak_set_rows(later_peaks, exclude=no_clock).query("window == 'methods_-1_0'")
.drop(columns=["units_ks_p<0.0001", "frac_p<0.0001", "pooled_ks_p"]).to_string(index=False))
print("\nWithout A 09172019 only:")
print(peak_set_rows(later_peaks, exclude=["09172019"]).query("window == 'methods_-1_0' and periodic_set == 'a2 (main)'")
.drop(columns=["units_ks_p<0.0001", "frac_p<0.0001", "pooled_ks_p"]).to_string(index=False))
periodic_set window monkey units_tested units_ks_p<0.001 median_second_s n_second median_third_s n_third median_interval_1_2_s median_interval_2_3_s second_in_-1.35_-1.25 third_in_-2.00_-1.90
a2 (main) code_-1_-0.1 A 55 45 -1.255 45 -1.735 45 0.525 0.456 True False
a2 (main) code_-1_-0.1 M 34 27 -1.115 27 -1.615 27 0.391 0.394 False False
a2 (main) methods_-1_0 A 55 45 -1.245 45 -1.445 45 0.550 0.407 False False
a2 (main) methods_-1_0 M 34 30 -1.070 30 -1.745 30 0.344 0.394 False False
a1 methods_-1_0 A 16 12 -1.310 12 -1.450 12 0.574 0.478 True False
b methods_-1_0 A 298 227 -1.115 227 -1.475 227 0.498 0.390 False False
b methods_-1_0 M 358 288 -1.015 288 -1.305 287 0.269 0.145 False False
a2, gp null methods_-1_0 A 36 28 -1.285 28 -1.445 28 0.561 0.477 True False
a2, gp null methods_-1_0 M 1 1 -1.145 1 -1.445 1 0.588 0.410 False False
Main reading, Methods window, both monkeys pooled: median second-peak mode -1.205 s; median third-peak mode -1.495 s
Per-unit intervals between peaks (median over bootstraps), main reading, Methods window, units with KS p < 0.001:
monkey A M
interval_1_2 count 45.000 30.000
mean 0.484 0.337
std 0.198 0.206
min 0.061 0.055
25% 0.394 0.135
50% 0.550 0.344
75% 0.622 0.520
max 0.737 0.692
interval_2_3 count 45.000 30.000
mean 0.398 0.361
std 0.235 0.207
min 0.058 0.043
25% 0.176 0.141
50% 0.407 0.394
75% 0.568 0.533
max 0.828 0.678
Without A 09172019, M 04132021 and M 04142021:
periodic_set window monkey units_tested units_ks_p<0.001 median_second_s n_second median_third_s n_third median_interval_1_2_s median_interval_2_3_s second_in_-1.35_-1.25 third_in_-2.00_-1.90
a2 (main) methods_-1_0 A 54 44 -1.240 44 -1.445 44 0.547 0.406 False False
a2 (main) methods_-1_0 M 34 30 -1.070 30 -1.745 30 0.344 0.394 False False
a1 methods_-1_0 A 16 12 -1.310 12 -1.450 12 0.574 0.478 True False
b methods_-1_0 A 266 210 -1.110 210 -1.475 210 0.497 0.391 False False
b methods_-1_0 M 346 280 -1.010 280 -1.265 279 0.264 0.137 False False
a2, gp null methods_-1_0 A 36 28 -1.285 28 -1.445 28 0.561 0.477 True False
a2, gp null methods_-1_0 M 1 1 -1.145 1 -1.445 1 0.588 0.410 False False
Without A 09172019 only:
periodic_set window monkey units_tested units_ks_p<0.001 median_second_s n_second median_third_s n_third median_interval_1_2_s median_interval_2_3_s second_in_-1.35_-1.25 third_in_-2.00_-1.90 a2 (main) methods_-1_0 A 54 44 -1.24 44 -1.445 44 0.547 0.406 False False a2 (main) methods_-1_0 M 34 30 -1.07 30 -1.745 30 0.344 0.394 False False
Result: not reproduced. Under the main reading and the Methods' window, the median second-peak mode is -1.245 s in monkey A and -1.070 s in monkey M, against -1.30 s, and the third -1.445 and -1.745 s, against -1.95 s. The median interval between the first and second peak within a unit is 0.550 s (A) and 0.344 s (M), and between the second and third 0.407 and 0.394 s, shorter than 0.65 s. Monkey A's second peak is within 0.05 s of -1.30 s under the code's window (-1.255 s, 45 units), rule a1 (-1.310 s, 12 units) and the GP-null reading (-1.285 s, 28 units), and 0.005 s outside it under the main reading; monkey M's second peak is outside it under every reading, and no reading puts the third peak near -1.95 s in either monkey. So the second peak holds at most in one monkey and the third in neither, which meets neither partial clause cleanly. Leaving out A 09172019 gives -1.240 s for A's second peak. Pooled over both monkeys the medians are -1.205 and -1.495 s.
ed3-no-eye-hand-periodicity¶
Eye position (horizontal and vertical) and joystick position during the navigation epoch show no periodicity; their PI at 650 ms is below chance.
In five monkey A sessions the eye and joystick timestamps span the trial times. In the other ten
(A 08152019 and every monkey M session) they are on another clock, whose offset the files do not
give; it is estimated here by cross-correlating the time the joystick is deflected with the
intervals between the trial table's joystick onsets and offsets (mnav.behaviour_clock_offset).
The same estimate in the five sessions on the trial clock gives -0.12 to -0.26 s, so its median
(-0.17 s) is taken as the estimate's bias and subtracted. Each shifted session is kept only if,
after the shift, the joystick deviates from rest at least twice as much 0.2 to 0.45 s after the
recorded onsets as before them, and more than 7 s later. The signals are treated as the rates are
for the PI (claims): PI trials of both directions pooled, joystick onset to 500 ms before offset,
the trace interpolated to 1 ms and smoothed with the rate kernel, each trial detrended, mean ACG,
PI (Methods lags). Bootstrap readings: (a) 100 resamples of the pooled trials with replacement,
full size; (b) 100 subsamples of a third without replacement. Null: Gaussian process surrogates
(SE kernel, length 100 ms) smoothed as the signals, without the Poisson step (the signals are
continuous), averaged over as many trials as a bootstrap sample, 1,000 times. Chance readings: (a)
the surrogate PI at 650 ms mean + 2 SD; (b) its mean. One-tailed one-sample t-test of the 100
bootstrap PIs at 650 ms against chance, per monkey.
from scipy.ndimage import gaussian_filter1d as _g1d
rng3 = np.random.default_rng(3003)
# Each session's eye and joystick traces. Where their timestamps do not span the trial times (A
# 08152019 and the monkey M sessions), the clock offset is estimated from the joystick trace
# (mnav.behaviour_clock_offset) and the traces are shifted by it; the estimate is also made, as a
# check, in the sessions whose traces are on the trial clock, where it should be near 0.
signals = {}
rows, ests = [], []
t = time.time()
# the sessions on the trial clock come first, so that the median of their estimates (the bias of
# the estimate) can be subtracted from the others
def _on_clock(s):
return (abs(s["info"]["eye_ts_first"] - s["trials"].start_time.min()) < 500
and abs(s["info"]["eye_ts_last"] - s["trials"].stop_time.max()) < 500)
calib = 0.0
for s in sorted(sessions, key=lambda s: not _on_clock(s)):
if not _on_clock(s) and ests and calib == 0.0:
calib = float(np.median([e for c, e in ests if c]))
print(f"median estimate in the sessions on the trial clock (bias subtracted from the others): {calib:.3f} s")
asset = [a for a in mnav.assets() if a.path == s["path"]][0]
f = mnav._open(asset)
on_clock = _on_clock(s)
hts = f["processing/behavior/hand_position/timestamps"][:]
hx = f["processing/behavior/hand_position/data"][:].ravel()
est, z, second = mnav.behaviour_clock_offset(hts, hx, s["trials"])
ests.append((on_clock, est))
shift = 0.0 if on_clock else est - calib
tr = mnav.pi_trials(s["trials"])
t0, t1 = mnav.nav_windows(tr)
ets = f["processing/behavior/eye_position/timestamps"][:] - shift
elo, ehi = np.searchsorted(ets, t0 - 0.01), np.searchsorted(ets, t1 + 0.01)
e0 = int(elo.min()) if len(elo) else 0
eye = f["processing/behavior/eye_position/data"][e0:int(ehi.max()) if len(ehi) else 0] # one read
for name, col in (("eye_x", 0), ("eye_y", 1), ("joystick", None)):
ts = hts - shift if col is None else ets
lo, hi = (np.searchsorted(ts, t0 - 0.01), np.searchsorted(ts, t1 + 0.01)) if col is None else (elo, ehi)
for a, b, x0, x1 in zip(lo, hi, t0, t1):
if b - a < 10:
continue
seg = hx[a:b] if col is None else eye[a - e0:b - e0, col]
signals.setdefault((s["monkey"], s["session"], name), []).append(np.interp(np.arange(x0, x1, 0.001), ts[a:b], seg))
# clock check: joystick deflection around the trial table's joystick onsets, after the shift
hts_s = hts - shift
on = s["trials"][(s["trials"].trial_type == 3) & (s["trials"].attempt == 1)].joy1onttl.values[:200]
pre = np.concatenate([hx[x:y] for x, y in zip(np.searchsorted(hts_s, on - 0.3), np.searchsorted(hts_s, on - 0.05))])
rest = float(np.median(pre)) if len(pre) else np.nan
def dev(times, a_, b_):
l_, h_ = np.searchsorted(hts_s, times + a_), np.searchsorted(hts_s, times + b_)
v = [np.abs(hx[x:y] - rest).mean() for x, y in zip(l_, h_) if y > x]
return float(np.mean(v)) if v else np.nan
rows.append({"monkey": s["monkey"], "session": s["session"], "on_trial_clock": on_clock,
"estimated_offset_s": round(est, 3), "peak_z": round(z, 1), "far_peak/peak": round(second, 3),
"shift_used_s": round(shift, 3), "PI_trials": len(tr), "traces": len(signals.get((s["monkey"], s["session"], "joystick"), [])),
"|dev| before onset (V)": round(dev(on, -0.3, -0.05), 3), "|dev| 0.2-0.45 s after onset": round(dev(on, 0.2, 0.45), 3),
"same, onset + 7 s": round(dev(on + 7.0, 0.2, 0.45), 3)})
beh = pd.DataFrame(rows)
print(f"streamed in {time.time() - t:.0f} s")
print(beh.to_string(index=False))
median estimate in the sessions on the trial clock (bias subtracted from the others): -0.169 s
streamed in 171 s
monkey session on_trial_clock estimated_offset_s peak_z far_peak/peak shift_used_s PI_trials traces |dev| before onset (V) |dev| 0.2-0.45 s after onset same, onset + 7 s
A 08212019 True -0.149 36.0 0.212 0.000 378 378 0.429 2.485 1.077
A 08232019 True -0.257 37.9 0.237 0.000 393 393 0.600 2.451 1.142
A 08272019 True -0.169 33.3 0.208 0.000 383 383 0.501 2.464 1.188
A 08292019 True -0.189 31.8 0.262 0.000 392 392 0.702 2.462 1.323
A 09172019 True -0.120 30.1 0.258 0.000 239 239 0.610 1.953 1.162
A 08152019 False 775445.202 30.8 0.295 775445.371 303 303 0.494 2.502 1.137
M 04122021 False 431609.942 40.5 0.261 431610.111 288 288 0.157 2.434 0.787
M 04132021 False 519842.944 18.6 0.358 519843.113 160 160 0.272 2.460 0.922
M 04142021 False 604086.332 24.3 0.292 604086.502 130 130 0.183 2.480 0.982
M 04212021 False 1207222.223 87.4 0.201 1207222.392 303 303 0.232 2.532 0.803
M 04222021 False 1298314.258 28.8 0.311 1298314.427 130 130 0.678 2.392 0.998
M 04252021 False 1561450.070 15.1 0.249 1561450.240 295 295 0.843 2.529 1.389
M 05042021 False 3326.229 83.0 0.213 3326.398 263 263 0.274 2.490 0.808
M 05062021 False 172434.218 74.4 0.204 172434.387 239 239 0.235 2.520 0.726
M 05072021 False 258441.419 38.5 0.226 258441.588 190 190 0.325 2.499 1.190
def smooth_trace(x):
return _g1d(x, mnav.SMOOTH_SD_MS, truncate=mnav.SMOOTH_TRUNCATE, mode="nearest")
def acgs_of(traces):
out = []
for x in traces:
a, _k = mnav.mean_acg([smooth_trace(x)], min_rate=-np.inf)
if a is not None:
out.append(a)
return np.array(out)
def gp_trace(rng, n):
sig = mnav.GP_LENGTH_MS / np.sqrt(2)
w = rng.standard_normal(n + 800)
return _g1d(w, sig, mode="wrap")[400:400 + n]
lags = mnav.LAGS_METHODS
i650 = int(np.nonzero(lags == 650)[0][0])
clock_ok = beh[(beh.on_trial_clock) | ((beh["|dev| 0.2-0.45 s after onset"] > 2 * beh["|dev| before onset (V)"])
& (beh["|dev| 0.2-0.45 s after onset"] > beh["same, onset + 7 s"]))]
rows = []
for m in ("A", "M"):
ses_m = clock_ok[clock_ok.monkey == m].session.tolist()
durations = [len(x) for ss in ses_m for x in signals.get((m, ss, "joystick"), [])]
if not durations:
continue
bank = acgs_of([gp_trace(rng3, n) for n in rng3.choice(durations, 3000)])
for name in ("eye_x", "eye_y", "joystick"):
A_ = acgs_of([x for ss in ses_m for x in signals.get((m, ss, name), [])])
n = len(A_)
for boot in ("a_full_with_replacement", "b_third_without"):
k = n if boot.startswith("a") else n // 3
w = np.zeros((100, n))
for j in range(100):
idx = rng3.integers(0, n, n) if boot.startswith("a") else rng3.choice(n, k, replace=False)
np.add.at(w[j], idx, 1.0 / k)
bpi = mnav.pi_curves(w @ A_, {"m": lags})["m"][:, i650]
wn = np.zeros((1000, len(bank)))
for j in range(1000):
np.add.at(wn[j], rng3.integers(0, len(bank), k), 1.0 / k)
null = mnav.pi_curves(wn @ bank, {"m": lags})["m"][:, i650]
for chance_name, chance in (("a_mean+2SD", null.mean() + 2 * null.std()), ("b_mean", null.mean())):
lo_ = stats.ttest_1samp(bpi, chance, alternative="less")
hi_ = stats.ttest_1samp(bpi, chance, alternative="greater")
rows.append({"monkey": m, "sessions": len(ses_m), "signal": name, "trials": n, "bootstrap": boot, "chance": chance_name,
"mean_PI650": round(bpi.mean(), 3), "chance_value": round(chance, 3), "t(99)": round(lo_.statistic, 1),
"p_below": lo_.pvalue, "p_above": hi_.pvalue,
"below_p<1e-4": bool(lo_.pvalue < 1e-4), "above_p<0.05": bool(hi_.pvalue < 0.05)})
ed3 = pd.DataFrame(rows)
print("Sessions used (on the trial clock, or shifted with the joystick deflection then aligned to the recorded onsets):")
print(clock_ok.groupby("monkey").session.apply(list).to_dict())
print(ed3.to_string(index=False))
print("\nPaper: mean PI at 650 ms, horizontal eye -0.25, vertical eye -0.44, joystick -0.53 (tolerance 10%).")
print("Readings (bootstrap x chance) with all three signals below chance at p < 1e-4, per monkey:")
print(ed3.groupby(["monkey", "bootstrap", "chance"])["below_p<1e-4"].all().to_string())
print("Any signal above chance at p < 0.05 under any reading:", ed3.groupby("monkey")["above_p<0.05"].any().to_dict())
Sessions used (on the trial clock, or shifted with the joystick deflection then aligned to the recorded onsets):
{'A': ['08212019', '08232019', '08272019', '08292019', '09172019', '08152019'], 'M': ['04122021', '04132021', '04142021', '04212021', '04222021', '04252021', '05042021', '05062021', '05072021']}
monkey sessions signal trials bootstrap chance mean_PI650 chance_value t(99) p_below p_above below_p<1e-4 above_p<0.05
A 6 eye_x 2088 a_full_with_replacement a_mean+2SD 0.155 0.270 -66.3 3.890710e-84 1.0 True False
A 6 eye_x 2088 a_full_with_replacement b_mean 0.155 0.231 -43.7 8.086464e-67 1.0 True False
A 6 eye_x 2088 b_third_without a_mean+2SD 0.156 0.298 -57.5 3.591335e-78 1.0 True False
A 6 eye_x 2088 b_third_without b_mean 0.156 0.230 -30.1 5.665076e-52 1.0 True False
A 6 eye_y 2088 a_full_with_replacement a_mean+2SD 0.117 0.274 -104.7 1.604867e-103 1.0 True False
A 6 eye_y 2088 a_full_with_replacement b_mean 0.117 0.233 -77.5 1.022332e-90 1.0 True False
A 6 eye_y 2088 b_third_without a_mean+2SD 0.121 0.301 -73.4 1.986588e-88 1.0 True False
A 6 eye_y 2088 b_third_without b_mean 0.121 0.232 -45.2 3.143917e-68 1.0 True False
A 6 joystick 2088 a_full_with_replacement a_mean+2SD -0.401 0.273 -472.2 4.360347e-168 1.0 True False
A 6 joystick 2088 a_full_with_replacement b_mean -0.401 0.232 -443.3 2.241296e-165 1.0 True False
A 6 joystick 2088 b_third_without a_mean+2SD -0.398 0.303 -300.0 1.342282e-148 1.0 True False
A 6 joystick 2088 b_third_without b_mean -0.398 0.232 -269.8 4.768551e-144 1.0 True False
M 9 eye_x 1998 a_full_with_replacement a_mean+2SD -0.370 0.350 -496.0 3.335624e-170 1.0 True False
M 9 eye_x 1998 a_full_with_replacement b_mean -0.370 0.310 -468.5 9.524408e-168 1.0 True False
M 9 eye_x 1998 b_third_without a_mean+2SD -0.374 0.381 -381.8 5.848472e-159 1.0 True False
M 9 eye_x 1998 b_third_without b_mean -0.374 0.311 -346.3 9.144491e-155 1.0 True False
M 9 eye_y 1998 a_full_with_replacement a_mean+2SD -0.159 0.348 -290.6 3.102508e-147 1.0 True False
M 9 eye_y 1998 a_full_with_replacement b_mean -0.159 0.310 -268.9 6.682622e-144 1.0 True False
M 9 eye_y 1998 b_third_without a_mean+2SD -0.156 0.381 -224.9 3.163314e-136 1.0 True False
M 9 eye_y 1998 b_third_without b_mean -0.156 0.311 -195.3 3.589552e-130 1.0 True False
M 9 joystick 1998 a_full_with_replacement a_mean+2SD -0.452 0.350 -678.4 1.153675e-183 1.0 True False
M 9 joystick 1998 a_full_with_replacement b_mean -0.452 0.310 -645.2 1.675104e-181 1.0 True False
M 9 joystick 1998 b_third_without a_mean+2SD -0.450 0.384 -451.8 3.407206e-166 1.0 True False
M 9 joystick 1998 b_third_without b_mean -0.450 0.312 -412.7 2.646602e-162 1.0 True False
Paper: mean PI at 650 ms, horizontal eye -0.25, vertical eye -0.44, joystick -0.53 (tolerance 10%).
Readings (bootstrap x chance) with all three signals below chance at p < 1e-4, per monkey:
monkey bootstrap chance
A a_full_with_replacement a_mean+2SD True
b_mean True
b_third_without a_mean+2SD True
b_mean True
M a_full_with_replacement a_mean+2SD True
b_mean True
b_third_without a_mean+2SD True
b_mean True
Any signal above chance at p < 0.05 under any reading: {'A': False, 'M': False}
Result: reproduced. In both monkeys the PI at 650 ms is below chance for all three signals under every bootstrap and chance reading (one-tailed t-test, p 5.7e-52 to 1.2e-183), and no signal is above chance (one-tailed p = 1.0). Monkey A (6 sessions, 2,088 trials): mean PI 0.155 (horizontal eye), 0.117 (vertical eye) and -0.401 (joystick), against chance values of 0.230 to 0.303. Monkey M (9 sessions, 1,998 trials): -0.370, -0.159 and -0.452, against 0.310 to 0.384. The mean PIs differ from the paper's -0.25, -0.44 and -0.53, by more than the 10% tolerance in every case, and the eye PIs of monkey A are positive. The clock offset of ten sessions is estimated, not given by the files: after the shift the joystick deviates from rest by 2.39 to 2.53 V after the recorded onsets, against 0.16 to 0.84 V before them and 0.73 to 1.39 V 7 s later, the pattern of the five sessions on the trial clock (0.43 to 0.70 V, 1.95 to 2.49 V and 1.08 to 1.32 V).
fig4g-reset-model-fits-better¶
A Bayesian observer model whose timing variability grows sublinearly, as with resets at multiples of 0.65 s, fits the monkeys' behaviour better than a scalar-variability model without resets.
Criterion (claims): per session, both observer models of the code's utilities (models 1 and 2:
Bayes least squares estimate of a noisy measurement under the code's prior, density 0.5846 -
0.1026 t on 0.65 to 3.25 s; measurement and production noise SD wm t and wp t_e without resets, wm
sqrt(t) and wp sqrt(t_e) with resets, as in utilities/bls_offset_*_nonuniformPrior_self.m;
offset b) fit by maximum likelihood (Nelder-Mead from [0.5, 0.5, 0], as fminsearch) to
|v_p| against |v_a| of correct MNAV trials (seqq < 3, as the prior covers 0.65 to 3.25 s). Lapses
excluded first by a Gaussian and uniform mixture per interval (Methods; the uniform's range "0 to
3 v_p" is read as 0 to 3 times the interval, an assumption); without exclusion (code) as a
variant. Model error per session: the absolute difference between the data's and the model's
summed squared bias plus variance across the five intervals (code line 118). Paired t-test of
no-reset against reset error across sessions, all sessions and, as the code does, without the
first; the code's two-sample test (line 1133) as a variant. The dandiset holds 6 and 9 sessions,
where the paper's degrees of freedom imply 79 and 103.
t = time.time()
rows = []
for s in sessions:
tr = s["trials"]
g = tr[(tr.trial_type == 3) & (tr.attempt == 1) & (tr.seqq < 3) & (tr.tp.abs() > 0.01) & (tr.tp.abs() < 99)]
ta, tp = np.round(g.ta.abs().values, 3), g.tp.abs().values
for lapse in ("Methods: lapses excluded", "code: no exclusion"):
keep = mnav.lapse_mask(ta, tp) if lapse.startswith("Methods") else np.ones(len(tp), bool)
a, p = ta[keep], tp[keep]
ts = np.unique(a)
data_bias = np.array([p[a == x].mean() - x for x in ts])
data_var = np.array([p[a == x].var() for x in ts])
row = {"monkey": s["monkey"], "session": s["session"], "lapses": lapse, "trials": int(keep.sum()), "dropped": int((~keep).sum())}
for model, reset in (("no_reset", False), ("reset", True)):
fit = mnav.bls_fit(a, p, reset)
mb, mv = mnav.bls_bias_var(fit.x, ts, reset)
row[f"{model}_wm"], row[f"{model}_wp"], row[f"{model}_b"] = np.round([abs(fit.x[0]), abs(fit.x[1]), fit.x[2]], 3)
row[f"{model}_nll"] = round(fit.fun, 1)
row[f"{model}_converged"] = bool(fit.success)
row[f"{model}_error"] = abs(np.sum(data_bias ** 2 + data_var) - np.sum(mb ** 2 + mv))
rows.append(row)
fig4g = pd.DataFrame(rows)
print(f"fits in {time.time() - t:.0f} s")
print(fig4g.round(4).to_string(index=False))
res = []
for (lapse, m), g in fig4g.groupby(["lapses", "monkey"]):
for test_name, sub in (("paired, all sessions", g), ("paired, first session left out (code)", g.sort_values("session").iloc[1:])):
tt = stats.ttest_rel(sub.no_reset_error, sub.reset_error)
res.append({"lapses": lapse, "monkey": m, "test": test_name, "sessions": len(sub),
"mean_error_no_reset": round(sub.no_reset_error.mean(), 4), "mean_error_reset": round(sub.reset_error.mean(), 4),
"reset_lower_in": f"{int((sub.reset_error < sub.no_reset_error).sum())}/{len(sub)}",
"t": round(tt.statistic, 2), "df": len(sub) - 1, "p_two_sided": tt.pvalue,
"meets(reset lower, p<1e-4)": bool(tt.statistic > 0 and tt.pvalue < 1e-4)})
tt2 = stats.ttest_ind(g.no_reset_error, g.reset_error)
res.append({"lapses": lapse, "monkey": m, "test": "two-sample (code line 1133)", "sessions": len(g), "t": round(tt2.statistic, 2),
"df": 2 * len(g) - 2, "p_two_sided": tt2.pvalue, "meets(reset lower, p<1e-4)": bool(tt2.statistic > 0 and tt2.pvalue < 1e-4)})
print(pd.DataFrame(res).to_string(index=False))
print("Log-likelihood difference (no reset minus reset NLL; positive favours reset), per monkey, median:",
(fig4g.no_reset_nll - fig4g.reset_nll).groupby([fig4g.lapses, fig4g.monkey]).median().round(1).to_dict())
# Effect size of the paired test (mean difference over its SD, d_z = t / sqrt(n)), beside the paper's
# from its t and degrees of freedom (claims file: A t(77) = 7.93, M t(101) = 16.56), and the t that a
# two-sided p < 0.0001 needs with the sessions held here.
paper_t = {"A": (7.93, 77), "M": (16.56, 101)}
for m, g in fig4g[fig4g.lapses == "Methods: lapses excluded"].groupby("monkey"):
d = g.no_reset_error - g.reset_error
t_need = stats.t.isf(0.5e-4, len(d) - 1)
print(f"monkey {m}: d_z = {d.mean() / d.std(ddof=1):.2f} over {len(d)} sessions; paper d_z = "
f"{paper_t[m][0] / np.sqrt(paper_t[m][1] + 1):.2f} over {paper_t[m][1] + 1}; "
f"t needed for p < 0.0001 with {len(d)} sessions {t_need:.1f}")
fits in 39 s
monkey session lapses trials dropped no_reset_wm no_reset_wp no_reset_b no_reset_nll no_reset_converged no_reset_error reset_wm reset_wp reset_b reset_nll reset_converged reset_error
A 08152019 Methods: lapses excluded 754 0 0.030 0.113 -0.019 -355.1 True 0.1293 0.046 0.118 -0.030 -395.3 True 0.0200
A 08152019 code: no exclusion 754 0 0.030 0.113 -0.019 -355.1 True 0.1293 0.046 0.118 -0.030 -395.3 True 0.0200
A 08212019 Methods: lapses excluded 936 0 0.022 0.118 -0.036 -406.3 True 0.1669 0.025 0.124 -0.038 -483.5 True 0.0085
A 08212019 code: no exclusion 936 0 0.022 0.118 -0.036 -406.3 True 0.1669 0.025 0.124 -0.038 -483.5 True 0.0085
A 08232019 Methods: lapses excluded 982 0 0.003 0.131 0.003 -341.8 True 0.1970 0.003 0.135 -0.027 -441.5 True 0.0211
A 08232019 code: no exclusion 982 0 0.003 0.131 0.003 -341.8 True 0.1970 0.003 0.135 -0.027 -441.5 True 0.0211
A 08272019 Methods: lapses excluded 953 0 0.004 0.120 -0.064 -411.6 True 0.1619 0.003 0.129 -0.057 -474.6 True 0.0168
A 08272019 code: no exclusion 953 0 0.004 0.120 -0.064 -411.6 True 0.1619 0.003 0.129 -0.057 -474.6 True 0.0168
A 08292019 Methods: lapses excluded 978 0 0.004 0.115 -0.046 -466.2 True 0.1254 0.003 0.126 -0.039 -505.4 True 0.0296
A 08292019 code: no exclusion 978 0 0.004 0.115 -0.046 -466.2 True 0.1254 0.003 0.126 -0.039 -505.4 True 0.0296
A 09172019 Methods: lapses excluded 610 0 0.003 0.108 -0.011 -330.7 True 0.1124 0.002 0.122 -0.005 -339.7 True 0.0167
A 09172019 code: no exclusion 610 0 0.003 0.108 -0.011 -330.7 True 0.1124 0.002 0.122 -0.005 -339.7 True 0.0167
M 04122021 Methods: lapses excluded 934 0 0.071 0.105 0.004 -508.7 True 0.1512 0.081 0.109 -0.002 -516.5 True 0.0211
M 04122021 code: no exclusion 934 0 0.071 0.105 0.004 -508.7 True 0.1512 0.081 0.109 -0.002 -516.5 True 0.0211
M 04132021 Methods: lapses excluded 665 0 0.053 0.127 -0.038 -338.0 True 0.2191 0.074 0.121 -0.055 -354.2 True 0.0004
M 04132021 code: no exclusion 665 0 0.053 0.127 -0.038 -338.0 True 0.2191 0.074 0.121 -0.055 -354.2 True 0.0004
M 04142021 Methods: lapses excluded 495 0 0.034 0.130 0.003 -240.5 True 0.1990 0.086 0.116 -0.044 -259.6 True 0.0032
M 04142021 code: no exclusion 495 0 0.034 0.130 0.003 -240.5 True 0.1990 0.086 0.116 -0.044 -259.6 True 0.0032
M 04212021 Methods: lapses excluded 830 0 0.022 0.115 -0.023 -426.0 True 0.0903 0.054 0.122 -0.040 -412.1 True 0.0423
M 04212021 code: no exclusion 830 0 0.022 0.115 -0.023 -426.0 True 0.0903 0.054 0.122 -0.040 -412.1 True 0.0423
M 04222021 Methods: lapses excluded 440 0 0.002 0.121 0.039 -264.9 True 0.1618 0.001 0.126 0.042 -267.9 True 0.0204
M 04222021 code: no exclusion 440 0 0.002 0.121 0.039 -264.9 True 0.1618 0.001 0.126 0.042 -267.9 True 0.0204
M 04252021 Methods: lapses excluded 744 0 0.064 0.101 -0.028 -328.3 True 0.0637 0.099 0.101 -0.057 -337.8 True 0.0445
M 04252021 code: no exclusion 744 0 0.064 0.101 -0.028 -328.3 True 0.0637 0.099 0.101 -0.057 -337.8 True 0.0445
M 05042021 Methods: lapses excluded 912 0 0.047 0.109 -0.040 -534.0 True 0.0968 0.081 0.104 -0.066 -536.8 True 0.0358
M 05042021 code: no exclusion 912 0 0.047 0.109 -0.040 -534.0 True 0.0968 0.081 0.104 -0.066 -536.8 True 0.0358
M 05062021 Methods: lapses excluded 731 0 0.063 0.102 -0.017 -414.0 True 0.1077 0.096 0.096 -0.046 -428.1 True 0.0240
M 05062021 code: no exclusion 731 0 0.063 0.102 -0.017 -414.0 True 0.1077 0.096 0.096 -0.046 -428.1 True 0.0240
M 05072021 Methods: lapses excluded 579 0 0.049 0.115 -0.009 -303.4 True 0.1338 0.072 0.116 -0.027 -307.7 True 0.0304
M 05072021 code: no exclusion 579 0 0.049 0.115 -0.009 -303.4 True 0.1338 0.072 0.116 -0.027 -307.7 True 0.0304
lapses monkey test sessions mean_error_no_reset mean_error_reset reset_lower_in t df p_two_sided meets(reset lower, p<1e-4)
Methods: lapses excluded A paired, all sessions 6 0.1488 0.0188 6/6 9.26 5 0.000247 False
Methods: lapses excluded A paired, first session left out (code) 5 0.1527 0.0185 5/5 8.17 4 0.001225 False
Methods: lapses excluded A two-sample (code line 1133) 6 NaN NaN NaN 9.77 10 0.000002 True
Methods: lapses excluded M paired, all sessions 9 0.1359 0.0247 9/9 4.99 8 0.001065 False
Methods: lapses excluded M paired, first session left out (code) 8 0.1340 0.0251 8/8 4.33 7 0.003427 False
Methods: lapses excluded M two-sample (code line 1133) 9 NaN NaN NaN 6.18 16 0.000013 True
code: no exclusion A paired, all sessions 6 0.1488 0.0188 6/6 9.26 5 0.000247 False
code: no exclusion A paired, first session left out (code) 5 0.1527 0.0185 5/5 8.17 4 0.001225 False
code: no exclusion A two-sample (code line 1133) 6 NaN NaN NaN 9.77 10 0.000002 True
code: no exclusion M paired, all sessions 9 0.1359 0.0247 9/9 4.99 8 0.001065 False
code: no exclusion M paired, first session left out (code) 8 0.1340 0.0251 8/8 4.33 7 0.003427 False
code: no exclusion M two-sample (code line 1133) 9 NaN NaN NaN 6.18 16 0.000013 True
Log-likelihood difference (no reset minus reset NLL; positive favours reset), per monkey, median: {('Methods: lapses excluded', 'A'): 51.6, ('Methods: lapses excluded', 'M'): 7.8, ('code: no exclusion', 'A'): 51.6, ('code: no exclusion', 'M'): 7.8}
monkey A: d_z = 3.78 over 6 sessions; paper d_z = 0.90 over 78; t needed for p < 0.0001 with 6 sessions 11.2
monkey M: d_z = 1.66 over 9 sessions; paper d_z = 1.64 over 102; t needed for p < 0.0001 with 9 sessions 7.1
Result: partially reproduced. The model with resets has the lower error in all 15 sessions (mean error 0.0188 against 0.1488 in A, 0.0247 against 0.1359 in M) and the higher likelihood in 14 (median log-likelihood difference 51.6 in A and 7.8 in M). The paired t-test across sessions does not reach p < 0.0001 in either monkey (A t(5) = 9.26, p = 2.5e-04; M t(8) = 4.99, p = 1.1e-03, two-sided), nor without the first session as the code does (p 1.2e-03 and 3.4e-03). The code's two-sample test gives p < 0.0001 in both (about 2e-06 and 1.3e-05). The mixture model excluded no trial (the trials are first-attempt trials), so the Methods' and the code's fits are the same. The dandiset holds 6 and 9 sessions, where the paper's degrees of freedom imply 79 and 103 (78 and 102 pairs in its test). The paired effect size is as large as the paper's or larger (d_z printed above), and the t that p < 0.0001 needs with 6 and 9 sessions is above the observed t. The criterion defines only the reproduced outcome. The effect is present and in the paper's direction in every session, and its test is weaker because of the number of sessions, so the outcome is partially reproduced under the protocol's definition (effect present but weaker); the verifier would give not reproduced, and that disagreement is recorded.
fig4d-model-subscalar-variability¶
The CAN model without landmarks shows scalar variability, its SD growing linearly with temporal distance, while the model with endogenous landmarks shows subscalar variability, with a power-law exponent below 1.
A simulation, not a test on recorded data. The model is ported from the paper's code
(mnav.can_*, from gc_CAN_model_simulation.m and utilities/gc1d_*.m at commit 4bee109, the code
the claims file cites). Weber fraction 0.08: the claims file writes 0.8 for the file Fig. 4c loads
(wm8), but the code names its files with 100 times the Weber fraction, and its Fig. 4c section
sets wmm = .08; a Weber fraction of 0.8 would make some velocities negative. This is asked of a
reader (questions/000897.yaml). 50 simulations with
landmarks at velocity 0.42 (1 + 0.08 N(0, 1)); the time to reach the network state that the
median landmark trajectory has at 0.65, 1.30, ..., 3.25 s (code, get_simul_variance, steps
counted as 1/20,000 s as the code does). Without landmarks the dynamics depend only on the
velocity, so a grid of noise-free velocities (0.28 to 0.465) is simulated and the times of 50
noisy velocities are interpolated from it. No-landmark velocity: per base interval by grid search
so that the mean time matches the landmark model's (Methods, main); the script's single 0.35 as a
variant. 1,000 bootstrap SDs and means per interval; SD = a mean^b + c fitted per bootstrap
(claims), and a mean^b (the code's lsqcurvefit, a variant); one-tailed t-test of b against 1.
from scipy.optimize import curve_fit
rng4d = np.random.default_rng(404)
t = time.time()
W_can = mnav.can_setup()
s0_can = mnav.can_init(W_can, rng4d)
WM, V_LM, V_NOLM = 0.08, 0.42, 0.35 # Weber fraction of Fig. 4c (file wm8) and the script's speeds
N_SIM = 50
v_lm = V_LM * (1 + WM * rng4d.standard_normal(N_SIM))
xi_nolm = rng4d.standard_normal(N_SIM)
v_grid = np.round(np.arange(0.28, 0.4651, 0.005), 3) # no-landmark velocities, without noise
traj = mnav.can_run(W_can, s0_can, np.r_[v_lm, v_grid], np.r_[np.ones(N_SIM, bool), np.zeros(len(v_grid), bool)])
print(f"{N_SIM} simulations with landmarks and {len(v_grid)} without in {time.time() - t:.0f} s")
STEP_S = 1 / 20000 # the code's conversion of steps to seconds (mnav_ec_fig4.m)
tr_lm = traj[:, :N_SIM].astype(float)
end_step = mnav.first_reach(traj[:, :N_SIM], 360)
for j in range(N_SIM):
tr_lm[int(end_step[j]):, j] = np.nan
median_traj = np.nanmedian(tr_lm, axis=1)
base_time = 0.65 * np.arange(1, 6)
tt_ = np.arange(1, len(median_traj) + 1) * STEP_S
states = np.array([median_traj[np.argmax(tt_ > b)] for b in base_time])
rt_lm = np.stack([mnav.first_reach(traj[:, :N_SIM], s_) * STEP_S for s_ in states], axis=1)
# no-landmark times as a function of velocity, interpolated for noisy velocities
grid_rt = np.stack([mnav.first_reach(traj[:, N_SIM:], s_) * STEP_S for s_ in states], axis=1) # grid x 5
def nolm_rt(v_base):
v = v_base * (1 + WM * xi_nolm)
return np.stack([np.interp(v, v_grid, grid_rt[:, i]) for i in range(5)], axis=1)
# Methods (CAN model): the no-landmark velocity set by a grid search for each base interval so that
# the mean time matches the landmark model's
v_fit = []
cand = np.arange(0.30, 0.42, 0.0005)
for i in range(5):
err = [abs(nolm_rt(c)[:, i].mean() - rt_lm[:, i].mean()) for c in cand]
v_fit.append(cand[int(np.argmin(err))])
rt_nolm_methods = np.stack([nolm_rt(v_fit[i])[:, i] for i in range(5)], axis=1)
rt_nolm_script = nolm_rt(V_NOLM)
print("Network state at the base intervals (median landmark trajectory):", states.astype(int).tolist())
print("No-landmark velocity per base interval (grid search):", np.round(v_fit, 4).tolist(),
"; velocities outside the simulated grid:", int(((V_NOLM * (1 + WM * xi_nolm) < v_grid[0]) | (V_NOLM * (1 + WM * xi_nolm) > v_grid[-1])).sum()))
def power_b(rt, with_c, n_boot=1000):
idx = rng4d.integers(0, len(rt), (n_boot, len(rt)))
sds = np.stack([rt[idx, i].std(1, ddof=1) for i in range(5)], axis=1)
means = np.stack([rt[idx, i].mean(1) for i in range(5)], axis=1)
bs = []
for j in range(n_boot):
try:
if with_c:
# started from the linear fit (b = 1), a standard choice; from [0.5, 0.5, 0] the fit
# can stop in a degenerate minimum with b near 0
k_, c_ = np.polyfit(means[j], sds[j], 1)
p, _ = curve_fit(lambda x, a, b, c: a * x ** b + c, means[j], sds[j], p0=[k_, 1.0, c_], maxfev=10000)
else:
p, _ = curve_fit(lambda x, a, b: a * x ** b, means[j], sds[j], p0=[0.5, 0.5], maxfev=10000)
bs.append(p[1])
except RuntimeError:
bs.append(np.nan)
return np.array(bs)
rows = []
for model, rt in (("landmarks", rt_lm), ("no landmarks, velocity per interval (Methods)", rt_nolm_methods),
("no landmarks, velocity 0.35 (script)", rt_nolm_script)):
for fit_name, with_c in (("a m^b + c (claims)", True), ("a m^b (code)", False)):
bs = power_b(rt, with_c)
ok = np.isfinite(bs)
tt2 = stats.ttest_1samp(bs[ok], 1, alternative="less")
rows.append({"model": model, "fit": fit_name, "mean_s": np.round(rt.mean(0), 3).tolist(), "sd_s": np.round(rt.std(0, ddof=1), 4).tolist(),
"cv_first_last": f"{rt[:, 0].std(ddof=1) / rt[:, 0].mean():.3f}, {rt[:, 4].std(ddof=1) / rt[:, 4].mean():.3f}",
"fits": int(ok.sum()), "b_mean": round(float(bs[ok].mean()), 3), "b_median": round(float(np.median(bs[ok])), 3),
"t(df)": f"{tt2.statistic:.2f} ({ok.sum() - 1})", "p_b<1": tt2.pvalue})
fig4d = pd.DataFrame(rows)
print(fig4d.to_string(index=False))
print(f"in {time.time() - t:.0f} s")
50 simulations with landmarks and 38 without in 152 s
Network state at the base intervals (median landmark trajectory): [91, 151, 213, 275, 336] No-landmark velocity per base interval (grid search): [0.327, 0.334, 0.34, 0.3425, 0.345] ; velocities outside the simulated grid: 0
model fit mean_s sd_s cv_first_last fits b_mean b_median t(df) p_b<1
landmarks a m^b + c (claims) [0.654, 1.308, 1.963, 2.619, 3.248] [0.0422, 0.0921, 0.135, 0.1781, 0.218] 0.065, 0.067 1000 0.848 0.848 -1188.85 (999) 0.000000e+00
landmarks a m^b (code) [0.654, 1.308, 1.963, 2.619, 3.248] [0.0422, 0.0921, 0.135, 0.1781, 0.218] 0.065, 0.067 1000 0.977 0.978 -292.38 (999) 0.000000e+00
no landmarks, velocity per interval (Methods) a m^b + c (claims) [0.654, 1.307, 1.963, 2.62, 3.246] [0.0472, 0.0963, 0.1464, 0.1952, 0.2414] 0.072, 0.074 1000 0.977 0.976 -79.98 (999) 0.000000e+00
no landmarks, velocity per interval (Methods) a m^b (code) [0.654, 1.307, 1.963, 2.62, 3.246] [0.0472, 0.0963, 0.1464, 0.1952, 0.2414] 0.072, 0.074 1000 1.009 1.010 30.31 (999) 1.000000e+00
no landmarks, velocity 0.35 (script) a m^b + c (claims) [0.612, 1.249, 1.907, 2.565, 3.2] [0.0453, 0.0929, 0.1421, 0.1909, 0.2378] 0.074, 0.074 1000 0.990 0.990 -611.45 (999) 0.000000e+00
no landmarks, velocity 0.35 (script) a m^b (code) [0.612, 1.249, 1.907, 2.565, 3.2] [0.0453, 0.0929, 0.1421, 0.1909, 0.2378] 0.074, 0.074 1000 1.000 1.000 -13.22 (999) 3.383583e-37
in 154 s
Result: not reproduced. With the claims' fit (SD = a mean^b + c), the model with landmarks has b = 0.848 (mean of 1,000 bootstrap fits; t(999) = -1188.9, p < 1e-300), but the model without landmarks also has b below 1: 0.977 with the velocity set per base interval (t(999) = -80.0) and 0.990 at the script's 0.35 (t(999) = -611.5), both p < 1e-300, so the second condition fails. With the code's fit (SD = a mean^b) the landmark model has b = 0.977 (t(999) = -292.4) and the no-landmark model 1.009 (not below 1, p = 1.0) with the Methods' velocities, which would meet the criterion, while at 0.35 it is 1.000 (t(999) = -13.2, p = 3.4e-37). In the simulations the coefficient of variation hardly changes with the interval in either model (0.065 to 0.067 with landmarks, 0.072 to 0.074 without), so the ported model shows little of the subscalar variability of Fig. 4d. The t-test treats 1,000 bootstrap fits of the same 50 simulations as independent, so any b that differs from 1 by a few thousandths is significant. The model is a port of the paper's MATLAB code, run here for the first time; Fig. 4d's own simulations are not available to compare.
Summary¶
| Claim | Outcome |
|---|---|
| fig1bc-produced-matches-actual | reproduced |
| fig1e-immediate-generalization | blocked (data missing) |
| fig2-periodic-proportion | not reproduced (rests on an assumption about the null and the rule) |
| fig2d-periodicity-at-065 | partially reproduced (rests on an assumption) |
| ed3-no-eye-hand-periodicity | reproduced (clock offset estimated in 10 sessions) |
| ed2b-pi-navigation-above-iti | blocked (metadata insufficient) |
| ed2c-pi-error-below-correct | not reproduced (rests on an assumption) |
| fig2e-first-peak-near-065 | partially reproduced |
| fig2f-earlier-peaks-multiples | not reproduced |
| ed4d-overshoot-longer-period | partially reproduced (rests on an assumption) |
| fig3b-onset-and-offset-coding | reproduced |
| fig3d-offset-coding-stronger | partially reproduced |
| ed4e-rate-codes-distance-over-slope | not reproduced |
| ed4fg-distance-over-landmark | partially reproduced |
| fig3e-xcorr-conserved-iti | blocked (metadata insufficient) |
| ed5-xcorr-conserved-other-epochs | reproduced |
| fig4h-periodic-fano-factor | partially reproduced (rests on an assumption) |
| ed8-fano-periodic-without-rate-periodicity | reproduced |
| ed10-motor-preparation-signal | partially reproduced (rests on an assumption) |
| fig4g-reset-model-fits-better | partially reproduced |
| fig4b-learned-landmark-input | blocked (not attempted in time) |
| fig4d-model-subscalar-variability | not reproduced (a port of the paper's model) |
| ed6d-model-fano-periodic | blocked (not attempted in time) |
| ed5efg-tighter-phase-clustering | blocked (the model's phase distribution was not simulated) |
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 dandiset holds 15 of the 32 sessions, so counts are compared as proportions. The PI null turns the Gaussian process into a firing rate by exponentiation at the group's median unit rate, which the Methods do not describe; it is built from a bank of surrogate trials and matched to each unit's number of kept trials; the code's GP null without the Poisson step is run as a variant. The classification rule is open and every reading is shown; fig2d, ed2c, ed4d, ed10 and the other claims that use periodic units are run under the main reading (the code's score against the surrogate scores) and under rules a1 and b and the GP null, and their outcomes rest on the main reading. The population regression of fig3d centres each unit across both windows, as the Methods state; centring per window and time bin, as the cited method's constant term does, is run as a variant and gives offset above onset in both monkeys, where the Methods' centring gives a reversal in monkey M. Monkey M's sequence 3 trials are left out wherever distance is used, as the code does for the PI. All units are used, including those labelled noise in monkey A. The eye and joystick clock of ten sessions is an estimate made here. The Bayesian observer models and the CAN model are ported from the paper's MATLAB code to Python.
What was found about the dataset. It holds 15 of the 32 sessions, with 1,093 units (364 and
729) against the paper's 1,478; 103 monkey A units are labelled noise. Monkey M's seqq 3 trials use vectors in
steps of 0.975 s, which the paper does not describe. The eye and joystick timestamps are on another
clock than the trials in 10 of the 15 sessions, offset by 3,326 to 1,561,450 s by the estimate
made here; the files do not give the offset. In three sessions (A 09172019, M 04132021, M
04142021) no unit's rate changes around joystick offset, and in M 04132021 the spikes stop halfway
through the trials, so a shared clock could not be established. The ITI is not marked, and the
ITI segments longer than 4 s the paper uses cannot be found. Rows with attempt 2 hold the first
attempt's joystick times, which the files do not describe. Unit labels differ between monkeys
(quality in A, only KSLabel in M), and most unit columns have no description.