Reproducing Kim, Lycke et al. (2026) from dandiset 001868¶
Paper. Kim R, Lycke R, Zolotavin P, Montes J, Xie C, Luan L (2026). Learning induces
activation-mechanism-dependent neural plasticity in an intracortical microstimulation task.
Science Advances 12, eaef0343. doi:10.1126/sciadv.aef0343.
The claims file (claims/001868/10.1126-sciadv.aef0343.yaml) was written from the version of
record (PMC13544205) with its supplement; this attempt read the same text.
Dandiset. DANDI:001868, version 0.260715.2016, pinned because the paper's data availability statement cites it (doi:10.48324/dandi.001868/0.260715.2016), as the claims file does. It is the only published version; the draft has the same 85 assets and size and was last modified at publication (2026-07-15), so there is no other version to compare.
What the dandiset holds. 85 NWB files of 12 mice. Six behavioural mice (ICMS92, 93, 98, 100, 101 and 83, the paper's Mice 1 to 6) have 55 sessions with sorted units (spike times, one mean waveform, cell type, position on the probe), trial and train tables, wheel position, and in most sessions volumetric and planar per-ROI dF/F with ROI masks. Six control mice have 30 files with survey imaging (ROI masks and binary activation per channel and current) and, in some, sorted units. The behavioural files hold no stimulation-site position in the imaging volume, no activation calls, and no identity of planar ROIs across sessions; no file holds spike waveforms in and out of stimulation or multi-channel templates.
Method details from the paper's code. Where the claims file leaves a detail open, it is taken from the code the paper links (XieLuanLab/icms-activation-plasticity at d1922d1), each cited where it is used: the threshold at 50% of the fitted function (processing/dataloader.py), the pulse-locking null drawn from 20% of the pulses (responses_v2/pulse_locked_response_metrics.py), the position of the stimulation contact, (32 - stim_channel) x 60 um, and the cortical depths of fig. S4 (longitudinal_data_utils.py), and the control pipeline (processing/control, python/fig_s7, processing/pop_coupling/pop_coupling_control.py). The code is a source of method only; no data come from it.
Random numbers. Every random draw uses its own numpy generator with a fixed seed, given in
the cell that draws.
Outside inputs. None.
| Claim | Status |
|---|---|
| fig2a-threshold-decrease | tested |
| fig2a-steepest-decline-early | tested |
| figs1e-response-times-decrease | tested |
| fig2c-count-increase-fixed-current, fig2b-example-counts, figs2ab-density-peak-200um, fig2e-density-200um-increase, fig2e-radius-increase, figs2de-behaviour-increase-comparison, fig2g-stable-at-threshold | blocked: no stimulation-site position in the imaging volume and no activation calls |
| figs2de-control-stable | blocked: no voxel size, plane spacing or channel-to-contact mapping for the control survey, so no density; the count half is run and reported |
| fig3b, fig3df, fig3i, fig3j (two), fig3lm, fig3n (two), figs6 (planar tracking) | blocked: no identity of planar ROIs across sessions |
| fig4-unit-counts | tested |
| fig4d-waveform-separation | blocked: no spike waveforms in and out of stimulation |
| figs3b, figs3c, figs3d, fig4g, fig4ij (two) | tested |
| fig5d, fig5e, figs4b, fig5fg, figs4de, fig5h | tested |
| figs5c-wheel-onset-decrease, figs5-restricted-window-consistent | tested |
| figs7-control-no-change | tested, with the paper's control code |
| figs8-tracked-units | blocked: no multi-channel templates or tracked-unit labels |
| fig6c, fig6d, fig6f, figs9a | tested |
| figs9-coupling-window-and-control | tested: the 120 ms part, and the control part with the paper's control code |
In [1]:
import os
for _v in ("OMP_NUM_THREADS", "OPENBLAS_NUM_THREADS", "MKL_NUM_THREADS"):
os.environ[_v] = "1" # many small fits; one linear-algebra thread each (protocol, Keep it fast)
import pickle
import time
import warnings
from collections import Counter, defaultdict
import numpy as np
import pandas as pd
from scipy import stats
import icms
warnings.filterwarnings("ignore", category=RuntimeWarning)
pd.set_option("display.width", 220)
pd.set_option("display.max_columns", 40)
pd.set_option("display.max_rows", 300)
started = time.time()
if os.environ.get("ICMS_DEV_PICKLE"): # development only: a cached copy of the same loaded data
sessions_all = pickle.load(open(os.environ["ICMS_DEV_PICKLE"], "rb"))
else:
sessions_all = icms.load_behaviour_sessions(threads=12)
print(f"loaded {len(sessions_all)} behavioural sessions from dandiset {icms.DANDISET} "
f"version {icms.VERSION} in {time.time() - started:.0f} s")
loaded 55 behavioural sessions from dandiset 001868 version 0.260715.2016 in 200 s
Data check¶
Each check of the description against the values is stated under its output.
In [2]:
assets = icms.list_assets()
kinds = Counter()
for a in assets:
subj = a["path"].split("/")[0][4:]
kinds[(subj in icms.BEHAVIOUR, a["path"].split("_")[-1][:-4])] += 1
print("files:", len(assets), "| behavioural mice:", sum(v for (b, _), v in kinds.items() if b),
"| control mice:", sum(v for (b, _), v in kinds.items() if not b))
for (b, k), v in sorted(kinds.items()):
print(f" {'behaviour' if b else 'control'} files of kind {k}: {v}")
tab = pd.DataFrame([{"subject": s["subject"], "mouse": icms.MOUSE[s["subject"]], "sex": s["sex"],
"date": s["date"], "day": s["day"], "week": s["week"],
"units": len(s["spikes"]), "trials": len(s["trials"]["trial_index"]),
"trains_ephys": len(icms.ephys_trains(s)["start_time"]),
"catch": int((s["trials"]["current_uA"] == 0).sum()),
"catch_hit": int(s["trials"]["is_hit"][s["trials"]["current_uA"] == 0].sum()),
"catch_rt": int((~np.isnan(s["trials"]["response_time"][s["trials"]["current_uA"] == 0])).sum()),
"currents": ",".join(str(int(c)) for c in np.unique(s["trials"]["current_uA"])),
"channels": ",".join(str(c) for c in np.unique(s["trials"]["stim_channel"]) if c),
"wheel": "wheel" in s} for s in sessions_all])
print(tab.to_string())
print("\nsessions per mouse and week:")
print(pd.crosstab(tab["mouse"], tab["week"]))
sessions = [s for s in sessions_all if s["week"] <= 4]
print(f"\nsessions in weeks 0 to 4: {len(sessions)} (the paper: 46 sessions); weeks above 4 dropped: "
f"{len(sessions_all) - len(sessions)}")
print("catch trials:", tab["catch"].sum(), "| catch trials marked hit:", tab["catch_hit"].sum(),
"| catch trials with a response time:", tab["catch_rt"].sum())
print("units by stored cell_type:", dict(Counter(ct for s in sessions for ct in s["cell_type"])))
files: 85 | behavioural mice: 55 | control mice: 30
control files of kind ecephys: 6
control files of kind ecephys+ophys: 9
control files of kind ophys: 15
behaviour files of kind behavior+ecephys: 10
behaviour files of kind behavior+ecephys+ophys: 45
subject mouse sex date day week units trials trains_ephys catch catch_hit catch_rt currents channels wheel
0 ICMS100 Mouse 4 M 2023-10-26 0 0 16 400 360 40 3 3 0,5,6,7,8,9 3,5,6 True
1 ICMS100 Mouse 4 M 2023-11-02 7 1 15 181 163 18 0 0 0,6,7,8,9 3,5,6 True
2 ICMS100 Mouse 4 M 2023-11-03 8 2 13 190 175 15 5 6 0,4,5,6,7,8,9 3,5,6 True
3 ICMS100 Mouse 4 M 2023-11-09 14 2 10 400 360 40 5 5 0,2,3,4,5,6,7,8,9 3,5,6 True
4 ICMS100 Mouse 4 M 2023-11-16 21 3 6 400 360 40 1 1 0,2,3,4,5,6,7,8,9,10 3,5,6 True
5 ICMS100 Mouse 4 M 2023-11-21 26 4 16 400 360 40 1 1 0,1,2,3,4,5,6,7,8,9,10,11,12,13 3,5,6 True
6 ICMS100 Mouse 4 M 2023-12-05 40 6 17 400 360 40 0 0 0,3,4,5,6,7,8,9,10,11,12,13,14 3,5,6 True
7 ICMS101 Mouse 5 F 2023-10-27 0 0 10 400 360 40 5 6 0,6,7,8,9,10,11,12,13 7,9,10 True
8 ICMS101 Mouse 5 F 2023-11-01 5 1 14 215 155 21 5 6 0,3,4,5,6,7,8 10,11,12 True
9 ICMS101 Mouse 5 F 2023-11-03 7 1 12 190 176 14 1 1 0,4,5,6,7,8,9 10,11,12 True
10 ICMS101 Mouse 5 F 2023-11-16 20 3 12 500 360 50 7 7 0,2,3,4,5,6 10,11,12 True
11 ICMS101 Mouse 5 F 2023-11-21 25 4 15 400 360 40 3 3 0,3,4,5,6 10,11,12 True
12 ICMS101 Mouse 5 F 2023-11-22 26 4 18 400 360 40 4 5 0,3,4,5,6,7 10,11,12 True
13 ICMS101 Mouse 5 F 2023-11-29 33 5 18 499 360 50 0 0 0,2,3,4,5,6,7,8,9,10 10,11,12 True
14 ICMS101 Mouse 5 F 2023-12-01 35 5 15 400 360 40 2 2 0,4,5,6,7,8,9 10,11,12 True
15 ICMS101 Mouse 5 F 2023-12-05 39 6 13 400 360 40 0 0 0,4,5,6,7,8,9 10,11,12 True
16 ICMS83 Mouse 6 M 2023-07-20 0 0 5 480 480 0 0 0 4,5,6,7,8,9,10 3,4,5 True
17 ICMS83 Mouse 6 M 2023-07-24 4 1 4 240 240 0 0 0 4,5,6 3,4 True
18 ICMS83 Mouse 6 M 2023-07-26 6 1 5 320 320 0 0 0 2,3,4,5,6 3,4 True
19 ICMS83 Mouse 6 M 2023-07-28 8 2 6 180 180 0 0 0 3,4,5,6 3,4 True
20 ICMS83 Mouse 6 M 2023-07-31 11 2 8 430 430 0 0 0 2,3,4,5 3,4,5 True
21 ICMS83 Mouse 6 M 2023-08-02 13 2 4 440 440 0 0 0 2,3,4,5,6 3,4,5 True
22 ICMS83 Mouse 6 M 2023-08-04 15 3 9 360 360 0 0 0 2,3,4,5 3,4,5 True
23 ICMS83 Mouse 6 M 2023-08-07 18 3 8 440 440 0 0 0 2,3,4,5 3,4,5 True
24 ICMS83 Mouse 6 M 2023-08-10 21 3 8 475 475 0 0 0 3,4,5,6,7 3,4,5 True
25 ICMS92 Mouse 1 M 2023-08-30 0 0 20 499 360 50 3 3 0,3,4,5,6,7 9,10,11,12,25,29 True
26 ICMS92 Mouse 1 M 2023-09-01 2 1 19 400 360 40 0 0 0,3,4,5,6 9,11,12 True
27 ICMS92 Mouse 1 M 2023-09-06 7 1 18 499 360 50 1 1 0,2,3,4,5,6 9,11,12 True
28 ICMS92 Mouse 1 M 2023-09-08 9 2 19 499 360 50 4 4 0,3,4,5 9,11,12 True
29 ICMS92 Mouse 1 M 2023-09-12 13 2 19 400 360 40 1 1 0,2,3,4,5 9,11,12 True
30 ICMS92 Mouse 1 M 2023-09-14 15 3 19 499 360 50 0 0 0,2,3,4,5 9,11,12 True
31 ICMS92 Mouse 1 M 2023-09-19 20 3 23 499 360 50 0 0 0,2,3,4,5 9,11,12 True
32 ICMS92 Mouse 1 M 2023-09-21 22 4 23 499 360 50 4 4 0,2,3,4,5 9,11,12 True
33 ICMS92 Mouse 1 M 2023-09-25 26 4 24 499 360 50 6 6 0,1,2,3,4 9,11,12 True
34 ICMS92 Mouse 1 M 2023-09-27 28 4 24 499 360 50 1 1 0,2,3,4 9,11,12 True
35 ICMS93 Mouse 2 M 2023-08-30 0 0 30 499 360 50 2 2 0,6,7,8,9 9,10,11,12,25,29 True
36 ICMS93 Mouse 2 M 2023-09-06 7 1 28 499 360 50 2 2 0,4,5,6,7 9,11,12 True
37 ICMS93 Mouse 2 M 2023-09-12 13 2 32 400 360 40 4 4 0,3,4,5,6 9,11,12 True
38 ICMS93 Mouse 2 M 2023-09-14 15 3 30 499 360 50 3 3 0,3,4,5,6 9,11,12 True
39 ICMS93 Mouse 2 M 2023-09-20 21 3 29 400 360 40 1 1 0,2,3,4,5,6 9,11,12 True
40 ICMS93 Mouse 2 M 2023-09-22 23 4 30 499 360 50 2 2 0,2,3,4,5,6 9,11,12 True
41 ICMS93 Mouse 2 M 2023-09-26 27 4 29 499 360 49 2 2 0,2,3,4,5,6 9,11,12 True
42 ICMS93 Mouse 2 M 2023-09-29 30 5 28 499 360 50 2 2 0,2,3,4,5,6,7,8 9,11,12 True
43 ICMS93 Mouse 2 M 2023-10-04 35 5 33 400 360 40 0 0 0,2,3,4,5,6 9,11,12 True
44 ICMS93 Mouse 2 M 2023-10-06 37 6 33 204 184 20 1 1 0,3,4,5,6,7 9,11,12 True
45 ICMS98 Mouse 3 F 2023-10-20 0 0 21 499 360 50 8 8 0,3,4,5,6,7 4,5,6 True
46 ICMS98 Mouse 3 F 2023-10-24 4 1 21 499 360 50 6 6 0,2,3,4,5,6,7 4,5,6 True
47 ICMS98 Mouse 3 F 2023-10-26 6 1 17 499 360 49 6 6 0,3,4,5,6,7,8 4,5,6 True
48 ICMS98 Mouse 3 F 2023-10-31 11 2 20 499 360 49 0 0 0,1,2,3,4,5,6 4,5,6 True
49 ICMS98 Mouse 3 F 2023-11-02 13 2 21 212 155 20 1 1 0,1,2,3,4,5,6,7,8 4,5,6 True
50 ICMS98 Mouse 3 F 2023-11-07 18 3 20 259 188 25 1 1 0,2,3,4,5,6 4,5,6 True
51 ICMS98 Mouse 3 F 2023-11-14 25 4 24 400 360 40 0 0 0,1,2,3,4,5,6 4,5,6 True
52 ICMS98 Mouse 3 F 2023-11-17 28 4 21 400 359 40 2 2 0,2,3,4,5,6,7,8,9 4,5,6 True
53 ICMS98 Mouse 3 F 2023-11-20 31 5 24 400 360 40 4 4 0,4,5,6,7,8,9,10,11 4,5,6 True
54 ICMS98 Mouse 3 F 2023-11-22 33 5 23 400 360 40 8 8 0,7,8,9,10,11 4,5,6 True
sessions per mouse and week:
week 0 1 2 3 4 5 6
mouse
Mouse 1 1 2 2 2 3 0 0
Mouse 2 1 1 1 2 2 2 1
Mouse 3 1 2 2 1 2 2 0
Mouse 4 1 1 2 1 1 0 1
Mouse 5 1 2 0 1 2 2 1
Mouse 6 1 2 3 3 0 0 0
sessions in weeks 0 to 4: 46 (the paper: 46 sessions); weeks above 4 dropped: 9
catch trials: 1890 | catch trials marked hit: 117 | catch trials with a response time: 121
units by stored cell_type: {np.str_('pyramidal'): 536, np.str_('interneuron'): 202, np.str_('Pyramidal Cell'): 48, np.str_('unknown'): 7, np.str_('Narrow Interneuron'): 2}
In [3]:
# Trains and trials: units, ranges and the shared clock
rows = []
for s in sessions_all:
st, T = s["stim"], s["trials"]
ti = dict(zip(T["trial_index"], range(len(T["trial_index"]))))
lag = np.array([st["start_time"][i] - T["start_time"][ti[k]] for i, k in enumerate(st["trial_index"])])
tr = icms.ephys_trains(s)
rows.append({"session": f'{s["subject"]} {s["date"]}', "trains": len(st["start_time"]),
"trains_70_in_ephys": len(tr["start_time"]),
"pulse_counts": ",".join(map(str, sorted(set(st["pulse_count"].tolist())))) if len(st["pulse_count"]) < 0 else
f'{st["pulse_count"].min()}-{st["pulse_count"].max()}',
"freq": ",".join(map(str, sorted(set(st["frequency_hz"].tolist())))),
"phase_us": ",".join(map(str, sorted(set(st["pulse_width_us"].tolist())))),
"stim_minus_trial_ms": f"{1000 * lag.min():.1f} to {1000 * lag.max():.1f}",
"rt_hit_s": f'{np.nanmin(T["response_time"]):.3f} to {np.nanmax(T["response_time"]):.3f}',
"ephys_end_s": round(s["ephys_end"], 1), "last_train_s": round(float(st["stop_time"].max()), 1)})
chk = pd.DataFrame(rows)
print(chk.to_string())
long_out = sum(((s["stim"]["pulse_count"] != 70) & (s["stim"]["stop_time"] < s["ephys_end"])).sum() for s in sessions_all)
print("trains that are not 70 pulses but lie inside the electrophysiology recording:", long_out)
print("hit trials with no response time:",
sum(int(np.isnan(s["trials"]["response_time"][s["trials"]["is_hit"]]).sum()) for s in sessions_all),
"| miss trials with a response time:",
sum(int((~np.isnan(s["trials"]["response_time"][~s["trials"]["is_hit"]])).sum()) for s in sessions_all))
dur = np.concatenate([1000 * (icms.ephys_trains(s)["stop_time"] - icms.ephys_trains(s)["start_time"]) for s in sessions_all])
print(f"70-pulse trains: stop - start from {dur.min():.2f} to {dur.max():.2f} ms (median {np.median(dur):.2f}); "
"70 pulses 10 ms apart span 690 ms plus the last pulse; no per-pulse times are stored")
session trains trains_70_in_ephys pulse_counts freq phase_us stim_minus_trial_ms rt_hit_s ephys_end_s last_train_s 0 ICMS100 2023-10-26 360 360 70-70 100.0 167.0 5.5 to 6.5 0.130 to 0.695 2523.9 2522.0 1 ICMS100 2023-11-02 163 163 70-70 100.0 167.0 5.5 to 8.5 0.044 to 0.584 1200.1 1196.1 2 ICMS100 2023-11-03 175 175 70-70 100.0 167.0 5.5 to 6.5 0.037 to 0.691 1200.1 1196.9 3 ICMS100 2023-11-09 360 360 70-70 100.0 167.0 5.5 to 6.8 0.030 to 0.692 2419.5 2417.4 4 ICMS100 2023-11-16 360 360 70-70 100.0 167.0 5.5 to 6.9 0.079 to 0.692 2588.1 2586.6 5 ICMS100 2023-11-21 360 360 70-70 100.0 167.0 5.5 to 10.1 0.134 to 0.695 2379.0 2371.6 6 ICMS100 2023-12-05 360 360 70-70 100.0 167.0 5.5 to 9.4 0.039 to 0.674 2464.1 2462.1 7 ICMS101 2023-10-27 360 360 70-70 100.0 167.0 56.0 to 84.9 0.091 to 0.692 2556.0 2555.1 8 ICMS101 2023-11-01 194 155 70-220 100.0 167.0 5.5 to 270.3 0.020 to 0.692 1200.0 3440.5 9 ICMS101 2023-11-03 176 176 70-70 100.0 167.0 5.5 to 9.2 0.093 to 0.679 1200.0 1196.2 10 ICMS101 2023-11-16 450 360 70-232 100.0 167.0 5.5 to 32.5 0.032 to 0.695 2642.2 3766.3 11 ICMS101 2023-11-21 360 360 70-70 100.0 167.0 5.5 to 9.6 0.035 to 0.672 2561.7 2559.6 12 ICMS101 2023-11-22 360 360 70-70 100.0 167.0 5.5 to 9.6 0.028 to 0.687 2501.0 2499.0 13 ICMS101 2023-11-29 449 360 70-208 100.0 167.0 5.5 to 32.5 0.038 to 0.688 2417.7 3359.5 14 ICMS101 2023-12-01 360 360 70-70 100.0 167.0 5.5 to 6.6 0.041 to 0.692 2538.8 2537.2 15 ICMS101 2023-12-05 360 360 70-70 100.0 167.0 5.5 to 6.7 0.082 to 0.695 2407.4 2405.4 16 ICMS83 2023-07-20 480 480 70-70 100.0 167.0 0.0 to 0.0 0.057 to 0.690 3227.7 3225.7 17 ICMS83 2023-07-24 240 240 70-70 100.0 167.0 0.0 to 0.0 0.155 to 0.594 1714.8 1712.8 18 ICMS83 2023-07-26 320 320 70-70 100.0 167.0 0.0 to 0.0 0.127 to 0.692 2347.8 2333.8 19 ICMS83 2023-07-28 180 180 70-70 100.0 167.0 0.0 to 0.0 0.145 to 0.690 1406.9 1399.5 20 ICMS83 2023-07-31 430 430 70-70 100.0 167.0 0.0 to 0.0 0.153 to 0.692 3405.1 3397.1 21 ICMS83 2023-08-02 440 440 70-70 100.0 167.0 0.0 to 0.0 0.108 to 0.689 3401.3 3399.7 22 ICMS83 2023-08-04 360 360 70-70 100.0 167.0 0.0 to 0.0 0.121 to 0.649 3015.6 3008.2 23 ICMS83 2023-08-07 440 440 70-70 100.0 167.0 0.0 to 0.0 0.057 to 0.681 3317.5 3315.9 24 ICMS83 2023-08-10 475 475 70-70 100.0 167.0 0.0 to 0.0 0.134 to 0.694 3667.7 3666.3 25 ICMS92 2023-08-30 449 360 70-242 100.0 167.0 5.5 to 32.5 0.065 to 0.664 2396.5 4379.3 26 ICMS92 2023-09-01 360 360 70-70 100.0 167.0 5.5 to 11.2 0.106 to 0.679 2268.1 2266.6 27 ICMS92 2023-09-06 449 360 70-183 100.0 167.0 5.5 to 32.5 0.036 to 0.637 2404.1 3203.6 28 ICMS92 2023-09-08 449 360 70-198 100.0 167.0 5.5 to 32.5 0.062 to 0.666 2372.1 3150.9 29 ICMS92 2023-09-12 360 360 70-70 100.0 167.0 5.5 to 37.4 0.082 to 0.694 2428.7 2427.2 30 ICMS92 2023-09-14 449 360 70-232 100.0 167.0 5.5 to 32.5 0.046 to 0.659 2384.2 3276.1 31 ICMS92 2023-09-19 449 360 70-155 100.0 167.0 5.5 to 32.5 0.092 to 0.658 2327.9 3320.8 32 ICMS92 2023-09-21 449 360 70-226 100.0 167.0 5.5 to 32.5 0.030 to 0.667 2305.1 3873.2 33 ICMS92 2023-09-25 449 360 70-238 100.0 167.0 5.5 to 32.5 0.018 to 0.694 2359.4 4122.5 34 ICMS92 2023-09-27 449 360 70-217 100.0 167.0 5.5 to 32.5 0.059 to 0.687 2291.5 3183.6 35 ICMS93 2023-08-30 449 360 70-223 100.0 167.0 5.5 to 32.5 0.047 to 0.694 2695.9 3774.5 36 ICMS93 2023-09-06 449 360 70-207 100.0 167.0 5.5 to 32.5 0.032 to 0.681 2416.9 3597.2 37 ICMS93 2023-09-12 360 360 70-70 100.0 167.0 5.5 to 6.4 0.064 to 0.685 2683.5 2681.5 38 ICMS93 2023-09-14 449 360 70-229 100.0 167.0 5.5 to 32.5 0.033 to 0.693 2385.5 3424.2 39 ICMS93 2023-09-20 360 360 70-70 100.0 167.0 5.5 to 14.2 0.028 to 0.618 2329.1 2327.0 40 ICMS93 2023-09-22 449 360 70-235 100.0 167.0 5.5 to 32.5 0.025 to 0.616 2406.4 3221.4 41 ICMS93 2023-09-26 450 360 70-232 100.0 167.0 5.5 to 32.5 0.072 to 0.666 2357.7 3255.1 42 ICMS93 2023-09-29 449 360 70-174 100.0 167.0 5.5 to 32.5 0.107 to 0.613 2324.6 3392.7 43 ICMS93 2023-10-04 360 360 70-70 100.0 167.0 5.5 to 7.1 0.116 to 0.674 2341.9 2339.8 44 ICMS93 2023-10-06 184 184 70-70 100.0 167.0 5.5 to 6.5 0.075 to 0.673 1200.1 1196.2 45 ICMS98 2023-10-20 449 360 70-229 100.0 167.0 5.5 to 32.5 0.026 to 0.631 2496.5 3352.7 46 ICMS98 2023-10-24 449 360 70-208 100.0 167.0 5.5 to 32.5 0.028 to 0.662 2404.2 3740.9 47 ICMS98 2023-10-26 450 360 70-220 100.0 167.0 5.5 to 32.5 0.073 to 0.694 2447.8 3319.8 48 ICMS98 2023-10-31 450 360 70-235 100.0 167.0 5.5 to 32.5 0.071 to 0.676 2369.8 3200.1 49 ICMS98 2023-11-02 192 155 70-229 100.0 167.0 5.5 to 32.5 0.074 to 0.654 1200.1 2340.0 50 ICMS98 2023-11-07 234 188 70-155 100.0 167.0 5.5 to 32.5 0.091 to 0.637 1200.1 2559.9 51 ICMS98 2023-11-14 360 360 70-70 100.0 167.0 5.5 to 6.7 0.078 to 0.634 2366.3 2364.8 52 ICMS98 2023-11-17 360 359 68-70 100.0 167.0 5.5 to 19.4 0.036 to 0.665 2509.1 2507.0 53 ICMS98 2023-11-20 360 360 70-70 100.0 167.0 5.5 to 6.6 0.056 to 0.651 2405.0 2403.0 54 ICMS98 2023-11-22 360 360 70-70 100.0 167.0 5.5 to 10.3 0.039 to 0.670 2457.0 2455.5 trains that are not 70 pulses but lie inside the electrophysiology recording: 1 hit trials with no response time: 3 | miss trials with a response time: 98 70-pulse trains: stop - start from 700.00 to 700.00 ms (median 700.00); 70 pulses 10 ms apart span 690 ms plus the last pulse; no per-pulse times are stored
In [4]:
# Shared clock: spikes and wheel against the train onsets, with the onsets shifted as a contrast.
rng = np.random.default_rng(100)
clock = []
for s in sessions:
tr = icms.ephys_trains(s)
on = np.sort(tr["start_time"][np.isin(tr["current_uA"], icms.CURRENTS)])
if len(on) < 10:
continue
sh = on + rng.uniform(1.5, 2.5, len(on)) # shifted into the intertrial period
allspk = np.sort(np.concatenate(s["spikes"]))
d_true = (icms.counts_in(allspk, on, 0, 0.7) - icms.counts_in(allspk, on, -0.7, 0)).mean()
d_shift = (icms.counts_in(allspk, sh, 0, 0.7) - icms.counts_in(allspk, sh, -0.7, 0)).mean()
T = s["trials"]
hit = T["is_hit"] & (T["current_uA"] > 0)
miss = ~T["is_hit"] & (T["current_uA"] > 0)
w_hit = np.nanmedian(np.abs(icms.wheel_displacement(s, T["start_time"][hit], 0, 0.7)[:, -1])) if "wheel" in s else np.nan
w_miss = np.nanmedian(np.abs(icms.wheel_displacement(s, T["start_time"][miss], 0, 0.7)[:, -1])) if "wheel" in s else np.nan
w_pre_hit = np.nanmedian(np.abs(icms.wheel_displacement(s, T["start_time"][hit] - 0.7, 0, 0.7)[:, -1])) if "wheel" in s else np.nan
clock.append({"session": f'{s["subject"]} {s["date"]}', "delta_spikes_true": d_true,
"delta_spikes_shifted": d_shift, "wheel_hit_0_700ms": w_hit, "wheel_miss_0_700ms": w_miss,
"wheel_hit_before_onset": w_pre_hit})
clock = pd.DataFrame(clock)
print(clock.round(2).to_string())
print("sessions where the train raises the summed spike count more than the shifted onsets:",
int((clock.delta_spikes_true > clock.delta_spikes_shifted).sum()), "of", len(clock))
print("sessions where hit trials move the wheel more than miss trials in the 700 ms after trial start:",
int((clock.wheel_hit_0_700ms > clock.wheel_miss_0_700ms).sum()), "of", len(clock))
session delta_spikes_true delta_spikes_shifted wheel_hit_0_700ms wheel_miss_0_700ms wheel_hit_before_onset 0 ICMS100 2023-10-26 58.51 -3.55 67.43 1.51 1.13 1 ICMS100 2023-11-02 68.05 -14.10 601.87 1.03 0.73 2 ICMS100 2023-11-03 39.50 -2.05 337.58 2.85 4.96 3 ICMS100 2023-11-09 61.86 -3.08 275.82 1.42 0.69 4 ICMS100 2023-11-16 39.49 5.98 66.80 1.67 1.15 5 ICMS100 2023-11-21 100.32 1.81 40.52 1.48 1.30 6 ICMS101 2023-10-27 68.30 2.56 37.54 1.70 1.14 7 ICMS101 2023-11-01 172.26 4.92 62.88 2.04 0.99 8 ICMS101 2023-11-03 181.50 1.51 21.67 1.64 0.80 9 ICMS101 2023-11-16 154.13 10.85 41.52 1.69 1.41 10 ICMS101 2023-11-21 210.77 1.98 29.88 1.62 1.35 11 ICMS101 2023-11-22 212.52 3.47 73.02 1.20 0.90 12 ICMS83 2023-07-20 21.59 0.04 5.95 3.38 5.03 13 ICMS83 2023-07-24 37.61 1.70 1.48 1.79 0.60 14 ICMS83 2023-07-26 46.67 -0.52 1.89 1.01 1.62 15 ICMS83 2023-07-28 60.67 -2.05 0.57 1.36 0.51 16 ICMS83 2023-07-31 67.45 -3.58 7.98 0.80 0.52 17 ICMS83 2023-08-02 45.89 -1.83 70.88 0.25 0.45 18 ICMS83 2023-08-04 69.59 -2.68 123.75 0.73 0.62 19 ICMS83 2023-08-07 58.12 -0.69 0.62 0.73 0.55 20 ICMS83 2023-08-10 57.56 -0.24 0.65 0.29 0.14 21 ICMS92 2023-08-30 74.79 -0.18 242.83 1.39 1.02 22 ICMS92 2023-09-01 86.09 -9.20 566.23 1.51 1.03 23 ICMS92 2023-09-06 74.40 -34.44 461.58 1.44 1.21 24 ICMS92 2023-09-08 79.46 -35.89 405.91 1.34 0.80 25 ICMS92 2023-09-12 104.39 -21.05 392.81 1.47 1.12 26 ICMS92 2023-09-14 95.34 -20.04 422.72 1.16 1.16 27 ICMS92 2023-09-19 157.31 -11.23 301.67 1.50 0.89 28 ICMS92 2023-09-21 121.63 -20.83 504.93 1.48 0.96 29 ICMS92 2023-09-25 93.62 -13.85 472.64 1.51 1.00 30 ICMS92 2023-09-27 138.71 -23.76 393.03 1.14 0.72 31 ICMS93 2023-08-30 52.12 14.59 105.12 2.47 1.27 32 ICMS93 2023-09-06 71.21 6.58 129.51 1.73 0.93 33 ICMS93 2023-09-12 150.88 13.72 135.55 1.51 1.12 34 ICMS93 2023-09-14 146.73 4.92 250.32 1.30 1.15 35 ICMS93 2023-09-20 152.87 2.05 399.02 0.90 0.97 36 ICMS93 2023-09-22 126.40 1.20 420.40 0.94 0.91 37 ICMS93 2023-09-26 121.64 -0.66 303.82 1.38 1.05 38 ICMS98 2023-10-20 59.48 11.80 58.94 1.69 1.29 39 ICMS98 2023-10-24 91.58 10.52 67.08 1.74 1.12 40 ICMS98 2023-10-26 69.90 7.39 101.75 1.31 0.85 41 ICMS98 2023-10-31 110.31 8.28 44.33 1.24 0.64 42 ICMS98 2023-11-02 126.73 3.40 62.20 1.21 0.97 43 ICMS98 2023-11-07 125.79 16.49 42.19 2.74 0.90 44 ICMS98 2023-11-14 113.31 -3.64 75.80 1.83 0.97 45 ICMS98 2023-11-17 93.33 0.32 35.66 1.26 1.15 sessions where the train raises the summed spike count more than the shifted onsets: 46 of 46 sessions where hit trials move the wheel more than miss trials in the 700 ms after trial start: 43 of 46
In [5]:
# Recordings deposited twice: identical trial tables or spike trains across files
sig_t, sig_u = defaultdict(list), defaultdict(list)
for s in sessions_all:
sig_t[hash(np.round(s["trials"]["start_time"], 4).tobytes())].append(f'{s["subject"]} {s["date"]}')
for i, spk in enumerate(s["spikes"]):
if len(spk) > 20:
sig_u[hash(np.round(spk[:200], 5).tobytes())].append(f'{s["subject"]} {s["date"]} unit {i}')
print("identical trial tables:", [v for v in sig_t.values() if len(v) > 1])
dup_u = [v for v in sig_u.values() if len(v) > 1]
print("identical spike trains (first 200 spikes):", len(dup_u), dup_u[:10])
identical trial tables: [] identical spike trains (first 200 spikes): 0 []
In [6]:
# Cell type: the stored label against the trough-to-peak of the stored mean waveform
# (below 0.45 ms is a narrow-waveform interneuron; Fig. 4H legend, claims common method)
def trough_to_peak_ms(w, fs=30000.0):
i = int(np.argmin(w))
j = i + int(np.argmax(w[i:]))
return (j - i) / fs * 1000.0
ct_rows = []
for s in sessions:
for i, w in enumerate(s["waveform"]):
ct_rows.append({"mouse": icms.MOUSE[s["subject"]], "stored": s["cell_type"][i],
"t2p_ms": trough_to_peak_ms(w)})
ct = pd.DataFrame(ct_rows)
ct["t2p_class"] = np.where(ct.t2p_ms < 0.45, "interneuron", "pyramidal")
print(pd.crosstab(ct.stored, ct.t2p_class))
print("trough-to-peak range (ms):", ct.t2p_ms.min().round(3), "to", ct.t2p_ms.max().round(3))
print("waveform_mean shape per unit:", sessions[0]["waveform"].shape[1], "samples (3 ms at 30 kHz), one channel")
print("unit_y_um range:", np.nanmin(np.concatenate([s["unit_y"] for s in sessions])).round(1), "to",
np.nanmax(np.concatenate([s["unit_y"] for s in sessions])).round(1), "um; electrode rel_y range:",
min(s["el_y"].min() for s in sessions), "to", max(s["el_y"].max() for s in sessions), "um")
_f = icms.open_asset(assets[[a["path"] for a in assets].index(sessions[0]["path"])]["asset_id"])
for _k in ["units/cell_type", "units/waveform_mean", "units/accepted", "units/unit_y_um",
"general/extracellular_ephys/electrodes/channel_name", "intervals/electrical_stimulation/stim_channel"]:
print(f"description of {_k}: {_f[_k].attrs.get('description', '')!r}")
_f.close()
acc = pd.DataFrame([{"session": f'{s["subject"]} {s["date"]}', "units": len(s["accepted"]),
"accepted": int(s["accepted"].sum())} for s in sessions])
print("units not accepted (weeks 0 to 4):", int((acc.units - acc.accepted).sum()), "of", int(acc.units.sum()),
"| sessions with any:", acc[acc.units > acc.accepted].to_dict("records"))
t2p_class interneuron pyramidal stored Narrow Interneuron 2 0 Pyramidal Cell 0 48 interneuron 173 29 pyramidal 43 493 unknown 0 7 trough-to-peak range (ms): 0.0 to 1.133 waveform_mean shape per unit: 90 samples (3 ms at 30 kHz), one channel unit_y_um range: -15.7 to 1899.1 um; electrode rel_y range: 0.0 to 1860.0 um
description of units/cell_type: 'Putative cell type from t2p classification: pyramidal, interneuron, or unknown' description of units/waveform_mean: 'the spike waveform mean for each spike unit' description of units/accepted: 'Whether this unit passed manual curation; True for all units if manual curation has not been performed on this session' description of units/unit_y_um: 'Estimated y position of the unit on the probe (microns)' description of general/extracellular_ephys/electrodes/channel_name: 'Original Ripple system channel number (1-32)' description of intervals/electrical_stimulation/stim_channel: 'Ripple stimulation channel number' units not accepted (weeks 0 to 4): 0 of 795 | sessions with any: []
In [7]:
# Imaging tables of the behavioural files: ROI counts, segmentation descriptions, reference images,
# mask layout (metadata and shapes only). The blocked imaging claims rest on these outputs.
from concurrent.futures import ThreadPoolExecutor
img_assets = [a for a in assets if "behavior" in a["path"] and "ophys" in a["path"]]
with ThreadPoolExecutor(12) as _ex:
IMG = pd.DataFrame(list(_ex.map(icms.imaging_metadata, img_assets))).sort_values(["subject", "date"])
IMG["mouse"] = IMG.subject.map(icms.MOUSE)
print(IMG[["mouse", "date", "Planar_rois", "Volumetric_rois", "Volumetric_mask_planes"]].to_string(index=False))
print("files with imaging:", len(IMG), "| with planar ROIs:", int((IMG.Planar_rois > 0).sum()))
print("planar ROIs per session, Mouse 1:", IMG.loc[(IMG.mouse == "Mouse 1") & (IMG.Planar_rois > 0), "Planar_rois"].tolist())
print("planar segmentation description (one file):", repr(IMG.Planar_description.dropna().iloc[0]))
print("volumetric segmentation description (one file):", repr(IMG.Volumetric_description.dropna().iloc[0]))
print("planar columns:", IMG.Planar_columns.dropna().unique().tolist(), "| volumetric columns:",
IMG.Volumetric_columns.dropna().unique().tolist())
_pl = IMG[IMG.Planar_rois > 0]
print("planar roi_id_source is the row index 0..N-1 in every session:", bool(_pl.Planar_roi_id_is_row_index.astype(bool).all()))
print("reference images stored in the segmentations: planar", sorted(_pl.Planar_reference_images.astype(int).unique().tolist()),
"| volumetric", sorted(IMG.Volumetric_reference_images.dropna().astype(int).unique().tolist()))
print("processing/ophys entries:", sorted({e for v in IMG.ophys_entries for e in v}))
print("imaging planes:", sorted({e for v in IMG.imaging_planes for e in v}))
mouse date Planar_rois Volumetric_rois Volumetric_mask_planes Mouse 4 2023-10-26 0 1117 13 Mouse 4 2023-11-02 0 982 13 Mouse 4 2023-11-09 0 1235 13 Mouse 4 2023-11-21 0 1312 13 Mouse 5 2023-10-27 0 1356 13 Mouse 5 2023-11-01 53 1398 13 Mouse 5 2023-11-03 0 1379 13 Mouse 5 2023-11-16 521 1134 13 Mouse 5 2023-11-22 0 1427 13 Mouse 5 2023-11-29 402 1091 11 Mouse 6 2023-07-20 0 1429 12 Mouse 6 2023-07-24 0 1397 13 Mouse 6 2023-07-26 0 1429 13 Mouse 6 2023-07-28 0 1761 13 Mouse 6 2023-07-31 0 1479 13 Mouse 6 2023-08-02 0 1717 13 Mouse 6 2023-08-04 0 1234 13 Mouse 6 2023-08-07 0 1657 13 Mouse 6 2023-08-10 0 1631 13 Mouse 1 2023-08-30 385 460 9 Mouse 1 2023-09-01 0 950 13 Mouse 1 2023-09-06 172 511 9 Mouse 1 2023-09-08 388 637 13 Mouse 1 2023-09-14 359 935 13 Mouse 1 2023-09-19 697 965 11 Mouse 1 2023-09-21 329 677 11 Mouse 1 2023-09-25 543 750 12 Mouse 1 2023-09-27 443 764 11 Mouse 2 2023-08-30 105 329 9 Mouse 2 2023-09-06 340 659 11 Mouse 2 2023-09-12 0 792 11 Mouse 2 2023-09-14 965 944 13 Mouse 2 2023-09-20 0 1136 13 Mouse 2 2023-09-22 663 1160 13 Mouse 2 2023-09-26 782 793 9 Mouse 2 2023-09-29 947 1228 13 Mouse 2 2023-10-04 0 1005 13 Mouse 2 2023-10-06 0 1101 13 Mouse 3 2023-10-20 508 1041 11 Mouse 3 2023-10-24 855 1102 13 Mouse 3 2023-10-26 763 1243 13 Mouse 3 2023-10-31 1019 1102 13 Mouse 3 2023-11-02 384 1094 13 Mouse 3 2023-11-07 590 1080 13 Mouse 3 2023-11-14 0 1050 13 files with imaging: 45 | with planar ROIs: 23 planar ROIs per session, Mouse 1: [385, 172, 388, 359, 697, 329, 543, 443] planar segmentation description (one file): 'Planar ROI segmentation (53 ROIs, 512x512 field of view). 2D image_mask per ROI. Note: planar ROIs are a separate segmentation from the volumetric A1 segmentation (different ROI IDs).' volumetric segmentation description (one file): 'Volumetric ROI segmentation (1117 ROIs across 13 Z-planes, 512x512 per plane). Each ROI stored as a voxel_mask; the associated ImagingPlane is the topmost plane (Z0) by convention. ROI coordinates are (x=col, y=row, z=plane_index) in 0-indexed pixel space. Source: BehaveAnalysis_03-06-2024_12-50_A1\\blockanalysis\\ROI.mat.' planar columns: ['roi_id_source,n_pixels,image_mask'] | volumetric columns: ['roi_id_source,n_voxels,voxel_mask'] planar roi_id_source is the row index 0..N-1 in every session: True reference images stored in the segmentations: planar [0] | volumetric [0] processing/ophys entries: ['DfOverF', 'ImageSegmentation_Planar', 'ImageSegmentation_Volumetric'] imaging planes: ['ImagingPlane_Planar', 'ImagingPlane_Vol_Z0', 'ImagingPlane_Vol_Z1', 'ImagingPlane_Vol_Z10', 'ImagingPlane_Vol_Z11', 'ImagingPlane_Vol_Z12', 'ImagingPlane_Vol_Z2', 'ImagingPlane_Vol_Z3', 'ImagingPlane_Vol_Z4', 'ImagingPlane_Vol_Z5', 'ImagingPlane_Vol_Z6', 'ImagingPlane_Vol_Z7', 'ImagingPlane_Vol_Z8', 'ImagingPlane_Vol_Z9']
fig2a-threshold-decrease¶
Behavioural detection thresholds fall with learning; thresholds in the late phase (weeks 2 to 4) are lower than in the early phase (weeks 0 to 1), median 5.5 against 4.2 uA.
Per channel and block, a cumulative Gaussian with free guess and lapse is fitted with psignifit 4
at its default grid to the hit proportion at each current, with the block's catch-trial false-alarm
rate as a 0 uA point; the session threshold of a channel is the mean of its block thresholds. The
threshold is the current at which the fitted function reaches 50% detection (ps_x50), as the
criterion says and as the paper's code takes it (ps.getThreshold(result, 0.5) in
processing/dataloader.py, whose default is the scaled function). psignifit's own threshold
parameter (the midpoint of the unscaled sigmoid) and a bounded maximum likelihood fit (its mu and
its x50) are shown as variants. The dataset stores no thresholds, so only refitted ones are shown.
Mouse 6 (ICMS83) has no catch trials in the files; its fits have no 0 uA point. Effect size
r = |Z| / sqrt(n1 + n2) (claims, common method, for thresholds).
In [8]:
def r_from_p(p, n1, n2):
return abs(stats.norm.ppf(p / 2)) / np.sqrt(n1 + n2)
from concurrent.futures import ProcessPoolExecutor
t0 = time.time()
items = []
for s in sessions:
for d in icms.block_data(s):
items.append(({"subject": s["subject"], "mouse": icms.MOUSE[s["subject"]], "date": s["date"],
"week": s["week"], "phase": s["phase"]}, d))
with ProcessPoolExecutor(4) as ex:
fits = list(ex.map(icms.fit_block, [d for _, d in items], chunksize=4))
thr_rows = [{**meta, **f} for (meta, _), f in zip(items, fits)]
print(f"{len(thr_rows)} block fits in {time.time() - t0:.0f} s")
blk = pd.DataFrame(thr_rows)
print("block fits:", len(blk), "| with a catch point:", int(blk.has_catch.sum()),
"| psignifit x50 range", blk.ps_x50.min().round(2), "to", blk.ps_x50.max().round(2),
"| x50 undefined (guess above 0.5):", int(blk.ps_x50.isna().sum()),
"| MLE mu at a bound (0 or 30 uA):", int(((blk.mu < 0.01) | (blk.mu > 29.99)).sum()),
"| guess range", blk.guess.min().round(3), "to", blk.guess.max().round(3),
"| lapse range", blk.lapse.min().round(3), "to", blk.lapse.max().round(3))
print("blocks per channel and session:", blk.groupby(["subject", "date", "channel"]).size().value_counts().to_dict())
thr = blk.groupby(["subject", "mouse", "date", "week", "phase", "channel"], as_index=False).agg(
ps=("ps_x50", "mean"), ps_u=("ps_threshold", "mean"), mu=("mu", "mean"), x50=("x50", "mean"),
n_blocks=("block", "size"))
THR_NAMES = {"ps": "psignifit x50 (main)", "ps_u": "psignifit unscaled threshold", "mu": "MLE mu", "x50": "MLE x50"}
print("channel-session thresholds:", len(thr), "| range (psignifit):", thr.ps.min().round(2), "to", thr.ps.max().round(2), "uA")
res_thr = {}
for name in THR_NAMES:
e, l = thr.loc[thr.phase == "early", name], thr.loc[thr.phase == "late", name]
t = icms.mwu(e, l)
res_thr[name] = t
print(f"{THR_NAMES[name]}: early median {t['med1']:.2f} uA (n={t['n1']}), late median {t['med2']:.2f} uA (n={t['n2']}), "
f"two-sided MWU p = {t['p']:.3g}, r = {r_from_p(t['p'], t['n1'], t['n2']):.2f}")
print("\nby mouse (psignifit x50):")
print(thr.groupby(["mouse", "phase"]).ps.agg(["median", "size"]).unstack().round(2))
# without Mouse 6, whose blocks lack catch trials
t6 = icms.mwu(thr.loc[(thr.phase == "early") & (thr.subject != "ICMS83"), "ps"],
thr.loc[(thr.phase == "late") & (thr.subject != "ICMS83"), "ps"])
print(f"\nwithout Mouse 6: early {t6['med1']:.2f} (n={t6['n1']}), late {t6['med2']:.2f} (n={t6['n2']}), p = {t6['p']:.3g}")
/home/runner/work/_temp/executed/.venv/lib/python3.13/site-packages/psignifit/psignifit.py:175: UserWarning: psignifit:boundWarning
The posterior marginal distribution of the parameter <threshold> can not be fully captured inside the estimation grid which is based on the prior
The probability at the lower bound (i.e. the lowest sampled value) is not near 0.
This indicates that your data is not sufficient to exclude lower thresholds than are included in the estimation space.
Either change the prior or ensure that the data is sufficient to constrain the posterior.
Refer to the paper or the manual for more info on this topic.
warnings.warn('psignifit:boundWarning\n'
/home/runner/work/_temp/executed/.venv/lib/python3.13/site-packages/psignifit/psignifit.py:182: UserWarning: psignifit:boundWarning
The posterior marginal distribution of the parameter <threshold> can not be fully captured inside the estimation grid which is based on the prior
The probability at the upper bound (i.e. the highest sampled value) is not near 0.
This indicates that your data is not sufficient to exclude higher thresholds than are included in the estimation space.
Either change the prior or ensure that the data is sufficient to constrain the posterior.
Refer to the paper or the manual for more info on this topic.
warnings.warn('psignifit:boundWarning\n'
/home/runner/work/_temp/executed/.venv/lib/python3.13/site-packages/psignifit/psignifit.py:175: UserWarning: psignifit:boundWarning
The posterior marginal distribution of the parameter <threshold> can not be fully captured inside the estimation grid which is based on the prior
The probability at the lower bound (i.e. the lowest sampled value) is not near 0.
This indicates that your data is not sufficient to exclude lower thresholds than are included in the estimation space.
Either change the prior or ensure that the data is sufficient to constrain the posterior.
Refer to the paper or the manual for more info on this topic.
warnings.warn('psignifit:boundWarning\n'
/home/runner/work/_temp/executed/.venv/lib/python3.13/site-packages/psignifit/psignifit.py:175: UserWarning: psignifit:boundWarning
The posterior marginal distribution of the parameter <threshold> can not be fully captured inside the estimation grid which is based on the prior
The probability at the lower bound (i.e. the lowest sampled value) is not near 0.
This indicates that your data is not sufficient to exclude lower thresholds than are included in the estimation space.
Either change the prior or ensure that the data is sufficient to constrain the posterior.
Refer to the paper or the manual for more info on this topic.
warnings.warn('psignifit:boundWarning\n'
/home/runner/work/_temp/executed/.venv/lib/python3.13/site-packages/psignifit/psignifit.py:182: UserWarning: psignifit:boundWarning
The posterior marginal distribution of the parameter <threshold> can not be fully captured inside the estimation grid which is based on the prior
The probability at the upper bound (i.e. the highest sampled value) is not near 0.
This indicates that your data is not sufficient to exclude higher thresholds than are included in the estimation space.
Either change the prior or ensure that the data is sufficient to constrain the posterior.
Refer to the paper or the manual for more info on this topic.
warnings.warn('psignifit:boundWarning\n'
/home/runner/work/_temp/executed/.venv/lib/python3.13/site-packages/psignifit/psignifit.py:161: UserWarning: All provided data blocks contain <= 5 trials.
Did you sample adaptively?
If so please specify a range which contains the whole psychometric function in
conf.stimulus_range.
An appropriate prior prior will be then chosen. For now we use the standard
heuristic, assuming that the psychometric function is covered by the stimulus
levels, which is frequently invalid for adaptive procedures!
warnings.warn("""All provided data blocks contain <= 5 trials.
/home/runner/work/_temp/executed/.venv/lib/python3.13/site-packages/psignifit/psignifit.py:182: UserWarning: psignifit:boundWarning
The posterior marginal distribution of the parameter <threshold> can not be fully captured inside the estimation grid which is based on the prior
The probability at the upper bound (i.e. the highest sampled value) is not near 0.
This indicates that your data is not sufficient to exclude higher thresholds than are included in the estimation space.
Either change the prior or ensure that the data is sufficient to constrain the posterior.
Refer to the paper or the manual for more info on this topic.
warnings.warn('psignifit:boundWarning\n'
456 block fits in 961 s
block fits: 456 | with a catch point: 408 | psignifit x50 range 0.15 to 12.83 | x50 undefined (guess above 0.5): 0 | MLE mu at a bound (0 or 30 uA): 3 | guess range 0.0 to 0.5 | lapse range 0.0 to 0.5
blocks per channel and session: {4: 93, 2: 38, 1: 2, 3: 2}
channel-session thresholds: 135 | range (psignifit): 0.83 to 10.8 uA
psignifit x50 (main): early median 5.39 uA (n=46), late median 3.88 uA (n=89), two-sided MWU p = 5.32e-06, r = 0.39
psignifit unscaled threshold: early median 5.41 uA (n=46), late median 3.84 uA (n=89), two-sided MWU p = 5.32e-06, r = 0.39
MLE mu: early median 5.40 uA (n=46), late median 3.71 uA (n=89), two-sided MWU p = 4.9e-07, r = 0.43
MLE x50: early median 5.46 uA (n=46), late median 3.70 uA (n=86), two-sided MWU p = 5.15e-06, r = 0.40
by mouse (psignifit x50):
median size
phase early late early late
mouse
Mouse 1 3.70 3.04 9 21
Mouse 2 5.82 4.56 6 15
Mouse 3 5.05 3.88 9 15
Mouse 4 6.94 5.45 6 12
Mouse 5 6.52 5.21 9 9
Mouse 6 5.59 3.87 7 17
without Mouse 6: early 5.17 (n=39), late 3.91 (n=72), p = 7.61e-05
Result: reproduced. With psignifit at its default grid, the median channel-session threshold (the current at 50% of the fitted function) is 5.39 uA early (n = 46) and 3.88 uA late (n = 89), both within tolerance (late 3.78 to 4.62), and late is lower (two-sided Mann-Whitney p = 5.32e-06, r = 0.39). The early n is within 43 +- 4; the late n of 89 is one above 80 +- 8, which the criterion does not judge. psignifit's unscaled threshold gives 5.41 and 3.84 uA. The bounded maximum likelihood variant, which is not the criterion's fit, puts the late median below the tolerance (mu 3.71, x50 3.70 uA) with the same direction (p = 4.9e-07, 5.15e-06). Every mouse is lower late. Without Mouse 6, whose blocks have no catch trials, early 5.17 and late 3.91 uA (p = 7.61e-05).
fig2a-steepest-decline-early¶
The steepest fall in threshold is between weeks 0 and 1, followed by smaller changes and relative stability in weeks 2 to 4.
In [9]:
wk = thr.groupby("week").ps.agg(["median", "size"])
wk["change_from_previous"] = wk["median"].diff()
print(wk.round(3))
for name in THR_NAMES:
m = thr.groupby("week")[name].median()
d = m.diff()
print(f"{THR_NAMES[name]}: week medians {m.round(3).to_dict()} | fall week 0 to 1 {-d.loc[1]:.3f} uA, "
f"largest later change {d.loc[2:].abs().max():.3f} uA (weeks {int(d.loc[2:].abs().idxmax()) - 1} to "
f"{int(d.loc[2:].abs().idxmax())})")
print("by mouse, median psignifit x50 per week:")
print(thr.pivot_table(index="mouse", columns="week", values="ps", aggfunc="median").round(2))
median size change_from_previous
week
0 6.495 18 NaN
1 4.658 28 -1.837
2 3.866 29 -0.792
3 4.310 30 0.444
4 3.808 30 -0.503
psignifit x50 (main): week medians {0: 6.495, 1: 4.658, 2: 3.866, 3: 4.31, 4: 3.808} | fall week 0 to 1 1.837 uA, largest later change 0.792 uA (weeks 1 to 2)
psignifit unscaled threshold: week medians {0: 6.407, 1: 4.632, 2: 3.785, 3: 4.323, 4: 3.826} | fall week 0 to 1 1.775 uA, largest later change 0.846 uA (weeks 1 to 2)
MLE mu: week medians {0: 5.913, 1: 4.897, 2: 3.711, 3: 3.815, 4: 3.599} | fall week 0 to 1 1.016 uA, largest later change 1.186 uA (weeks 1 to 2)
MLE x50: week medians {0: 6.081, 1: 4.942, 2: 3.685, 3: 3.976, 4: 3.633} | fall week 0 to 1 1.139 uA, largest later change 1.257 uA (weeks 1 to 2)
by mouse, median psignifit x50 per week:
week 0 1 2 3 4
mouse
Mouse 1 4.89 3.60 3.37 3.14 2.69
Mouse 2 7.22 5.08 4.56 4.68 4.29
Mouse 3 6.01 4.47 3.81 4.44 3.60
Mouse 4 7.48 5.17 5.68 5.41 5.36
Mouse 5 8.75 5.79 NaN 5.21 5.23
Mouse 6 5.59 5.07 3.81 4.57 NaN
Result: reproduced, on the criterion's reading of the figure. With the psignifit thresholds of fig2a, the median falls by 1.837 uA from week 0 to week 1; the largest later week-to-week change is 0.792 uA (week 1 to 2), and weeks 2 to 4 stay between 3.808 and 4.310 uA. psignifit's unscaled threshold gives 1.775 against 0.846 uA. The maximum likelihood variant reverses the order (mu 1.016 against 1.186 uA, x50 1.139 against 1.257 uA, both for week 1 to 2); it is not the criterion's fit, but the ordering depends on the fitting method.
figs1e-response-times-decrease¶
Response times follow the same time course as thresholds, falling over training.
Response time of hit trials at 4 to 6 uA. The dataset's response_time is measured from trial
start, which precedes the train onset by 0 to 85 ms (data check); the time from train onset is
shown as a variant.
In [10]:
rt_rows = []
for s in sessions:
T = s["trials"]
st = s["stim"]
on = dict(zip(st["trial_index"], st["start_time"]))
m = T["is_hit"] & np.isin(T["current_uA"], icms.CURRENTS) & ~np.isnan(T["response_time"])
for k, t0, rt in zip(T["trial_index"][m], T["start_time"][m], T["response_time"][m]):
rt_rows.append({"subject": s["subject"], "mouse": icms.MOUSE[s["subject"]], "week": s["week"],
"phase": s["phase"], "rt_ms": 1000 * rt,
"rt_from_train_ms": 1000 * (t0 + rt - on.get(k, np.nan))})
rt = pd.DataFrame(rt_rows)
print("hit trials at 4 to 6 uA:", len(rt), "| range", rt.rt_ms.min().round(1), "to", rt.rt_ms.max().round(1), "ms")
for col in ["rt_ms", "rt_from_train_ms"]:
per = rt.pivot_table(index="mouse", columns="week", values=col, aggfunc="median")
late = rt[rt.phase == "late"].groupby("mouse")[col].median()
med = per.median()
print(f"\n{col}: median per mouse and week")
print(per.round(1))
print("median over mice per week:", med.round(1).to_dict(), "| late phase:", round(late.median(), 1))
fall = med[0] - late.median()
print(f"fall week 0 to late: {fall:.1f} ms; fall week 0 to 1: {med[0] - med[1]:.1f} ms "
f"({100 * (med[0] - med[1]) / fall:.0f}% of the fall)" if fall > 0 else f"no fall: {fall:.1f} ms")
t = icms.mwu(rt.loc[rt.phase == "early", col], rt.loc[rt.phase == "late", col])
print(f"variant, MWU early vs late hit trials: early {t['med1']:.1f} (n={t['n1']}), late {t['med2']:.1f} "
f"(n={t['n2']}), p = {t['p']:.3g}, rank-biserial r = {t['r']:.2f}")
hit trials at 4 to 6 uA: 5976 | range 17.7 to 665.1 ms
rt_ms: median per mouse and week
week 0 1 2 3 4
mouse
Mouse 1 202.3 154.4 136.3 130.8 121.8
Mouse 2 272.4 177.5 187.9 164.9 166.9
Mouse 3 200.9 186.0 165.6 193.4 141.0
Mouse 4 408.9 212.5 215.8 195.4 309.2
Mouse 5 349.5 268.2 NaN 198.1 229.8
Mouse 6 277.4 249.3 359.9 347.5 NaN
median over mice per week: {0: 274.9, 1: 199.2, 2: 187.9, 3: 194.4, 4: 166.9} | late phase: 193.9
fall week 0 to late: 81.0 ms; fall week 0 to 1: 75.7 ms (93% of the fall)
variant, MWU early vs late hit trials: early 190.1 (n=1808), late 173.8 (n=4168), p = 7.19e-17, rank-biserial r = -0.14
rt_from_train_ms: median per mouse and week
week 0 1 2 3 4
mouse
Mouse 1 193.5 146.6 128.1 121.7 110.9
Mouse 2 253.3 169.2 182.4 153.7 157.4
Mouse 3 191.5 177.7 157.7 183.8 135.2
Mouse 4 403.0 206.8 210.1 189.7 303.1
Mouse 5 288.6 262.5 NaN 190.7 224.0
Mouse 6 277.5 249.3 359.9 347.5 NaN
median over mice per week: {0: 265.4, 1: 192.3, 2: 182.4, 3: 186.7, 4: 157.4} | late phase: 186.8
fall week 0 to late: 78.6 ms; fall week 0 to 1: 73.1 ms (93% of the fall)
variant, MWU early vs late hit trials: early 183.0 (n=1808), late 166.3 (n=4168), p = 6.91e-16, rank-biserial r = -0.13
Result: reproduced. The median over mice of the per-mouse median hit response time is 274.9 ms in week 0 and 193.9 ms over the late phase; 75.7 of the 81.0 ms fall (93%) is between weeks 0 and 1. Measured from the train onset the values are 265.4 and 186.8 ms (93% again). Over hit trials, early 190.1 ms against late 173.8 ms (p = 7.19e-17).
figs5c-wheel-onset-decrease¶
Wheel-movement onset decreases from early to late training; median onset is 122 ms in the fastest animal (Mouse 1), 135 ms over all hit trials and 129 ms in the late phase.
The stored wheel_position_processed is already low-passed at 20 Hz forward and backward (its
description). The response threshold is the 25th percentile of absolute displacement at response
time; the dataset's miss trials have no response time (98 of them do, data check), so the hit
trials define it, and the miss trials with a response time are added as a variant. Onset is
measured from the train onset, the first time |displacement| exceeds 50% of the threshold for 50 ms.
Mice 1 to 5, hit trials at 4 to 6 uA. Readings: threshold per animal and pooled.
In [11]:
wrows = []
for s in sessions:
if s["subject"] == "ICMS83":
continue
T = s["trials"]
on = dict(zip(s["stim"]["trial_index"], s["stim"]["start_time"]))
m = np.isin(T["current_uA"], icms.CURRENTS) & ~np.isnan(T["response_time"])
for k, t0, rt_, h in zip(T["trial_index"][m], T["start_time"][m], T["response_time"][m], T["is_hit"][m]):
if k not in on:
continue
wrows.append({"si": id(s), "s": s, "subject": s["subject"], "mouse": icms.MOUSE[s["subject"]],
"week": s["week"], "phase": s["phase"], "hit": bool(h), "onset_t": on[k],
"resp_after_onset": t0 + rt_ - on[k]})
W = pd.DataFrame(wrows)
disp = np.vstack([icms.wheel_displacement(r.s, [r.onset_t], 0, 1.0)[0] for r in W.itertuples()])
ridx = np.clip(np.round(W.resp_after_onset.values * 1000).astype(int), 0, disp.shape[1] - 1)
W["abs_at_resp"] = np.abs(disp[np.arange(len(W)), ridx])
print("trials whose wheel trace is missing (outside the stored signal):", int(np.isnan(disp).all(1).sum()))
print("trials:", len(W), "| hit:", int(W.hit.sum()), "| |displacement| at response, range:",
W.abs_at_resp.min().round(1), "to", W.abs_at_resp.max().round(1), "a.u.")
wheel_res = {}
for reading in ["per animal", "pooled"]:
for which in ["hit only", "hit and miss with a response time"]:
base = W[W.hit] if which == "hit only" else W
if reading == "pooled":
thr_w = {sub: np.nanpercentile(base.abs_at_resp, 25) for sub in W.subject.unique()}
else:
thr_w = base.groupby("subject").abs_at_resp.quantile(0.25).to_dict()
hits = W[W.hit].copy()
hd = disp[W.hit.values]
hits["onset_ms"] = [1000 * icms.movement_onset(hd[i:i + 1], thr_w[sub])[0]
for i, sub in enumerate(hits.subject)]
m1 = hits.loc[hits.subject == "ICMS92", "onset_ms"].median()
allm = hits.onset_ms.median()
latem = hits.loc[hits.phase == "late", "onset_ms"].median()
t = icms.mwu(hits.loc[hits.phase == "early", "onset_ms"], hits.loc[hits.phase == "late", "onset_ms"])
wheel_res[(reading, which)] = (m1, allm, latem, t)
print(f"\nthreshold {reading}, from {which}: thresholds {dict((k, round(v, 1)) for k, v in thr_w.items())}")
print(f" median onset Mouse 1 {m1:.1f} ms, all hit trials {allm:.1f} ms, late {latem:.1f} ms; "
f"no onset found in {int(hits.onset_ms.isna().sum())} of {len(hits)}")
print(f" early {t['med1']:.1f} ms (n={t['n1']}) vs late {t['med2']:.1f} ms (n={t['n2']}), "
f"p = {t['p']:.3g}, rank-biserial r = {t['r']:.2f}")
no_on = hits.onset_ms.isna().values
print(f" of the hit trials with no onset, {int((no_on & np.isnan(hd).all(1)).sum())} have no stored wheel trace "
f"and {int((no_on & ~np.isnan(hd).all(1)).sum())} have one")
print(" median per mouse:", hits.groupby("mouse").onset_ms.median().round(1).to_dict())
trials whose wheel trace is missing (outside the stored signal): 568 trials: 4789 | hit: 4772 | |displacement| at response, range: 1.6 to 8111.3 a.u.
threshold per animal, from hit only: thresholds {'ICMS100': 11.9, 'ICMS101': 11.8, 'ICMS92': 14.5, 'ICMS93': 14.0, 'ICMS98': 12.6}
median onset Mouse 1 119.0 ms, all hit trials 133.0 ms, late 128.0 ms; no onset found in 568 of 4772
early 144.0 ms (n=1307) vs late 128.0 ms (n=2897), p = 6.09e-32, rank-biserial r = -0.23
of the hit trials with no onset, 568 have no stored wheel trace and 0 have one
median per mouse: {'Mouse 1': 119.0, 'Mouse 2': 141.0, 'Mouse 3': 131.0, 'Mouse 4': 165.5, 'Mouse 5': 171.0}
threshold per animal, from hit and miss with a response time: thresholds {'ICMS100': 11.9, 'ICMS101': 11.8, 'ICMS92': 14.5, 'ICMS93': 14.0, 'ICMS98': 12.6}
median onset Mouse 1 119.0 ms, all hit trials 133.0 ms, late 128.0 ms; no onset found in 568 of 4772
early 144.0 ms (n=1307) vs late 128.0 ms (n=2897), p = 6.25e-32, rank-biserial r = -0.23
of the hit trials with no onset, 568 have no stored wheel trace and 0 have one
median per mouse: {'Mouse 1': 119.0, 'Mouse 2': 141.0, 'Mouse 3': 131.0, 'Mouse 4': 165.5, 'Mouse 5': 171.0}
threshold pooled, from hit only: thresholds {'ICMS100': np.float64(13.7), 'ICMS101': np.float64(13.7), 'ICMS92': np.float64(13.7), 'ICMS93': np.float64(13.7), 'ICMS98': np.float64(13.7)}
median onset Mouse 1 118.0 ms, all hit trials 134.0 ms, late 129.0 ms; no onset found in 568 of 4772
early 144.0 ms (n=1307) vs late 129.0 ms (n=2897), p = 1.94e-29, rank-biserial r = -0.22
of the hit trials with no onset, 568 have no stored wheel trace and 0 have one
median per mouse: {'Mouse 1': 118.0, 'Mouse 2': 141.0, 'Mouse 3': 134.0, 'Mouse 4': 174.5, 'Mouse 5': 177.0}
threshold pooled, from hit and miss with a response time: thresholds {'ICMS100': np.float64(13.7), 'ICMS101': np.float64(13.7), 'ICMS92': np.float64(13.7), 'ICMS93': np.float64(13.7), 'ICMS98': np.float64(13.7)}
median onset Mouse 1 118.0 ms, all hit trials 134.0 ms, late 129.0 ms; no onset found in 568 of 4772
early 144.0 ms (n=1307) vs late 129.0 ms (n=2897), p = 1.94e-29, rank-biserial r = -0.22
of the hit trials with no onset, 568 have no stored wheel trace and 0 have one
median per mouse: {'Mouse 1': 118.0, 'Mouse 2': 141.0, 'Mouse 3': 134.0, 'Mouse 4': 174.5, 'Mouse 5': 177.0}
Result: reproduced. With the response threshold per animal, the median onset is 119.0 ms in Mouse 1, 133.0 ms over all hit trials and 128.0 ms in the late phase, all within tolerance, and onset is earlier late than early (144.0 against 128.0 ms, n = 1307 and 2897 against the paper's 1305 and 2922; p = 6.09e-32, r = -0.23). The pooled threshold gives 118.0, 134.0 and 129.0 ms (p = 1.94e-29). The 568 of 4772 hit trials with no onset are the 568 trials with no stored wheel trace; every hit trial with a trace has an onset.
Electrophysiology: responses¶
One response is one unit at one stimulation channel and current in one session (claims, common
method). All currents are computed, since fig. S3D needs the threshold current; the other claims
use 4, 5 and 6 uA. mod50 and mod100 are the modulated calls with the code's filters and a
z cap of 50 or 100; mod_methods uses the Methods' filters alone. Pulse-locking is computed for
every response modulated under mod100 (a superset of mod50), as the paper's code computes it
(icms.pulse_locking, from processing/batch_process/postprocessing/responses_v2/
pulse_locked_response_metrics.py at d1922d1): 500 draws of 20% of the pulses for the spike
probability curve, and a null of 100 uniform placements of the spikes of one 20% draw. Pulses are
taken 10 ms apart from the train onset, as the files store train onsets only.
In [12]:
t0 = time.time()
ch_days = icms.channel_days(sessions)
resp_rows = icms.response_table(sessions, rng_seed=0, with_trials=True)
for r in resp_rows:
r["mod50"] = icms.modulated(r, 50, True, ch_days)
r["mod100"] = icms.modulated(r, 100, True, ch_days)
r["mod_methods"] = icms.modulated(r, code_filters=False)
r["mod"] = r["mod100"]
r["ctype"] = {"pyramidal": "pyramidal", "Pyramidal Cell": "pyramidal", "interneuron": "interneuron",
"Narrow Interneuron": "interneuron"}.get(str(r["cell_type"]), "unknown")
icms.add_pulse_locking(resp_rows, sessions, rng_seed=1)
R = pd.DataFrame([{k: v for k, v in r.items() if k not in ("trial_delta", "trial_hit", "starts")} for r in resp_rows])
R["mouse"] = R.subject.map(icms.MOUSE)
# the criterion's rule applied to the stored mean waveform (variant; see the data check)
R["t2p_type"] = ["interneuron" if trough_to_peak_ms(sessions[si]["waveform"][u]) < 0.45 else "pyramidal"
for si, u in zip(R.si, R.unit)]
print(f"responses: {len(R)} in {time.time() - t0:.0f} s; modulated mod50 {int(R.mod50.sum())}, mod100 "
f"{int(R.mod100.sum())}, Methods filters only {int(R.mod_methods.sum())}")
print("channels dropped by the 5-day filter:", sorted(k for k, v in ch_days.items() if v < 5),
"| counting days over all 55 sessions instead:", sorted(k for k, v in icms.channel_days(sessions_all).items() if v < 5))
print("responses by current:", R.current.value_counts().sort_index().to_dict())
print("z range (mod50):", R.loc[R.mod50, "z"].min().round(2), "to", R.loc[R.mod50, "z"].max().round(2),
"| time to max range (ms):", R.loc[R.mod50, "tmax_ms"].min(), "to", R.loc[R.mod50, "tmax_ms"].max())
print("responses with 5 trials or fewer:", int((R.n_trials <= 5).sum()), "| z >= 50 among Methods-modulated:",
int((R.mod_methods & (R.z >= 50)).sum()))
print("PL among mod50 at 4-6 uA:", int(R.loc[R.mod50 & R.current.isin(icms.CURRENTS), "pl"].sum()), "of",
int((R.mod50 & R.current.isin(icms.CURRENTS)).sum()))
print("modulated (mod50, 4-6 uA) by mouse and phase:")
print(pd.crosstab(R.loc[R.mod50 & R.current.isin(icms.CURRENTS), "mouse"], R.loc[R.mod50 & R.current.isin(icms.CURRENTS), "phase"]))
responses: 8879 in 101 s; modulated mod50 3417, mod100 3419, Methods filters only 3522
channels dropped by the 5-day filter: [('ICMS101', 7), ('ICMS101', 9)] | counting days over all 55 sessions instead: [('ICMS101', 7), ('ICMS101', 9)]
responses by current: {1.0: 105, 2.0: 832, 3.0: 1629, 4.0: 1926, 5.0: 1645, 6.0: 1291, 7.0: 693, 8.0: 408, 9.0: 205, 10.0: 47, 11.0: 36, 12.0: 36, 13.0: 26}
z range (mod50): 0.36 to 49.06 | time to max range (ms): 5.0 to 695.0000000000001
responses with 5 trials or fewer: 197 | z >= 50 among Methods-modulated: 3
PL among mod50 at 4-6 uA: 1413 of 2267
modulated (mod50, 4-6 uA) by mouse and phase:
phase early late
mouse
Mouse 1 279 381
Mouse 2 65 294
Mouse 3 197 275
Mouse 4 47 105
Mouse 5 130 295
Mouse 6 57 142
In [13]:
# The pulse-locking call: how it treats synthetic responses, and where the PL calls peak
rng_syn = np.random.default_rng(50)
def synth_pulses(shape, rate, n=2100):
"""n pulses (30 trains of 70) with Poisson spikes in the unblanked 1.5 to 9.5 ms, shaped by `shape`."""
c = (np.arange(20) + 0.5) * 0.5
valid = (c >= 1.5) & (c <= 9.5)
w = np.zeros(20)
w[valid] = shape(c[valid])
w /= w.sum()
tot = rng_syn.poisson(rate * n)
M = np.zeros((n, 20), np.int32)
np.add.at(M, (rng_syn.integers(0, n, tot), rng_syn.choice(20, tot, p=w)), 1)
return M
shapes = {"uniform (no locking)": lambda c: np.ones_like(c), "monotone ramp, no peak": lambda c: c,
"peak at 4 ms (sd 0.5 ms)": lambda c: 1 + 3 * np.exp(-0.5 * ((c - 4) / 0.5) ** 2)}
syn = []
for name, f in shapes.items():
for rate in [0.02, 0.1, 0.3]:
calls = [icms.pulse_locking(synth_pulses(f, rate), rng_syn)["pl"] for _ in range(200)]
syn.append({"spike train": name, "spikes per pulse": rate, "fraction called PL": np.mean(calls)})
print("fraction of 200 synthetic responses called PL (nominal false-positive rate of the call: 1%):")
print(pd.DataFrame(syn).pivot(index="spike train", columns="spikes per pulse", values="fraction called PL"))
Q = R[R.mod50 & R.current.isin(icms.CURRENTS) & R.phase.notna()].copy()
Q["spikes_per_pulse"] = Q.n_pl_spikes / (Q.n_trials * 70)
print("\nmodulated responses at 4-6 uA, unblanked spikes per pulse: quartiles",
Q.spikes_per_pulse.quantile([0.25, 0.5, 0.75]).round(3).tolist())
print(f"PL: {int(Q.pl.sum())} of {len(Q)} ({Q.pl.mean():.3f})")
print("PL calls by the time of the peak of the mean curve (ms; 1.5 is the inserted point, 9.25 the last unblanked bin):")
print(Q[Q.pl.astype(bool)].peak_ms.value_counts().sort_index().to_dict())
late_peak = Q.pl.astype(bool) & (Q.peak_ms >= 8.75)
print(f"PL calls peaking in the last two unblanked bins (8.75, 9.25 ms): {int(late_peak.sum())}; "
f"PL fraction without them: {(Q.pl.astype(bool) & ~late_peak).sum() / len(Q):.3f}")
print("PL latency (mean bootstrap peak time) range:", Q[Q.pl.astype(bool)].latency_ms.min().round(2), "to",
Q[Q.pl.astype(bool)].latency_ms.max().round(2), "ms")
R["pl_strict"] = R.pl.fillna(False).astype(bool) & (R.peak_ms < 8.75) # variant: PL calls that do not peak at the end
fraction of 200 synthetic responses called PL (nominal false-positive rate of the call: 1%):
spikes per pulse 0.02 0.10 0.30
spike train
monotone ramp, no peak 0.0 0.045 0.705
peak at 4 ms (sd 0.5 ms) 0.0 0.775 1.000
uniform (no locking) 0.0 0.000 0.000
modulated responses at 4-6 uA, unblanked spikes per pulse: quartiles [0.09, 0.167, 0.298]
PL: 1413 of 2267 (0.623)
PL calls by the time of the peak of the mean curve (ms; 1.5 is the inserted point, 9.25 the last unblanked bin):
{1.5: 15, 2.25: 66, 2.75: 240, 3.25: 244, 3.75: 246, 4.25: 143, 4.75: 158, 5.25: 99, 5.75: 48, 6.25: 53, 6.75: 26, 7.25: 37, 7.75: 19, 8.25: 12, 8.75: 6, 9.25: 1}
PL calls peaking in the last two unblanked bins (8.75, 9.25 ms): 7; PL fraction without them: 0.620
PL latency (mean bootstrap peak time) range: 1.5 to 8.69 ms
The call is conservative on uniform spikes: none of 600 synthetic responses without locking is called PL, below the nominal 1%. It still calls a monotone ramp PL when spikes are dense (0.705 of ramps at 0.3 spikes per pulse), but only 7 of the 1413 PL calls on the data peak in the last two unblanked bins, and without them the PL fraction is 0.620 instead of 0.623. The PL calls peak mostly at 2.75 to 4.75 ms. Round 1 built the null from the spikes of all pulses, which narrowed it and called 80% of responses PL; the paper's n imply about 62%.
fig4-unit-counts¶
Electrophysiology recorded 973 units in six animals (672 putative pyramidal, 301 interneurons), of which 484 were significantly modulated by the stimulation train (311 pyramidal, 173 interneurons).
Units are counted per session (reading a); the files carry no identity of a unit across
sessions, so reading (b), unique units per animal, cannot be formed. Cell type is the stored
cell_type (ICMS83 uses other spellings, which are the same labels), with the trough-to-peak of
the stored mean waveform as a variant.
In [14]:
def unit_counts(sel_sessions, mod_col, currents):
rows = []
for si, s in enumerate(sessions):
if s not in sel_sessions:
continue
sub = R[(R.si == si) & R.current.isin(currents)]
modu = set(sub.loc[sub[mod_col], "unit"])
for u in range(len(s["spikes"])):
ct_ = {"pyramidal": "pyramidal", "Pyramidal Cell": "pyramidal", "interneuron": "interneuron",
"Narrow Interneuron": "interneuron"}.get(str(s["cell_type"][u]), "unknown")
rows.append({"mouse": icms.MOUSE[s["subject"]], "ctype": ct_, "mod": u in modu,
"accepted": bool(s["accepted"][u]),
"t2p": "interneuron" if trough_to_peak_ms(s["waveform"][u]) < 0.45 else "pyramidal"})
return pd.DataFrame(rows)
uc_main = None
for label, sel, mod_col, cur in [("weeks 0-4, mod50, 4-6 uA (main)", sessions, "mod50", icms.CURRENTS),
("weeks 0-4, Methods filters, 4-6 uA", sessions, "mod_methods", icms.CURRENTS),
("weeks 0-4, mod50, all currents", sessions, "mod50", tuple(sorted(R.current.unique())))]:
U = unit_counts(sel, mod_col, cur)
if uc_main is None:
uc_main = U
print(f"{label}: units {len(U)}, pyramidal {int((U.ctype == 'pyramidal').sum())}, interneurons "
f"{int((U.ctype == 'interneuron').sum())}, unknown {int((U.ctype == 'unknown').sum())}; modulated "
f"{int(U['mod'].sum())} (pyramidal {int((U['mod'] & (U.ctype == 'pyramidal')).sum())}, interneurons "
f"{int((U['mod'] & (U.ctype == 'interneuron')).sum())})")
U = uc_main
print("variant, cell type by the trough-to-peak rule on the stored mean waveform: all", U.t2p.value_counts().to_dict(),
"| modulated", U.loc[U["mod"], "t2p"].value_counts().to_dict())
UA_ = U[U.accepted]
print(f"variant, accepted units only: units {len(UA_)}, pyramidal {int((UA_.ctype == 'pyramidal').sum())}, interneurons "
f"{int((UA_.ctype == 'interneuron').sum())}; modulated {int(UA_['mod'].sum())}")
# variant: all 55 sessions, weeks 5 and 6 included
extra = [s for s in sessions_all if s["week"] > 4]
rx = icms.response_table(extra, rng_seed=5, currents=icms.CURRENTS)
ch_days_all = icms.channel_days(sessions_all)
mod_x = {(id(extra[r["si"]]), r["unit"]) for r in rx if icms.modulated(r, 50, True, ch_days_all)}
ux = pd.DataFrame([{"ctype": {"pyramidal": "pyramidal", "Pyramidal Cell": "pyramidal", "interneuron": "interneuron",
"Narrow Interneuron": "interneuron"}.get(str(s["cell_type"][u]), "unknown"),
"mod": (id(s), u) in mod_x} for s in extra for u in range(len(s["spikes"]))])
UA = pd.concat([U[["ctype", "mod"]], ux])
print(f"variant, all 55 sessions: units {len(UA)}, pyramidal {int((UA.ctype == 'pyramidal').sum())}, interneurons "
f"{int((UA.ctype == 'interneuron').sum())}; modulated {int(UA['mod'].sum())} (pyramidal "
f"{int((UA['mod'] & (UA.ctype == 'pyramidal')).sum())}, interneurons {int((UA['mod'] & (UA.ctype == 'interneuron')).sum())})")
print("per mouse:", U.groupby("mouse").size().to_dict())
t2p_all = np.array([trough_to_peak_ms(w) for s in sessions for w in s["waveform"]])
h, e = np.histogram(t2p_all, bins=np.arange(0, 1.2, 0.05))
print("trough-to-peak histogram (0.05 ms bins from 0):", h.tolist())
# trough-to-peak takes whole samples at 30 kHz (1/30 ms), so a 0.05 ms bin holds one or two possible values
n_vals, _ = np.histogram(np.arange(40) / 30000.0 * 1000.0, bins=e)
print("possible values per bin:", n_vals.tolist())
print("units per possible value, by bin (count / values):",
[round(a / b, 1) if b else None for a, b in zip(h.tolist(), n_vals.tolist())])
_v, _c = np.unique(np.round(t2p_all, 4), return_counts=True)
print("units at each value from 0.30 to 0.70 ms:", dict(zip(np.round(_v[(_v >= 0.3) & (_v <= 0.7)], 3).tolist(),
_c[(_v >= 0.3) & (_v <= 0.7)].tolist())))
weeks 0-4, mod50, 4-6 uA (main): units 795, pyramidal 584, interneurons 204, unknown 7; modulated 488 (pyramidal 339, interneurons 146)
weeks 0-4, Methods filters, 4-6 uA: units 795, pyramidal 584, interneurons 204, unknown 7; modulated 488 (pyramidal 339, interneurons 146)
weeks 0-4, mod50, all currents: units 795, pyramidal 584, interneurons 204, unknown 7; modulated 516 (pyramidal 363, interneurons 150)
variant, cell type by the trough-to-peak rule on the stored mean waveform: all {'pyramidal': 577, 'interneuron': 218} | modulated {'pyramidal': 334, 'interneuron': 154}
variant, accepted units only: units 795, pyramidal 584, interneurons 204; modulated 488
variant, all 55 sessions: units 999, pyramidal 732, interneurons 260; modulated 585 (pyramidal 405, interneurons 177)
per mouse: {'Mouse 1': 208, 'Mouse 2': 208, 'Mouse 3': 165, 'Mouse 4': 76, 'Mouse 5': 81, 'Mouse 6': 57}
trough-to-peak histogram (0.05 ms bins from 0): [1, 0, 0, 0, 30, 45, 39, 44, 59, 22, 44, 128, 61, 150, 71, 47, 37, 9, 2, 0, 4, 1, 1]
possible values per bin: [2, 1, 2, 1, 2, 2, 1, 1, 2, 1, 2, 2, 1, 2, 1, 1, 2, 1, 2, 1, 2, 1, 2]
units per possible value, by bin (count / values): [0.5, 0.0, 0.0, 0.0, 15.0, 22.5, 39.0, 44.0, 29.5, 22.0, 22.0, 64.0, 61.0, 75.0, 71.0, 47.0, 18.5, 9.0, 1.0, 0.0, 2.0, 1.0, 0.5]
units at each value from 0.30 to 0.70 ms: {0.3: 30, 0.333: 39, 0.367: 44, 0.4: 41, 0.433: 18, 0.467: 22, 0.5: 24, 0.533: 20, 0.567: 54, 0.6: 74, 0.633: 61, 0.667: 66, 0.7: 84}
Result: not reproduced. In weeks 0 to 4 the files hold 795 unit-sessions (584 pyramidal, 204 interneurons, 7 unknown), of which 488 are modulated (339 pyramidal, 146 interneurons): the modulated total and the modulated pyramidal count are within tolerance, the other four counts are below it, 2 of 6 where the criterion asks 5. Over all 55 sessions there are 999 units (732 and 260) and 585 modulated (405 and 177), 3 of 6. With the trough-to-peak rule applied to the stored mean waveforms, 577 pyramidal and 218 interneurons, 334 and 154 modulated. Every unit is marked accepted, so restricting to accepted units changes nothing. Reading (b) cannot be formed. The trough-to-peak histogram in 0.05 ms bins alternates (59, 22, 44 at 0.40 to 0.55 ms; 128, 61, 150 at 0.55 to 0.70 ms), but at 30 kHz a bin holds one or two possible values, and per value the counts at 0.40 to 0.55 ms are 29.5, 22 and 22, so the dip in the 0.45 to 0.50 ms bin comes from the binning. No bimodality is judged.
figs3b-modulation-increase-all¶
The modulation z-score of all modulated responses at fixed currents increases from the early to the late phase at 4, 5 and 6 uA.
In [15]:
def phase_test(df, col, label, alternative="two-sided"):
out = {}
for c in icms.CURRENTS:
d = df[df.current == c]
t = icms.mwu(d.loc[d.phase == "early", col], d.loc[d.phase == "late", col], alternative)
out[c] = t
print(f"{label} {int(c)} uA: early median {t['med1']:.3f} (n={t['n1']}), late median {t['med2']:.3f} "
f"(n={t['n2']}), p = {t['p']:.3g}, r = {t['r']:.2f}")
return out
M50 = R[R.mod50 & R.current.isin(icms.CURRENTS) & R.phase.notna()]
res_s3b = phase_test(M50, "z", "z")
print("variant, Methods filters only:")
_ = phase_test(R[R.mod_methods & R.current.isin(icms.CURRENTS) & R.phase.notna()], "z", "z")
print("by mouse (median z, early / late, all three currents):")
print(M50.pivot_table(index="mouse", columns="phase", values="z", aggfunc="median").round(2))
z 4 uA: early median 2.706 (n=191), late median 3.765 (n=738), p = 9.57e-05, r = 0.18 z 5 uA: early median 3.399 (n=284), late median 3.875 (n=498), p = 0.00656, r = 0.12 z 6 uA: early median 3.559 (n=300), late median 4.570 (n=256), p = 0.000404, r = 0.17 variant, Methods filters only: z 4 uA: early median 2.824 (n=195), late median 3.765 (n=738), p = 0.00032, r = 0.17 z 5 uA: early median 3.402 (n=285), late median 3.875 (n=498), p = 0.00853, r = 0.11 z 6 uA: early median 3.592 (n=309), late median 4.570 (n=256), p = 0.00055, r = 0.17 by mouse (median z, early / late, all three currents): phase early late mouse Mouse 1 2.76 3.52 Mouse 2 3.32 3.11 Mouse 3 2.84 3.63 Mouse 4 3.23 3.52 Mouse 5 5.68 4.70 Mouse 6 4.21 6.64
Result: reproduced. The modulation z-score of modulated responses is higher late than early at 4 uA (2.706 against 3.765, p = 9.57e-05, r = 0.18), 5 uA (3.399 against 3.875, p = 0.00656, r = 0.12) and 6 uA (3.559 against 4.570, p = 0.000404, r = 0.17). The response counts (191/738, 284/498, 300/256) are close to the paper's (189/729, 283/489, 303/250).
figs3c-latency-all¶
Time to maximum firing rate of all modulated responses decreases at 5 uA but not significantly at 4 or 6 uA.
In [16]:
res_s3c = phase_test(M50, "tmax_ms", "time to max (ms)")
time to max (ms) 4 uA: early median 75.000 (n=191), late median 55.000 (n=738), p = 0.0465, r = -0.09 time to max (ms) 5 uA: early median 80.000 (n=284), late median 45.000 (n=498), p = 0.00522, r = -0.12 time to max (ms) 6 uA: early median 105.000 (n=300), late median 75.000 (n=256), p = 0.168, r = -0.07
Result: partially reproduced. Time to maximum falls at 5 uA (80 to 45 ms, p = 0.00522) and does not change significantly at 6 uA (p = 0.168), as the paper reports, but it also falls significantly at 4 uA (75 to 55 ms, p = 0.0465, r = -0.09), where the paper reports p = 0.18.
figs3d-modulation-stable-at-threshold¶
Modulation at each session's detection threshold does not change over training.
The session threshold of a channel (psignifit, fig2a) rounded to the nearest uA picks the current.
In [17]:
thr_cur = {(r.subject, r.date, r.channel): int(np.round(r.ps)) for r in thr.itertuples()}
R["thr_cur"] = [thr_cur.get((a, b, c), -1) for a, b, c in zip(R.subject, R.date, R.channel)]
AT = R[R.mod50 & (R.current == R.thr_cur) & R.phase.notna()]
per_cs = AT.groupby(["subject", "date", "channel", "phase"], as_index=False).z.median()
ta = icms.mwu(per_cs.loc[per_cs.phase == "early", "z"], per_cs.loc[per_cs.phase == "late", "z"])
tb = icms.mwu(AT.loc[AT.phase == "early", "z"], AT.loc[AT.phase == "late", "z"])
print("threshold currents used:", AT.current.value_counts().sort_index().to_dict())
print(f"(a) channel-session medians: early {ta['med1']:.2f} (n={ta['n1']}), late {ta['med2']:.2f} (n={ta['n2']}), "
f"p = {ta['p']:.3g}, r = {ta['r']:.2f}")
print(f"(b) all responses: early {tb['med1']:.2f} (n={tb['n1']}), late {tb['med2']:.2f} (n={tb['n2']}), "
f"p = {tb['p']:.3g}, r = {tb['r']:.2f}")
# x50 threshold reading
thr_cur2 = {(r.subject, r.date, r.channel): int(np.round(r.x50)) if not np.isnan(r.x50) else -1 for r in thr.itertuples()}
AT2 = R[R.mod50 & (R.current == [thr_cur2.get((a, b, c), -1) for a, b, c in zip(R.subject, R.date, R.channel)]) & R.phase.notna()]
p2 = AT2.groupby(["subject", "date", "channel", "phase"], as_index=False).z.median()
t2 = icms.mwu(p2.loc[p2.phase == "early", "z"], p2.loc[p2.phase == "late", "z"])
print(f"(a) with the MLE x50 threshold: early n={t2['n1']}, late n={t2['n2']}, p = {t2['p']:.3g}, r = {t2['r']:.2f}")
threshold currents used: {2.0: 52, 3.0: 245, 4.0: 243, 5.0: 227, 6.0: 112, 7.0: 64, 8.0: 18, 9.0: 14, 11.0: 5}
(a) channel-session medians: early 3.93 (n=42), late 3.30 (n=87), p = 0.16, r = -0.15
(b) all responses: early 3.02 (n=299), late 3.12 (n=681), p = 0.789, r = 0.01
(a) with the MLE x50 threshold: early n=40, late n=81, p = 0.702, r = -0.04
Result: reproduced. At the current nearest each session's psignifit threshold, the channel-session median z-scores do not differ between phases (reading a: early 3.93, n = 42; late 3.30, n = 87; p = 0.16, r = -0.15), nor do all modulated responses (reading b, p = 0.789); with the maximum likelihood x50 threshold, p = 0.702. The n of reading (a) are about twice the paper's 21 and 34, so the criterion's choice of (a) as the reading whose n matches does not hold here.
fig4g-modulated-count-mouse1¶
The number of modulated units rose moderately in some animals (Mouse 1 at 4 uA), with no significant change at group level.
In [18]:
cnt = (R[(R.current == 4) & R.phase.notna()].groupby(["subject", "mouse", "date", "channel", "phase"], as_index=False)
.mod50.sum().rename(columns={"mod50": "n_mod"}))
m1 = cnt[cnt.subject == "ICMS92"]
t = icms.mwu(m1.loc[m1.phase == "early", "n_mod"], m1.loc[m1.phase == "late", "n_mod"], "less")
print(f"Mouse 1, 4 uA: early median {t['med1']} (n={t['n1']}), late median {t['med2']} (n={t['n2']}), "
f"one-sided p = {t['p']:.3g}, r = {t['r']:.2f}")
res_fig4g = t
tg = icms.mwu(cnt.loc[cnt.phase == "early", "n_mod"], cnt.loc[cnt.phase == "late", "n_mod"], "less")
print(f"all six mice, 4 uA: early median {tg['med1']} (n={tg['n1']}), late median {tg['med2']} (n={tg['n2']}), "
f"one-sided p = {tg['p']:.3g}, r = {tg['r']:.2f}")
for sub in icms.BEHAVIOUR:
d = cnt[cnt.subject == sub]
tt = icms.mwu(d.loc[d.phase == "early", "n_mod"], d.loc[d.phase == "late", "n_mod"], "less")
print(f" {icms.MOUSE[sub]}: early mean {d.loc[d.phase == 'early', 'n_mod'].mean():.2f} (n={tt['n1']}), "
f"late mean {d.loc[d.phase == 'late', 'n_mod'].mean():.2f} (n={tt['n2']}), one-sided p = {tt['p']:.3g}")
Mouse 1, 4 uA: early median 12.0 (n=8), late median 14.0 (n=21), one-sided p = 0.0443, r = 0.42 all six mice, 4 uA: early median 7.0 (n=26), late median 9.0 (n=84), one-sided p = 0.0728, r = 0.19 Mouse 1: early mean 12.00 (n=8), late mean 13.76 (n=21), one-sided p = 0.0443 Mouse 2: early mean 6.00 (n=1), late mean 7.87 (n=15), one-sided p = 0.255 Mouse 3: early mean 8.17 (n=6), late mean 9.92 (n=13), one-sided p = 0.221 Mouse 4: early mean nan (n=0), late mean 3.78 (n=9), one-sided p = nan Mouse 5: early mean 4.40 (n=5), late mean 10.89 (n=9), one-sided p = 0.00232 Mouse 6: early mean 3.00 (n=6), late mean 4.12 (n=17), one-sided p = 0.0165
Result: reproduced. In Mouse 1 at 4 uA the number of modulated responses per channel and session is higher late than early (median 12 against 14, n = 8 and 21 as in the paper, one-sided p = 0.0443, r = 0.42). Over all six mice the same test is not significant (p = 0.0728).
fig4ij-pyramidal-modulation-increase¶
Increases in modulation over training are carried by pyramidal cells, at 4, 5 and 6 uA, with latency decreasing at 5 and 6 uA.
z cap 100 (claims). The code's one-sided tests are a variant.
In [19]:
M100 = R[R.mod100 & R.current.isin(icms.CURRENTS) & R.phase.notna()]
PYR = M100[M100.ctype == "pyramidal"]
res_pyr_z = phase_test(PYR, "z", "pyramidal z")
res_pyr_t = phase_test(PYR, "tmax_ms", "pyramidal time to max (ms)")
print("variant, one-sided (late greater for z, late smaller for time to max):")
_ = phase_test(PYR, "z", "pyramidal z", "less")
_ = phase_test(PYR, "tmax_ms", "pyramidal time to max", "greater")
print("variant, cell type by the trough-to-peak rule on the stored mean waveform:")
_ = phase_test(M100[M100.t2p_type == "pyramidal"], "z", "pyramidal (t2p) z")
_ = phase_test(M100[M100.t2p_type == "pyramidal"], "tmax_ms", "pyramidal (t2p) time to max")
pyramidal z 4 uA: early median 2.443 (n=134), late median 3.913 (n=511), p = 4.79e-05, r = 0.23 pyramidal z 5 uA: early median 3.157 (n=194), late median 3.945 (n=334), p = 0.000587, r = 0.18 pyramidal z 6 uA: early median 3.528 (n=208), late median 4.444 (n=176), p = 0.000878, r = 0.20 pyramidal time to max (ms) 4 uA: early median 75.000 (n=134), late median 35.000 (n=511), p = 0.0208, r = -0.13 pyramidal time to max (ms) 5 uA: early median 65.000 (n=194), late median 25.000 (n=334), p = 0.000329, r = -0.18 pyramidal time to max (ms) 6 uA: early median 85.000 (n=208), late median 35.000 (n=176), p = 0.025, r = -0.13 variant, one-sided (late greater for z, late smaller for time to max): pyramidal z 4 uA: early median 2.443 (n=134), late median 3.913 (n=511), p = 2.4e-05, r = 0.23 pyramidal z 5 uA: early median 3.157 (n=194), late median 3.945 (n=334), p = 0.000294, r = 0.18 pyramidal z 6 uA: early median 3.528 (n=208), late median 4.444 (n=176), p = 0.000439, r = 0.20 pyramidal time to max 4 uA: early median 75.000 (n=134), late median 35.000 (n=511), p = 0.0104, r = -0.13 pyramidal time to max 5 uA: early median 65.000 (n=194), late median 25.000 (n=334), p = 0.000165, r = -0.18 pyramidal time to max 6 uA: early median 85.000 (n=208), late median 35.000 (n=176), p = 0.0125, r = -0.13 variant, cell type by the trough-to-peak rule on the stored mean waveform: pyramidal (t2p) z 4 uA: early median 2.519 (n=128), late median 3.784 (n=510), p = 0.000212, r = 0.21 pyramidal (t2p) z 5 uA: early median 3.133 (n=191), late median 4.103 (n=334), p = 0.000435, r = 0.18 pyramidal (t2p) z 6 uA: early median 3.451 (n=208), late median 4.547 (n=180), p = 0.000492, r = 0.21 pyramidal (t2p) time to max 4 uA: early median 45.000 (n=128), late median 35.000 (n=510), p = 0.0434, r = -0.11 pyramidal (t2p) time to max 5 uA: early median 45.000 (n=191), late median 25.000 (n=334), p = 0.000362, r = -0.18 pyramidal (t2p) time to max 6 uA: early median 85.000 (n=208), late median 35.000 (n=180), p = 0.0587, r = -0.11
Result: reproduced. Pyramidal modulation rises at 4, 5 and 6 uA (p = 4.79e-05, 0.000587, 0.000878) and time to maximum falls at 5 uA (p = 0.000329) and 6 uA (p = 0.025). It also falls at 4 uA (p = 0.0208), where the paper reports p = 0.32. With cell types from the trough-to-peak rule on the stored mean waveform, modulation rises at all three currents (p = 0.000212, 0.000435, 0.000492) and time to maximum falls at 5 uA (p = 0.000362) but not significantly at 6 uA (p = 0.0587).
fig4ij-interneurons-little-change¶
Interneurons show no significant change in modulation at any current; their latency changes significantly only at 6 uA.
In [20]:
INT = M100[M100.ctype == "interneuron"]
res_int_z = phase_test(INT, "z", "interneuron z")
res_int_t = phase_test(INT, "tmax_ms", "interneuron time to max (ms)")
print("interneuron responses by mouse:", INT.mouse.value_counts().to_dict())
print("variant, cell type by the trough-to-peak rule on the stored mean waveform:")
_ = phase_test(M100[M100.t2p_type == "interneuron"], "z", "interneuron (t2p) z")
_ = phase_test(M100[M100.t2p_type == "interneuron"], "tmax_ms", "interneuron (t2p) time to max")
interneuron z 4 uA: early median 3.383 (n=55), late median 3.519 (n=221), p = 0.331, r = 0.08
interneuron z 5 uA: early median 3.746 (n=90), late median 3.378 (n=158), p = 0.884, r = -0.01
interneuron z 6 uA: early median 3.540 (n=90), late median 5.071 (n=76), p = 0.0885, r = 0.15
interneuron time to max (ms) 4 uA: early median 105.000 (n=55), late median 85.000 (n=221), p = 0.929, r = -0.01
interneuron time to max (ms) 5 uA: early median 115.000 (n=90), late median 125.000 (n=158), p = 0.361, r = 0.07
interneuron time to max (ms) 6 uA: early median 130.000 (n=90), late median 190.000 (n=76), p = 0.0335, r = 0.19
interneuron responses by mouse: {'Mouse 3': 198, 'Mouse 1': 191, 'Mouse 2': 162, 'Mouse 5': 108, 'Mouse 4': 31}
variant, cell type by the trough-to-peak rule on the stored mean waveform:
interneuron (t2p) z 4 uA: early median 3.049 (n=63), late median 3.657 (n=228), p = 0.117, r = 0.13
interneuron (t2p) z 5 uA: early median 3.761 (n=94), late median 3.285 (n=164), p = 0.706, r = -0.03
interneuron (t2p) z 6 uA: early median 3.900 (n=92), late median 4.599 (n=76), p = 0.229, r = 0.11
interneuron (t2p) time to max 4 uA: early median 115.000 (n=63), late median 85.000 (n=228), p = 0.486, r = -0.06
interneuron (t2p) time to max 5 uA: early median 115.000 (n=94), late median 125.000 (n=164), p = 0.684, r = 0.03
interneuron (t2p) time to max 6 uA: early median 135.000 (n=92), late median 155.000 (n=76), p = 0.419, r = 0.07
Result: reproduced. Interneuron modulation does not change at 4, 5 or 6 uA (p = 0.331, 0.884, 0.0885; r = 0.08, -0.01, 0.15). Their time to maximum changes only at 6 uA (130 to 190 ms, p = 0.0335, r = 0.19; paper p = 0.010, r = 0.21), not at 4 or 5 uA (p = 0.929, 0.361). With trough-to-peak cell types, modulation does not change either (p = 0.117, 0.706, 0.229), and time to maximum changes at no current (6 uA p = 0.419).
fig5d-pulse-locking-stable-across-currents¶
Over 80% of units keep their PL or NPL identity at a given stimulation contact as the current rises from 4 to 6 uA.
In [21]:
def identity_patterns(df):
rows_ = []
for (si, ch, u), g in df.groupby(["si", "channel", "unit"]):
if len(g) < 2:
continue
seq = g.sort_values("current").pl.astype(bool).tolist()
p = ("all PL" if all(seq) else "all NPL" if not any(seq) else "NPL to PL" if seq == sorted(seq)
else "PL to NPL" if seq == sorted(seq, reverse=True) else "mixed")
rows_.append({"mouse": g.mouse.iloc[0], "pattern": p, "n_currents": len(g)})
return pd.DataFrame(rows_)
pat_rows = []
for (si, ch, u), g in M50.groupby(["si", "channel", "unit"]):
if len(g) < 2:
continue
seq = g.sort_values("current").pl.astype(bool).tolist()
if all(seq):
p = "all PL"
elif not any(seq):
p = "all NPL"
elif seq == sorted(seq):
p = "NPL to PL"
elif seq == sorted(seq, reverse=True):
p = "PL to NPL"
else:
p = "mixed"
pat_rows.append({"mouse": g.mouse.iloc[0], "pattern": p, "n_currents": len(g)})
PAT = pd.DataFrame(pat_rows)
vc = PAT.pattern.value_counts()
frac_keep = (vc.get("all PL", 0) + vc.get("all NPL", 0)) / len(PAT)
print("unit-channel pairs seen at two or more currents:", len(PAT), "| patterns:", vc.to_dict())
print(f"fraction keeping PL or NPL identity: {frac_keep:.3f}")
print("by mouse:", PAT.groupby("mouse").pattern.apply(lambda x: round(x.isin(["all PL", "all NPL"]).mean(), 3)).to_dict())
# The PL call is random (the bootstrap draws and the null). The same responses called again with three
# other seeds, for this claim and for fig5e and fig5h:
pl_reseeded = {}
t0 = time.time()
for seed in [11, 12, 13]:
rng_s = np.random.default_rng(seed)
labs = []
for r in M50.itertuples():
s = sessions[r.si]
tr = icms.ephys_trains(s)
m = (tr["stim_channel"] == r.channel) & (tr["current_uA"] == r.current)
Mx = icms.pulse_matrix(s["spikes"][r.unit], icms.pulse_starts(np.sort(tr["start_time"][m])))
labs.append(icms.pulse_locking(Mx, rng_s)["pl"])
pl_reseeded[seed] = np.array(labs)
P2 = identity_patterns(M50.assign(pl=pl_reseeded[seed]))
print(f"seed {seed}: PL {int(pl_reseeded[seed].sum())} of {len(M50)}, pairs {len(P2)}, fraction keeping identity "
f"{P2.pattern.isin(['all PL', 'all NPL']).mean():.3f}")
print(f"(in {time.time() - t0:.0f} s)")
unit-channel pairs seen at two or more currents: 740 | patterns: {'all PL': 389, 'all NPL': 197, 'PL to NPL': 89, 'NPL to PL': 49, 'mixed': 16}
fraction keeping PL or NPL identity: 0.792
by mouse: {'Mouse 1': 0.796, 'Mouse 2': 0.72, 'Mouse 3': 0.788, 'Mouse 4': 0.804, 'Mouse 5': 0.784, 'Mouse 6': 0.908}
seed 11: PL 1404 of 2267, pairs 740, fraction keeping identity 0.774
seed 12: PL 1401 of 2267, pairs 740, fraction keeping identity 0.793
seed 13: PL 1415 of 2267, pairs 740, fraction keeping identity 0.789 (in 174 s)
Result: not reproduced. Of 740 unit-channel pairs seen at two or more currents (paper 729), 389 stay PL and 197 stay NPL, a fraction of 0.792, below the 0.80 the claim states; by mouse it ranges from 0.72 to 0.908. Called again with three other seeds, the fraction is 0.774, 0.793 and 0.789, below 0.80 each time.
fig5e-pyramidal-fraction-pl¶
Pyramidal cells are more common among PL units (76.2%) than among NPL units (49.8%).
Two-step majority: over stimulation channels at each current, then over currents, ties NPL. Units of unknown type are left out of the fractions.
In [22]:
def unit_classes(df):
out_ = []
for (si, u), g in df.groupby(["si", "unit"]):
per_cur = []
for c, gc in g.groupby("current"):
k = int(gc.pl.astype(bool).sum())
per_cur.append(k > len(gc) - k)
k = sum(per_cur)
out_.append({"si": si, "unit": u, "ctype": g.ctype.iloc[0], "PL": k > len(per_cur) - k})
return pd.DataFrame(out_)
ucls = []
for (si, u), g in M50.groupby(["si", "unit"]):
per_cur = []
for c, gc in g.groupby("current"):
k = int(gc.pl.astype(bool).sum())
per_cur.append(k > len(gc) - k)
k = sum(per_cur)
ucls.append({"si": si, "unit": u, "ctype": g.ctype.iloc[0], "mouse": g.mouse.iloc[0],
"PL": k > len(per_cur) - k})
UC = pd.DataFrame(ucls)
tab5e = pd.crosstab(UC.PL, UC.ctype)
print(tab5e)
known = UC[UC.ctype != "unknown"]
f_pl = (known[known.PL].ctype == "pyramidal").mean()
f_npl = (known[~known.PL].ctype == "pyramidal").mean()
odds, p_f = stats.fisher_exact(pd.crosstab(known.PL, known.ctype).values)
print(f"pyramidal fraction: PL units {f_pl:.3f} (n={int(known.PL.sum())}), NPL units {f_npl:.3f} "
f"(n={int((~known.PL).sum())}); Fisher exact p = {p_f:.3g}, odds ratio {odds:.2f}")
# variant: cell type from the trough-to-peak of the stored waveform
t2p_of = {(si, u): ("interneuron" if trough_to_peak_ms(sessions[si]["waveform"][u]) < 0.45 else "pyramidal")
for si, u in zip(UC.si, UC.unit)}
UC["t2p"] = [t2p_of[(a, b)] for a, b in zip(UC.si, UC.unit)]
fp2, fn2 = (UC[UC.PL].t2p == "pyramidal").mean(), (UC[~UC.PL].t2p == "pyramidal").mean()
print(f"variant, trough-to-peak types: PL {fp2:.3f}, NPL {fn2:.3f}, Fisher p = "
f"{stats.fisher_exact(pd.crosstab(UC.PL, UC.t2p).values)[1]:.3g}")
for seed, labs in pl_reseeded.items():
U2 = unit_classes(M50.assign(pl=labs))
U2 = U2[U2.ctype != "unknown"]
print(f"seed {seed}: pyramidal fraction PL {(U2[U2.PL].ctype == 'pyramidal').mean():.3f} (n={int(U2.PL.sum())}), "
f"NPL {(U2[~U2.PL].ctype == 'pyramidal').mean():.3f} (n={int((~U2.PL).sum())}), Fisher p = "
f"{stats.fisher_exact(pd.crosstab(U2.PL, U2.ctype).values)[1]:.3g}")
ctype interneuron pyramidal unknown PL False 96 119 0 True 50 220 3 pyramidal fraction: PL units 0.815 (n=270), NPL units 0.553 (n=215); Fisher exact p = 7.73e-10, odds ratio 3.55 variant, trough-to-peak types: PL 0.791, NPL 0.549, Fisher p = 1.67e-08
seed 11: pyramidal fraction PL 0.813 (n=262), NPL 0.565 (n=223), Fisher p = 3.39e-09
seed 12: pyramidal fraction PL 0.817 (n=263), NPL 0.559 (n=222), Fisher p = 8.74e-10
seed 13: pyramidal fraction PL 0.826 (n=264), NPL 0.548 (n=221), Fisher p = 3.15e-11
Result: partially reproduced. Pyramidal cells are 0.815 of PL units (n = 270), within tolerance, and 0.553 of NPL units (n = 215), 0.005 above the tolerance of 0.448 to 0.548; PL units are more often pyramidal (Fisher exact p = 7.73e-10, odds ratio 3.55). With three other seeds of the PL call, the NPL fraction is 0.565, 0.559 and 0.548, and the PL fraction 0.813 to 0.826. With cell types from the trough-to-peak rule on the stored mean waveform, 0.791 and 0.549 (p = 1.67e-08).
figs4b-pl-latency-distance¶
PL spike latency increases linearly with distance from the stimulation site, giving conduction velocities of about 0.2 to 0.27 m/s; latency jitter does not vary with distance.
Distance follows the paper's code (processing/batch_process/postprocessing/longitudinal_data_utils.py,
add_data_to_df and correct_df_unit_location; python/fig_s4/generate_figure.py): the stimulation
contact sits at (32 - stim_channel) x 60 um up the probe, the unit at its estimated position
(unit_y_um); both are turned into cortical depth, (1860 - height - the animal's length of probe
outside the brain) x cos 30 degrees, responses with either depth at or above the surface are left
out, and the distance is the difference in depth. The code has these lengths for Mice 1 to 5 only, so
Mouse 6 is left out, as in the code. PL responses are the modulated ones with the z cap of 100 that
the code's figure uses (filter_modulated(max_z_score=100)). Latency readings over responses at 120 to
840 um, in 120 um bins from 60 um: (a) binned means (legend), (b) binned medians (code); jitter (c)
binned means and (d) unbinned responses.
Round 1 placed the contact at the rel_y of the electrode whose channel_name equals
stim_channel; that mapping is kept as a variant. The two disagree: the electrode table puts
channel 32 at 1860 um and the code puts stimulation channel 32 at 0 um. A check of which is right
is printed first: the depth of each animal's stimulation contacts against the Methods' 200 to 400 um,
and the fraction of modulated responses against distance under each mapping.
In [23]:
# lengths of probe outside the brain (um), from the paper's code (correct_df_unit_location)
OUTSIDE_UM = {"ICMS92": 300, "ICMS93": 180, "ICMS98": 60, "ICMS100": -300, "ICMS101": 360}
def cortical_depth(height_um, subject):
return (1860.0 - height_um - OUTSIDE_UM[subject]) * np.cos(np.deg2rad(30))
def stim_height_code(ch):
return (32 - ch) * 60.0
def stim_height_table(si, ch):
s = sessions[si]
k = np.flatnonzero(s["el_channel"] == ch)
return s["el_y"][k[0]] if len(k) else np.nan
# check of the mapping 1: depth of the stimulation contacts (Methods: 200 to 400 um)
dep = []
for sub in OUTSIDE_UM:
si0 = [i for i, s in enumerate(sessions) if s["subject"] == sub]
chs = sorted(set(R.loc[R.si.isin(si0) & R.current.isin(icms.CURRENTS), "channel"]))
for ch in chs:
dep.append({"mouse": icms.MOUSE[sub], "stim_channel": ch,
"depth_code_um": float(np.round(cortical_depth(stim_height_code(ch), sub))),
"depth_electrode_table_um": float(np.round(cortical_depth(stim_height_table(si0[0], ch), sub)))})
print(pd.DataFrame(dep).to_string())
# check of the mapping 2: modulated fraction against distance, all responses at 4-6 uA
A = R[R.current.isin(icms.CURRENTS) & R.phase.notna()].copy()
A["dist_code"] = [abs(sessions[si]["unit_y"][u] - stim_height_code(ch)) for si, u, ch in zip(A.si, A.unit, A.channel)]
A["dist_table"] = [abs(sessions[si]["unit_y"][u] - stim_height_table(si, ch)) for si, u, ch in zip(A.si, A.unit, A.channel)]
for col in ["dist_code", "dist_table"]:
print(f"\nall responses at 4-6 uA by distance along the probe, {col} (um):")
print(A.groupby(pd.cut(A[col], [0, 60, 120, 240, 480, 960, 2000]), observed=False).agg(
fraction_modulated=("mod50", "mean"), n=("mod50", "size"), median_z=("z", "median")).round(3).T)
print(" fraction modulated within 240 um and beyond 960 um, by mouse:",
{m: (round(g.loc[g[col] <= 240, "mod50"].mean(), 2), round(g.loc[g[col] > 960, "mod50"].mean(), 2))
for m, g in A.groupby("mouse")})
print("\nfraction of modulated responses by distance along the probe (code mapping, um), by mouse:")
_g = A.groupby(["mouse", pd.cut(A.dist_code, [0, 60, 120, 240, 480, 960, 2000])], observed=False).mod50
print(pd.concat({"fraction": _g.mean().round(2), "n": _g.size()}, axis=1).unstack(1).to_string())
def s4b_fits(PLR, label, verbose=True):
edges = np.arange(60, 901, 120)
centres = (edges[:-1] + edges[1:]) / 2
out = {}
for c in icms.CURRENTS:
d = PLR[(PLR.current == c) & PLR.dist.notna()]
dr = d[(d.dist >= 120) & (d.dist <= 840)]
b = np.digitize(dr.dist, edges) - 1
o = {}
for agg in ["mean", "median"]:
y = np.array([getattr(np, agg)(dr.latency_ms[b == i]) if (b == i).sum() else np.nan for i in range(len(centres))])
use = ~np.isnan(y)
lr = stats.linregress(centres[use], y[use])
o[agg] = (lr.slope, lr.pvalue, lr.rvalue ** 2, 1e-3 / lr.slope if lr.slope > 0 else np.nan)
if verbose:
print(f"{label} {int(c)} uA latency, binned {agg}s 120-840 um: slope {lr.slope:.5f} ms/um, p = {lr.pvalue:.3g}, "
f"r2 = {lr.rvalue ** 2:.2f}, velocity {o[agg][3]:.3f} m/s; bin means/medians "
f"{np.round(y, 2).tolist()}, n {[int((b == i).sum()) for i in range(len(centres))]}")
bj = np.digitize(d.dist, edges) - 1
yj = np.array([np.mean(d.jitter_ms2[bj == i]) if (bj == i).sum() else np.nan for i in range(len(centres))])
use = ~np.isnan(yj)
lj = stats.linregress(centres[use], yj[use])
lj2 = stats.linregress(d.dist, d.jitter_ms2)
o["jitter_binned"], o["jitter_unbinned"] = (lj.slope, lj.pvalue), (lj2.slope, lj2.pvalue)
if verbose:
print(f" jitter: binned means slope {lj.slope:.2e} ms2/um p = {lj.pvalue:.3g}; unbinned slope {lj2.slope:.2e} "
f"p = {lj2.pvalue:.3g}; PL responses {len(d)}, at 120-840 um {len(dr)}; latency range "
f"{d.latency_ms.min():.2f} to {d.latency_ms.max():.2f} ms")
out[c] = o
return out
P100 = M100[M100.pl.astype(bool) & (M100.subject != "ICMS83")].copy()
P100["d_unit"] = [cortical_depth(sessions[si]["unit_y"][u], sub) for si, u, sub in zip(P100.si, P100.unit, P100.subject)]
P100["d_stim"] = [cortical_depth(stim_height_code(ch), sub) for ch, sub in zip(P100.channel, P100.subject)]
inside = (P100.d_unit > 0) & (P100.d_stim > 0)
print(f"\nPL responses (z cap 100, Mice 1 to 5): {len(P100)}; left out with a depth at or above the surface: {int((~inside).sum())}")
PLR = P100[inside].copy()
PLR["dist"] = (PLR.d_unit - PLR.d_stim).abs()
print("distance range", PLR.dist.min().round(1), "to", PLR.dist.max().round(1), "um")
res_s4b = s4b_fits(PLR, "code mapping")
print("\nvariant: electrode-table mapping (round 1), distance along the probe, PL responses of all six mice")
PLT = M100[M100.pl.astype(bool)].copy()
PLT["dist"] = [abs(sessions[si]["unit_y"][u] - stim_height_table(si, ch)) for si, u, ch in zip(PLT.si, PLT.unit, PLT.channel)]
print("PL responses without a distance (stimulation contact not in the electrode table):", int(PLT.dist.isna().sum()))
res_s4b_table = s4b_fits(PLT, "table mapping")
print("\nvariant: code mapping without Mouse 5, whose modulated fraction does not fall with distance")
res_s4b_no5 = s4b_fits(PLR[PLR.subject != "ICMS101"], "code mapping, Mice 1 to 4")
print("\ncode mapping, leaving out one mouse at a time: velocity (m/s) from binned means / binned medians")
for sub in OUTSIDE_UM:
o_ = s4b_fits(PLR[PLR.subject != sub], "", verbose=False)
print(f" without {icms.MOUSE[sub]} (PL responses {int((PLR.subject != sub).sum())}): " + "; ".join(
f"{int(c)} uA {o_[c]['mean'][3]:.3f} / {o_[c]['median'][3]:.3f} (p {o_[c]['mean'][1]:.2g} / {o_[c]['median'][1]:.2g})"
for c in icms.CURRENTS))
print("PL responses per mouse at 120-840 um:", PLR[(PLR.dist >= 120) & (PLR.dist <= 840)].groupby("mouse").size().to_dict())
mouse stim_channel depth_code_um depth_electrode_table_um
0 Mouse 1 9 156.0 727.0
1 Mouse 1 11 260.0 -104.0
2 Mouse 1 12 312.0 1143.0
3 Mouse 2 9 260.0 831.0
4 Mouse 2 11 364.0 0.0
5 Mouse 2 12 416.0 1247.0
6 Mouse 3 4 104.0 1559.0
7 Mouse 3 5 156.0 1039.0
8 Mouse 3 6 208.0 624.0
9 Mouse 4 3 364.0 NaN
10 Mouse 4 5 468.0 1351.0
11 Mouse 4 6 520.0 935.0
12 Mouse 5 7 0.0 -52.0
13 Mouse 5 9 104.0 NaN
14 Mouse 5 10 156.0 260.0
15 Mouse 5 11 208.0 -156.0
16 Mouse 5 12 260.0 1091.0
all responses at 4-6 uA by distance along the probe, dist_code (um):
dist_code (0, 60] (60, 120] (120, 240] (240, 480] (480, 960] (960, 2000]
fraction_modulated 0.315 0.711 0.623 0.564 0.444 0.390
n 181.000 201.000 390.000 855.000 1533.000 1697.000
median_z -0.090 3.090 2.494 3.523 2.274 0.593
fraction modulated within 240 um and beyond 960 um, by mouse: {'Mouse 1': (np.float64(0.61), np.float64(0.51)), 'Mouse 2': (np.float64(0.6), np.float64(0.22)), 'Mouse 3': (np.float64(0.6), np.float64(0.31)), 'Mouse 4': (np.float64(0.42), np.float64(0.06)), 'Mouse 5': (np.float64(0.31), np.float64(0.48)), 'Mouse 6': (np.float64(0.08), np.float64(0.88))}
all responses at 4-6 uA by distance along the probe, dist_table (um):
dist_table (0, 60] (60, 120] (120, 240] (240, 480] (480, 960] (960, 2000]
fraction_modulated 0.426 0.299 0.431 0.503 0.513 0.465
n 338.000 261.000 645.000 876.000 1567.000 1027.000
median_z 0.905 1.283 1.278 2.372 2.077 1.757
fraction modulated within 240 um and beyond 960 um, by mouse: {'Mouse 1': (np.float64(0.59), np.float64(0.52)), 'Mouse 2': (np.float64(0.13), np.float64(0.23)), 'Mouse 3': (np.float64(0.39), np.float64(0.5)), 'Mouse 4': (np.float64(0.36), np.float64(0.06)), 'Mouse 5': (np.float64(0.68), np.float64(0.56)), 'Mouse 6': (np.float64(0.73), np.float64(0.83))}
fraction of modulated responses by distance along the probe (code mapping, um), by mouse:
fraction n
dist_code (0, 60] (60, 120] (120, 240] (240, 480] (480, 960] (960, 2000] (0, 60] (60, 120] (120, 240] (240, 480] (480, 960] (960, 2000]
mouse
Mouse 1 0.22 0.73 0.65 0.72 0.84 0.51 27 48 124 227 195 411
Mouse 2 0.39 0.76 0.62 0.30 0.03 0.22 64 68 136 220 296 556
Mouse 3 0.37 0.72 0.64 0.44 0.44 0.31 54 72 76 168 408 312
Mouse 4 0.20 0.50 0.58 0.77 0.19 0.06 25 8 31 111 181 86
Mouse 5 0.00 0.00 0.43 0.81 0.76 0.48 6 3 23 114 323 159
Mouse 6 0.20 0.00 NaN 0.00 0.35 0.88 5 2 0 15 130 173
PL responses (z cap 100, Mice 1 to 5): 1225; left out with a depth at or above the surface: 1
distance range 0.0 to 1195.1 um
code mapping 4 uA latency, binned means 120-840 um: slope 0.00591 ms/um, p = 0.000384, r2 = 0.93, velocity 0.169 m/s; bin means/medians [3.02, 3.13, 3.56, 4.24, 4.95, 5.96, 7.29], n [42, 96, 98, 75, 28, 26, 2]
code mapping 4 uA latency, binned medians 120-840 um: slope 0.00621 ms/um, p = 0.000147, r2 = 0.96, velocity 0.161 m/s; bin means/medians [2.76, 3.03, 3.57, 4.28, 4.9, 5.99, 7.29], n [42, 96, 98, 75, 28, 26, 2]
jitter: binned means slope 3.75e-04 ms2/um p = 0.148; unbinned slope -2.13e-05 p = 0.79; PL responses 524, at 120-840 um 367; latency range 1.50 to 7.99 ms
code mapping 5 uA latency, binned means 120-840 um: slope 0.00445 ms/um, p = 0.000314, r2 = 0.94, velocity 0.225 m/s; bin means/medians [3.3, 3.48, 3.72, 4.04, 5.21, 5.74, 6.27], n [31, 67, 82, 65, 29, 29, 8]
code mapping 5 uA latency, binned medians 120-840 um: slope 0.00533 ms/um, p = 0.000356, r2 = 0.94, velocity 0.188 m/s; bin means/medians [2.78, 3.18, 3.42, 3.85, 5.0, 5.36, 6.76], n [31, 67, 82, 65, 29, 29, 8]
jitter: binned means slope 2.68e-04 ms2/um p = 0.33; unbinned slope -2.31e-05 p = 0.819; PL responses 415, at 120-840 um 311; latency range 1.50 to 8.52 ms
code mapping 6 uA latency, binned means 120-840 um: slope 0.00442 ms/um, p = 0.00079, r2 = 0.91, velocity 0.226 m/s; bin means/medians [3.32, 3.78, 3.68, 4.1, 5.05, 5.58, 6.62], n [20, 45, 60, 55, 24, 30, 4]
code mapping 6 uA latency, binned medians 120-840 um: slope 0.00502 ms/um, p = 0.000486, r2 = 0.93, velocity 0.199 m/s; bin means/medians [2.85, 3.19, 3.4, 3.89, 4.82, 5.19, 6.65], n [20, 45, 60, 55, 24, 30, 4]
jitter: binned means slope -1.98e-04 ms2/um p = 0.377; unbinned slope -1.85e-04 p = 0.143; PL responses 285, at 120-840 um 238; latency range 1.60 to 7.41 ms
variant: electrode-table mapping (round 1), distance along the probe, PL responses of all six mice
PL responses without a distance (stimulation contact not in the electrode table): 32
table mapping 4 uA latency, binned means 120-840 um: slope -0.00049 ms/um, p = 0.294, r2 = 0.22, velocity nan m/s; bin means/medians [3.82, 4.16, 3.83, 3.94, 3.92, 4.06, 3.3], n [18, 60, 59, 50, 86, 73, 31]
table mapping 4 uA latency, binned medians 120-840 um: slope -0.00036 ms/um, p = 0.468, r2 = 0.11, velocity nan m/s; bin means/medians [3.44, 3.91, 3.44, 4.01, 3.66, 3.51, 3.23], n [18, 60, 59, 50, 86, 73, 31]
jitter: binned means slope 7.03e-05 ms2/um p = 0.459; unbinned slope -6.09e-05 p = 0.252; PL responses 596, at 120-840 um 377; latency range 1.50 to 7.99 ms
table mapping 5 uA latency, binned means 120-840 um: slope 0.00086 ms/um, p = 0.334, r2 = 0.19, velocity 1.161 m/s; bin means/medians [3.03, 4.43, 3.89, 4.37, 4.18, 4.52, 3.84], n [11, 64, 46, 50, 51, 68, 21]
table mapping 5 uA latency, binned medians 120-840 um: slope -0.00004 ms/um, p = 0.962, r2 = 0.00, velocity nan m/s; bin means/medians [2.99, 4.57, 3.71, 4.1, 3.83, 3.85, 3.38], n [11, 64, 46, 50, 51, 68, 21]
jitter: binned means slope 4.32e-04 ms2/um p = 0.103; unbinned slope 1.94e-06 p = 0.978; PL responses 476, at 120-840 um 311; latency range 1.50 to 8.69 ms
table mapping 6 uA latency, binned means 120-840 um: slope 0.00027 ms/um, p = 0.783, r2 = 0.02, velocity 3.770 m/s; bin means/medians [3.28, 4.49, 4.12, 3.83, 4.07, 4.76, 3.42], n [10, 54, 31, 35, 30, 28, 6]
table mapping 6 uA latency, binned medians 120-840 um: slope 0.00077 ms/um, p = 0.519, r2 = 0.09, velocity 1.294 m/s; bin means/medians [2.79, 4.55, 4.01, 3.76, 3.88, 4.85, 3.5], n [10, 54, 31, 35, 30, 28, 6]
jitter: binned means slope 2.60e-04 ms2/um p = 0.487; unbinned slope 1.48e-04 p = 0.0632; PL responses 309, at 120-840 um 194; latency range 1.60 to 8.59 ms
variant: code mapping without Mouse 5, whose modulated fraction does not fall with distance
code mapping, Mice 1 to 4 4 uA latency, binned means 120-840 um: slope 0.00558 ms/um, p = 0.00068, r2 = 0.92, velocity 0.179 m/s; bin means/medians [3.02, 3.19, 3.62, 4.28, 5.03, 5.46, 7.29], n [42, 86, 83, 59, 13, 15, 2]
code mapping, Mice 1 to 4 4 uA latency, binned medians 120-840 um: slope 0.00558 ms/um, p = 0.00134, r2 = 0.89, velocity 0.179 m/s; bin means/medians [2.76, 3.07, 3.67, 4.28, 5.02, 4.97, 7.29], n [42, 86, 83, 59, 13, 15, 2]
jitter: binned means slope 2.24e-04 ms2/um p = 0.275; unbinned slope -6.85e-05 p = 0.412; PL responses 457, at 120-840 um 300; latency range 1.50 to 7.99 ms
code mapping, Mice 1 to 4 5 uA latency, binned means 120-840 um: slope 0.00458 ms/um, p = 0.0037, r2 = 0.84, velocity 0.218 m/s; bin means/medians [3.3, 3.41, 3.66, 4.18, 4.97, 4.9, 7.0], n [31, 56, 62, 45, 8, 10, 6]
code mapping, Mice 1 to 4 5 uA latency, binned medians 120-840 um: slope 0.00525 ms/um, p = 0.00265, r2 = 0.86, velocity 0.190 m/s; bin means/medians [2.78, 3.18, 3.42, 3.95, 4.87, 4.77, 7.12], n [31, 56, 62, 45, 8, 10, 6]
jitter: binned means slope 1.01e-04 ms2/um p = 0.597; unbinned slope -6.02e-05 p = 0.536; PL responses 321, at 120-840 um 218; latency range 1.50 to 8.52 ms
code mapping, Mice 1 to 4 6 uA latency, binned means 120-840 um: slope 0.00437 ms/um, p = 0.000882, r2 = 0.91, velocity 0.229 m/s; bin means/medians [3.32, 3.45, 3.6, 4.33, 4.94, 5.17, 6.62], n [20, 38, 47, 36, 2, 9, 4]
code mapping, Mice 1 to 4 6 uA latency, binned medians 120-840 um: slope 0.00492 ms/um, p = 0.000759, r2 = 0.91, velocity 0.203 m/s; bin means/medians [2.85, 3.13, 3.39, 4.09, 4.94, 4.9, 6.65], n [20, 38, 47, 36, 2, 9, 4]
jitter: binned means slope -3.40e-05 ms2/um p = 0.882; unbinned slope -1.30e-04 p = 0.387; PL responses 203, at 120-840 um 156; latency range 1.60 to 7.25 ms
code mapping, leaving out one mouse at a time: velocity (m/s) from binned means / binned medians
without Mouse 1 (PL responses 814): 4 uA 0.189 / 0.180 (p 0.00047 / 0.00016); 5 uA 0.227 / 0.188 (p 0.00016 / 0.0003); 6 uA 0.218 / 0.192 (p 0.00087 / 0.00067)
without Mouse 2 (PL responses 1050): 4 uA 0.160 / 0.160 (p 0.00015 / 0.00015); 5 uA 0.206 / 0.186 (p 6.6e-05 / 0.00033); 6 uA 0.208 / 0.194 (p 0.00037 / 0.00035)
without Mouse 3 (PL responses 941): 4 uA 0.152 / 0.146 (p 0.00086 / 0.00042); 5 uA 0.394 / 0.334 (p 0.12 / 0.037); 6 uA 0.255 / 0.231 (p 0.0021 / 0.00082)
without Mouse 4 (PL responses 1110): 4 uA 0.169 / 0.161 (p 0.00043 / 0.00015); 5 uA 0.231 / 0.191 (p 0.00052 / 0.0004); 6 uA 0.248 / 0.216 (p 0.0022 / 0.0014)
without Mouse 5 (PL responses 981): 4 uA 0.179 / 0.179 (p 0.00068 / 0.0013); 5 uA 0.218 / 0.190 (p 0.0037 / 0.0027); 6 uA 0.229 / 0.203 (p 0.00088 / 0.00076)
PL responses per mouse at 120-840 um: {'Mouse 1': 246, 'Mouse 2': 107, 'Mouse 3': 221, 'Mouse 4': 100, 'Mouse 5': 242}
Result: partially reproduced. With the stimulation contact placed as the paper's code places it, PL latency rises linearly with distance at every current (binned means over 120 to 840 um: slopes 0.00591, 0.00445, 0.00442 ms/um, p = 0.000384, 0.000314, 0.00079, r2 0.91 to 0.94; binned medians: p = 0.000147, 0.000356, 0.000486), and jitter shows no slope (binned p = 0.148, 0.33, 0.377; unbinned 0.79, 0.819, 0.143). The velocities are 0.169, 0.225 and 0.226 m/s from binned means and 0.161, 0.188 and 0.199 m/s from binned medians: at 4 uA both are below the 0.18 to 0.297 m/s the criterion allows, so no reading has all three inside. Pooled over mice, the outputs favour the code's mapping: the contacts lie at 0 to 520 um depth, 12 of 17 between 150 and 420 um (Methods: 200 to 400 um), against -156 to 1559 um with the electrode table's rel_y, and the fraction of modulated responses falls from 0.711 at 60 to 120 um to 0.390 beyond 960 um, where under the electrode-table mapping it does not fall (0.299 to 0.513). The check is not clean: the 0 to 60 um bin is the lowest (0.315), and per mouse the fall reverses in Mouse 5 (0.31 within 240 um, 0.48 beyond 960 um) and Mouse 6 (0.08 and 0.88; Mouse 6 is not in the fit), while Mouse 5 looks local under the electrode-table mapping (0.68 and 0.56). Without Mouse 5 the 4 uA velocity is 0.179 m/s from both binned means and medians (5 and 6 uA 0.190 to 0.229), still below 0.18; leaving out one mouse at a time, the 4 uA velocity is 0.146 to 0.189 m/s, and reaches 0.18 under both aggregations only without Mouse 1 (0.189 and 0.180). Under the electrode-table mapping no latency slope is significant (p 0.294 to 0.962).
fig5fg-pl-modulation-increase¶
PL modulation over the 700 ms train increases from early to late at 4, 5 and 6 uA, and PL time to maximum firing rate decreases at 5 and 6 uA.
In [24]:
PLm = M50[M50.pl.astype(bool)]
NPLm = M50[~M50.pl.astype(bool)]
res_pl_z = phase_test(PLm, "z", "PL z")
res_pl_t = phase_test(PLm, "tmax_ms", "PL time to max (ms)")
print("variant, PL calls peaking in the last two unblanked bins counted as NPL:")
PLs = M50[M50.pl_strict]
_ = phase_test(PLs, "z", "PL (strict) z")
_ = phase_test(PLs, "tmax_ms", "PL (strict) time to max (ms)")
PL z 4 uA: early median 2.644 (n=137), late median 4.129 (n=468), p = 2.64e-05, r = 0.24 PL z 5 uA: early median 3.566 (n=187), late median 4.572 (n=300), p = 0.000333, r = 0.19 PL z 6 uA: early median 3.412 (n=180), late median 5.037 (n=141), p = 0.000345, r = 0.23 PL time to max (ms) 4 uA: early median 75.000 (n=137), late median 45.000 (n=468), p = 0.0858, r = -0.09 PL time to max (ms) 5 uA: early median 75.000 (n=187), late median 30.000 (n=300), p = 0.0031, r = -0.16 PL time to max (ms) 6 uA: early median 75.000 (n=180), late median 35.000 (n=141), p = 0.041, r = -0.13 variant, PL calls peaking in the last two unblanked bins counted as NPL: PL (strict) z 4 uA: early median 2.644 (n=137), late median 4.129 (n=468), p = 2.64e-05, r = 0.24 PL (strict) z 5 uA: early median 3.566 (n=187), late median 4.572 (n=296), p = 0.000335, r = 0.19 PL (strict) z 6 uA: early median 3.412 (n=180), late median 5.075 (n=138), p = 0.000306, r = 0.24 PL (strict) time to max (ms) 4 uA: early median 75.000 (n=137), late median 45.000 (n=468), p = 0.0858, r = -0.09 PL (strict) time to max (ms) 5 uA: early median 75.000 (n=187), late median 35.000 (n=296), p = 0.00415, r = -0.15 PL (strict) time to max (ms) 6 uA: early median 75.000 (n=180), late median 35.000 (n=138), p = 0.0458, r = -0.13
Result: reproduced. PL modulation rises at 4, 5 and 6 uA (p = 2.64e-05, 0.000333, 0.000345; r = 0.24, 0.19, 0.23) and PL time to maximum falls at 5 uA (75 to 30 ms, p = 0.0031) and 6 uA (75 to 35 ms, p = 0.041), not significantly at 4 uA (p = 0.0858), as the paper reports. The PL n (137/468, 187/300, 180/141) are close to the paper's (136/469, 188/285, 179/135). Counting the PL calls that peak in the last two unblanked bins as NPL gives the same six results.
figs4de-npl-no-change¶
NPL responses show no significant change in modulation or time to maximum firing rate over training at any current.
In [25]:
res_npl_z = phase_test(NPLm, "z", "NPL z")
res_npl_t = phase_test(NPLm, "tmax_ms", "NPL time to max (ms)")
print("variant, PL calls peaking in the last two unblanked bins counted as NPL:")
NPLs = M50[~M50.pl_strict]
_ = phase_test(NPLs, "z", "NPL (strict) z")
_ = phase_test(NPLs, "tmax_ms", "NPL (strict) time to max (ms)")
NPL z 4 uA: early median 2.998 (n=54), late median 3.391 (n=270), p = 0.185, r = 0.11 NPL z 5 uA: early median 2.894 (n=97), late median 2.943 (n=198), p = 0.756, r = 0.02 NPL z 6 uA: early median 3.642 (n=120), late median 4.149 (n=115), p = 0.174, r = 0.10 NPL time to max (ms) 4 uA: early median 90.000 (n=54), late median 60.000 (n=270), p = 0.0824, r = -0.15 NPL time to max (ms) 5 uA: early median 105.000 (n=97), late median 85.000 (n=198), p = 0.386, r = -0.06 NPL time to max (ms) 6 uA: early median 125.000 (n=120), late median 125.000 (n=115), p = 0.942, r = 0.01 variant, PL calls peaking in the last two unblanked bins counted as NPL: NPL (strict) z 4 uA: early median 2.998 (n=54), late median 3.391 (n=270), p = 0.185, r = 0.11 NPL (strict) z 5 uA: early median 2.894 (n=97), late median 3.029 (n=202), p = 0.712, r = 0.03 NPL (strict) z 6 uA: early median 3.642 (n=120), late median 4.193 (n=118), p = 0.172, r = 0.10 NPL (strict) time to max (ms) 4 uA: early median 90.000 (n=54), late median 60.000 (n=270), p = 0.0824, r = -0.15 NPL (strict) time to max (ms) 5 uA: early median 105.000 (n=97), late median 80.000 (n=202), p = 0.3, r = -0.07 NPL (strict) time to max (ms) 6 uA: early median 125.000 (n=120), late median 120.000 (n=118), p = 0.933, r = -0.01
Result: reproduced. NPL modulation does not change at 4, 5 or 6 uA (p = 0.185, 0.756, 0.174; r = 0.11, 0.02, 0.10), nor does NPL time to maximum (p = 0.0824, 0.386, 0.942). The NPL n (54/270, 97/198, 120/115) are close to the paper's (53/260, 95/204, 124/115). The variant that counts late-peaking PL calls as NPL gives the same.
fig5h-pl-spike-probability-increase¶
PL responses show increased peak spike probability within the 10 ms interpulse interval over training at 5 and 6 uA, not at 4 uA.
In [26]:
res_pl_pk = phase_test(PLm, "peak_prob", "PL peak spike probability")
print("peak probability range:", PLm.peak_prob.min().round(4), "to", PLm.peak_prob.max().round(4))
print("variant, PL calls peaking in the last two unblanked bins counted as NPL:")
_ = phase_test(PLs, "peak_prob", "PL (strict) peak spike probability")
print("the PL calls under three other seeds (see fig5d):")
for seed, labs in pl_reseeded.items():
_ = phase_test(M50[labs], "peak_prob", f"seed {seed}, PL peak spike probability")
PL peak spike probability 4 uA: early median 0.040 (n=137), late median 0.042 (n=468), p = 0.105, r = 0.09 PL peak spike probability 5 uA: early median 0.045 (n=187), late median 0.051 (n=300), p = 0.0613, r = 0.10 PL peak spike probability 6 uA: early median 0.047 (n=180), late median 0.061 (n=141), p = 0.0113, r = 0.16 peak probability range: 0.0067 to 0.304 variant, PL calls peaking in the last two unblanked bins counted as NPL: PL (strict) peak spike probability 4 uA: early median 0.040 (n=137), late median 0.042 (n=468), p = 0.105, r = 0.09 PL (strict) peak spike probability 5 uA: early median 0.045 (n=187), late median 0.052 (n=296), p = 0.0499, r = 0.11 PL (strict) peak spike probability 6 uA: early median 0.047 (n=180), late median 0.061 (n=138), p = 0.00695, r = 0.18 the PL calls under three other seeds (see fig5d): seed 11, PL peak spike probability 4 uA: early median 0.040 (n=140), late median 0.042 (n=462), p = 0.0515, r = 0.11 seed 11, PL peak spike probability 5 uA: early median 0.044 (n=190), late median 0.052 (n=294), p = 0.0197, r = 0.13 seed 11, PL peak spike probability 6 uA: early median 0.047 (n=179), late median 0.061 (n=139), p = 0.0086, r = 0.17 seed 12, PL peak spike probability 4 uA: early median 0.040 (n=137), late median 0.042 (n=464), p = 0.0938, r = 0.09 seed 12, PL peak spike probability 5 uA: early median 0.044 (n=186), late median 0.053 (n=296), p = 0.0275, r = 0.12 seed 12, PL peak spike probability 6 uA: early median 0.047 (n=180), late median 0.061 (n=138), p = 0.0102, r = 0.17 seed 13, PL peak spike probability 4 uA: early median 0.040 (n=138), late median 0.041 (n=471), p = 0.137, r = 0.08 seed 13, PL peak spike probability 5 uA: early median 0.044 (n=191), late median 0.051 (n=296), p = 0.0485, r = 0.11 seed 13, PL peak spike probability 6 uA: early median 0.047 (n=179), late median 0.061 (n=140), p = 0.0108, r = 0.17
Result: partially reproduced. The PL peak spike probability rises at 6 uA (0.047 to 0.061, p = 0.0113, r = 0.16) but not significantly at 5 uA (0.045 to 0.051, p = 0.0613, r = 0.10); at 4 uA p = 0.105, as the paper reports no change there. The 5 uA test depends on the random PL call: with three other seeds it gives p = 0.0197, 0.0275 and 0.0485, and the variant that counts late-peaking PL calls as NPL gives p = 0.0499.
figs5-restricted-window-consistent¶
In a 120 ms window after stimulation onset, before movement, PL modulation still increases at all currents while NPL modulation does not change, and PL time to maximum falls at 4 and 5 uA.
Spikes in the first 120 ms of the train against the 120 ms before it; time to maximum in 2 ms bins smoothed by a 5 ms Gaussian. Readings: (keep) modulated set and PL labels from the 700 ms window; (recompute) modulated by the 120 ms test and PL from the first 12 pulses.
In [27]:
t0 = time.time()
r120 = icms.response_table(sessions, rng_seed=3, currents=icms.CURRENTS, window=0.12, tmax_bin=0.002, tmax_sigma=0.005)
key700 = {(r["si"], r["unit"], r["channel"], r["current"]): r for r in resp_rows}
rng120 = np.random.default_rng(4)
for r in r120:
r["mod120"] = icms.modulated(r, 50, True, ch_days)
k = key700[(r["si"], r["unit"], r["channel"], r["current"])]
r["mod700"], r["pl700"] = k["mod50"], bool(k.get("pl", False))
if r["mod120"]:
s = sessions[r["si"]]
tr = icms.ephys_trains(s)
m = (tr["stim_channel"] == r["channel"]) & (tr["current_uA"] == r["current"])
Mx = icms.pulse_matrix(s["spikes"][r["unit"]], icms.pulse_starts(np.sort(tr["start_time"][m]), 12))
r["pl120"] = icms.pulse_locking(Mx, rng120)["pl"]
R120 = pd.DataFrame(r120)
R120 = R120[R120.phase.notna()]
print(f"120 ms responses: {len(R120)} in {time.time() - t0:.0f} s; modulated in 120 ms: {int(R120.mod120.sum())}")
res_s5 = {}
for reading, sel, lab in [("keep", R120.mod700, "pl700"), ("recompute", R120.mod120, "pl120")]:
d = R120[sel]
pl = d[d[lab].astype(bool)]
npl = d[~d[lab].astype(bool)]
print(f"\nreading {reading}:")
res_s5[reading] = {"pl_z": phase_test(pl, "z", "PL z"), "npl_z": phase_test(npl, "z", "NPL z"),
"pl_t": phase_test(pl, "tmax_ms", "PL time to max"),
"npl_t": phase_test(npl, "tmax_ms", "NPL time to max")}
120 ms responses: 4862 in 30 s; modulated in 120 ms: 1528 reading keep: PL z 4 uA: early median 1.860 (n=136), late median 2.652 (n=466), p = 0.000508, r = 0.20 PL z 5 uA: early median 2.401 (n=187), late median 3.028 (n=299), p = 0.000218, r = 0.20 PL z 6 uA: early median 2.235 (n=178), late median 3.587 (n=141), p = 0.000116, r = 0.25 NPL z 4 uA: early median 2.065 (n=53), late median 2.530 (n=269), p = 0.254, r = 0.10 NPL z 5 uA: early median 1.950 (n=94), late median 2.525 (n=194), p = 0.633, r = 0.03 NPL z 6 uA: early median 2.501 (n=118), late median 2.520 (n=113), p = 0.776, r = 0.02 PL time to max 4 uA: early median 35.000 (n=137), late median 23.000 (n=468), p = 0.0107, r = -0.14 PL time to max 5 uA: early median 29.000 (n=187), late median 17.000 (n=300), p = 0.0119, r = -0.14 PL time to max 6 uA: early median 25.000 (n=180), late median 27.000 (n=141), p = 0.994, r = -0.00 NPL time to max 4 uA: early median 43.000 (n=54), late median 43.000 (n=270), p = 0.992, r = -0.00 NPL time to max 5 uA: early median 33.000 (n=97), late median 39.000 (n=197), p = 0.964, r = 0.00 NPL time to max 6 uA: early median 43.000 (n=120), late median 55.000 (n=114), p = 0.607, r = -0.04 reading recompute: PL z 4 uA: early median 2.062 (n=70), late median 2.955 (n=219), p = 0.000476, r = 0.28 PL z 5 uA: early median 2.798 (n=97), late median 3.779 (n=142), p = 0.00113, r = 0.25 PL z 6 uA: early median 3.595 (n=79), late median 4.221 (n=59), p = 0.14, r = 0.15 NPL z 4 uA: early median 3.094 (n=63), late median 3.786 (n=318), p = 0.0368, r = 0.17 NPL z 5 uA: early median 2.677 (n=92), late median 3.675 (n=197), p = 0.000146, r = 0.28 NPL z 6 uA: early median 3.000 (n=94), late median 4.320 (n=98), p = 5.56e-05, r = 0.34 PL time to max 4 uA: early median 33.000 (n=70), late median 21.000 (n=219), p = 0.016, r = -0.19 PL time to max 5 uA: early median 31.000 (n=97), late median 17.000 (n=142), p = 0.0157, r = -0.18 PL time to max 6 uA: early median 19.000 (n=79), late median 19.000 (n=59), p = 0.612, r = 0.05 NPL time to max 4 uA: early median 27.000 (n=63), late median 29.000 (n=318), p = 0.931, r = 0.01 NPL time to max 5 uA: early median 26.000 (n=92), late median 23.000 (n=197), p = 0.19, r = -0.10 NPL time to max 6 uA: early median 31.000 (n=94), late median 36.000 (n=98), p = 0.851, r = 0.02
Result: reproduced, under the reading that keeps the 700 ms labels. In the first 120 ms, PL modulation rises at 4, 5 and 6 uA (p = 0.000508, 0.000218, 0.000116), NPL modulation does not change (p = 0.254, 0.633, 0.776), and PL time to maximum falls at 4 and 5 uA (p = 0.0107, 0.0119) and not at 6 uA (p = 0.994), as the paper reports; NPL time to maximum does not change (p = 0.992, 0.964, 0.607). With the modulated set and labels recomputed in the 120 ms window, NPL modulation rises at all three currents (p = 0.0368, 0.000146, 5.56e-05) and PL modulation does not rise at 6 uA (p = 0.14).
fig6c-kl-divergence¶
Hit-miss separation of delta-spikes distributions (KL divergence) increases over time for both PL and NPL responses, and NPL responses separate hit from miss more than PL responses in both early and late phases.
Trials of the modulated responses (mod50) at 4 to 6 uA, pooled per week and label. Readings (k1) winsorizing 1%, fraction 0.4 and (k2) 2.5%, fraction 0.8; 1000 draws.
In [28]:
lab_of = {(r.si, r.unit, r.channel, r.current): bool(r.pl) for r in M50.itertuples()}
trial_pool = defaultdict(lambda: {"hit": [], "miss": []})
for r in resp_rows:
k = (r["si"], r["unit"], r["channel"], r["current"])
if k not in lab_of or r["phase"] is None:
continue
key = (r["week"], "PL" if lab_of[k] else "NPL")
trial_pool[key]["hit"].extend(r["trial_delta"][r["trial_hit"]].tolist())
trial_pool[key]["miss"].extend(r["trial_delta"][~r["trial_hit"]].tolist())
print("trials per week and label (hit, miss):", {k: (len(v["hit"]), len(v["miss"])) for k, v in sorted(trial_pool.items())})
res_kl = {}
for reading, wins, frac in [("k1", 0.01, 0.4), ("k2", 0.025, 0.8)]:
rngk = np.random.default_rng(10)
kl = {k: icms.kl_bootstrap(v["hit"], v["miss"], rngk, wins, frac) for k, v in sorted(trial_pool.items())}
def pool(phase, lab):
return np.concatenate([kl[(w, lab)] for w in (icms.EARLY_WEEKS if phase == "early" else icms.LATE_WEEKS)
if (w, lab) in kl])
tests = {"PL early vs late": icms.mwu(pool("early", "PL"), pool("late", "PL")),
"NPL early vs late": icms.mwu(pool("early", "NPL"), pool("late", "NPL")),
"NPL vs PL early": icms.mwu(pool("early", "PL"), pool("early", "NPL")),
"NPL vs PL late": icms.mwu(pool("late", "PL"), pool("late", "NPL"))}
q = icms.fdr_bh([t["p"] for t in tests.values()])
print(f"\nreading {reading}: median KL per week and label:",
{f"w{w} {l}": round(float(np.median(v)), 4) for (w, l), v in kl.items()})
for (name, t), qq in zip(tests.items(), q):
print(f" {name}: medians {t['med1']:.4f} vs {t['med2']:.4f} (n={t['n1']}, {t['n2']}), FDR p = {qq:.3g}, r = {t['r']:.2f}")
res_kl[reading] = (tests, q)
trials per week and label (hit, miss): {(0, 'NPL'): (1057, 1293), (0, 'PL'): (2655, 2965), (1, 'NPL'): (3894, 1766), (1, 'PL'): (7473, 3660), (2, 'NPL'): (3082, 982), (2, 'PL'): (6449, 2196), (3, 'NPL'): (3459, 1749), (3, 'PL'): (6834, 4597), (4, 'NPL'): (6484, 3246), (4, 'PL'): (7334, 4536)}
reading k1: median KL per week and label: {'w0 NPL': 0.1504, 'w0 PL': 0.0433, 'w1 NPL': 0.2526, 'w1 PL': 0.0383, 'w2 NPL': 0.2268, 'w2 PL': 0.0445, 'w3 NPL': 0.2277, 'w3 PL': 0.0136, 'w4 NPL': 0.1999, 'w4 PL': 0.0289}
PL early vs late: medians 0.0403 vs 0.0279 (n=2000, 3000), FDR p = 3e-192, r = -0.49
NPL early vs late: medians 0.2089 vs 0.2204 (n=2000, 3000), FDR p = 5.26e-08, r = 0.09
NPL vs PL early: medians 0.0403 vs 0.2089 (n=2000, 2000), FDR p = 0, r = 1.00
NPL vs PL late: medians 0.0279 vs 0.2204 (n=3000, 3000), FDR p = 0, r = 0.99
reading k2: median KL per week and label: {'w0 NPL': 0.9182, 'w0 PL': 0.0475, 'w1 NPL': 0.2498, 'w1 PL': 0.0659, 'w2 NPL': 0.2407, 'w2 PL': 0.0348, 'w3 NPL': 0.1969, 'w3 PL': 0.0082, 'w4 NPL': 0.2224, 'w4 PL': 0.859}
PL early vs late: medians 0.0586 vs 0.0348 (n=2000, 3000), FDR p = 1.48e-59, r = -0.27
NPL early vs late: medians 0.4185 vs 0.2192 (n=2000, 3000), FDR p = 0, r = -0.75
NPL vs PL early: medians 0.0586 vs 0.4185 (n=2000, 2000), FDR p = 0, r = 1.00
NPL vs PL late: medians 0.0348 vs 0.2192 (n=3000, 3000), FDR p = 1.81e-115, r = 0.34
Result: not reproduced. NPL responses separate hit from miss more than PL in both phases (reading k1, median KL 0.2089 against 0.0403 early and 0.2204 against 0.0279 late, FDR p = 0), and NPL separation rises over training (0.2089 to 0.2204, FDR p = 5.26e-08, r = 0.09), but PL separation falls (0.0403 to 0.0279, FDR p = 3e-192, r = -0.49). Under reading k2 both fall (PL 0.0586 to 0.0348, NPL 0.4185 to 0.2192), and NPL exceeds PL in both phases.
fig6d-npl-hit-miss-divergence¶
PL responses show no significant hit-miss divergence in firing rate at any time, while NPL responses diverge from 65 ms after ICMS onset early in training and from 35 ms late in training, both before movement onset.
Each row is one modulated response's trial-averaged rate for hit and for miss trials (responses with at least one of each), from train onset to 700 ms. Labels are permuted by swapping a response's hit and miss rows, which permutes the labels over rows while keeping each response's pair (the code's permutation of hit and miss labels over rows; the exact scheme is not in the paper). 1000 permutations. Readings (b1) 10 ms bins and (b2) 5 ms bins, onset at the bin centre; (o1) first and (o2) last significant cluster.
In [29]:
t0 = time.time()
res_6d = {}
for bin_s, bl in [(0.010, "b1"), (0.005, "b2")]:
rows_h = defaultdict(list)
rows_m = defaultdict(list)
for r in resp_rows:
k = (r["si"], r["unit"], r["channel"], r["current"])
if k not in lab_of or r["phase"] is None or r["trial_hit"].sum() == 0 or (~r["trial_hit"]).sum() == 0:
continue
spk = sessions[r["si"]]["spikes"][r["unit"]]
key = (r["phase"], "PL" if lab_of[k] else "NPL")
rows_h[key].append(icms.psth(spk, r["starts"][r["trial_hit"]], bin_s))
rows_m[key].append(icms.psth(spk, r["starts"][~r["trial_hit"]], bin_s))
for key in sorted(rows_h):
rngc = np.random.default_rng(20)
H, Mm = np.array(rows_h[key]), np.array(rows_m[key])
cl, tstat = icms.cluster_test(H, Mm, rngc, n_perm=1000)
sig = [c for c in cl if c["p"] < 0.05]
on1 = (sig[0]["start_bin"] + 0.5) * bin_s * 1000 if sig else np.nan
on2 = (sig[-1]["start_bin"] + 0.5) * bin_s * 1000 if sig else np.nan
res_6d[(bl, key)] = (len(H), sig, on1, on2)
print(f"{bl} {key[0]} {key[1]}: rows {len(H)}, clusters {len(cl)}, significant "
f"{[(round((c['start_bin'] + 0.5) * bin_s * 1000, 1), round((c['end_bin'] + 0.5) * bin_s * 1000, 1), round(c['mass'], 1), c['p']) for c in sig]}; "
f"onset first {on1} ms, last {on2} ms")
print(f"cluster tests in {time.time() - t0:.0f} s")
b1 early NPL: rows 223, clusters 7, significant [(75.0, 125.0, 18.3, 0.0), (145.0, 155.0, 4.3, 0.006), (175.0, 275.0, 29.2, 0.0), (295.0, 315.0, 7.3, 0.0), (335.0, 365.0, 10.1, 0.0), (385.0, 695.0, 122.4, 0.0)]; onset first 75.0 ms, last 385.0 ms
b1 early PL: rows 436, clusters 9, significant [(365.0, 375.0, 3.6, 0.002), (605.0, 615.0, 2.9, 0.008), (645.0, 665.0, 5.3, 0.0), (685.0, 685.0, 2.0, 0.037)]; onset first 365.0 ms, last 685.0 ms
b1 late NPL: rows 506, clusters 4, significant [(5.0, 5.0, 3.8, 0.018), (45.0, 695.0, 727.0, 0.0)]; onset first 5.0 ms, last 45.0 ms
b1 late PL: rows 826, clusters 24, significant [(5.0, 5.0, 1.9, 0.013), (25.0, 45.0, 5.6, 0.0), (65.0, 115.0, 12.0, 0.0), (135.0, 145.0, -3.0, 0.002), (165.0, 165.0, 1.9, 0.008), (175.0, 175.0, -2.2, 0.004), (195.0, 195.0, 2.4, 0.003), (205.0, 205.0, -2.1, 0.004), (255.0, 255.0, 1.9, 0.012), (285.0, 285.0, -1.6, 0.043), (305.0, 315.0, -3.7, 0.0), (335.0, 335.0, 2.0, 0.004), (345.0, 385.0, -10.0, 0.0), (395.0, 395.0, 3.1, 0.001), (405.0, 505.0, -28.3, 0.0), (515.0, 515.0, 3.1, 0.001), (525.0, 545.0, -7.8, 0.0), (555.0, 555.0, 2.9, 0.002), (565.0, 585.0, -10.1, 0.0), (595.0, 595.0, 4.4, 0.0), (605.0, 635.0, -15.1, 0.0), (645.0, 685.0, 29.4, 0.0), (695.0, 695.0, -3.3, 0.001)]; onset first 5.0 ms, last 695.0000000000001 ms
b2 early NPL: rows 223, clusters 19, significant [(87.5, 117.5, 23.2, 0.0), (127.5, 132.5, 5.2, 0.009), (142.5, 147.5, 4.3, 0.02), (172.5, 187.5, 8.8, 0.001), (197.5, 222.5, 15.3, 0.0), (232.5, 267.5, 19.8, 0.0), (277.5, 277.5, 3.8, 0.039), (287.5, 302.5, 9.3, 0.001), (312.5, 322.5, 5.2, 0.01), (337.5, 362.5, 17.9, 0.0), (372.5, 447.5, 40.7, 0.0), (457.5, 457.5, 3.7, 0.042), (467.5, 567.5, 74.9, 0.0), (577.5, 627.5, 40.3, 0.0), (637.5, 697.5, 41.3, 0.0)]; onset first 87.50000000000001 ms, last 637.5000000000001 ms
b2 early PL: rows 436, clusters 27, significant [(107.5, 112.5, 3.7, 0.002), (347.5, 347.5, 3.2, 0.013), (357.5, 372.5, 7.3, 0.0), (527.5, 532.5, 3.9, 0.002), (587.5, 587.5, 2.9, 0.015), (602.5, 607.5, 3.8, 0.002), (647.5, 647.5, 3.3, 0.013), (657.5, 667.5, 6.0, 0.0)]; onset first 107.5 ms, last 657.5 ms
b2 late NPL: rows 506, clusters 20, significant [(7.5, 7.5, 4.1, 0.031), (12.5, 17.5, -4.6, 0.021), (22.5, 32.5, 9.7, 0.0), (37.5, 42.5, -4.6, 0.021), (47.5, 57.5, 12.8, 0.0), (62.5, 62.5, -4.2, 0.031), (67.5, 167.5, 181.5, 0.0), (172.5, 172.5, -7.4, 0.003), (177.5, 252.5, 161.7, 0.0), (257.5, 257.5, -4.0, 0.033), (262.5, 277.5, 41.0, 0.0), (282.5, 282.5, -5.2, 0.014), (287.5, 307.5, 48.2, 0.0), (312.5, 312.5, -6.3, 0.007), (317.5, 347.5, 64.7, 0.0), (352.5, 352.5, -8.2, 0.001), (357.5, 427.5, 143.4, 0.0), (432.5, 432.5, -7.0, 0.005), (437.5, 697.5, 664.9, 0.0)]; onset first 7.5 ms, last 437.5 ms
b2 late PL: rows 826, clusters 71, significant [(2.5, 7.5, 8.4, 0.0), (12.5, 12.5, -2.3, 0.014), (17.5, 32.5, 13.1, 0.0), (37.5, 42.5, -3.5, 0.002), (47.5, 52.5, 9.1, 0.0), (57.5, 57.5, -2.0, 0.041), (62.5, 77.5, 17.6, 0.0), (82.5, 82.5, -3.1, 0.003), (87.5, 87.5, 3.5, 0.002), (92.5, 92.5, -2.1, 0.025), (97.5, 97.5, 7.9, 0.0), (102.5, 102.5, -2.6, 0.009), (107.5, 117.5, 13.0, 0.0), (127.5, 127.5, 4.4, 0.0), (137.5, 137.5, 4.9, 0.0), (147.5, 147.5, 5.0, 0.0), (157.5, 157.5, -5.3, 0.0), (162.5, 172.5, 10.3, 0.0), (177.5, 187.5, -11.2, 0.0), (192.5, 192.5, 4.3, 0.0), (197.5, 202.5, -5.0, 0.0), (207.5, 207.5, 6.2, 0.0), (217.5, 217.5, 4.7, 0.0), (222.5, 247.5, -19.6, 0.0), (252.5, 252.5, 3.6, 0.001), (257.5, 287.5, -20.0, 0.0), (297.5, 297.5, 5.2, 0.0), (302.5, 302.5, -2.6, 0.009), (307.5, 307.5, 5.2, 0.0), (312.5, 332.5, -13.8, 0.0), (337.5, 337.5, 4.6, 0.0), (342.5, 372.5, -23.3, 0.0), (377.5, 382.5, 8.0, 0.0), (387.5, 387.5, -4.8, 0.0), (392.5, 392.5, 4.1, 0.0), (397.5, 412.5, -12.6, 0.0), (417.5, 417.5, 3.5, 0.002), (422.5, 422.5, -2.6, 0.01), (427.5, 427.5, 4.7, 0.0), (432.5, 437.5, -5.8, 0.0), (442.5, 442.5, 4.3, 0.0), (447.5, 452.5, -6.3, 0.0), (457.5, 457.5, 5.5, 0.0), (462.5, 512.5, -47.6, 0.0), (517.5, 517.5, 5.2, 0.0), (522.5, 527.5, -6.6, 0.0), (532.5, 532.5, 6.2, 0.0), (537.5, 552.5, -13.8, 0.0), (557.5, 557.5, 7.3, 0.0), (562.5, 562.5, -3.7, 0.001), (567.5, 567.5, 6.2, 0.0), (572.5, 572.5, -3.9, 0.0), (577.5, 577.5, 8.7, 0.0), (582.5, 592.5, -14.1, 0.0), (597.5, 597.5, 6.2, 0.0), (602.5, 612.5, -14.0, 0.0), (617.5, 617.5, 8.1, 0.0), (622.5, 622.5, -4.1, 0.0), (627.5, 632.5, 10.0, 0.0), (637.5, 637.5, -4.2, 0.0), (642.5, 647.5, 12.3, 0.0), (652.5, 652.5, -5.4, 0.0), (657.5, 677.5, 38.0, 0.0), (682.5, 682.5, -5.1, 0.0), (687.5, 687.5, 8.6, 0.0), (692.5, 692.5, -3.1, 0.003), (697.5, 697.5, 7.3, 0.0)]; onset first 2.5 ms, last 697.5 ms cluster tests in 79 s
Result: not reproduced. PL responses diverge between hit and miss in both phases (with 10 ms bins, significant clusters from 365 ms early and from 5 ms late; with 5 ms bins from 107.5 and 2.5 ms). With 10 ms bins the first NPL cluster starts at 75 ms early, at the edge of 65 +- 10 ms, and at 5 ms late (a single bin), outside 35 +- 10 ms; the last NPL cluster starts at 385 ms early and 45 ms late. With 5 ms bins, 87.5 and 7.5 ms (first) and 637.5 and 437.5 ms (last). No reading has both NPL onsets within tolerance, though late is earlier than early under the first-cluster readings.
fig6f-population-coupling¶
PL population coupling increases from early to late while NPL coupling stays stable, and NPL units have higher population coupling than PL units in both phases.
In [30]:
def coupling_tests(rows_c, label):
C = pd.DataFrame(rows_c)
C = C[C.phase.notna() & C.coupling.notna()]
C["lab"] = np.where(C.pl.astype(bool), "PL", "NPL")
tests = {"PL early vs late": icms.mwu(C[(C.lab == "PL") & (C.phase == "early")].coupling,
C[(C.lab == "PL") & (C.phase == "late")].coupling),
"NPL early vs late": icms.mwu(C[(C.lab == "NPL") & (C.phase == "early")].coupling,
C[(C.lab == "NPL") & (C.phase == "late")].coupling),
"NPL vs PL early": icms.mwu(C[(C.lab == "PL") & (C.phase == "early")].coupling,
C[(C.lab == "NPL") & (C.phase == "early")].coupling),
"NPL vs PL late": icms.mwu(C[(C.lab == "PL") & (C.phase == "late")].coupling,
C[(C.lab == "NPL") & (C.phase == "late")].coupling)}
q = icms.fdr_bh([t["p"] for t in tests.values()])
print(f"{label}: values {len(C)}, range {C.coupling.min():.2f} to {C.coupling.max():.2f}")
for (name, t), qq in zip(tests.items(), q):
print(f" {name}: medians {t['med1']:.3f} vs {t['med2']:.3f} (n={t['n1']}, {t['n2']}), FDR p = {qq:.3g}, r = {t['r']:.2f}")
print(" median by mouse and label:", C.groupby([C.subject.map(icms.MOUSE), "lab"]).coupling.median().round(2).to_dict())
return tests, q
for r in resp_rows:
r["mod"] = r["mod50"] # population and units: the mod50 responses
t0 = time.time()
pc_main = icms.population_coupling(sessions, resp_rows, rng_seed=30)
res_6f = coupling_tests(pc_main, "700 ms, all six mice")
print(f"in {time.time() - t0:.0f} s")
700 ms, all six mice: values 2260, range -12.03 to 9.63
PL early vs late: medians -0.134 vs 0.481 (n=503, 906), FDR p = 0.000671, r = 0.11
NPL early vs late: medians 1.593 vs 1.468 (n=270, 581), FDR p = 0.303, r = -0.04
NPL vs PL early: medians -0.134 vs 1.593 (n=503, 270), FDR p = 6.08e-22, r = 0.42
NPL vs PL late: medians 0.481 vs 1.468 (n=906, 581), FDR p = 3.83e-22, r = 0.30
median by mouse and label: {('Mouse 1', 'NPL'): 2.69, ('Mouse 1', 'PL'): 0.56, ('Mouse 2', 'NPL'): 0.23, ('Mouse 2', 'PL'): -0.9, ('Mouse 3', 'NPL'): 2.03, ('Mouse 3', 'PL'): 1.05, ('Mouse 4', 'NPL'): 1.45, ('Mouse 4', 'PL'): 0.7, ('Mouse 5', 'NPL'): 1.0, ('Mouse 5', 'PL'): -0.39, ('Mouse 6', 'NPL'): 2.58, ('Mouse 6', 'PL'): 0.4}
in 7 s
Result: reproduced. PL coupling rises from early to late (median -0.134 to 0.481, FDR p = 0.000671, r = 0.11), NPL coupling does not change (1.593 to 1.468, FDR p = 0.303), and NPL exceeds PL early (FDR p = 6.08e-22, r = 0.42) and late (FDR p = 3.83e-22, r = 0.30). The n (PL 503 and 906, NPL 270 and 581) are close to the paper's (494 and 899, 263 and 553). NPL is above PL in every mouse.
figs9a-coupling-700ms¶
In the 700 ms window of fig. S9A, PL population coupling increases from early to late while NPL coupling does not change, and PL coupling is lower than NPL coupling in both phases.
Readings (a) as fig6f without Mouse 6, (b) the population built from the other modulated PL units only.
In [31]:
pc_a = [r for r in pc_main if r["subject"] != "ICMS83"]
res_s9a_a = coupling_tests(pc_a, "(a) without Mouse 6")
pc_b = icms.population_coupling(sessions, resp_rows, rng_seed=31, pop_label=True)
res_s9a_b = coupling_tests(pc_b, "(b) population of PL units only")
(a) without Mouse 6: values 2062, range -10.82 to 9.63
PL early vs late: medians -0.206 vs 0.513 (n=453, 769), FDR p = 7.76e-06, r = 0.16
NPL early vs late: medians 1.582 vs 1.440 (n=264, 576), FDR p = 0.366, r = -0.04
NPL vs PL early: medians -0.206 vs 1.582 (n=453, 264), FDR p = 2.48e-22, r = 0.44
NPL vs PL late: medians 0.513 vs 1.440 (n=769, 576), FDR p = 5.6e-19, r = 0.29
median by mouse and label: {('Mouse 1', 'NPL'): 2.69, ('Mouse 1', 'PL'): 0.56, ('Mouse 2', 'NPL'): 0.23, ('Mouse 2', 'PL'): -0.9, ('Mouse 3', 'NPL'): 2.03, ('Mouse 3', 'PL'): 1.05, ('Mouse 4', 'NPL'): 1.45, ('Mouse 4', 'PL'): 0.7, ('Mouse 5', 'NPL'): 1.0, ('Mouse 5', 'PL'): -0.39}
(b) population of PL units only: values 2245, range -16.04 to 7.00
PL early vs late: medians -0.585 vs 0.014 (n=501, 903), FDR p = 4.84e-05, r = 0.13
NPL early vs late: medians 1.261 vs 1.113 (n=270, 571), FDR p = 0.24, r = -0.05
NPL vs PL early: medians -0.585 vs 1.261 (n=501, 270), FDR p = 1.49e-22, r = 0.43
NPL vs PL late: medians 0.014 vs 1.113 (n=903, 571), FDR p = 2.75e-19, r = 0.28
median by mouse and label: {('Mouse 1', 'NPL'): 2.08, ('Mouse 1', 'PL'): -0.2, ('Mouse 2', 'NPL'): -0.33, ('Mouse 2', 'PL'): -0.98, ('Mouse 3', 'NPL'): 1.8, ('Mouse 3', 'PL'): 0.57, ('Mouse 4', 'NPL'): 1.42, ('Mouse 4', 'PL'): 0.35, ('Mouse 5', 'NPL'): 0.67, ('Mouse 5', 'PL'): -0.97, ('Mouse 6', 'NPL'): 2.51, ('Mouse 6', 'PL'): 0.08}
Result: reproduced. Without Mouse 6 (reading a), PL coupling rises (-0.206 to 0.513, FDR p = 7.76e-06), NPL does not change (FDR p = 0.366), and PL is below NPL early (FDR p = 2.48e-22) and late (FDR p = 5.6e-19). The PL n (453, 769) are near the legend's 446 and 763. Reading (b) gives the same four results (FDR p = 4.84e-05, 0.24, 1.49e-22, 2.75e-19).
figs9-coupling-window-and-control (120 ms part)¶
NPL coupling exceeds PL coupling also within the 120 ms premotion window, and the difference is absent in control animals without behavioural learning.
The control part is run at the end, with the control animals. Here, spikes in the first 120 ms of each train, a 10-spike minimum, labels from the 700 ms window.
In [32]:
pc_120 = icms.population_coupling(sessions, resp_rows, rng_seed=32, window=0.12, min_spikes=10)
res_s9_120 = coupling_tests(pc_120, "120 ms window")
120 ms window: values 2186, range -7.77 to 6.21
PL early vs late: medians -0.520 vs -0.349 (n=497, 892), FDR p = 0.408, r = 0.03
NPL early vs late: medians 0.751 vs 0.420 (n=257, 540), FDR p = 0.0954, r = -0.08
NPL vs PL early: medians -0.520 vs 0.751 (n=497, 257), FDR p = 1.12e-14, r = 0.35
NPL vs PL late: medians -0.349 vs 0.420 (n=892, 540), FDR p = 1.04e-15, r = 0.26
median by mouse and label: {('Mouse 1', 'NPL'): 0.77, ('Mouse 1', 'PL'): -0.24, ('Mouse 2', 'NPL'): -0.21, ('Mouse 2', 'PL'): -1.7, ('Mouse 3', 'NPL'): 0.61, ('Mouse 3', 'PL'): -0.24, ('Mouse 4', 'NPL'): 0.3, ('Mouse 4', 'PL'): -0.13, ('Mouse 5', 'NPL'): 1.04, ('Mouse 5', 'PL'): 0.29, ('Mouse 6', 'NPL'): -0.41, ('Mouse 6', 'PL'): -0.83}
Result (120 ms part): holds. In the first 120 ms, NPL coupling exceeds PL early (0.751 against -0.520, FDR p = 1.12e-14, r = 0.35) and late (0.420 against -0.349, FDR p = 1.04e-15, r = 0.26), as the paper reports. The outcome of the claim is given with its control part, below.
Control animals¶
The control files are first described: their stimulation protocol, their imaging tables and their units. The control claims are then tested as far as the files allow, with the paper's control code (processing/control at d1922d1) for what the paper does not say.
In [33]:
import re
ctrl = [a for a in assets if a["path"].split("/")[0][4:] not in icms.BEHAVIOUR]
crow = []
stim_desc = set()
trial_desc = set()
for a in ctrl:
f = icms.open_asset(a["asset_id"])
d = {"file": a["path"].split("/")[1], "date": f["session_start_time"][()].decode()[:10],
"subject": f["general/subject/subject_id"][()].decode()}
if "intervals/electrical_stimulation" in f:
st = f["intervals/electrical_stimulation"]
stim_desc.add(str(st.attrs["description"]))
d["pulses_per_burst"] = ",".join(map(str, np.unique(st["pulse_count"][()])))
d["freq_hz"] = ",".join(map(str, np.unique(st["frequency_hz"][()])))
d["bursts_per_trial"] = int(np.max(st["burst_index_in_trial"][()]) + 1) if "burst_index_in_trial" in st else None
d["burst_ms"] = float(np.round(np.median(st["stop_time"][()] - st["start_time"][()]) * 1000, 1))
_s, _t = st["start_time"][()], icms.control_trial_ids({k: st[k][()] for k in ["start_time", "burst_index_in_trial"]})
_o = np.argsort(_s)
_s, _t = _s[_o], _t[_o]
_same = _t[1:] == _t[:-1]
d["burst_onset_gap_s"] = f"{np.diff(_s)[_same].min():.3f}-{np.diff(_s)[_same].max():.3f}"
_first = np.array([_s[_t == k].min() for k in np.unique(_t)])
d["trial_onset_gap_s"] = f"{np.diff(_first).min():.2f}-{np.diff(_first).max():.2f}" if len(_first) > 1 else ""
d["trials_rows"] = len(f["intervals/trials/id"]) if "intervals/trials" in f else None
d["stim_trials"] = len(np.unique(_t))
d["trial_index_values"] = len(np.unique(st["trial_index"][()]))
if "intervals/trials" in f:
trial_desc.add(str(f["intervals/trials/trial_index"].attrs.get("description", "")))
d["trial_pulses"] = ",".join(map(str, np.unique(f["intervals/trials/n_pulses"][()])))
d["reps"] = len(np.unique(f["intervals/trials/rep_num"][()]))
d["units"] = len(f["units/id"]) if "units" in f else 0
d["cell_types"] = ",".join(np.unique(f["units/cell_type"][()].astype(str))) if "units" in f else ""
o = "processing/ophys"
if o + "/stimulus_conditions" in f:
g = f[o + "/ImageSegmentation_ControlSurvey/PlaneSegmentation_ControlSurvey"]
act = g["activation"][()]
cur = f[o + "/stimulus_conditions/current_uA"][()]
d["survey_rois"] = act.shape[0]
d["currents_imaged"] = ",".join(str(int(c)) for c in np.unique(cur))
d["activated_at_5uA_per_channel"] = ",".join(str(int(v)) for v in act[:, cur == 5].sum(0))
ip = f["general/optophysiology/ImagingPlane_ControlSurvey"]
d["grid_spacing_given"] = "grid_spacing" in ip
_pl = re.search(r"(\d+) planes", ip["description"][()].decode())
d["planes_in_description"] = int(_pl.group(1)) if _pl else None
d["stim_contacts_at_origin"] = int(((f[o + "/stim_electrode_positions/x"][()] == 0) &
(f[o + "/stim_electrode_positions/y"][()] == 0)).sum())
f.close()
crow.append(d)
CT = pd.DataFrame(crow).sort_values(["subject", "date"])
print(CT.to_string())
print("\ndescriptions of the control electrical_stimulation tables:", sorted(stim_desc))
print("descriptions of the control trials tables' trial_index:", sorted(trial_desc),
"| trials per distinct trial_index value: trials_rows / trial_index_values above")
print("\ncontrol mice:", CT.subject.nunique(), "| control files:", len(CT), "| with survey imaging:",
int(CT.survey_rois.notna().sum()), "| with units:", int((CT.get("units", pd.Series(dtype=float)).fillna(0) > 0).sum()))
file date subject survey_rois currents_imaged activated_at_5uA_per_channel grid_spacing_given planes_in_description stim_contacts_at_origin pulses_per_burst freq_hz bursts_per_trial burst_ms burst_onset_gap_s trial_onset_gap_s trials_rows stim_trials trial_index_values trial_pulses reps units cell_types 0 sub-ICMS43_ses-2022-01-11_ophys.nwb 2022-01-11 ICMS43 627.0 2,5,7,10 42,22,3,19,12,20,22 False 171.0 6.0 NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN 1 sub-ICMS43_ses-2022-01-21_ophys.nwb 2022-01-21 ICMS43 403.0 2,5,7,10 51,16,16,40,52,11,16 False 156.0 6.0 NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN 3 sub-ICMS43_ses-2022-01-27_ophys.nwb 2022-01-27 ICMS43 448.0 2,5,7,10 7,23,1,55,6,2,1 False 153.0 6.0 NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN 2 sub-ICMS43_ses-2022-02-05_ophys.nwb 2022-02-05 ICMS43 416.0 2,5,7,10 101,30,0,0,0,50,19 False 156.0 6.0 NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN 4 sub-ICMS43_ses-2022-02-12_ophys.nwb 2022-02-12 ICMS43 488.0 2,5,7,10 68,37,26,49,21 False 156.0 6.0 NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN 5 sub-ICMS43_ses-2022-02-19_ophys.nwb 2022-02-19 ICMS43 541.0 2,5,7,10 101,29,0,37,87,14 False 151.0 6.0 NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN 6 sub-ICMS43_ses-2022-02-25_ophys.nwb 2022-02-25 ICMS43 395.0 2,5,7,10 75,20,22,24,49,14,0 False 161.0 6.0 NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN 8 sub-ICMS43_ses-2022-03-03_ophys.nwb 2022-03-03 ICMS43 250.0 2,5,7,10 11,28,64,1 False 161.0 6.0 NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN 7 sub-ICMS45_ses-2022-02-05_ecephys.nwb 2022-02-05 ICMS45 NaN NaN NaN NaN NaN NaN 10 50.0 3.0 200.0 0.500-0.501 10.00-10.00 17.0 17.0 1.0 100 17.0 8.0 unknown 9 sub-ICMS45_ses-2022-02-12_ecephys+ophys.nwb 2022-02-12 ICMS45 53.0 2,5,7,10 3,5,20,1,3,19,4,12 False 131.0 5.0 10 50.0 3.0 200.0 0.500-0.502 10.00-10.00 96.0 96.0 32.0 30 3.0 7.0 unknown 10 sub-ICMS45_ses-2022-02-21_ophys.nwb 2022-02-21 ICMS45 31.0 2,5,7,10 0,1,1,0,0,2,0,0 False 116.0 5.0 NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN 12 sub-ICMS45_ses-2022-02-26_ecephys+ophys.nwb 2022-02-26 ICMS45 60.0 2,5,7,10 4,8,1,1,0,9,0,3 False 136.0 5.0 10 50.0 3.0 200.0 0.500-0.504 9.94-10.00 160.0 160.0 32.0 30 5.0 8.0 unknown 11 sub-ICMS45_ses-2022-03-04_ecephys.nwb 2022-03-04 ICMS45 NaN NaN NaN NaN NaN NaN 10 50.0 3.0 200.0 0.500-0.503 9.95-10.00 224.0 224.0 32.0 30 7.0 8.0 unknown 15 sub-ICMS48_ses-2022-04-08_ecephys+ophys.nwb 2022-04-08 ICMS48 721.0 2,5,7,10 53,42,79,114,0,67,71,66,115 False 151.0 2.0 10 50.0 3.0 200.0 0.500-0.515 10.00-10.00 180.0 180.0 36.0 30 5.0 3.0 unknown 13 sub-ICMS48_ses-2022-04-15_ecephys+ophys.nwb 2022-04-15 ICMS48 689.0 2,5,7,10 46,57,46,28,45,58,112 False 151.0 2.0 10 50.0 3.0 200.0 0.500-0.517 9.94-10.00 196.0 196.0 28.0 30 7.0 2.0 unknown 14 sub-ICMS48_ses-2022-04-21_ecephys.nwb 2022-04-21 ICMS48 NaN NaN NaN NaN NaN NaN 10 50.0 3.0 200.0 0.500-0.515 9.94-10.00 140.0 140.0 28.0 30 5.0 1.0 unknown 16 sub-ICMS54_ses-2022-05-20_ecephys.nwb 2022-05-20 ICMS54 NaN NaN NaN NaN NaN NaN 10 50.0 3.0 200.0 0.500-0.515 10.00-10.00 126.0 126.0 10.0 210,340 13.0 13.0 unknown 17 sub-ICMS54_ses-2022-05-25_ecephys+ophys.nwb 2022-05-25 ICMS54 115.0 2,5,7,10 1,7,9,16,13,4,3,3,19,10 False 125.0 9.0 10 50.0 3.0 200.0 0.161-0.515 9.96-10.00 185.0 185.0 58.0 30,132,404 5.0 0.0 18 sub-ICMS54_ses-2022-06-01_ecephys+ophys.nwb 2022-06-01 ICMS54 229.0 2,5,7,10 0,10,9,11,2,2,1,2,7 False 63.0 8.0 10 50.0 3.0 200.0 0.500-0.515 9.95-10.00 360.0 360.0 36.0 30 10.0 8.0 unknown 21 sub-ICMS54_ses-2022-06-07_ophys.nwb 2022-06-07 ICMS54 29.0 2,5,7,10 0,1,3,0,0,1,0,1 False 63.0 7.0 NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN 20 sub-ICMS56_ses-2022-06-14_ophys.nwb 2022-06-14 ICMS56 1.0 2,5,7,10 0,0,0,0,0,0,0,0 False 51.0 5.0 NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN 19 sub-ICMS56_ses-2022-06-21_ecephys.nwb 2022-06-21 ICMS56 NaN NaN NaN NaN NaN NaN 10 50.0 3.0 200.0 0.500-0.515 9.96-10.00 352.0 352.0 32.0 30 11.0 10.0 unknown 23 sub-ICMS56_ses-2022-06-30_ecephys+ophys.nwb 2022-06-30 ICMS56 491.0 2,5,7,10 8,4,0,35,47,20,0,1,33,9 False 76.0 5.0 10 50.0 3.0 200.0 0.500-0.514 9.68-10.00 360.0 360.0 40.0 29,30 9.0 6.0 unknown 22 sub-ICMS56_ses-2022-07-06_ecephys.nwb 2022-07-06 ICMS56 NaN NaN NaN NaN NaN NaN 10 50.0 3.0 200.0 0.500-0.514 9.94-10.00 144.0 144.0 36.0 30 4.0 8.0 unknown 27 sub-ICMS56_ses-2022-07-11_ecephys+ophys.nwb 2022-07-11 ICMS56 127.0 2,5,7,10 3,3,0,0,11,1,0,1,0 False 61.0 5.0 10 50.0 3.0 200.0 0.161-0.514 9.95-10.00 324.0 323.0 53.0 2,30,224 9.0 10.0 unknown 24 sub-ICMS56_ses-2022-07-15_ecephys+ophys.nwb 2022-07-15 ICMS56 76.0 2,5,7,10 3,1,2,0,4,2,0,0,0 False 63.0 5.0 10 50.0 3.0 200.0 0.162-0.515 10.00-10.00 362.0 362.0 55.0 30,62,224 10.0 6.0 unknown 25 sub-ICMS56_ses-2022-07-23_ophys.nwb 2022-07-23 ICMS56 23.0 2,5,7,10 0,2,2,0,2,0,0,0,0 False 64.0 5.0 NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN 28 sub-ICMS57_ses-2022-07-27_ophys.nwb 2022-07-27 ICMS57 196.0 2,5,7,10 32,0,83 False 81.0 6.0 NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN 26 sub-ICMS57_ses-2022-08-02_ophys.nwb 2022-08-02 ICMS57 16.0 2,5,7,10 3,0 False 63.0 6.0 NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN 29 sub-ICMS57_ses-2022-08-10_ophys.nwb 2022-08-10 ICMS57 74.0 2,5,7,10 9,2,19 False 63.0 6.0 NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN descriptions of the control electrical_stimulation tables: ['Control passive-stim protocol: 3 bursts per trial, 10 biphasic pulses per burst at 50 Hz, 167 us per phase. One row per burst.'] descriptions of the control trials tables' trial_index: ['Trial index within session'] | trials per distinct trial_index value: trials_rows / trial_index_values above control mice: 6 | control files: 30 | with survey imaging: 24 | with units: 14
figs2de-control-stable (count half)¶
In control animals stimulated without behavioural training, activation density within 200 um and the number of activated neurons do not change from early to late weeks.
The control survey files store each ROI's binary activation per stimulation channel and current; the number of activated ROIs per channel and current is their sum. Density within 200 um needs ROI distances from the stimulation site, which need the voxel size and the contact positions that the files do not give (data check above), so only the count half is run. The 800 um limit of the count cannot be applied for the same reason; all activated ROIs are counted. Weeks are counted from each control animal's first imaging session (common method). The controls were imaged at 2, 5, 7 and 10 uA, so of the 4 to 6 uA that reading (b) pools only 5 uA exists. Readings: (a) 5 uA, late weeks 2 and 3, (b) 5 uA, late weeks 2 to 4; each on channel x session values and on channel x week medians (the units of fig2c); all currents pooled is a further variant.
In [34]:
crows = []
for a in ctrl:
f = icms.open_asset(a["asset_id"])
o = "processing/ophys"
if o + "/stimulus_conditions" in f:
act = f[o + "/ImageSegmentation_ControlSurvey/PlaneSegmentation_ControlSurvey/activation"][()]
cur = f[o + "/stimulus_conditions/current_uA"][()]
chn = f[o + "/stimulus_conditions/stim_channel"][()].astype(str).astype(int)
sub, date = f["general/subject/subject_id"][()].decode(), f["session_start_time"][()].decode()[:10]
for j in range(len(cur)):
crows.append({"subject": sub, "date": date, "channel": chn[j], "current": cur[j], "count": int(act[:, j].sum())})
f.close()
CC = pd.DataFrame(crows)
_cs = [{"subject": a, "date": b} for a, b in CC[["subject", "date"]].drop_duplicates().values]
icms.add_weeks(_cs)
CC["week"] = CC.set_index(["subject", "date"]).index.map({(s["subject"], s["date"]): s["week"] for s in _cs})
print("control imaging sessions per mouse:", CC.groupby("subject").date.nunique().to_dict(),
"| channel-current values per week:", CC.week.value_counts().sort_index().to_dict())
print("count range:", CC["count"].min(), "to", CC["count"].max(), "| currents:", sorted(CC.current.unique()))
res_s2de = {}
for reading, late, curs in [("(a) 5 uA, late weeks 2-3", (2, 3), (5.0,)), ("(b) 5 uA, late weeks 2-4", (2, 3, 4), (5.0,)),
("variant, all currents, late weeks 2-4", (2, 3, 4), (2.0, 5.0, 7.0, 10.0))]:
d = CC[CC.current.isin(curs) & (CC.week <= 4)]
t1 = icms.mwu(d.loc[d.week.isin(icms.EARLY_WEEKS), "count"], d.loc[d.week.isin(late), "count"])
dm = d.groupby(["subject", "channel", "current", "week"], as_index=False)["count"].median()
t2 = icms.mwu(dm.loc[dm.week.isin(icms.EARLY_WEEKS), "count"], dm.loc[dm.week.isin(late), "count"])
res_s2de[reading] = (t1, t2)
for unit, t in [("channel x session", t1), ("channel x week medians", t2)]:
print(f"{reading}, {unit}: early median {t['med1']:.1f} (n={t['n1']}), late median {t['med2']:.1f} (n={t['n2']}), "
f"p = {t['p']:.3g}, r = {r_from_p(t['p'], t['n1'], t['n2']):.2f}")
print("median count at 5 uA by mouse and week:")
print(CC[CC.current == 5].pivot_table(index="subject", columns="week", values="count", aggfunc="median"))
control imaging sessions per mouse: {'ICMS43': 8, 'ICMS45': 3, 'ICMS48': 2, 'ICMS54': 3, 'ICMS56': 5, 'ICMS57': 3} | channel-current values per week: {0: 180, 1: 72, 2: 136, 3: 68, 4: 64, 5: 56, 6: 60, 7: 28, 8: 16}
count range: 0 to 311 | currents: [np.float64(2.0), np.float64(5.0), np.float64(7.0), np.float64(10.0)]
(a) 5 uA, late weeks 2-3, channel x session: early median 10.0 (n=63), late median 2.0 (n=51), p = 0.0146, r = 0.23
(a) 5 uA, late weeks 2-3, channel x week medians: early median 10.0 (n=63), late median 4.0 (n=43), p = 0.0816, r = 0.17
(b) 5 uA, late weeks 2-4, channel x session: early median 10.0 (n=63), late median 2.0 (n=67), p = 0.00573, r = 0.24
(b) 5 uA, late weeks 2-4, channel x week medians: early median 10.0 (n=63), late median 3.0 (n=59), p = 0.0261, r = 0.20
variant, all currents, late weeks 2-4, channel x session: early median 10.0 (n=252), late median 2.0 (n=268), p = 0.000173, r = 0.16
variant, all currents, late weeks 2-4, channel x week medians: early median 10.0 (n=252), late median 2.8 (n=236), p = 0.00221, r = 0.14
median count at 5 uA by mouse and week:
week 0 1 2 3 4 5 6 7 8
subject
ICMS43 20.0 NaN 16.0 6.0 19.0 37.0 33.0 22.0 19.5
ICMS45 4.5 NaN 1.0 NaN NaN NaN NaN NaN NaN
ICMS48 67.0 46.0 NaN NaN NaN NaN NaN NaN NaN
ICMS54 8.0 2.0 0.5 NaN NaN NaN NaN NaN NaN
ICMS56 0.0 NaN NaN 8.5 1.0 1.0 0.0 NaN NaN
ICMS57 32.0 1.5 9.0 NaN NaN NaN NaN NaN NaN
Result: blocked (density cannot be formed); the count half does not show stability. At 5 uA the number of activated ROIs per channel and session falls from a median of 10 early to 2 late (reading a, late weeks 2 and 3: p = 0.0146, r = 0.23; reading b, weeks 2 to 4: p = 0.00573), and with channel x week medians p = 0.0816 (a) and 0.0261 (b). So the count test gives P >= 0.05 under one of the four combinations, where the paper reports p = 0.066 (n 42 and 26; here 63 and 43 to 67). The density half, which the criterion also requires, cannot be formed.
figs7-control-no-change¶
In control animals stimulated without behavioural training, neither PL nor NPL responses change in modulation, time to maximum firing rate or PL peak spike probability.
The paper says only that the controls had similar stimulation parameters, but its code runs the
control units through the behavioural response classes with a 20 ms pulse window
(processing/control/postprocessing/create_session_responses.py) and tests them at 5 uA with its own
filters (python/fig_s7/generate_figure.py: P < 0.05, z > 0, baseline at least 0.5 Hz, more than 5
trains, more than 15 pulse-window spikes, weeks below 5). That code is followed (icms.control_responses):
each 200 ms burst is a train, the baseline is the 700 ms before it, pulse-locking uses 20 ms windows.
The code's timing file (control_timing_params.json) is not in the repository, so the 20 ms window
comes from the code's docstring and the blanking (1.5 ms and 0.5 ms) is the behavioural one; time to
maximum is taken over the 200 ms burst. These are assumptions. Weeks count from each animal's first
electrophysiology session, the first dates the code lists (python/utils/config.py, CONTROL_SESSIONS).
Readings: (a) 5 uA (code) and (b) all currents pooled; the controls have no 4 or 6 uA.
The bursts of a trial start 0.5 s apart (control data check above), so the 700 ms baseline of the second and third burst holds the burst before it. A baseline with stimulation in it can hide a change, and this is a claim of absence, so a variant takes every burst's baseline from the 700 ms before the first burst of its trial (trials are 10 s apart), with the same pulse-locking draws. The claim is judged on the code's test, as followed above (protocol, claims of absence); the variant is shown.
In [35]:
cephys = []
for a in ctrl:
f = icms.open_asset(a["asset_id"])
has_units = "units" in f
f.close()
if has_units:
cephys.append(icms.load_control_ephys(a))
cephys.sort(key=lambda s: (s["subject"], s["date"]))
icms.add_weeks(cephys)
print("control ephys sessions with units:", len(cephys), "| mice:", sorted({s["subject"] for s in cephys}))
print([(s["subject"], s["date"], s["week"], len(s["spikes"])) for s in cephys])
t0 = time.time()
CR = pd.DataFrame(icms.control_responses(cephys, rng_seed=40, currents=None))
CR["mod"] = ((CR.p < 0.05) & (CR.z > 0) & (CR.base_rate >= 0.5) & (CR.num_trials > 5) & (CR.num_spikes > 15)
& (CR.week < 5))
print(f"control responses: {len(CR)} in {time.time() - t0:.0f} s; bursts per response {CR.num_trials.min()} to "
f"{CR.num_trials.max()}; modulated {int(CR['mod'].sum())}")
res_s7 = {}
for reading, curs in [("(a) 5 uA", (5.0,)), ("(b) all currents", tuple(sorted(CR.current.unique())))]:
d = CR[CR["mod"] & CR.current.isin(curs)]
print(f"\nreading {reading}: modulated {len(d)}, PL {int(d.pl.sum())}; by phase and label:",
pd.crosstab(d.phase, d.pl.map({True: "PL", False: "NPL"})).to_dict())
out = {}
for lab, g, col in [("PL", d[d.pl], "z"), ("NPL", d[~d.pl], "z"), ("PL", d[d.pl], "tmax_ms"),
("NPL", d[~d.pl], "tmax_ms"), ("PL", d[d.pl], "peak_prob")]:
t = icms.mwu(g.loc[g.phase == "early", col], g.loc[g.phase == "late", col])
out[(lab, col)] = t
print(f" {lab} {col}: early median {t['med1']:.3f} (n={t['n1']}), late median {t['med2']:.3f} (n={t['n2']}), "
f"p = {t['p']:.3g}, r = {t['r']:.2f}")
res_s7[reading] = out
# How many bursts have an earlier burst inside their 700 ms baseline (bursts 0.5 s apart, cell above)
_n_in, _n_all = 0, 0
for s in cephys:
st = s["stim"]
keep = st["stop_time"] < s["ephys_end"]
on, off = st["start_time"][keep], st["stop_time"][keep]
for a in on:
_n_all += 1
_n_in += bool(np.any((off > a - icms.CONTROL_PRE_S) & (on < a)))
print(f"\nbursts whose 700 ms baseline holds an earlier burst: {_n_in} of {_n_all}")
# Variant: the baseline of every burst is the 700 ms before the first burst of its trial (no stimulation
# in it). Same seed, so the pulse-locking draws and PL labels are those of the main reading.
CRt = pd.DataFrame(icms.control_responses(cephys, rng_seed=40, currents=None, baseline="trial"))
CRt["mod"] = ((CRt.p < 0.05) & (CRt.z > 0) & (CRt.base_rate >= 0.5) & (CRt.num_trials > 5) & (CRt.num_spikes > 15)
& (CRt.week < 5))
print("PL labels identical to the main reading:", bool((CRt.pl.values == CR.pl.values).all()))
dt = CRt[CRt["mod"] & (CRt.current == 5)]
print(f"variant, baseline before the trial's first burst, reading (a) 5 uA: modulated {len(dt)}; by phase and label:",
pd.crosstab(dt.phase, dt.pl.map({True: "PL", False: "NPL"})).to_dict())
res_s7_trial = {}
for lab, g, col in [("PL", dt[dt.pl], "z"), ("NPL", dt[~dt.pl], "z"), ("PL", dt[dt.pl], "tmax_ms"),
("NPL", dt[~dt.pl], "tmax_ms"), ("PL", dt[dt.pl], "peak_prob")]:
t = icms.mwu(g.loc[g.phase == "early", col], g.loc[g.phase == "late", col])
res_s7_trial[(lab, col)] = t
print(f" {lab} {col}: early median {t['med1']:.3f} (n={t['n1']}), late median {t['med2']:.3f} (n={t['n2']}), "
f"p = {t['p']:.3g}, r = {t['r']:.2f}")
d = CR[CR["mod"] & (CR.current == 5)]
print("\n5 uA modulated responses by mouse and phase:", pd.crosstab(d.subject, d.phase).to_dict())
print("z range", d.z.min().round(2), "to", d.z.max().round(2), "| time to max range", d.tmax_ms.min(), "to", d.tmax_ms.max(), "ms")
control ephys sessions with units: 14 | mice: ['ICMS45', 'ICMS48', 'ICMS54', 'ICMS56']
[('ICMS45', '2022-02-05', 0, 8), ('ICMS45', '2022-02-12', 1, 7), ('ICMS45', '2022-02-26', 3, 8), ('ICMS45', '2022-03-04', 4, 8), ('ICMS48', '2022-04-08', 0, 3), ('ICMS48', '2022-04-15', 1, 2), ('ICMS48', '2022-04-21', 2, 1), ('ICMS54', '2022-05-20', 0, 13), ('ICMS54', '2022-06-01', 2, 8), ('ICMS56', '2022-06-21', 0, 10), ('ICMS56', '2022-06-30', 2, 6), ('ICMS56', '2022-07-06', 3, 8), ('ICMS56', '2022-07-11', 3, 10), ('ICMS56', '2022-07-15', 4, 6)]
control responses: 2958 in 39 s; bursts per response 3 to 36; modulated 309
reading (a) 5 uA: modulated 111, PL 51; by phase and label: {'NPL': {'early': 27, 'late': 33}, 'PL': {'early': 32, 'late': 19}}
PL z: early median 1.456 (n=32), late median 1.878 (n=19), p = 0.946, r = -0.01
NPL z: early median 0.324 (n=27), late median 0.529 (n=33), p = 0.122, r = 0.23
PL tmax_ms: early median 5.000 (n=32), late median 55.000 (n=19), p = 0.112, r = 0.25
NPL tmax_ms: early median 65.000 (n=27), late median 75.000 (n=33), p = 0.399, r = 0.13
PL peak_prob: early median 0.097 (n=32), late median 0.075 (n=19), p = 0.182, r = -0.23
reading (b) all currents: modulated 309, PL 132; by phase and label: {'NPL': {'early': 74, 'late': 103}, 'PL': {'early': 59, 'late': 73}}
PL z: early median 1.450 (n=59), late median 2.050 (n=73), p = 0.0804, r = 0.18
NPL z: early median 0.260 (n=74), late median 0.446 (n=103), p = 0.00809, r = 0.23
PL tmax_ms: early median 5.000 (n=59), late median 55.000 (n=73), p = 0.0352, r = 0.20
NPL tmax_ms: early median 65.000 (n=74), late median 85.000 (n=103), p = 0.36, r = 0.08
PL peak_prob: early median 0.089 (n=59), late median 0.086 (n=73), p = 0.815, r = -0.02
bursts whose 700 ms baseline holds an earlier burst: 5863 of 8795
PL labels identical to the main reading: True
variant, baseline before the trial's first burst, reading (a) 5 uA: modulated 107; by phase and label: {'NPL': {'early': 29, 'late': 28}, 'PL': {'early': 30, 'late': 20}}
PL z: early median 2.456 (n=30), late median 3.379 (n=20), p = 0.744, r = 0.06
NPL z: early median 0.334 (n=29), late median 0.988 (n=28), p = 0.0059, r = 0.43
PL tmax_ms: early median 10.000 (n=30), late median 5.000 (n=20), p = 0.957, r = 0.01
NPL tmax_ms: early median 55.000 (n=29), late median 85.000 (n=28), p = 0.44, r = 0.12
PL peak_prob: early median 0.104 (n=30), late median 0.081 (n=20), p = 0.07, r = -0.31
5 uA modulated responses by mouse and phase: {'early': {'ICMS45': 0, 'ICMS48': 5, 'ICMS54': 29, 'ICMS56': 25}, 'late': {'ICMS45': 6, 'ICMS48': 1, 'ICMS54': 28, 'ICMS56': 17}}
z range 0.01 to 5.2 | time to max range 5.0 to 195.0 ms
Result: reproduced, resting on the assumptions about the control code's timing. At 5 uA (reading a, the code's current), 111 control responses are modulated, PL 32 early and 19 late, NPL 27 and 33 (paper 29 and 12, 30 and 29). None of the five tests is significant: PL modulation p = 0.946, NPL modulation p = 0.122, PL time to maximum p = 0.112, NPL time to maximum p = 0.399, PL peak spike probability p = 0.182. Pooling all currents (reading b, n 59/73 PL and 74/103 NPL, further from the paper's) gives NPL modulation p = 0.00809 and PL time to maximum p = 0.0352. The criterion judges under the reading whose n is closest to the paper's, (a). The 700 ms baseline holds an earlier burst for 5863 of 8795 bursts. With the baseline taken before each trial's first burst (variant), NPL modulation rises (0.334 to 0.988, p = 0.0059, r = 0.43; n 29 and 28) and the other four tests stay above 0.05 (PL modulation p = 0.744, PL time to maximum p = 0.957, NPL time to maximum p = 0.44, PL peak spike probability p = 0.07). So the absence of a change in NPL modulation holds with the code's baseline and not with a baseline free of stimulation. PL time to maximum also moves the other way from the paper's: here it rises (5 to 55 ms, r = 0.25), where the paper reports r = -0.39.
figs9-coupling-window-and-control (control part)¶
Coupling in the control animals at 5 uA (fig. S9C legend), early against late, PL against NPL, FDR
over the two tests. The paper's control code (processing/pop_coupling/pop_coupling_control.py)
concatenates, per channel and current, segments from the first to the last pulse of each row of its
trial table, builds the population from every unit with a response at that condition (not only the
modulated ones) and labels each unit PL or NPL by its call at that condition; units need 50 spikes.
Its trial table is not public, so whether a row is a trial (3 bursts) or a burst is open: readings
(s1) trial segments and (s2) burst segments. The stored trial_index restarts in each recording block
(the trials table's rep_num), so trials are numbered from the bursts' burst_index_in_trial in time
order (icms.control_trial_ids); the previous revision grouped by trial_index and so joined trials of different
blocks into one segment. Populations: (p1) all units (code), (p2) the modulated
units only (as fig6f).
In [36]:
res_s9c = {}
for seg in ["s1 trial", "s2 burst"]:
for popr in ["p1 all units", "p2 modulated"]:
rngc9 = np.random.default_rng(41)
out = []
for (si, ch, cur), g in CR[(CR.current == 5) & (CR.week < 5)].groupby(["si", "channel", "current"]):
s = cephys[si]
st = s["stim"]
m = (st["stim_channel"] == ch) & (st["current_uA"] == cur) & (st["stop_time"] < s["ephys_end"])
tid = icms.control_trial_ids(st) # trial_index restarts in each block (cell above)
if seg.startswith("s2"):
segs = [(a, a + 0.18) for a in np.sort(st["start_time"][m])]
else:
segs = [(st["start_time"][m & (tid == k)].min(),
st["start_time"][m & (tid == k)].max() + 0.18) for k in np.unique(tid[m])]
use = g if popr.startswith("p1") else g[g["mod"]]
if len(segs) < 3 or len(use) < 2:
continue
pc = icms.coupling_on_segments({u: s["spikes"][u] for u in use.unit}, sorted(segs), rngc9)
lab = dict(zip(use.unit, use.pl))
out += [{"phase": g.phase.iloc[0], "subject": s["subject"], "pl": lab[u], "c": v} for u, v in pc.items()]
O = pd.DataFrame(out)
O = O[np.isfinite(O.c)]
tests = [icms.mwu(O[O.pl & (O.phase == ph)].c, O[~O.pl & (O.phase == ph)].c) for ph in ["early", "late"]]
q = icms.fdr_bh([t["p"] for t in tests])
res_s9c[(seg, popr)] = (tests, q)
print(f"{seg}, {popr}: " + "; ".join(
f"{ph} PL {t['med1']:.2f} (n={t['n1']}) vs NPL {t['med2']:.2f} (n={t['n2']}), FDR p = {qq:.3g}, r = {t['r']:.2f}"
for ph, t, qq in zip(["early", "late"], tests, q)))
s1 trial, p1 all units: early PL 2.03 (n=45) vs NPL 1.56 (n=48), FDR p = 0.431, r = -0.10; late PL 2.46 (n=31) vs NPL 1.58 (n=90), FDR p = 0.134, r = -0.22 s1 trial, p2 modulated: early PL 2.12 (n=30) vs NPL 1.91 (n=15), FDR p = 0.99, r = -0.00; late PL 2.37 (n=13) vs NPL 3.96 (n=15), FDR p = 0.99, r = 0.04
s2 burst, p1 all units: early PL -1.09 (n=41) vs NPL 0.59 (n=18), FDR p = 0.00428, r = 0.51; late PL 0.47 (n=28) vs NPL 0.75 (n=46), FDR p = 0.0831, r = 0.24 s2 burst, p2 modulated: early PL -1.18 (n=27) vs NPL 1.03 (n=6), FDR p = 0.173, r = 0.37; late PL -0.15 (n=13) vs NPL 0.10 (n=10), FDR p = 0.173, r = 0.35
Result: reproduced, resting on the assumption about the control trial segments. The 120 ms part holds (above). In the controls at 5 uA, with the code's population of all units and segments from the first to the last pulse of each trial (s1, p1), PL and NPL coupling do not differ early (2.03 against 1.56, n 45 and 48, FDR p = 0.431) or late (2.46 against 1.58, n 31 and 90, FDR p = 0.134); the paper reports n 49 and 64, 25 and 128, FDR p = 0.48. With the modulated units as population (s1, p2), FDR p = 0.99 in both phases. With burst segments (s2), the early difference is significant with all units (FDR p = 0.00428, n 41 and 18) and not with modulated units (0.173). The trial reading is taken because the criterion says the full train, the dataset's own trials table has one row per trial of 3 bursts (as many rows as such trials in 14 of the 15 files with stimulation, and n_pulses 30, the pulses of 3 bursts, in most; control data check), and its n are closer to the paper's in both phases (45/48 and 31/90 against 41/18 and 28/46 for bursts).
Summary¶
A failure to reproduce here may be a failure of this attempt, not of the paper: see the reasoning and the differences from the paper's method.
| Claim | Outcome |
|---|---|
| fig2a-threshold-decrease | reproduced |
| fig2a-steepest-decline-early | reproduced (rests on the criterion's reading of the figure) |
| figs1e-response-times-decrease | reproduced |
| figs5c-wheel-onset-decrease | reproduced |
| fig4-unit-counts | not reproduced |
| figs3b-modulation-increase-all | reproduced |
| figs3c-latency-all | partially reproduced |
| figs3d-modulation-stable-at-threshold | reproduced |
| fig4g-modulated-count-mouse1 | reproduced |
| fig4ij-pyramidal-modulation-increase | reproduced |
| fig4ij-interneurons-little-change | reproduced |
| fig5d-pulse-locking-stable-across-currents | not reproduced |
| fig5e-pyramidal-fraction-pl | partially reproduced |
| figs4b-pl-latency-distance | partially reproduced |
| fig5fg-pl-modulation-increase | reproduced |
| figs4de-npl-no-change | reproduced |
| fig5h-pl-spike-probability-increase | partially reproduced |
| figs5-restricted-window-consistent | reproduced |
| fig6c-kl-divergence | not reproduced |
| fig6d-npl-hit-miss-divergence | not reproduced |
| fig6f-population-coupling | reproduced |
| figs9a-coupling-700ms | reproduced |
| figs9-coupling-window-and-control | reproduced (rests on the control trial segments) |
| figs7-control-no-change | reproduced (rests on the control code's timing) |
| fig2c, fig2b, figs2ab, fig2e (two), figs2de-behaviour, fig2g | blocked (metadata insufficient) |
| figs2de-control-stable | blocked (metadata insufficient; the count half is reported) |
| fig3b, fig3df, fig3i, fig3j (two), fig3lm, fig3n (two), figs6 | blocked (metadata insufficient) |
| fig4d-waveform-separation, figs8-tracked-units | blocked (data missing) |
How the analysis differs from the paper's. Thresholds are refitted with the Python port of
psignifit; the dataset stores none. The movement-onset threshold is taken from hit trials, as miss
trials have no response time. Pulse-locking follows the paper's code, with pulses placed 10 ms
apart from the train onsets, as no pulse times are stored; the call is random, and the claims
that rest on it near a threshold (fig5d, fig5e, fig5h) are shown under three other seeds. The
hit-miss PSTH bin width and smoothing are not in the paper or its code; two widths were run. The
stimulation contact is placed as the paper's code places it, not by the electrode table. The
control analyses follow the paper's control code, whose timing file is not public; its 700 ms
baseline holds the previous burst for most bursts, and a baseline free of stimulation is shown as a
variant. Control trials are numbered from burst_index_in_trial, as the stored trial_index
restarts in each recording block. Cell types are
the stored labels, which the files describe as the trough-to-peak classification; the rule applied
to the stored mean waveform is shown as a variant.
What was found about the dataset. The 46 sessions of weeks 0 to 4 match the paper's count, and the
ephys response counts are close to the paper's, but the unit total (795 in weeks 0 to 4) is below
the paper's 973. Behavioural files give no stimulation-site position in the imaging volume and no
activation calls, and the planar ROIs are segmented per session with no identity across sessions,
so no behavioural imaging claim can be tested. The stored stim_channel does not index the
electrode table's channel_name, though both are described as Ripple channel numbers. The control
files use 3 bursts of 10 pulses at 50 Hz, where the paper says the controls had similar stimulation
parameters. In the control files trial_index, described as the trial index within the session, restarts in
each recording block. Most miss trials carry no response time, Mouse 6 has no catch trials, and cell types
are spelled differently in Mouse 6.