Hires et al. (2015) from dandiset 000013¶
Paper. Hires, Gutnisky, Yu, O'Connor & Svoboda (2015), "Low-noise encoding of active touch by layer 4 in the somatosensory cortex", eLife 4, e06619, doi:10.7554/eLife.06619.
Dataset. DANDI:000013, version
0.220126.2143 (doi:10.48324/dandi.000013/0.220126.2143).
The paper (2015) cites no version of the dandiset, which was created in 2020. This is its only
published version, and so the most recent one; the draft holds the same 52 files (the same
blobs, checked through the API when the attempt was prepared), so no other version could match
the claims file better.
Claims. claims/000013/10.7554-elife.06619.yaml, written by a reader from the version of
record (PMC4525079, as served by Europe PMC), its figure legends and figure supplements, and
checked by a second reader. The claims file decides what is tested and how; where it leaves a
method detail open, every reading it lists is run, and where it lists none for a detail that
can change a result, the question is asked in questions/000013.yaml and the readings this
notebook can see are all shown.
The data are streamed from the archive; nothing is written to disk. The helper module
hires.py holds the loading and every reusable measure, with the source of each parameter.
| Claim | Tested here |
|---|---|
| recordings-counts, recordings-mice | yes |
| behavior-performance, behavior-sample-epoch, behavior-whisking-bouts, behavior-touches-reaction | yes |
| example-cell-rates, example-cell-sample-ff, example-cell-touch-window, example-cell-ff-around-touch, example-cell-phase-modulation | yes |
| text-spikes-per-touch, text-whisking-rates | yes |
| table1-l4-c2, table1-l4-outside, table1-l5, latency-determinable | yes |
| l4-outside-weaker-touch, l5-higher-rates, phase-mod-l4-inside-outside, phase-mod-l4-vs-l5, whisking-rate-l4-no-change, whisking-rate-l5-increase | yes |
| fig3c-touch-spikes-l4c2, fig3c-touch-spikes-outside, fig3c-touch-spikes-l5 | yes |
| fig3d-touch-phase-l4c2, fig3d-touch-phase-outside, fig3d-touch-phase-l5, fig3d-phase-only | yes |
| fig4a-touch-detection-c2, fig4a-touch-detection-outside, fig4b-touch-time-precision, fig4c-phase-decoding | yes |
| fig5c-modulation-index, fig5-modulation-correlations | yes |
| fig6s1b-clusters-less-variable, fig6cd-ff-near-binomial, fig6d-not-refractory | yes |
| fig6e-sample-epoch-ff, fig6e-random-38ms-ff, fig6e-touch-38ms-ff, fig6e-touch-10ms-ff, fig6e-similar-touch-ff | yes |
| whisk-cycle-spikes | yes |
Every claim's data are in the dandiset, so none is blocked before the analysis; the outcomes,
some of them blocked as methods_underspecified where the readings disagree, are in the
summary at the end.
Answers. A first run asked seven questions (questions/000013.yaml). A reader answered them
from the paper, and the claims file now carries the answers. This run follows the updated
criteria: the touch-window cutoff is the transition where the curve levels off, the literal
search from 1 ms being a variant only (fig3c-touch-spikes-l4c2, fig3d-touch-phase-l4c2, and the
claims that use the touch window); the longest window of the cutoff curve is an open detail
with readings of 100 and 200 ms; touch detection adds the reading with every pre-touch window
(fig4a-touch-detection-c2); the time decoder's classes and smoothing are readings
(fig4b-touch-time-precision); the phase decoder uses the cycles of phase tuning, with two
readings for empty bins (fig4c-phase-decoding); the refractory model's rate is raised so that
its output matches the cell's touch response (fig6d-not-refractory); and the L5 whisking
increase is the median of the per-cell differences (whisking-rate-l5-increase).
Verification, round 1. An independent verifier found three problems that bear on conclusions, and this run corrects them. Whisking periods now leave out licking cycles and the 70 ms after each touch, as the claims file's definitions say (they were left in before). Phase tuning now uses every phase-analysis cycle of the trial, as the criteria say, and not only the cycles inside the sample epoch, a restriction the criteria do not give; the sample-epoch cycles are kept for example-cell-phase-modulation, whose criterion names exploration periods, and as a variant elsewhere. The touch detector's fraction of touches detected is now read at the time point of best percent correct, with ties between thresholds broken at the middle one; the earlier rule took the highest hit rate over time points, which rewards a threshold that calls everything a touch. Touch detection now runs 50 repetitions per population size (20 before). The outputs the verifier asked for are added: per-file touch peaks against the clock offset, per-cell touch windows for the L4 outside and L5 cells, the matching rows of example-cell-ff-around-touch, the detector at N = 1, and phase depths with at least 20 spikes.
import os
import time
import warnings
os.environ.setdefault("OMP_NUM_THREADS", "1")
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy import stats
import hires
warnings.filterwarnings("ignore", category=RuntimeWarning)
pd.set_option("display.width", 200)
pd.set_option("display.max_columns", 40)
pd.set_option("display.max_rows", 200)
T_START = time.time()
t0 = time.time()
cells = hires.load_all(max_workers=8)
print(f"loaded {len(cells)} files from dandiset {hires.DANDISET_ID} version {hires.DANDISET_VERSION} in {time.time() - t0:.0f} s")
GROUP_ORDER = ["L4C2", "L4out", "L5"]
by_no = {c.cell_no: c for c in cells}
by_path = {c.path: c.cell_no for c in cells}
loaded 52 files from dandiset 000013 version 0.220126.2143 in 37 s
Data check¶
What was loaded, and each description checked against the values.
meta = pd.DataFrame([dict(cell=c.cell_no, identifier=c.identifier, subject=c.subject,
genotype=c.genotype or "(blank)", region=c.region, layer=c.layer,
depth_um=c.depth, group=c.group, trials=len(c.trials),
spikes=len(c.spike_t), behaviour_ms=len(c.vt),
spike_block=c.spike_block_n[0], spike_offset_ms=1000 * c.extra["spike_offset"],
volt_minus_spike_ms=round(1000 * c.volt_offset, 3),
spot_check_max_err_s=c.extra["spot_check_err"])
for c in cells]).set_index("cell")
print(meta.drop(columns=["identifier"]).to_string())
print("\nfiles", len(meta), "| subjects", meta.subject.nunique(), "| cells per group", meta.group.value_counts().to_dict())
print("trials per file: mean %.1f, s.d. %.1f, min %d, max %d" % (meta.trials.mean(), meta.trials.std(), meta.trials.min(), meta.trials.max()))
print("cell numbers in the identifiers:", sorted(meta.index) == list(range(1, 53)))
subject genotype region layer depth_um group trials spikes behaviour_ms spike_block spike_offset_ms volt_minus_spike_ms spot_check_max_err_s
cell
1 anm244028 PV-Cre C2 4 472.0 L4C2 48 210 240000 50000 0.0 0.0 0.000000e+00
2 anm266951 PV-Cre C2 4 501.0 L4C2 19 283 95000 50000 0.0 0.0 0.000000e+00
3 anm244028 PV-Cre C2 4 559.0 L4C2 64 829 320000 50000 0.0 0.0 5.684342e-14
4 anm244025 PV-Cre C2 4 447.0 L4C2 55 168 275000 50000 0.0 0.0 0.000000e+00
5 anm266951 PV-Cre C2 4 444.0 L4C2 34 102 170000 50000 0.0 0.0 0.000000e+00
6 anm244024 PV-Cre C2 4 532.0 L4C2 48 95 240000 50000 0.0 0.0 4.440892e-16
7 anm244024 PV-Cre C2 4 465.0 L4C2 26 35 130000 50000 0.0 0.0 0.000000e+00
8 anm266945 PV-Cre C2 4 552.0 L4C2 30 202 150000 50000 0.0 0.0 0.000000e+00
9 anm215592 (blank) C2 4 471.0 L4C2 122 9796 610000 50000 -10.0 10.0 1.136868e-13
10 anm215592 (blank) C2 4 560.0 L4C2 241 2009 1205000 50000 -10.0 10.0 2.273737e-13
11 anm234232 PV-Cre C2 4 470.0 L4C2 80 1052 400000 50000 -10.0 10.0 1.136868e-13
12 anm226694 Ai32 x PV-Cre C2 4 444.0 L4C2 145 242 725000 50000 -10.0 10.0 2.273737e-13
13 anm266945 PV-Cre C2 4 539.0 L4C2 63 1305 315000 50000 0.0 0.0 4.440892e-16
14 jf35892 C57BL6 C2 4 565.0 L4C2 146 1097 657000 45000 -10.0 10.0 1.136868e-13
15 anm244024 PV-Cre C2 4 555.0 L4C2 133 781 665000 50000 0.0 0.0 2.220446e-16
16 jf37166 C57BL6 C2 4 555.0 L4C2 94 1581 423000 45000 -10.0 10.0 1.136868e-13
17 anm106211 C57BL6 C2 4 528.0 L4C2 166 1894 747166 45010 -10.0 10.0 2.273737e-13
18 anm226695 Ai32 x PV-Cre C2 4 457.0 L4C2 63 2505 315000 50000 -10.0 10.0 1.136868e-13
19 anm226695 Ai32 x PV-Cre C2 4 512.0 L4C2 98 130 490000 50000 -10.0 10.0 1.136868e-13
20 anm226694 Ai32 x PV-Cre C2 4 542.0 L4C2 158 478 790000 50000 0.0 0.0 2.273737e-13
21 anm226695 Ai32 x PV-Cre C2 4 500.0 L4C2 114 1486 570000 50000 -10.0 10.0 1.136868e-13
22 anm244024 PV-Cre C2 4 448.0 L4C2 24 107 120000 50000 0.0 0.0 5.551115e-17
23 jf35892 C57BL6 C2 4 442.0 L4C2 121 664 544500 45000 -10.0 10.0 1.136868e-13
24 anm244025 PV-Cre C2 4 450.0 L4C2 87 141 435000 50000 0.0 0.0 0.000000e+00
25 anm266945 PV-Cre C2 4 486.0 L4C2 103 725 515000 50000 0.0 0.0 0.000000e+00
26 jf35892 C57BL6 C2 4 483.0 L4C2 198 368 891000 45000 -10.0 10.0 2.273737e-13
27 anm17702 C57BL6 C2 4 410.0 L4C2 224 901 1008000 45000 -10.0 10.0 2.273737e-13
28 anm215592 (blank) C2 4 439.0 L4C2 88 187 440000 50000 -10.0 10.0 1.136868e-13
29 anm244028 PV-Cre C2 4 502.0 L4C2 50 128 250000 50000 0.0 0.0 0.000000e+00
30 jf35892 C57BL6 C2 4 536.0 L4C2 98 330 441000 45000 -10.0 10.0 1.136868e-13
31 anm226694 Ai32 x PV-Cre C2 4 547.0 L4C2 96 4371 480000 50000 -10.0 10.0 1.136868e-13
32 anm106213 C57BL6 surround 4 405.0 L4out 192 158 864192 45010 -10.0 10.0 2.273737e-13
33 anm106213 C57BL6 surround 4 475.0 L4out 187 111 841687 45010 -10.0 10.0 2.273737e-13
34 anm131970 C57BL6 surround 4 453.0 L4out 106 4964 477106 45010 -10.0 10.0 2.273737e-13
35 anm134333 C57BL6 surround 4 536.0 L4out 101 195 505101 50010 -10.0 10.0 1.136868e-13
36 anm211295 (blank) surround 4 549.0 L4out 172 1001 860000 50000 -10.0 10.0 2.273737e-13
37 anm215592 (blank) surround 4 510.0 L4out 205 719 1025000 50000 -10.0 10.0 2.273737e-13
38 anm217330 (blank) surround 4 560.0 L4out 148 796 740000 50000 -10.0 10.0 2.273737e-13
39 anm234230 PV-Cre surround 4 528.0 L4out 98 2232 490000 50000 -10.0 10.0 1.136868e-13
40 anm234231 PV-Cre surround 4 451.0 L4out 103 769 515000 50000 0.0 0.0 1.136868e-13
41 anm226692 Ai32 x PV-Cre surround 4 537.0 L4out 240 3298 1200000 50000 -10.0 10.0 2.273737e-13
42 anm17702 C57BL6 C2 5 844.0 L5 96 22298 432000 45000 -10.0 10.0 1.136868e-13
43 jf42400 C57BL6 C2 5 804.0 L5 119 3515 535619 45010 -10.0 10.0 2.273737e-13
44 jf42400 C57BL6 C2 5 830.0 L5 92 10921 414092 45010 -10.0 10.0 1.136868e-13
45 jf42400 C57BL6 C2 5 707.0 L5 93 656 418593 45010 -10.0 10.0 1.136868e-13
46 anm101105 C57BL6 C2 5 675.0 L5 117 7341 526617 45010 -10.0 10.0 1.136868e-13
47 anm101105 C57BL6 C2 5 632.0 L5 143 10749 643643 45010 -10.0 10.0 1.136868e-13
48 anm101105 C57BL6 C2 5 627.0 L5 104 16282 468104 45010 -10.0 10.0 1.136868e-13
49 anm106213 C57BL6 C2 5 637.0 L5 99 3360 445599 45010 -10.0 10.0 2.273737e-13
50 anm106213 C57BL6 C2 5 683.0 L5 118 316 531118 45010 -10.0 10.0 1.136868e-13
51 anm134333 C57BL6 C2 5 815.0 L5 107 6428 535107 50010 -10.0 10.0 2.273737e-13
52 anm134333 C57BL6 C2 5 784.0 L5 302 10574 1510302 50010 -10.0 10.0 4.547474e-13
files 52 | subjects 23 | cells per group {'L4C2': 31, 'L5': 11, 'L4out': 10}
trials per file: mean 115.0, s.d. 60.6, min 19, max 302
cell numbers in the identifiers: True
# Units and ranges of every series used, pooled over files
def rng_of(x):
x = np.asarray(x, float)
return f"{np.nanmin(x):.3g} to {np.nanmax(x):.3g}"
print("theta_at_base (deg):", rng_of(np.concatenate([c.theta for c in cells])))
print("amplitude (described a.u., read as deg):", rng_of(np.concatenate([c.amp for c in cells])))
print("phase (rad):", rng_of(np.concatenate([c.phase for c in cells])))
print("delta_kappa (a.u.):", rng_of(np.concatenate([c.dkappa for c in cells])))
seglen = np.concatenate([np.diff(np.r_[c.seg_start, len(c.vt)]) for c in cells])
print("behaviour samples per trial:", np.unique(seglen, return_counts=True))
print("spike samples per trial block:", meta.spike_block.value_counts().to_dict(), "(10 kHz)")
pp = pd.concat([c.trials.pole_position for c in cells], ignore_index=True)
print("pole_position range per file:", sorted({(float(c.trials.pole_position.min()), float(c.trials.pole_position.max())) for c in cells}))
print("pole_position, described as mm:", sorted(pp.unique())[:6], "... values per file:",
sorted({len(c.trials.pole_position.unique()) for c in cells}))
resp = pd.concat([c.trials.assign(cell=c.cell_no) for c in cells], ignore_index=True)
print("trial outcomes:", resp.response.value_counts().to_dict(), "| types:", resp.type.value_counts().to_dict())
print("stim_present values:", resp.stim_present.value_counts().to_dict())
print("outcome against type (rows: type):\n", pd.crosstab(resp.type, resp.response))
theta_at_base (deg): -113 to 105 amplitude (described a.u., read as deg): 7.05e-09 to 139 phase (rad): -3.14 to 3.14
delta_kappa (a.u.): -0.372 to 0.457
behaviour samples per trial: (array([4500, 4501, 5000, 5001]), array([ 977, 1536, 2955, 510]))
spike samples per trial block: {50000: 30, 45010: 12, 45000: 7, 50010: 3} (10 kHz)
pole_position range per file: [(0.0, 5952.0), (0.1984, 0.8928), (0.1984, 0.992), (0.5456, 1.3888), (0.5456, 1.488), (0.5952, 1.3888), (0.5952, 1.488), (0.6448, 1.3888), (0.7936, 1.5872), (992.0, 11904.0), (1984.0, 13888.0), (3968.0, 8928.0), (4960.0, 9920.0), (4960.0, 10912.0), (4960.0, 11904.0), (5456.0, 11408.0), (5952.0, 10912.0), (5952.0, 11904.0), (49600000.0, 119040000.0), (49600000.0, 128960000.0), (59520000.0, 128960000.0), (59520000.0, 138880000.0), (74400000.0, 148800000.0), (95232000.0, 148800000.0), (95232000.0, 163680000.0), (95232000.0, 168640000.0), (595200000000.0, 1388800000000.0)]
pole_position, described as mm: [np.float64(0.0), np.float64(0.1984), np.float64(0.5456), np.float64(0.5952), np.float64(0.6448), np.float64(0.7936)] ... values per file: [2, 3, 4, 5, 6]
trial outcomes: {'CR': 2232, 'Hit': 2205, 'Miss': 784, 'FA': 757} | types: {'Go': 2989, 'Nogo': 2989}
stim_present values: {0: 5978}
outcome against type (rows: type):
response CR FA Hit Miss
type
Go 0 0 2205 784
Nogo 2232 757 0 0
# Touches: every onset should have an offset in its trial; licks against first_lick_time
TCH = {c.cell_no: hires.touches(c) for c in cells}
n_on = sum(len(T.on) for T in TCH.values())
n_inc = sum((~T.complete).sum() for T in TCH.values())
dur = np.concatenate([T.dur for T in TCH.values()])
print(f"touch onsets {n_on}, without an offset in their trial {n_inc}; touch duration ms: median {np.median(dur):.0f}, range {dur.min()} to {dur.max()}, > 100 ms {np.mean(dur > 100):.3f}")
lick_err = []
for c in cells:
for k, fl in enumerate(c.trials.first_lick_time):
a = c.seg_start[k]
b = c.seg_start[k + 1] if k + 1 < len(c.seg_start) else len(c.vt)
li = np.flatnonzero(c.lick[a:b])
if np.isnan(fl):
lick_err.append(np.nan if len(li) == 0 else -1.0)
else:
lick_err.append((c.seg_t0[k] + li[0] / 1000 - fl) * 1000 if len(li) else np.inf)
lick_err = np.array(lick_err)
fin = lick_err[np.isfinite(lick_err)]
print(f"first_lick_time against the first lick onset of the lick series: {len(fin)} trials, |difference| median {np.median(np.abs(fin)):.2f} ms, max {np.abs(fin).max():.2f} ms;"
f" lick series without first_lick_time {np.sum(lick_err == -1)}, first_lick_time without a lick in the series {np.sum(np.isinf(lick_err))}")
pa_first = []
for c in cells:
for k in range(len(c.trials)):
a = c.seg_start[k]
b = c.seg_start[k + 1] if k + 1 < len(c.seg_start) else len(c.vt)
r = hires.runs(c.pole[a:b])
pa_first.append(np.nan if len(r) == 0 else (c.seg_t0[k] + r[0, 0] / 1000 - c.trials.pole_in_time.iloc[k]))
pa_first = np.array(pa_first)
print(f"pole_available onset minus pole_in_time (s): median {np.nanmedian(pa_first):.3f}, range {np.nanmin(pa_first):.3f} to {np.nanmax(pa_first):.3f}; trials without pole_available {np.isnan(pa_first).sum()}")
touch onsets 41428, without an offset in their trial 0; touch duration ms: median 21, range 1 to 3042, > 100 ms 0.080 first_lick_time against the first lick onset of the lick series: 3473 trials, |difference| median 0.76 ms, max 3142.16 ms; lick series without first_lick_time 47, first_lick_time without a lick in the series 136
pole_available onset minus pole_in_time (s): median 0.370, range 0.117 to 0.709; trials without pole_available 0
# The stored amplitude and phase against the Methods' filter (6-60 Hz, 4th-order Butterworth, Hilbert),
# and the stored theta_filt and set_point against their descriptions, on the first 3 trials of every file
from scipy.signal import butter, filtfilt, hilbert
bb, aa = butter(2, [6 / 500, 60 / 500], "bandpass")
chk = []
for c in cells:
for k in range(3):
a, b = c.seg_start[k], c.seg_start[k + 1]
y = filtfilt(bb, aa, c.theta[a:b].astype(float))
h = hilbert(y)
chk.append(dict(cell=c.cell_no, amp_r=np.corrcoef(np.abs(h), c.amp[a:b])[0, 1],
amp_maxdiff=np.abs(np.abs(h) - c.amp[a:b]).max(),
phase_maxdiff=np.abs(np.angle(np.exp(1j * (np.angle(h) - c.phase[a:b])))).max()))
chk = pd.DataFrame(chk)
print("recomputed Hilbert amplitude vs stored, %d trials: r median %.5f, min %.4f; max |diff| median %.3f deg, max %.2f deg" %
(len(chk), chk.amp_r.median(), chk.amp_r.min(), chk.amp_maxdiff.median(), chk.amp_maxdiff.max()))
# theta_filt, set_point and the voltage, read for the first seconds of every file
SER = hires.check_series_all()
SER["cell"] = [by_path[p] for p in SER.path]
SER = SER.set_index("cell").sort_index()
SER["volt_ptp_mV"] = SER.volt_ptp_stored * SER.volt_conversion * 1000
SER["spike_offset_ms"] = meta.spike_offset_ms
print(SER[["theta_filt_r", "theta_filt_maxdiff", "set_point_r", "volt_ptp_mV", "volt_unit", "spike_offset_ms", "file_name_ogen", "stimulus_groups"]].round(4).to_string())
print(SER.groupby("spike_offset_ms").agg(files=("path", "size"), theta_filt_r_min=("theta_filt_r", "min"), theta_filt_maxdiff=("theta_filt_maxdiff", "max"),
volt_mV_min=("volt_ptp_mV", "min"), volt_mV_max=("volt_ptp_mV", "max")).to_string())
SER["theta_filt_is_bandpass"] = SER.theta_filt_maxdiff < 0.01
SER["volt_in_mV_range"] = SER.volt_ptp_mV > 1
print("files by spike offset (rows) and by whether theta_filt equals amplitude*cos(phase) and the voltage spans more than 1 mV (columns):")
print(pd.crosstab(SER.spike_offset_ms, [SER.theta_filt_is_bandpass, SER.volt_in_mV_range]))
print("files whose set_point holds missing values (r undefined):", int(SER.set_point_r.isna().sum()), "| r where defined: %.3f to %.3f" % (SER.set_point_r.min(), SER.set_point_r.max()))
print("files named 'ogen': %d; with any stimulus series: %d; trials with stim_present = 1: %d" %
(SER.file_name_ogen.sum(), (SER.stimulus_groups > 0).sum(), int((resp.stim_present == 1).sum())))
recomputed Hilbert amplitude vs stored, 156 trials: r median 0.99996, min 0.9839; max |diff| median 0.576 deg, max 4.82 deg
theta_filt_r theta_filt_maxdiff set_point_r volt_ptp_mV volt_unit spike_offset_ms file_name_ogen stimulus_groups
cell
1 1.0000 0.0000 0.9842 9.9322 volts 0.0 True 0
2 1.0000 0.0000 0.9866 12.6230 volts 0.0 True 0
3 1.0000 0.0000 0.9930 12.1394 volts 0.0 True 0
4 1.0000 0.0000 0.9982 3.1777 volts 0.0 True 0
5 1.0000 0.0000 0.9968 14.3952 volts 0.0 True 0
6 1.0000 0.0000 0.9822 9.2798 volts 0.0 True 0
7 1.0000 0.0000 0.9832 6.0697 volts 0.0 True 0
8 1.0000 0.0000 0.9872 7.9750 volts 0.0 True 0
9 0.9916 1.0884 NaN 0.0023 volts -10.0 True 0
10 0.9833 33.7196 0.9966 0.0113 volts -10.0 True 0
11 0.9873 26.5349 NaN 0.0048 volts -10.0 True 0
12 0.9873 11.1904 0.9952 0.0076 volts -10.0 True 0
13 1.0000 0.0000 0.9913 28.3620 volts 0.0 True 0
14 0.9840 21.1456 0.9923 0.0046 volts -10.0 False 0
15 1.0000 0.0000 0.9574 7.8971 volts 0.0 True 0
16 0.9892 18.5044 0.9955 0.0052 volts -10.0 False 0
17 0.9861 22.3218 NaN 0.0090 volts -10.0 False 0
18 0.9894 8.8093 0.9958 0.0059 volts -10.0 True 0
19 0.9754 5.8178 NaN 0.0086 volts -10.0 True 0
20 0.8935 33.7075 0.9984 0.0072 volts 0.0 True 0
21 0.9893 4.3747 NaN 0.0046 volts -10.0 True 0
22 1.0000 0.0000 0.9756 6.7156 volts 0.0 True 0
23 0.9768 31.8798 NaN 0.0039 volts -10.0 False 0
24 1.0000 0.0000 0.9720 21.9969 volts 0.0 True 0
25 1.0000 0.0000 0.9980 9.1922 volts 0.0 True 0
26 0.9916 11.2635 0.9968 0.0079 volts -10.0 False 0
27 0.9866 18.7602 0.9931 0.0123 volts -10.0 False 0
28 0.9904 9.6222 0.9990 0.0047 volts -10.0 True 0
29 1.0000 0.0000 0.9891 7.9474 volts 0.0 True 0
30 0.9958 2.7767 0.9957 0.0041 volts -10.0 False 0
31 0.9796 15.7750 NaN 0.0123 volts -10.0 True 0
32 0.9867 16.8916 NaN 0.0042 volts -10.0 False 0
33 0.9941 5.3349 NaN 0.0028 volts -10.0 False 0
34 0.9920 5.3390 NaN 0.0096 volts -10.0 False 0
35 0.9932 5.7710 NaN 0.0073 volts -10.0 False 0
36 0.9944 3.9676 0.9974 0.0085 volts -10.0 True 0
37 0.9557 46.1250 NaN 0.0028 volts -10.0 True 0
38 0.9903 15.4984 0.9973 0.0108 volts -10.0 True 0
39 0.9822 21.6269 0.9956 0.0260 volts -10.0 True 0
40 0.9889 10.4071 NaN 0.0020 volts 0.0 True 0
41 0.9679 23.6640 0.9924 0.0039 volts -10.0 True 0
42 0.9924 3.7868 0.9926 0.0056 volts -10.0 False 0
43 0.9935 7.2988 NaN 0.0089 volts -10.0 False 0
44 0.9948 6.1587 NaN 0.0037 volts -10.0 False 0
45 0.9710 23.5843 NaN 0.0056 volts -10.0 False 0
46 0.9935 3.0854 NaN 0.0119 volts -10.0 False 0
47 0.9788 17.5848 NaN 0.0035 volts -10.0 False 0
48 0.9755 13.8123 NaN 0.0036 volts -10.0 False 0
49 0.9903 11.5763 NaN 0.0164 volts -10.0 False 0
50 0.9883 9.1587 NaN 0.0059 volts -10.0 False 0
51 0.9749 15.5144 NaN 0.0039 volts -10.0 False 0
52 0.9883 9.1633 NaN 0.0019 volts -10.0 False 0
files theta_filt_r_min theta_filt_maxdiff volt_mV_min volt_mV_max
spike_offset_ms
-10.0 36 0.955731 46.124981 0.001876 0.026002
0.0 16 0.893471 33.707486 0.002048 28.362004
files by spike offset (rows) and by whether theta_filt equals amplitude*cos(phase) and the voltage spans more than 1 mV (columns):
theta_filt_is_bandpass False True
volt_in_mV_range False True
spike_offset_ms
-10.0 36 0
0.0 2 14
files whose set_point holds missing values (r undefined): 23 | r where defined: 0.957 to 0.999
files named 'ogen': 30; with any stimulus series: 0; trials with stim_present = 1: 0
Clocks. The trials table, the behaviour series and the spike train each carry timestamps.
Behaviour segments start within 5 ms of each trial's start_time. In 36 files the
SpikeTrain timestamps start 10 ms before the voltage (CurrentClampSeries) timestamps, though
the spike samples line up with the voltage samples one for one; in 16 files the two agree. The
cell below asks which clock puts the spikes on the behaviour's time: the touch response is the
contrast. If the SpikeTrain timestamps are right, the touch-response peak should fall at the
same latency in both kinds of file; if the voltage timestamps are right, the peaks in the 36
files should sit 10 ms earlier than the others under the SpikeTrain clock.
rows = []
for c in cells:
T = TCH[c.cell_no]
for clock in hires.CLOCKS:
s = c.spikes_vt(clock)
t, h, n = hires.pcth(c, s, T.on)
post = t >= 0
rows.append(dict(cell=c.cell_no, group=c.group, offset_ms=1000 * c.extra["spike_offset"], clock=clock,
peak_ms=int(t[post][np.argmax(h[post])]), peak_rate=h[post].max() * 1000,
pre_rate=h[~post].mean() * 1000))
clk = pd.DataFrame(rows)
l4 = clk[clk.group == "L4C2"]
print(l4.groupby(["clock", "offset_ms"]).peak_ms.describe()[["count", "mean", "50%", "min", "max"]])
for clock in hires.CLOCKS:
x = l4[(l4.clock == clock) & (l4.offset_ms == 0)].peak_ms
y = l4[(l4.clock == clock) & (l4.offset_ms == -10)].peak_ms
print(f"{clock} clock: L4 C2 touch peak, files without offset median {x.median():.0f} ms (n={len(x)}), files with the 10 ms offset median {y.median():.0f} ms (n={len(y)}); rank sum p = {stats.mannwhitneyu(x, y).pvalue:.3g}")
print("peak over pre-touch rate, L4 C2 cells, spike clock: median %.1f-fold, min %.1f" %
((l4[l4.clock == "spike"].peak_rate / l4[l4.clock == "spike"].pre_rate.clip(lower=1e-9)).median(),
(l4[l4.clock == "spike"].peak_rate / l4[l4.clock == "spike"].pre_rate.clip(lower=1e-9)).min()))
# Shifted-label control: the same contrast with touch onsets moved by a random amount (200 to 1000 ms) within the trial
rng_clk = np.random.default_rng(11)
ratio_true, ratio_shift = [], []
for c in cells:
if c.group != "L4C2":
continue
T = TCH[c.cell_no]
s = c.spikes_vt("spike")
_, h, _ = hires.pcth(c, s, T.on)
ratio_true.append(h[50:80].mean() / max(h[:50].mean(), 1e-9))
shifted = T.on + rng_clk.integers(200, 1000, len(T.on)) * rng_clk.choice([-1, 1], len(T.on))
lo = c.seg_start[c.trial_of[T.on]]
hi = np.r_[c.seg_start[1:], len(c.vt)][c.trial_of[T.on]] - 1
shifted = np.clip(shifted, lo, hi)
_, h2, _ = hires.pcth(c, s, shifted)
ratio_shift.append(h2[50:80].mean() / max(h2[:50].mean(), 1e-9))
print("L4 C2, rate 0-30 ms after over 50 ms before touch: true onsets median %.2f, shifted onsets median %.2f" %
(np.median(ratio_true), np.median(ratio_shift)))
print("L4 C2 touch-peak latency per file under the spike clock, by the file's spike-clock offset (ms):")
for off, x in l4[l4.clock == "spike"].groupby("offset_ms").peak_ms:
print(f" offset {off:+.0f} ms, {len(x)} files:", " ".join(str(v) for v in sorted(x)))
CLOCK = "spike"
count mean 50% min max
clock offset_ms
spike -10.0 16.0 13.000000 11.5 7.0 23.0
0.0 15.0 10.466667 8.0 7.0 25.0
voltage -10.0 16.0 23.000000 21.5 17.0 33.0
0.0 15.0 10.466667 8.0 7.0 25.0
spike clock: L4 C2 touch peak, files without offset median 8 ms (n=15), files with the 10 ms offset median 12 ms (n=16); rank sum p = 0.0157
voltage clock: L4 C2 touch peak, files without offset median 8 ms (n=15), files with the 10 ms offset median 22 ms (n=16); rank sum p = 2.14e-05
peak over pre-touch rate, L4 C2 cells, spike clock: median 15.3-fold, min 2.9
L4 C2, rate 0-30 ms after over 50 ms before touch: true onsets median 3.72, shifted onsets median 0.99 L4 C2 touch-peak latency per file under the spike clock, by the file's spike-clock offset (ms): offset -10 ms, 16 files: 7 10 10 10 11 11 11 11 12 13 14 14 17 17 17 23 offset +0 ms, 15 files: 7 7 7 8 8 8 8 8 10 10 10 11 13 17 25
Recordings deposited twice. Pairs of files from the same subject are compared trial by trial on the whisker angle and the spike train.
dups = []
sig = {}
for c in cells:
th = [c.theta[a:b][:200].tobytes() for a, b in zip(c.seg_start, np.r_[c.seg_start[1:], len(c.vt)])]
sig[c.cell_no] = (set(th), np.round(c.spike_t, 4))
for i, a in enumerate(cells):
for b in cells[i + 1:]:
if a.subject != b.subject:
continue
same_theta = len(sig[a.cell_no][0] & sig[b.cell_no][0])
same_spk = len(np.intersect1d(sig[a.cell_no][1], sig[b.cell_no][1]))
dups.append(dict(a=a.cell_no, b=b.cell_no, subject=a.subject, trials_a=len(a.trials), trials_b=len(b.trials),
identical_theta_trials=same_theta, shared_spike_times=same_spk,
spikes_a=len(a.spike_t), spikes_b=len(b.spike_t)))
dups = pd.DataFrame(dups)
print(dups.to_string(index=False))
print("pairs with any identical theta trial:", int((dups.identical_theta_trials > 0).sum()), "of", len(dups))
a b subject trials_a trials_b identical_theta_trials shared_spike_times spikes_a spikes_b 1 3 anm244028 48 64 0 0 210 829 1 29 anm244028 48 50 0 0 210 128 2 5 anm266951 19 34 0 0 283 102 3 29 anm244028 64 50 0 0 829 128 4 24 anm244025 55 87 0 0 168 141 6 7 anm244024 48 26 0 0 95 35 6 15 anm244024 48 133 0 0 95 781 6 22 anm244024 48 24 0 0 95 107 7 15 anm244024 26 133 0 0 35 781 7 22 anm244024 26 24 0 0 35 107 8 13 anm266945 30 63 0 0 202 1305 8 25 anm266945 30 103 0 0 202 725 9 10 anm215592 122 241 0 1 9796 2009 9 28 anm215592 122 88 0 0 9796 187 9 37 anm215592 122 205 0 0 9796 719 10 28 anm215592 241 88 0 0 2009 187 10 37 anm215592 241 205 0 0 2009 719 12 20 anm226694 145 158 0 0 242 478 12 31 anm226694 145 96 0 0 242 4371 13 25 anm266945 63 103 0 0 1305 725 14 23 jf35892 146 121 0 0 1097 664 14 26 jf35892 146 198 0 0 1097 368 14 30 jf35892 146 98 0 1 1097 330 15 22 anm244024 133 24 0 0 781 107 18 19 anm226695 63 98 0 0 2505 130 18 21 anm226695 63 114 0 1 2505 1486 19 21 anm226695 98 114 0 0 130 1486 20 31 anm226694 158 96 0 0 478 4371 23 26 jf35892 121 198 0 0 664 368 23 30 jf35892 121 98 0 0 664 330 26 30 jf35892 198 98 0 0 368 330 27 42 anm17702 224 96 0 2 901 22298 28 37 anm215592 88 205 0 0 187 719 32 33 anm106213 192 187 0 0 158 111 32 49 anm106213 192 99 0 0 158 3360 32 50 anm106213 192 118 0 0 158 316 33 49 anm106213 187 99 0 0 111 3360 33 50 anm106213 187 118 0 0 111 316 35 51 anm134333 101 107 0 1 195 6428 35 52 anm134333 101 302 0 0 195 10574 43 44 jf42400 119 92 0 3 3515 10921 43 45 jf42400 119 93 0 0 3515 656 44 45 jf42400 92 93 0 1 10921 656 46 47 anm101105 117 143 0 7 7341 10749 46 48 anm101105 117 104 0 23 7341 16282 47 48 anm101105 143 104 0 17 10749 16282 49 50 anm106213 99 118 0 0 3360 316 51 52 anm134333 107 302 0 1 6428 10574 pairs with any identical theta trial: 0 of 48
Measures shared by the claims¶
Computed once per cell, on the SpikeTrain clock, with every reading the claims file lists.
- Touches: onset and offset series paired within each trial. Exploration time: whisking (amplitude > 2.5 deg) or touch, leaving out touches > 100 ms, while the pole is in reach.
- Onset latency: first 1 ms bin of the PCTH 0 to 50 ms after onset above the pre-touch mean
(-50 to 0 ms) plus two s.d.; from first touches (
lat_first) and from all touches (lat_all). - Phase tuning: spikes in the phase-analysis cycles of the whole trial (cycles with amplitude
2.5 deg throughout, monotonic phase, no lick, not within 70 ms after a touch), 12-bin rate histogram, cosine fit, depth A / (A + B). The same cycles counted only in exploration time (
depth_mod_x) serve example-cell-phase-modulation and a variant elsewhere. - Rates: whisking (amplitude > 2.5 deg, leaving out the whisk cycles that hold a lick onset and
the 70 ms after each touch, read with and without touch epochs) and non-whisking (>= 100 ms
with amplitude < 1.25 deg, no touch, no lick), over all tracked time. The whisking rate without
the licking and after-touch exclusions, the earlier run's, is kept as a variant
(
_noexcl). The files give lick onsets only, so a cycle counts as "during licking" when a lick onset falls in it, a standard choice. - Minimum ISI: the shortest interval between two spikes of one trial.
t0 = time.time()
D = {}
base_rows = []
for c in cells:
T = TCH[c.cell_no]
s = c.spikes_vt(CLOCK)
E = hires.exploration_mask(c, T)
touch = hires.touch_mask(c, T)
cyc = hires.cycle_table(c, T)
first = T.on[T.order == 1]
lat_first = hires.onset_latency(c, s, first)
lat_all = hires.onset_latency(c, s, T.on)
phase_cyc = cyc[cyc.phase_ok]
cen, rate, cnt, occ = hires.phase_tuning(c, s, phase_cyc, ~touch)
A, B, phi, depth = hires.cos_fit(cen, np.nan_to_num(rate))
if cnt.sum() == 0:
phi, depth = np.nan, np.nan
# the same cycles, only their samples in exploration time (pole in reach): example-cell-phase-modulation's
# "exploration periods", and a variant elsewhere
_, rate_x, cnt_x, occ_x = hires.phase_tuning(c, s, phase_cyc, E & ~touch)
depth_x = hires.cos_fit(cen, np.nan_to_num(rate_x))[3] if cnt_x.sum() else np.nan
nw = hires.nonwhisk_mask(c, T)
D[c.cell_no] = dict(T=T, s=s, E=E, touch=touch, cyc=cyc, phi=phi, phase_rate=rate, phase_cnt=cnt, phase_occ=occ,
phase_rate_x=rate_x, phase_cnt_x=cnt_x)
base_rows.append(dict(
cell=c.cell_no, group=c.group, subject=c.subject, n_touch=len(T.on), n_first=len(first),
lat_first=lat_first, lat_all=lat_all, depth_mod=depth, phi_pref=phi, phase_spikes=int(cnt.sum()),
phase_time_s=float(occ.sum()), depth_mod_x=depth_x, phase_spikes_x=int(cnt_x.sum()), phase_time_x_s=float(occ_x.sum()),
rate_whisk_no_touch=hires.rate_in(c, s, hires.whisk_mask(c, T, False)),
rate_whisk_with_touch=hires.rate_in(c, s, hires.whisk_mask(c, T, True)),
rate_whisk_no_touch_noexcl=hires.rate_in(c, s, hires.whisk_mask(c, T, False, exclude=False)),
rate_whisk_with_touch_noexcl=hires.rate_in(c, s, hires.whisk_mask(c, T, True, exclude=False)),
rate_nonwhisk=hires.rate_in(c, s, nw), nonwhisk_s=nw.sum() / 1000,
min_isi_ms=hires.isi_min(c), explore_s=E.sum() / 1000, n_spikes=len(s)))
base = pd.DataFrame(base_rows).set_index("cell")
print(f"{time.time() - t0:.0f} s")
print(base.round(3).to_string())
print("phase-analysis spikes per cell, median by group:", base.groupby("group").phase_spikes.median().to_dict(),
"| restricted to exploration time (variant):", base.groupby("group").phase_spikes_x.median().to_dict())
print("cells with a depth of exactly 1 (> 0.999), of those with a fit:", {g: f"{int((x.dropna() > 0.999).sum())} of {x.notna().sum()}" for g, x in base.groupby("group").depth_mod})
print("cells with fewer than 10 spikes in the phase-analysis cycles:", base[base.phase_spikes < 10].groupby("group").size().to_dict(),
"| with none (no depth):", base[base.phase_spikes == 0].index.tolist())
13 s
group subject n_touch n_first lat_first lat_all depth_mod phi_pref phase_spikes phase_time_s depth_mod_x phase_spikes_x phase_time_x_s rate_whisk_no_touch rate_whisk_with_touch rate_whisk_no_touch_noexcl rate_whisk_with_touch_noexcl rate_nonwhisk nonwhisk_s min_isi_ms explore_s n_spikes
cell
1 L4C2 anm244028 251 24 3.0 4.0 1.000 2.372 8 15.490 1.000 4 10.188 0.339 2.270 0.778 2.284 0.076 119.134 1.1 41.380 210
2 L4C2 anm266951 118 12 5.0 6.0 0.859 -0.951 24 7.824 0.929 10 3.328 3.123 4.836 3.219 4.520 1.379 30.450 3.3 23.681 283
3 L4C2 anm244028 193 28 2.0 4.0 0.631 -1.503 23 7.201 1.000 5 2.717 2.656 12.626 3.444 11.708 0.546 227.147 1.2 24.612 829
4 L4C2 anm244025 332 26 5.0 6.0 NaN NaN 0 5.943 NaN 0 3.375 0.033 1.490 0.148 0.998 0.385 129.975 2.6 42.313 168
5 L4C2 anm266951 148 17 6.0 6.0 NaN NaN 0 14.077 NaN 0 1.567 0.024 1.601 0.133 1.384 0.171 58.556 2.0 26.067 102
6 L4C2 anm244024 142 27 8.0 2.0 1.000 -2.880 1 11.697 NaN 0 3.715 0.349 0.646 0.527 0.758 0.275 141.728 2.6 20.371 95
7 L4C2 anm244024 88 13 4.0 4.0 1.000 1.324 5 10.994 1.000 5 6.874 0.297 0.701 0.262 0.547 0.022 46.130 2.2 25.693 35
8 L4C2 anm266945 186 17 4.0 5.0 NaN NaN 0 6.646 NaN 0 2.165 0.376 2.307 0.886 2.293 1.130 58.385 2.3 21.656 202
9 L4C2 anm215592 679 75 1.0 9.0 0.123 -1.663 694 75.947 0.112 138 13.405 10.625 17.120 15.096 19.789 13.677 280.617 0.9 94.421 9775
10 L4C2 anm215592 1369 179 7.0 9.0 0.824 2.901 14 125.375 0.470 4 26.732 0.157 2.345 2.067 3.533 0.719 535.837 2.3 191.142 2005
11 L4C2 anm234232 684 56 6.0 8.0 0.480 0.385 180 70.844 0.602 28 14.665 2.294 3.193 2.358 3.091 2.144 118.479 11.0 89.503 1051
12 L4C2 anm226694 1104 93 7.0 7.0 0.247 -2.873 14 65.104 1.000 1 13.652 0.241 0.583 0.437 0.671 0.149 321.330 2.3 125.738 242
13 L4C2 anm266945 382 39 2.0 8.0 0.444 0.503 83 12.111 0.389 39 6.772 4.102 11.878 4.767 10.536 0.337 83.107 2.0 52.236 1305
14 L4C2 jf35892 1516 115 8.0 9.0 1.000 1.813 17 65.696 1.000 6 17.445 0.201 3.778 0.781 3.468 0.148 189.037 2.9 158.175 1095
15 L4C2 anm244024 356 51 8.0 9.0 1.000 -1.887 17 36.897 1.000 3 10.510 0.617 0.949 0.823 1.083 1.333 359.419 2.3 62.788 781
16 L4C2 jf37166 728 81 6.0 8.0 0.179 -0.908 407 57.092 0.084 91 14.721 6.121 6.507 5.353 5.734 1.065 112.708 4.9 106.443 1575
17 L4C2 anm106211 1429 156 10.0 10.0 0.401 -1.014 185 107.638 0.370 79 58.801 1.596 2.800 3.896 4.312 0.611 256.994 2.4 227.721 1889
18 L4C2 anm226695 314 40 7.0 6.0 0.314 -2.948 64 18.867 0.239 12 3.903 4.723 8.016 6.579 8.732 8.706 148.296 1.2 48.887 2497
19 L4C2 anm226695 514 72 9.0 9.0 1.000 -2.614 2 38.782 1.000 1 4.462 0.037 0.523 0.228 0.607 0.075 227.036 2.1 60.353 130
20 L4C2 anm226694 1299 119 16.0 18.0 0.265 -2.081 8 66.115 0.435 3 24.742 0.175 0.896 0.561 1.036 0.276 337.530 3.0 171.755 478
21 L4C2 anm226695 608 79 5.0 11.0 0.365 -0.985 19 36.100 0.429 8 5.152 0.639 5.067 3.718 6.559 0.632 272.194 1.5 72.541 1486
22 L4C2 anm244024 98 14 4.0 5.0 1.000 0.785 1 7.252 NaN 0 3.304 0.161 2.954 0.291 2.090 0.043 46.825 1.7 22.581 107
23 L4C2 jf35892 1023 97 8.0 9.0 1.000 2.646 92 64.808 1.000 14 17.191 0.777 2.195 1.257 2.373 0.227 131.909 2.7 121.641 664
24 L4C2 anm244025 503 42 6.0 7.0 0.899 1.382 3 13.477 1.000 1 2.973 0.388 0.569 0.757 0.783 0.123 211.690 2.2 42.918 141
25 L4C2 anm266945 533 46 8.0 8.0 1.000 1.648 5 10.857 1.000 1 3.524 0.514 3.470 1.897 3.697 0.192 218.868 2.3 63.824 725
26 L4C2 jf35892 1316 121 9.0 11.0 NaN NaN 0 69.206 NaN 0 16.547 0.016 0.820 0.359 0.940 0.164 371.426 2.9 156.712 367
27 L4C2 anm17702 2207 157 7.0 9.0 1.000 -1.521 3 102.400 1.000 1 45.096 0.089 1.372 0.660 1.574 0.487 314.195 2.2 292.995 898
28 L4C2 anm215592 521 57 9.0 9.0 1.000 1.907 53 53.337 1.000 17 9.300 0.638 0.688 0.860 0.876 0.103 184.223 2.5 71.129 187
29 L4C2 anm244028 201 21 13.0 13.0 1.000 2.356 1 11.137 1.000 1 3.528 0.138 0.897 0.554 1.109 0.315 155.400 8.1 25.184 128
30 L4C2 jf35892 482 81 18.0 18.0 0.329 0.313 10 35.690 0.377 7 16.406 0.247 0.400 1.935 1.900 0.050 199.974 2.9 102.513 330
31 L4C2 anm226694 1072 78 15.0 14.0 0.132 2.738 112 43.912 0.506 53 11.257 3.642 11.136 9.841 14.300 4.838 184.999 2.2 113.809 4368
32 L4out anm106213 1471 149 NaN 23.0 1.000 0.394 84 88.302 1.000 28 44.913 0.641 0.470 0.545 0.422 0.025 285.236 2.3 187.838 156
33 L4out anm106213 642 83 16.0 16.0 1.000 1.740 14 77.487 1.000 5 32.922 0.187 0.218 0.375 0.371 0.040 452.041 2.5 131.509 111
34 L4out anm131970 555 59 NaN 16.0 0.077 -3.042 544 44.209 0.140 341 27.841 12.138 11.435 12.300 11.693 8.922 203.550 1.5 102.814 4941
35 L4out anm134333 847 77 20.0 3.0 1.000 -2.738 10 57.124 1.000 9 43.721 0.268 0.286 0.285 0.297 0.533 150.009 1.6 186.346 195
36 L4out anm211295 733 83 5.0 0.0 0.215 -1.051 82 61.502 0.514 16 11.131 1.465 1.464 2.019 1.956 0.853 506.583 3.4 79.426 997
37 L4out anm215592 948 105 3.0 0.0 1.000 -0.995 220 107.246 1.000 42 21.154 1.827 1.777 1.658 1.635 0.266 566.860 2.5 125.030 716
38 L4out anm217330 800 104 4.0 0.0 0.503 -1.055 75 98.987 0.496 18 21.148 0.859 0.915 0.851 0.894 1.413 276.006 3.5 140.484 792
39 L4out anm234230 774 72 8.0 7.0 0.593 -1.248 162 47.219 0.961 37 7.967 4.067 4.033 3.758 3.772 5.584 202.721 4.4 91.706 2224
40 L4out anm234231 385 46 1.0 33.0 0.488 -1.448 50 35.317 1.000 9 6.363 1.462 1.625 1.837 1.927 1.334 304.377 2.7 47.022 769
41 L4out anm226692 1381 161 NaN 2.0 0.162 0.578 164 111.225 0.180 29 20.014 2.095 2.024 2.476 2.369 2.859 598.023 6.7 195.354 3288
42 L5 anm17702 1232 72 3.0 12.0 0.049 1.053 1913 29.068 0.038 667 10.284 59.310 69.098 69.686 74.266 29.846 137.540 0.1 98.970 22218
43 L5 jf42400 741 90 11.0 15.0 0.385 3.065 1136 86.831 0.389 398 39.611 11.215 10.901 10.507 10.304 1.869 156.779 0.3 121.894 3491
44 L5 jf42400 1468 64 10.0 10.0 0.026 1.625 1756 62.850 0.044 648 22.167 26.976 28.320 26.972 28.039 23.101 59.738 2.4 137.577 10879
45 L5 jf42400 1120 71 4.0 9.0 0.126 -0.067 44 41.200 0.463 22 20.079 1.204 1.540 1.329 1.561 1.029 84.526 8.3 122.757 655
46 L5 anm101105 1079 80 0.0 9.0 0.106 -2.574 944 61.011 0.018 347 22.692 15.067 15.401 14.323 14.763 12.971 157.654 0.1 115.514 7304
47 L5 anm101105 1167 94 4.0 24.0 0.099 -0.506 923 56.374 0.141 299 16.812 17.009 16.332 17.743 16.984 15.858 276.950 1.5 103.146 10705
48 L5 anm101105 1096 62 4.0 7.0 0.098 -0.789 1106 44.464 0.102 476 13.905 24.623 45.087 34.934 50.101 21.836 139.955 1.6 93.632 16246
49 L5 anm106213 760 47 3.0 NaN 0.166 1.476 188 18.570 0.302 64 6.592 9.497 11.522 9.549 11.326 5.689 267.906 9.1 54.191 3345
50 L5 anm106213 800 66 0.0 0.0 0.615 -0.431 40 53.713 0.529 21 26.612 0.978 0.741 0.862 0.690 0.240 242.014 1.3 110.472 313
51 L5 anm134333 915 51 NaN NaN 0.262 -1.947 530 40.331 0.251 236 19.969 12.819 11.357 10.962 10.171 12.903 234.669 2.9 130.465 6401
52 L5 anm134333 2118 206 10.0 12.0 0.065 -0.420 511 114.759 0.052 399 85.559 4.565 6.326 6.532 7.499 7.082 605.337 1.3 448.907 10539
phase-analysis spikes per cell, median by group: {'L4C2': 14.0, 'L4out': 83.0, 'L5': 923.0} | restricted to exploration time (variant): {'L4C2': 4.0, 'L4out': 23.0, 'L5': 347.0}
cells with a depth of exactly 1 (> 0.999), of those with a fit: {'L4C2': '12 of 27', 'L4out': '4 of 10', 'L5': '0 of 11'}
cells with fewer than 10 spikes in the phase-analysis cycles: {'L4C2': 14} | with none (no depth): [4, 5, 8, 26]
Touch spike windows¶
The claims file (fig3c-touch-spikes-l4c2) leaves four details open and lists their readings:
the shuffle (uniform over exploration time, or circular shifts within each trial), the
closest touch onset (before the spike, or on either side), the 1 ms step, and the longest
window over which the curve is built and the degree-11 polynomial fitted, 100 or 200 ms
(tmax, readings r100 and r200). Where the search for the cutoff starts (start) was asked,
and the answer is now in the criterion: cross, the first increment below chance after the
increment has been above it, the transition where "the curve levels off" (Results, paragraph on
Figure 3). The literal search from t = 1 ms (first) is computed as well and shown only as a
variant; the fitted polynomial dips at the start of a steep curve, so it often stops at 1 ms.
Outcomes rest on cross.
Cells without a first-touch latency take reading (a), latency 0 ms, or (b), the all-touch latency where there is one and else 0 ms (claim fig3c-touch-spikes-outside). Cells with a first-touch latency use it under both. Each cell and reading has its own random generator.
SEED_WIN = 3000
WIN = [(sh, nr, tm) for sh in ("uniform", "circular") for nr in ("before", "either") for tm in hires.TMAX_READINGS]
t0 = time.time()
win_rows = []
for c in cells:
d = D[c.cell_no]
lf, la = base.loc[c.cell_no, "lat_first"], base.loc[c.cell_no, "lat_all"]
if pd.notna(lf):
lats = {("a", "b"): float(lf)}
else:
lats = {("a",): 0.0, ("b",): float(la) if pd.notna(la) else 0.0}
for li, (lr, L) in enumerate(lats.items()):
for k, (sh, nr, tm) in enumerate(WIN):
rng = np.random.default_rng([SEED_WIN, c.cell_no, k, li])
res = hires.touch_window(c, d["s"], d["T"], d["E"], L, sh, nr, tm, rng)
for st, r in res.items():
for lr1 in lr:
win_rows.append(dict(cell=c.cell_no, group=c.group, lat_reading=lr1, shuffle=sh, nearest=nr,
tmax=tm, start=st, latency=L, t_on=r["t_on"], frac_spikes=r["frac_spikes"],
frac_time=r["frac_time"], crossed=r["crossed"], n_explore_spikes=r["n_spikes"]))
W = pd.DataFrame(win_rows)
WKEY = ["lat_reading", "shuffle", "nearest", "tmax", "start"]
print(f"touch windows: {len(W)} rows in {time.time() - t0:.0f} s")
print("cells whose increment never rose above chance under `cross` (the literal cutoff is used):")
print(W[(W.start == "cross") & ~W.crossed].groupby(["group"] + WKEY[1:4]).cell.apply(list).to_string())
touch windows: 1664 rows in 40 s
cells whose increment never rose above chance under `cross` (the literal cutoff is used):
group shuffle nearest tmax
L4out circular before 100 [32]
200 [32, 32]
either 100 [32]
200 [32]
uniform before 100 [32, 38, 38]
200 [32, 32]
either 100 [32, 39, 39]
200 [39, 39]
L5 circular before 100 [43, 43, 51, 51]
200 [43, 43, 49, 49]
either 100 [43, 43, 51, 51]
200 [43, 43, 47, 47]
uniform before 100 [43, 43, 51, 51]
200 [43, 43, 50, 50]
either 100 [43, 43, 50, 50, 51, 51]
200 [43, 43, 47, 47, 50, 50, 51, 51]
# Spikes per touch in each cell's window, first touches of trials and later touches apart
t0 = time.time()
spt_rows = []
for (no, lr, sh, nr, tm, st), g in W.set_index(["cell"] + WKEY).groupby(level=[0, 1, 2, 3, 4, 5]):
r = g.iloc[0]
c = by_no[no]
d = D[no]
T = d["T"]
n = hires.spikes_per_touch(c, d["s"], T, r.latency, r.t_on)
spt_rows.append(dict(cell=no, lat_reading=lr, shuffle=sh, nearest=nr, tmax=tm, start=st,
spt_first=n[T.order == 1].mean() if (T.order == 1).any() else np.nan,
spt_later=n[T.order > 1].mean() if (T.order > 1).any() else np.nan))
SPT = pd.DataFrame(spt_rows)
W = W.merge(SPT, on=["cell"] + WKEY)
print(f"{time.time() - t0:.0f} s")
1 s
Touch and phase windows¶
After the touch spikes are taken out, the remaining exploration spikes are indexed by the time
to the closest preferred-phase time of a whisking cycle (cycles with amplitude > 2.5 deg
throughout, overlapping exploration time, monotonic or not), and a cutoff t_w is found by the
same rule, under the same readings as the touch window (claim fig3d-touch-phase-l4c2). Its open
detail (4), the phase window read as +-t_w (pm) or t_w wide in total (total), doubles the
readings. The shuffle runs over exploration time outside the touch windows.
SEED_PH = 4000
t0 = time.time()
ph_rows = []
for c in cells:
d = D[c.cell_no]
cyc = d["cyc"]
Ec = np.r_[0, np.cumsum(d["E"])]
wc = cyc[cyc.whisking & ((Ec[cyc.stop.to_numpy()] - Ec[cyc.start.to_numpy()]) > 0)]
sub = W[W.cell == c.cell_no]
for li, (L, g0) in enumerate(sub.groupby("latency")):
lrs = sorted(g0.lat_reading.unique())
g0 = g0[g0.lat_reading == lrs[0]]
for k, (sh, nr, tm) in enumerate(WIN):
g = g0[(g0.shuffle == sh) & (g0.nearest == nr) & (g0.tmax == tm)].set_index("start")
t_on = {st: int(g.loc[st, "t_on"]) for st in ("first", "cross")}
for wi, width in enumerate(("pm", "total")):
rng = np.random.default_rng([SEED_PH, c.cell_no, k, li, wi])
res = hires.phase_window(c, d["s"], d["T"], d["E"], L, t_on, nr, d["phi"], wc, sh, tm, width, rng)
for st, r in res.items():
for lr in lrs:
ph_rows.append(dict(cell=c.cell_no, group=c.group, lat_reading=lr, shuffle=sh, nearest=nr, tmax=tm,
start=st, width=width, t_w=r["t_w"], frac_tp=r["frac_spikes"], time_tp=r["frac_time"]))
PH = pd.DataFrame(ph_rows).merge(W[["cell"] + WKEY + ["frac_spikes", "frac_time"]], on=["cell"] + WKEY)
PH["frac_phase_only"] = PH.frac_tp - PH.frac_spikes
PKEY = WKEY + ["width"]
print(f"phase windows: {len(PH)} rows in {time.time() - t0:.0f} s")
phase windows: 3328 rows in 107 s
def within(v, exp, tol):
return bool(np.isfinite(v) and abs(v - exp) <= tol + 1e-9)
def msm(x):
x = np.asarray(x, float)
x = x[np.isfinite(x)]
return dict(n=len(x), mean=x.mean() if len(x) else np.nan, sd=x.std(ddof=1) if len(x) > 1 else np.nan,
median=np.median(x) if len(x) else np.nan)
def ranksum(a, b):
a, b = np.asarray(a, float), np.asarray(b, float)
a, b = a[np.isfinite(a)], b[np.isfinite(b)]
return stats.ranksums(a, b).pvalue
G = {g: base.index[base.group == g].tolist() for g in GROUP_ORDER}
recordings-counts¶
The study holds 52 single-unit recordings, 31 L4 cells in the C2 barrel, 10 L4 cells outside C2 and 11 L5 cells near C2, with 115 +- 60 (mean +- s.d.) trials per neuron.
Groups come from the files' electrode location (brain_subregion C2 or surround,
cortical_layer 4 or 5). The files mark no subset of trials as analysed, so readings (a), all
trials, and (b), the performing trials analysed, are the same here.
print(meta.groupby(["region", "layer"]).size())
print("L4 C2 %d, L4 outside C2 %d, L5 %d" % (len(G["L4C2"]), len(G["L4out"]), len(G["L5"])))
print("trials per neuron, all trials (a) = (b): mean %.1f, s.d. %.1f, n = %d" % (meta.trials.mean(), meta.trials.std(), len(meta)))
print("depth by group (um):", meta.groupby("group").depth_um.agg(["min", "max"]).to_dict("index"))
print("within tolerance: counts", [len(G["L4C2"]) == 31, len(G["L4out"]) == 10, len(G["L5"]) == 11],
"trials", within(meta.trials.mean(), 115, 60))
region layer
C2 4 31
5 11
surround 4 10
dtype: int64
L4 C2 31, L4 outside C2 10, L5 11
trials per neuron, all trials (a) = (b): mean 115.0, s.d. 60.6, n = 52
depth by group (um): {'L4C2': {'min': 410.0, 'max': 565.0}, 'L4out': {'min': 405.0, 'max': 560.0}, 'L5': {'min': 627.0, 'max': 844.0}}
within tolerance: counts [True, True, True] trials True
Result: reproduced. The files hold 31 L4 C2, 10 L4 surround and 11 L5 C2 recordings, one per file, with 115.0 +- 60.6 trials per neuron (mean +- s.d.), the paper's 115 +- 60. The depths fit the Methods' ranges (L4 405 to 565 um, L5 627 to 844 um).
recordings-mice¶
The 52 recordings come from 21 mice of four lines, 21 neurons from 7 C57BL/6 mice, 6 from 3 VGAT-ChR2(H134R) mice, 17 from 8 PV-ires-cre mice and 8 from 3 PV-ires-cre x Ai32 mice.
gt = meta.groupby("genotype").agg(mice=("subject", "nunique"), neurons=("subject", "size"))
print(gt)
print("subjects:", meta.subject.nunique(), "| recordings:", len(meta))
print("subjects whose identifier starts with 'jf':", sorted(meta.subject[meta.subject.str.startswith("jf")].unique()))
print("cells per subject:", meta.groupby("subject").size().to_dict())
mice neurons
genotype
(blank) 3 6
Ai32 x PV-Cre 3 7
C57BL6 9 22
PV-Cre 8 17
subjects: 23 | recordings: 52
subjects whose identifier starts with 'jf': ['jf35892', 'jf37166', 'jf42400']
cells per subject: {'anm101105': 3, 'anm106211': 1, 'anm106213': 4, 'anm131970': 1, 'anm134333': 3, 'anm17702': 2, 'anm211295': 1, 'anm215592': 4, 'anm217330': 1, 'anm226692': 1, 'anm226694': 3, 'anm226695': 3, 'anm234230': 1, 'anm234231': 1, 'anm234232': 1, 'anm244024': 4, 'anm244025': 2, 'anm244028': 3, 'anm266945': 3, 'anm266951': 2, 'jf35892': 4, 'jf37166': 1, 'jf42400': 3}
Result: not reproduced. The dandiset has 52 recordings, as the paper, but 23 subjects, not
21, and its genotypes do not give the paper's counts per line: C57BL6 9 mice and 22 neurons
(paper 7 and 21), Ai32 x PV-Cre 3 mice and 7 neurons (paper 3 and 8), PV-Cre 8 mice and 17
neurons (as the paper), and 3 mice and 6 neurons with a blank genotype, which match the
paper's VGAT-ChR2 counts but carry no label. Two subject identifiers more than the paper's mice
may be one animal under two names (the jf and anm series), but nothing in the files says so.
behavior-performance¶
Across the 52 sessions the fraction of trials correct was 0.740 +- 0.086 (mean +- s.d.).
Correct trials are hits and correct rejections. Reading (a), all trials, and reading (b), the dataset's own marking of the trials used, are the same: the files mark none.
perf = pd.Series({c.cell_no: c.trials.response.isin(["Hit", "CR"]).mean() for c in cells})
print("fraction correct per session: mean %.3f, s.d. %.3f, range %.3f to %.3f, n = %d" % (perf.mean(), perf.std(), perf.min(), perf.max(), len(perf)))
print("within 0.740 +- 0.086:", within(perf.mean(), 0.740, 0.086))
print("by group:", perf.groupby(base.group).mean().round(3).to_dict())
fraction correct per session: mean 0.741, s.d. 0.085, range 0.579 to 0.935, n = 52
within 0.740 +- 0.086: True
by group: {'L4C2': 0.727, 'L4out': 0.721, 'L5': 0.8}
Result: reproduced. The fraction of trials correct is 0.741 +- 0.085 over the 52 sessions, the paper's 0.740 +- 0.086.
behavior-sample-epoch¶
The sample epoch, when the pole is presented, lasted 1.54 to 4.05 s, 2.39 s on average.
The sample epoch is the period pole_available is true in a trial (first sample on to last
sample on). Readings (a) over all trials pooled and (b) over per-session means. Trials whose
pole is still in reach when the 4.5 to 5 s recording ends are cut there; they are counted, and
pole_out_time - pole_in_time from the trials table is shown as a variant.
ep_rows = []
for c in cells:
for k in range(len(c.trials)):
a = c.seg_start[k]
b = c.seg_start[k + 1] if k + 1 < len(c.seg_start) else len(c.vt)
r = hires.runs(c.pole[a:b])
if len(r) == 0:
continue
ep_rows.append(dict(cell=c.cell_no, trial=k, dur=(r[-1, 1] - r[0, 0]) / 1000, n_runs=len(r),
cut=bool(c.pole[a] or c.pole[b - 1]), on_ms=r[0, 0], off_ms=r[-1, 1],
table=c.trials.pole_out_time.iloc[k] - c.trials.pole_in_time.iloc[k]))
EP = pd.DataFrame(ep_rows)
sess = EP.groupby("cell").dur.mean()
print("trials %d, with more than one run of pole_available %d, cut at the end of the recording %d" % (len(EP), (EP.n_runs > 1).sum(), EP.cut.sum()))
for name, x in [("(a) pooled trials", EP.dur), ("(b) session means", sess)]:
print(f"{name}: mean {x.mean():.3f} s, min {x.min():.3f}, max {x.max():.3f}; within tolerance:",
[within(x.mean(), 2.39, 0.239), within(x.min(), 1.54, 0.154), within(x.max(), 4.05, 0.405)])
x = EP.dur[~EP.cut]
print(f"(a) without the cut trials: mean {x.mean():.3f}, min {x.min():.3f}, max {x.max():.3f}")
print(f"variant, pole_out_time - pole_in_time: mean {EP.table.mean():.3f}, min {EP.table.min():.3f}, max {EP.table.max():.3f}")
trials 5978, with more than one run of pole_available 0, cut at the end of the recording 617 (a) pooled trials: mean 2.394 s, min 1.543, max 4.049; within tolerance: [True, True, True] (b) session means: mean 2.403 s, min 1.930, max 3.125; within tolerance: [True, False, False] (a) without the cut trials: mean 2.235, min 1.543, max 3.990 variant, pole_out_time - pole_in_time: mean 2.777, min 1.602, max 22.479
Result: reproduced. Over all trials pooled (reading a) the pole is in reach for 2.394 s on average, 1.543 s at least and 4.049 s at most, the paper's 2.39 s and 1.54 to 4.05 s. Over session means (reading b) the mean is 2.403 s but the range is 1.930 to 3.125 s, outside the tolerance; the claim holds under reading (a). 617 of 5978 trials are cut by the end of the recording, and the paper's range is reproduced with them in.
behavior-whisking-bouts¶
Mice whisked in bouts with a mean duration of 261 ms, a peak-to-peak amplitude of 15.7 deg and a frequency of 15.4 Hz, separated by rest.
A bout is a contiguous run of amplitude > 2.5 deg in one trial. Open details from the claims
file: (1) peak-to-peak as twice the bout's maximum Hilbert amplitude (2amp) or the largest
range of band-passed theta within a cycle of the bout (range); (2) whether touch breaks a
bout. Frequency is 1/cycle length over every whisking cycle (amplitude > 2.5 deg throughout),
pooled. As a variant not in the claims file, bouts holding at least one complete cycle.
bout_rows = []
freqs = []
for c in cells:
d = D[c.cell_no]
wc = d["cyc"][d["cyc"].whisking]
freqs.append(1000.0 / wc.dur.to_numpy())
for tb in (False, True):
bt = hires.bouts(c, d["T"], tb, d["cyc"])
bout_rows.append(bt.assign(cell=c.cell_no, touch_breaks=tb))
BT = pd.concat(bout_rows, ignore_index=True)
freq = np.concatenate(freqs)
print(f"whisking cycles {len(freq)}: frequency mean {freq.mean():.2f} Hz, median {np.median(freq):.2f}; within 15.4 +- 1.54: {within(freq.mean(), 15.4, 1.54)}")
rows = []
for tb in (False, True):
for only_cycles in (False, True):
b = BT[(BT.touch_breaks == tb) & ((BT.n_cycles > 0) | (not only_cycles))]
for p2p in ("2amp", "range"):
amp = 2 * b.amp_max if p2p == "2amp" else b.p2p_cycle
rows.append(dict(touch_breaks=tb, bouts=">=1 cycle (variant)" if only_cycles else "all runs", p2p=p2p, n_bouts=len(b),
dur_ms=b.dur.mean(), median_dur=b.dur.median(), amp_deg=amp.mean(),
dur_ok=within(b.dur.mean(), 261, 26.1), amp_ok=within(amp.mean(), 15.7, 1.57)))
BTS = pd.DataFrame(rows)
print(BTS.round(2).to_string(index=False))
print("bout duration range %d to %d ms; bouts of 1 to 10 ms: %.2f of all runs" % (BT.dur.min(), BT.dur.max(), np.mean(BT[~BT.touch_breaks].dur <= 10)))
whisking cycles 67406: frequency mean 15.76 Hz, median 14.29; within 15.4 +- 1.54: True
touch_breaks bouts p2p n_bouts dur_ms median_dur amp_deg dur_ok amp_ok
False all runs 2amp 102986 91.64 49.0 15.74 False True
False all runs range 102986 91.64 49.0 10.92 False False
False >=1 cycle (variant) 2amp 26372 238.31 178.0 33.91 True False
False >=1 cycle (variant) range 26372 238.31 178.0 25.62 True False
True all runs 2amp 124654 69.71 38.0 16.82 False True
True all runs range 124654 69.71 38.0 10.77 False False
True >=1 cycle (variant) 2amp 23108 213.78 165.0 35.66 False False
True >=1 cycle (variant) range 23108 213.78 165.0 25.56 False False
bout duration range 1 to 2477 ms; bouts of 1 to 10 ms: 0.10 of all runs
Result: not reproduced. The amplitude and the frequency match: twice the maximum Hilbert amplitude averages 15.74 deg over bouts that touch does not break (15.7 in the paper) and the whisking frequency 15.76 Hz (15.4). The mean bout duration does not: 91.6 ms with touch not breaking bouts and 69.7 ms with it, against 261 +- 26.1 ms, so no reading puts all three values in tolerance, which the criterion requires. Most runs of amplitude above 2.5 deg are short; counting only bouts that hold a complete cycle (not a reading of the claims file) gives 213.8 to 238.3 ms, but then the amplitude is 33.9 to 35.7 deg.
behavior-touches-reaction¶
Mice touched the pole 2.33 times on average before responding, and the reaction time from first touch to first lick was 367 +- 234 ms (mean +- s.d.).
Reaction time: first lick (first_lick_time) minus the first touch onset, every trial with a
lick and a touch, pooled. Touches before the response: (a) trials with a lick, touches before
the first lick; (b) all trials with a touch, touches before the first lick or, without a lick,
touches in the sample epoch. Variants: only trials whose first lick follows the first touch,
and the first lick taken from the lick series.
rt_rows = []
for c in cells:
T = D[c.cell_no]["T"]
on_t = c.seg_t0[T.trial] + (T.on - c.seg_start[T.trial]) / 1000
for k in range(len(c.trials)):
fl = c.trials.first_lick_time.iloc[k]
ons = on_t[T.trial == k]
in_pole = c.pole[T.on[T.trial == k]]
a = c.seg_start[k]
b = c.seg_start[k + 1] if k + 1 < len(c.seg_start) else len(c.vt)
li = np.flatnonzero(c.lick[a:b])
fl2 = c.seg_t0[k] + li[0] / 1000 if len(li) else np.nan
rt_rows.append(dict(cell=c.cell_no, trial=k, lick=not np.isnan(fl), n_touch=len(ons),
rt=(fl - ons[0]) * 1000 if len(ons) and not np.isnan(fl) else np.nan,
rt_series=(fl2 - ons[0]) * 1000 if len(ons) and not np.isnan(fl2) else np.nan,
before_lick=int(np.sum(ons < fl)) if not np.isnan(fl) else np.nan,
in_sample=int(in_pole.sum())))
RT = pd.DataFrame(rt_rows)
rt = RT.rt.dropna()
print(f"reaction time, trials with a lick and a touch: n {len(rt)}, mean {rt.mean():.0f} ms, s.d. {rt.std():.0f}; negative (lick before first touch) {np.sum(rt < 0)}; within 367 +- 234: {within(rt.mean(), 367, 234)}")
rp = rt[rt > 0]
print(f"variant, first lick after the first touch: n {len(rp)}, mean {rp.mean():.0f} ms, s.d. {rp.std():.0f}")
rs = RT.rt_series.dropna()
print(f"variant, first lick from the lick series: n {len(rs)}, mean {rs.mean():.0f} ms, s.d. {rs.std():.0f}")
na = RT[RT.lick].before_lick
nb_ = np.where(RT.lick, RT.before_lick, RT.in_sample)[RT.n_touch > 0]
print(f"touches before the response: (a) {na.mean():.2f} over {len(na)} trials with a lick; (b) {np.mean(nb_):.2f} over {len(nb_)} trials with a touch; within 2.33 +- 0.233: {within(na.mean(), 2.33, 0.233)}, {within(np.mean(nb_), 2.33, 0.233)}")
nap = RT[RT.lick & (RT.rt > 0)].before_lick
print(f"variant, trials whose first lick follows the first touch: {nap.mean():.2f} touches over {len(nap)} trials")
reaction time, trials with a lick and a touch: n 3093, mean 190 ms, s.d. 700; negative (lick before first touch) 544; within 367 +- 234: True variant, first lick after the first touch: n 2549, mean 420 ms, s.d. 411 variant, first lick from the lick series: n 3108, mean 168 ms, s.d. 644 touches before the response: (a) 3.00 over 3562 trials with a lick; (b) 3.51 over 3875 trials with a touch; within 2.33 +- 0.233: False, False variant, trials whose first lick follows the first touch: 4.20 touches over 2549 trials
Result: not reproduced. The number of touches before the response is 3.00 under reading (a) and 3.51 under (b), outside 2.33 +- 0.233 under both. The reaction time is 190 ms on average over 3093 trials with a lick and a touch, inside the claim's tolerance only because the tolerance is the paper's s.d. (367 +- 234 ms); its s.d. here is 700 ms, as 544 trials have their first lick before the first touch. Counting only trials whose first lick follows the first touch (a variant) gives 420 +- 411 ms.
example-cell-rates¶
The example L4 C2 neuron of Figure 1 fired at 0.73 spikes/s in the 500 ms before the sample epoch and 1.91 spikes/s during it, a 2.6-fold increase (Wilcoxon rank sum, p = 3.38e-7).
For every L4 C2 cell: per-trial rates in the 500 ms before pole_available turns on and over
the time it is on; the cell's rates are the means of the per-trial rates (rates from pooled
spikes and time are shown too); the rank-sum test compares the per-trial rates.
ex_rows = []
EP_by = {no: g for no, g in EP.groupby("cell")}
for no in G["L4C2"]:
c = by_no[no]
s = np.sort(D[no]["s"])
e = EP_by[no]
on_vt = e.trial.to_numpy() * hires.TRIAL_STRIDE_MS + e.on_ms.to_numpy()
off_vt = e.trial.to_numpy() * hires.TRIAL_STRIDE_MS + e.off_ms.to_numpy()
ok = e.on_ms.to_numpy() >= 500
nb0 = hires.counts_in(s, on_vt[ok] - 500, 500)
ns = np.searchsorted(s, off_vt[ok]) - np.searchsorted(s, on_vt[ok])
dur = (off_vt[ok] - on_vt[ok]) / 1000
rb, rs_ = nb0 / 0.5, ns / dur
ex_rows.append(dict(cell=no, trials=int(ok.sum()), base=rb.mean(), samp=rs_.mean(), fold=rs_.mean() / rb.mean() if rb.mean() > 0 else np.inf,
base_pooled=nb0.sum() / (0.5 * ok.sum()), sample_pooled=ns.sum() / dur.sum(),
p=stats.ranksums(rb, rs_).pvalue))
EX = pd.DataFrame(ex_rows).set_index("cell")
EX["match"] = [within(r.base, 0.73, 0.073) and within(r.samp, 1.91, 0.191) and within(r.fold, 2.6, 0.26) and r.p < 0.05 for r in EX.itertuples()]
EX["distance"] = np.sqrt(((EX.base - 0.73) / 0.073) ** 2 + ((EX.samp - 1.91) / 0.191) ** 2 + ((EX.fold - 2.6) / 0.26) ** 2)
EX["fold_pooled"] = EX.sample_pooled / EX.base_pooled
print(EX.sort_values("distance").to_string(formatters={"p": "{:.3g}".format}, float_format="{:.3f}".format))
print("cell 17 has no trial whose pole comes into reach 500 ms or more after the start of the recording, so it has no baseline")
EXC = int(EX.distance.idxmin())
print("cells matching all three and p < 0.05:", EX.index[EX.match].tolist(), "| nearest:", EXC)
trials base samp fold base_pooled sample_pooled p match distance fold_pooled cell 23 121 0.727 1.754 2.411 0.727 1.909 1.36e-07 True 1.095 2.625 8 30 1.067 2.024 1.898 1.067 1.890 0.0085 False 5.378 1.772 15 133 1.008 1.336 1.326 1.008 1.335 3.5e-06 False 6.891 1.325 12 145 0.234 0.532 2.271 0.234 0.565 1.14e-06 False 9.985 2.408 6 48 0.625 0.526 0.842 0.625 0.530 0.00108 False 10.016 0.848 28 88 0.205 0.653 3.194 0.205 0.680 0.000101 False 10.016 3.324 24 87 0.230 0.491 2.135 0.230 0.493 7.32e-05 False 10.264 2.145 4 55 0.218 0.939 4.303 0.218 0.905 5.45e-07 False 10.858 4.146 20 158 0.190 0.957 5.039 0.190 0.929 3.67e-13 False 12.949 4.894 26 198 0.121 0.591 4.873 0.121 0.686 3.61e-14 False 13.918 5.657 25 103 0.466 2.853 6.122 0.466 2.440 5.39e-13 False 14.866 5.236 10 241 0.481 3.364 6.989 0.481 3.272 4.14e-37 False 18.827 6.798 30 98 0.163 1.151 7.048 0.163 1.101 3.61e-15 False 19.204 6.741 11 80 2.100 3.104 1.478 2.100 3.176 8.58e-06 False 20.246 1.512 27 224 0.152 1.219 8.032 0.152 1.338 7.77e-24 False 22.634 8.814 14 146 0.247 2.516 10.204 0.247 2.914 8.8e-23 False 30.155 11.817 3 64 0.438 4.735 10.823 0.438 4.995 2.8e-10 False 35.142 11.416 13 63 2.413 7.147 2.962 2.413 6.024 8.31e-06 False 35.847 2.497 2 19 3.684 3.492 0.948 3.684 3.225 0.907 False 41.794 0.875 18 63 4.635 9.158 1.976 4.635 8.917 9.53e-08 False 65.630 1.924 5 34 0.059 1.158 19.691 0.059 1.050 3.13e-05 False 66.489 17.850 29 50 0.040 0.829 20.732 0.040 0.871 1.92e-05 False 70.603 21.776 21 114 0.246 5.049 20.557 0.246 5.021 1.16e-22 False 71.302 20.441 31 96 3.146 14.100 4.482 3.146 14.069 1.63e-18 False 72.255 4.472 16 94 6.191 3.332 0.538 6.191 3.189 8.21e-08 False 75.602 0.515 1 48 0.042 1.579 37.896 0.042 1.634 2.62e-09 False 136.091 39.206 9 122 13.475 22.113 1.641 13.475 22.189 3.3e-10 False 204.170 1.647 7 26 0.000 0.603 inf 0.000 0.444 0.000866 False inf inf 22 24 0.000 2.452 inf 0.000 1.952 0.000206 False inf inf 19 98 0.000 0.504 inf 0.000 0.517 6.63e-09 False inf inf 17 0 NaN NaN inf NaN NaN NaN False NaN NaN cell 17 has no trial whose pole comes into reach 500 ms or more after the start of the recording, so it has no baseline cells matching all three and p < 0.05: [23] | nearest: 23
Result: reproduced. One L4 C2 cell matches, cell 23 (subject jf35892, 2010-02-19): 0.727 spikes/s in the 500 ms before the pole is in reach and 1.754 spikes/s during it (mean of per-trial rates; 1.909 spikes/s from pooled spikes and time), 2.41-fold (2.63 pooled), rank-sum p = 1.36e-7 over 121 trials. The paper gives 0.73, 1.91, 2.6-fold and p = 3.38e-7. The next nearest cell is 4.9 times farther in the scaled distance (5.38 against 1.10). Cell 23 is the example cell of the claims that follow.
example-cell-sample-ff¶
The Fano factor of spike counts over the whole sample epoch was large, 7.51, in the paragraph on the example neuron, which does not say whose value it is.
Per cell, the variance over the mean of the per-trial spike count while the pole is in reach. Readings: (a) the cell of example-cell-rates; (b) the mean over the 41 L4 cells; (c) the mean over all 52 cells.
ff_sample = {}
for c in cells:
s = np.sort(D[c.cell_no]["s"])
e = EP_by[c.cell_no]
on_vt = e.trial.to_numpy() * hires.TRIAL_STRIDE_MS + e.on_ms.to_numpy()
off_vt = e.trial.to_numpy() * hires.TRIAL_STRIDE_MS + e.off_ms.to_numpy()
ff_sample[c.cell_no] = hires.fano(np.searchsorted(s, off_vt) - np.searchsorted(s, on_vt))
ffs = pd.Series(ff_sample)
vals = {"(a) example cell %d" % EXC: ffs[EXC], "(b) mean over 41 L4 cells": ffs[G["L4C2"] + G["L4out"]].mean(),
"(c) mean over 52 cells": ffs.mean()}
for k, v in vals.items():
print(f"{k}: {v:.2f}; within 7.51 +- 0.751: {within(v, 7.51, 0.751)}")
print("L4 C2 cells within 7.51 +- 0.751:", ffs[G["L4C2"]][(ffs[G["L4C2"]] - 7.51).abs() <= 0.751].round(2).to_dict())
(a) example cell 23: 7.51; within 7.51 +- 0.751: True
(b) mean over 41 L4 cells: 6.31; within 7.51 +- 0.751: False
(c) mean over 52 cells: 7.89; within 7.51 +- 0.751: True
L4 C2 cells within 7.51 +- 0.751: {10: 7.72, 23: 7.51, 25: 7.57}
Result: reproduced. The sample-epoch Fano factor of cell 23 is 7.51 (reading a), the paper's value; the mean over all 52 cells is 7.89 (reading c, also within tolerance) and over the 41 L4 cells 6.31 (reading b, outside it).
example-cell-touch-window¶
For the example neuron, 74.4% of spikes fell in a window 8 to 40 ms after touch onset, which covers 15.3% of exploration time.
As fig3c-touch-spikes-l4c2 (below), for the cell of example-cell-rates, under every reading of the touch window. The window runs from the latency to latency + t_on.
ew = W[(W.cell == EXC) & (W.lat_reading == "a")].copy()
ew["end_ms"] = ew.latency + ew.t_on
ew["ok_frac"] = [within(100 * v, 74.4, 7.44) for v in ew.frac_spikes]
ew["ok_start"] = [within(v, 8, 1) for v in ew.latency]
ew["ok_end"] = [within(v, 40, 4) for v in ew.end_ms]
ew["ok_time"] = [within(100 * v, 15.3, 1.53) for v in ew.frac_time]
ew["all_ok"] = ew[["ok_frac", "ok_start", "ok_end", "ok_time"]].all(axis=1)
print(ew[WKEY[1:] + ["latency", "t_on", "end_ms", "frac_spikes", "frac_time", "all_ok"]].round(3).to_string(index=False))
print("readings with all four within tolerance: %d of %d" % (ew.all_ok.sum(), len(ew)))
for st, g in ew.groupby("start"):
print(f"start = {st}{' (the criterion)' if st == 'cross' else ' (variant)'}: spikes {100 * g.frac_spikes.min():.1f} to {100 * g.frac_spikes.max():.1f} %, "
f"window {g.latency.min():.0f} to {g.end_ms.min():.0f}-{g.end_ms.max():.0f} ms, time {100 * g.frac_time.min():.1f} to {100 * g.frac_time.max():.1f} %, "
f"all four within tolerance in {g.all_ok.sum()} of {len(g)}")
shuffle nearest tmax start latency t_on end_ms frac_spikes frac_time all_ok uniform before 100 first 8.0 29 37.0 0.846 0.196 False uniform before 100 cross 8.0 29 37.0 0.846 0.196 False uniform before 200 first 8.0 29 37.0 0.846 0.196 False uniform before 200 cross 8.0 29 37.0 0.846 0.196 False uniform either 100 first 8.0 28 36.0 0.784 0.190 False uniform either 100 cross 8.0 28 36.0 0.784 0.190 False uniform either 200 first 8.0 30 38.0 0.791 0.202 False uniform either 200 cross 8.0 30 38.0 0.791 0.202 False circular before 100 first 8.0 28 36.0 0.836 0.190 False circular before 100 cross 8.0 28 36.0 0.836 0.190 False circular before 200 first 8.0 28 36.0 0.836 0.190 False circular before 200 cross 8.0 28 36.0 0.836 0.190 False circular either 100 first 8.0 28 36.0 0.784 0.190 False circular either 100 cross 8.0 28 36.0 0.784 0.190 False circular either 200 first 8.0 28 36.0 0.784 0.190 False circular either 200 cross 8.0 28 36.0 0.784 0.190 False readings with all four within tolerance: 0 of 16 start = cross (the criterion): spikes 78.4 to 84.6 %, window 8 to 36-38 ms, time 19.0 to 20.2 %, all four within tolerance in 0 of 8 start = first (variant): spikes 78.4 to 84.6 %, window 8 to 36-38 ms, time 19.0 to 20.2 %, all four within tolerance in 0 of 8
Result: not reproduced. The window covers 19.0 to 20.2% of exploration time under every reading, against 15.3 +- 1.53%, so no reading puts the values in tolerance. For cell 23 the transition cutoff and the literal one give the same window under every reading: it starts at the 8 ms latency and ends at 36 to 38 ms (paper 8 to 40 ms, within tolerance). The fraction of exploration spikes in it is 78.4 to 79.1% when the closest onset is read on either side (within 74.4 +- 7.44%) and 83.6 to 84.6% when it is read as the onset before the spike (outside).
text-spikes-per-touch¶
Each touch evoked on the order of one spike (1.51 spikes per first touch, 0.31 per later touch) with an onset latency of 8 ms.
Readings: (a) the example cell, (b) the L4 C2 mean, (c) the L4 C2 median; spikes per touch in the cell's touch window under every window reading; latency from first touches (and, as for table1, from all touches).
rows = []
wl = W[(W.lat_reading == "a")]
for key, g in wl[wl.cell.isin(G["L4C2"])].groupby(WKEY[1:]):
g = g.set_index("cell")
for lat_src in ("lat_first", "lat_all"):
lat = base.loc[G["L4C2"], lat_src]
for rd, f in [("(a) example", lambda x: x.loc[EXC]), ("(b) mean", np.nanmean), ("(c) median", np.nanmedian)]:
v1, v2, v3 = f(g.spt_first), f(g.spt_later), f(lat)
rows.append(dict(zip(WKEY[1:], key), latency_from=lat_src, reading=rd, first=v1, later=v2, latency=v3,
ok=within(v1, 1.51, 0.151) and within(v2, 0.31, 0.031) and within(v3, 8, 1)))
TS = pd.DataFrame(rows)
print(TS.groupby(["start", "reading", "latency_from"]).agg(first_min=("first", "min"), first_max=("first", "max"), later_min=("later", "min"),
later_max=("later", "max"), latency=("latency", "mean"), n_ok=("ok", "sum"), n=("ok", "size")).round(3))
print(TS[TS.ok].to_string(index=False) if TS.ok.any() else "no reading puts all three values in tolerance")
first_min first_max later_min later_max latency n_ok n
start reading latency_from
cross (a) example lat_all 1.474 1.505 0.302 0.307 9.000 8 8
lat_first 1.474 1.505 0.302 0.307 8.000 8 8
(b) mean lat_all 1.164 1.282 0.623 0.685 8.419 0 8
lat_first 1.164 1.282 0.623 0.685 7.290 0 8
(c) median lat_all 0.704 0.732 0.314 0.331 8.000 0 8
lat_first 0.704 0.732 0.314 0.331 7.000 0 8
first (a) example lat_all 1.474 1.505 0.302 0.307 9.000 8 8
lat_first 1.474 1.505 0.302 0.307 8.000 8 8
(b) mean lat_all 0.529 0.843 0.279 0.434 8.419 0 8
lat_first 0.529 0.843 0.279 0.434 7.290 0 8
(c) median lat_all 0.175 0.691 0.100 0.290 8.000 0 8
lat_first 0.175 0.691 0.100 0.290 7.000 0 8
shuffle nearest tmax start latency_from reading first later latency ok
circular before 100 cross lat_first (a) example 1.474227 0.302376 8.0 True
circular before 100 cross lat_all (a) example 1.474227 0.302376 9.0 True
circular before 100 first lat_first (a) example 1.474227 0.302376 8.0 True
circular before 100 first lat_all (a) example 1.474227 0.302376 9.0 True
circular before 200 cross lat_first (a) example 1.474227 0.302376 8.0 True
circular before 200 cross lat_all (a) example 1.474227 0.302376 9.0 True
circular before 200 first lat_first (a) example 1.474227 0.302376 8.0 True
circular before 200 first lat_all (a) example 1.474227 0.302376 9.0 True
circular either 100 cross lat_first (a) example 1.474227 0.302376 8.0 True
circular either 100 cross lat_all (a) example 1.474227 0.302376 9.0 True
circular either 100 first lat_first (a) example 1.474227 0.302376 8.0 True
circular either 100 first lat_all (a) example 1.474227 0.302376 9.0 True
circular either 200 cross lat_first (a) example 1.474227 0.302376 8.0 True
circular either 200 cross lat_all (a) example 1.474227 0.302376 9.0 True
circular either 200 first lat_first (a) example 1.474227 0.302376 8.0 True
circular either 200 first lat_all (a) example 1.474227 0.302376 9.0 True
uniform before 100 cross lat_first (a) example 1.505155 0.304536 8.0 True
uniform before 100 cross lat_all (a) example 1.505155 0.304536 9.0 True
uniform before 100 first lat_first (a) example 1.505155 0.304536 8.0 True
uniform before 100 first lat_all (a) example 1.505155 0.304536 9.0 True
uniform before 200 cross lat_first (a) example 1.505155 0.304536 8.0 True
uniform before 200 cross lat_all (a) example 1.505155 0.304536 9.0 True
uniform before 200 first lat_first (a) example 1.505155 0.304536 8.0 True
uniform before 200 first lat_all (a) example 1.505155 0.304536 9.0 True
uniform either 100 cross lat_first (a) example 1.474227 0.302376 8.0 True
uniform either 100 cross lat_all (a) example 1.474227 0.302376 9.0 True
uniform either 100 first lat_first (a) example 1.474227 0.302376 8.0 True
uniform either 100 first lat_all (a) example 1.474227 0.302376 9.0 True
uniform either 200 cross lat_first (a) example 1.505155 0.306695 8.0 True
uniform either 200 cross lat_all (a) example 1.505155 0.306695 9.0 True
uniform either 200 first lat_first (a) example 1.505155 0.306695 8.0 True
uniform either 200 first lat_all (a) example 1.505155 0.306695 9.0 True
Result: reproduced. Under reading (a), the example cell, there are 1.474 to 1.505 spikes per first touch and 0.302 to 0.307 per later touch with an 8 ms latency (9 ms from all touches), within tolerance under all 8 window readings of the transition cutoff (and the 8 of the literal variant, which give the same window for this cell). The L4 C2 mean (1.16 to 1.28 and 0.62 to 0.69 under the transition cutoff) and median (0.70 to 0.73 and 0.31 to 0.33), readings (b) and (c), do not match, which agrees with the claims file's note that the sentence is about the example cell.
text-whisking-rates¶
Spike rates were low in both whisking (0.67 spikes/s) and non-whisking (0.23 spikes/s) periods.
Readings: (a) the example cell; (b) the L4 C2 mean; (c) the L4 C2 median; (d) all L4 C2 spikes over all L4 C2 time pooled; whisking read without and with touch epochs. Whisking leaves out licking cycles and the 70 ms after each touch; the earlier run's mask without these exclusions is shown as a variant.
rows = []
tot = {}
for no in G["L4C2"]:
c = by_no[no]
T = D[no]["T"]
s = D[no]["s"]
si = hires.spike_samples(c, s)
for wt in (False, True):
for ex in (True, False):
m = hires.whisk_mask(c, T, wt, exclude=ex)
tot.setdefault(("w", wt, ex), [0, 0])
tot[("w", wt, ex)][0] += m[si].sum()
tot[("w", wt, ex)][1] += m.sum() / 1000
m = hires.nonwhisk_mask(c, T)
tot.setdefault(("n",), [0, 0])
tot[("n",)][0] += m[si].sum()
tot[("n",)][1] += m.sum() / 1000
b4 = base.loc[G["L4C2"]]
for ex in (True, False):
for wt in (False, True):
col = ("rate_whisk_with_touch" if wt else "rate_whisk_no_touch") + ("" if ex else "_noexcl")
for rd, w_, n_ in [("(a) example", b4.loc[EXC, col], b4.loc[EXC, "rate_nonwhisk"]),
("(b) mean", b4[col].mean(), b4.rate_nonwhisk.mean()),
("(c) median", b4[col].median(), b4.rate_nonwhisk.median()),
("(d) pooled", tot[("w", wt, ex)][0] / tot[("w", wt, ex)][1], tot[("n",)][0] / tot[("n",)][1])]:
rows.append(dict(exclusions=("licking and 70 ms after touch left out" if ex else "none (earlier run, variant)"),
whisking_with_touch=wt, reading=rd, whisking=w_, nonwhisking=n_,
ok=within(w_, 0.67, 0.067) and within(n_, 0.23, 0.023)))
TW = pd.DataFrame(rows)
print(TW.round(3).to_string(index=False))
print("L4 C2 cells with both rates within tolerance (whisking without / with touch):",
b4.index[[within(r.rate_whisk_no_touch, 0.67, 0.067) and within(r.rate_nonwhisk, 0.23, 0.023) for r in b4.itertuples()]].tolist(),
b4.index[[within(r.rate_whisk_with_touch, 0.67, 0.067) and within(r.rate_nonwhisk, 0.23, 0.023) for r in b4.itertuples()]].tolist())
exclusions whisking_with_touch reading whisking nonwhisking ok
licking and 70 ms after touch left out False (a) example 0.777 0.227 False
licking and 70 ms after touch left out False (b) mean 1.462 1.303 False
licking and 70 ms after touch left out False (c) median 0.376 0.315 False
licking and 70 ms after touch left out False (d) pooled 1.475 1.406 False
licking and 70 ms after touch left out True (a) example 2.195 0.227 False
licking and 70 ms after touch left out True (b) mean 3.698 1.303 False
licking and 70 ms after touch left out True (c) median 2.270 0.315 False
licking and 70 ms after touch left out True (d) pooled 3.436 1.406 False
none (earlier run, variant) False (a) example 1.257 0.227 False
none (earlier run, variant) False (b) mean 2.402 1.303 False
none (earlier run, variant) False (c) median 0.860 0.315 False
none (earlier run, variant) False (d) pooled 2.634 1.406 False
none (earlier run, variant) True (a) example 2.373 0.227 False
none (earlier run, variant) True (b) mean 3.977 1.303 False
none (earlier run, variant) True (c) median 2.284 0.315 False
none (earlier run, variant) True (d) pooled 3.970 1.406 False
L4 C2 cells with both rates within tolerance (whisking without / with touch): [] []
Result: not reproduced. No reading puts both rates in tolerance. The example cell's non-whisking rate is 0.227 spikes/s (paper 0.23), but its whisking rate is 0.777 spikes/s with touch epochs left out and 2.195 with them in, against 0.67 +- 0.067. The L4 C2 means (1.462 or 3.698 and 1.303), medians (0.376 or 2.270 and 0.315) and pooled rates (1.475 or 3.436 and 1.406) are further off. Leaving out the 70 ms after touch and the licking cycles lowers the example cell's whisking rate from 1.257 (the earlier mask, a variant) to 0.777.
example-cell-phase-modulation¶
During exploration the example neuron's spikes are coupled to whisking phase, with a modulation depth close to 1 and near zero spiking at phase 0.
The criterion takes spikes in whisking of exploration periods outside touch, so the phase-analysis cycles are counted in exploration time (pole in reach). The fit over the same cycles in the whole trial, as Table 1 uses, is shown too.
r = base.loc[EXC]
cen = np.linspace(-np.pi, np.pi, 13)[:-1] + np.pi / 12
print(f"cell {EXC}, phase-analysis cycles in exploration time (pole in reach): modulation depth {r.depth_mod_x:.3f}, from {int(r.phase_spikes_x)} spikes")
print(f"cell {EXC}, phase-analysis cycles over the whole trial (as Table 1): modulation depth {r.depth_mod:.3f}, preferred phase {np.degrees(r.phi_pref):.0f} deg, "
f"from {int(r.phase_spikes)} spikes in {r.phase_time_s:.1f} s")
print("rate per 30 deg bin (spikes/s), bin centres in deg:")
print(pd.DataFrame({"exploration": D[EXC]["phase_rate_x"], "whole trial": D[EXC]["phase_rate"]}, index=np.round(np.degrees(cen)).astype(int)).round(2).to_string())
print("spike counts per bin, exploration:", D[EXC]["phase_cnt_x"].tolist(), "| whole trial:", D[EXC]["phase_cnt"].tolist())
print("within 1 +- 0.1: exploration", within(r.depth_mod_x, 1, 0.1), "| whole trial", within(r.depth_mod, 1, 0.1))
fig, ax = plt.subplots(figsize=(4, 2.5))
ax.bar(np.degrees(cen), D[EXC]["phase_rate_x"], width=28, color="0.6", label="exploration")
ax.step(np.degrees(cen), D[EXC]["phase_rate"], where="mid", color="k", lw=0.8, label="whole trial")
ax.set_xlabel("whisking phase (deg)"); ax.set_ylabel("spikes/s"); ax.set_title(f"cell {EXC}"); ax.legend(fontsize=7)
plt.tight_layout(); plt.show()
both = np.r_[D[EXC]["phase_rate_x"], D[EXC]["phase_rate"]]
print("plotted range: %.2f to %.2f spikes/s, none outside the axes" % (np.nanmin(both), np.nanmax(both)))
cell 23, phase-analysis cycles in exploration time (pole in reach): modulation depth 1.000, from 14 spikes
cell 23, phase-analysis cycles over the whole trial (as Table 1): modulation depth 1.000, preferred phase 152 deg, from 92 spikes in 64.8 s
rate per 30 deg bin (spikes/s), bin centres in deg:
exploration whole trial
-165 1.22 2.00
-135 0.65 0.70
-105 1.32 0.70
-75 0.00 0.00
-45 0.00 0.00
-15 0.00 0.19
15 0.00 0.00
45 0.00 0.00
75 0.85 0.23
105 0.74 1.80
135 1.93 6.80
165 2.46 3.38
spike counts per bin, exploration: [2, 1, 2, 0, 0, 0, 0, 0, 1, 1, 3, 4] | whole trial: [12, 4, 4, 0, 0, 1, 0, 0, 1, 9, 40, 21]
within 1 +- 0.1: exploration True | whole trial True
plotted range: 0.00 to 6.80 spikes/s, none outside the axes
Result: reproduced. Cell 23's modulation depth is 1.000 in exploration time and 1.000 over the whole trial. In exploration time it fired no spike in the five 30 deg bins centred from -75 to 45 deg, around phase 0, as the claim states; over the whole trial one spike falls in them (0.19 spikes/s in the bin centred at -15 deg). The exploration-time fit rests on 14 spikes, so it is not precise; the whole-trial fit rests on 92.
Table 1: table1-l4-c2, table1-l4-outside, table1-l5¶
For the 31 L4 C2 cells Table 1 gives, as mean +- s.d. (median), spikes per first touch 1.36 +- 1.32 (0.87), per later touch 0.74 +- 0.87 (0.33), touch-evoked spikes 70.6 +- 20.9% (74.7%), touch and whisking spikes 75.2 +- 19.0% (80.3%), phase modulation depth 0.67 +- 0.32 (0.70), non-whisking rate 1.31 +- 2.88 (0.32) spikes/s, whisking rate 1.39 +- 2.33 (0.33) spikes/s, onset latency 7.8 +- 3.0 (8) ms and minimum ISI 2.9 +- 1.9 (2.3) ms.
(and the same measures for the 10 L4 cells outside C2 and the 11 L5 cells). Each value is computed under every combination of readings: the touch window (8 with the transition cutoff, the criterion's, and 8 more with the literal cutoff as a variant, and the two latency readings for cells without a first-touch latency), the phase window width, whisking with or without touch epochs, and the latency from first touches or from all touches. The criterion counts the values within tolerance (15 of 18 for L4 C2 and L5, 13 of 17 for L4 outside C2).
Whisking periods are read here as the claims file defines them, outside touch epochs ("Touch epochs are periods of contact. For other epochs theta is band-passed ... whisking epochs have amplitude >2.5 deg"; the answer about fig4c-phase-decoding quotes the Methods to the same effect), leaving out licking cycles and the 70 ms after touch. Whisking with touch epochs, which the whisking-rate claims list as a reading of their own, is shown here as a variant. Modulation depth uses every phase-analysis cycle of the trial. For L4 outside C2 the two medians of spikes per touch (0 and 0) are checked as exact and counted among the 17 values.
EXPECT = {
"L4C2": [("spt_first", "mean", 1.36, 1.32), ("spt_first", "median", 0.87, 0.087), ("spt_later", "mean", 0.74, 0.87),
("spt_later", "median", 0.33, 0.033), ("touch_pct", "mean", 70.6, 20.9), ("touch_pct", "median", 74.7, 7.47),
("tp_pct", "mean", 75.2, 19.0), ("tp_pct", "median", 80.3, 8.03), ("depth", "mean", 0.67, 0.32),
("depth", "median", 0.70, 0.07), ("nonwhisk", "mean", 1.31, 2.88), ("nonwhisk", "median", 0.32, 0.032),
("whisk", "mean", 1.39, 2.33), ("whisk", "median", 0.33, 0.033), ("latency", "mean", 7.8, 3.0),
("latency", "median", 8, 1), ("min_isi", "mean", 2.9, 1.9), ("min_isi", "median", 2.3, 0.23)],
"L4out": [("spt_first", "mean", 0.08, 0.19), ("spt_later", "mean", 0.02, 0.04), ("touch_pct", "mean", 9.2, 13.7),
("touch_pct", "median", 3.6, 0.36), ("tp_pct", "mean", 33.0, 30.2), ("tp_pct", "median", 23.8, 2.38),
("depth", "mean", 0.58, 0.36), ("depth", "median", 0.61, 0.061), ("nonwhisk", "mean", 2.20, 2.91),
("nonwhisk", "median", 1.10, 0.11), ("whisk", "mean", 2.41, 3.48), ("whisk", "median", 1.31, 0.131),
("latency", "mean", 18.5, 6.0), ("min_isi", "mean", 3.3, 1.3), ("min_isi", "median", 2.8, 0.28)],
"L5": [("spt_first", "mean", 1.95, 2.96), ("spt_first", "median", 0.98, 0.098), ("spt_later", "mean", 1.28, 1.88),
("spt_later", "median", 0.42, 0.042), ("touch_pct", "mean", 24.0, 18.6), ("touch_pct", "median", 24.4, 2.44),
("tp_pct", "mean", 31.5, 16.4), ("tp_pct", "median", 36.0, 3.6), ("depth", "mean", 0.18, 0.12),
("depth", "median", 0.13, 0.013), ("nonwhisk", "mean", 12.0, 9.87), ("nonwhisk", "median", 12.8, 1.28),
("whisk", "mean", 16.2, 15.3), ("whisk", "median", 13.5, 1.35), ("latency", "mean", 9.7, 5.2),
("latency", "median", 8, 1), ("min_isi", "mean", 3.7, 2.5), ("min_isi", "median", 2.7, 0.27)],
}
# the L4 outside row of the claims file lists 17 values; its expected list above has the 15 of them that
# are numbers with a tolerance, and its two medians of spikes per touch (0 and 0) are checked as exact and
# counted as the other two
ZERO_MEDIANS = {"L4out": ["spt_first", "spt_later"]}
THRESH = {"L4C2": (15, 10, 18), "L4out": (13, 9, 17), "L5": (15, 10, 18)}
per_cell = W.merge(PH, on=["cell", "group"] + WKEY, suffixes=("", "_ph"))
t1_rows = []
for g in GROUP_ORDER:
sub = per_cell[per_cell.group == g]
for key, kk in sub.groupby(PKEY):
kk = kk.set_index("cell")
for wt in (False, True):
for lat_src in ("lat_first", "lat_all"):
b = base.loc[kk.index]
vals = dict(spt_first=kk.spt_first, spt_later=kk.spt_later, touch_pct=100 * kk.frac_spikes, tp_pct=100 * kk.frac_tp,
depth=b.depth_mod, nonwhisk=b.rate_nonwhisk,
whisk=b.rate_whisk_with_touch if wt else b.rate_whisk_no_touch, latency=b[lat_src], min_isi=b.min_isi_ms)
row = dict(group=g, **dict(zip(PKEY, key)), whisk_with_touch=wt, latency_from=lat_src)
n_ok = 0
for m, stat, e, tol in EXPECT[g]:
v = msm(vals[m])[stat]
row[f"{m}_{stat}"] = v
ok = within(v, e, tol)
row[f"ok_{m}_{stat}"] = ok
n_ok += ok
for m in ZERO_MEDIANS.get(g, []):
v = msm(vals[m])["median"]
row[f"{m}_median"] = v
row[f"ok_{m}_median"] = bool(v == 0)
n_ok += bool(v == 0)
if g == "L4out":
row["spt_medians_zero"] = row["ok_spt_first_median"] and row["ok_spt_later_median"]
row["n_ok"] = n_ok
hi, lo, n = THRESH[g]
row["verdict"] = "reproduced" if n_ok >= hi else ("partially" if n_ok >= lo else "not reproduced")
t1_rows.append(row)
T1 = pd.DataFrame(t1_rows)
# the criterion's readings: the transition cutoff (answer in questions/000013.yaml) and whisking periods outside touch
# epochs (the claims file's definitions, "For other epochs ... whisking epochs have amplitude >2.5 deg")
T1["criterion"] = (T1.start == "cross") & ~T1.whisk_with_touch
FIRST = dict(lat_reading="a", shuffle="uniform", nearest="before", tmax=100, start="cross", width="pm", whisk_with_touch=False, latency_from="lat_first")
def first_reading(t):
"""The values under the first reading of every open detail, in the claims file's order."""
r = t[np.logical_and.reduce([t[k] == v for k, v in FIRST.items()])].iloc[0]
return [(f"{m}_{stat}", round(float(r[f"{m}_{stat}"]), 3), e, tol, bool(r[f"ok_{m}_{stat}"])) for m, stat, e, tol in EXPECT[r["group"]]]
for g in GROUP_ORDER:
t = T1[T1.group == g]
print(f"\n{g}: {len(t)} reading combinations; values within tolerance: min {t.n_ok.min()}, max {t.n_ok.max()} (of {THRESH[g][2]} values)")
print(t.verdict.value_counts().to_string())
print(t.groupby("start").n_ok.agg(["min", "max"]).to_string())
cols = [c for c in t.columns if c.startswith("ok_")]
print("fraction of combinations in which each value is within tolerance:")
print(t[cols].mean().round(2).to_string())
num = [c for c in t.columns if c.endswith("_mean") or c.endswith("_median")]
print("range of each value over the combinations:")
print(pd.DataFrame({"min": t[num].min(), "max": t[num].max()}).round(3).to_string())
if g == "L4out":
print("medians of spikes per first and later touch both 0:", t.spt_medians_zero.value_counts().to_dict())
print("\nmodulation depth, variants (the criterion's values use every phase-analysis cycle of the trial):")
for g in GROUP_ORDER:
b = base.loc[G[g]]
v20 = b.depth_mod[b.phase_spikes >= 20]
print(f" {g}: all cells with a fit mean {b.depth_mod.mean():.3f} median {b.depth_mod.median():.3f} (n {b.depth_mod.notna().sum()});"
f" cells with >= 20 phase spikes mean {v20.mean():.3f} median {v20.median():.3f} (n {len(v20)});"
f" exploration time only mean {b.depth_mod_x.mean():.3f} median {b.depth_mod_x.median():.3f} (n {b.depth_mod_x.notna().sum()})")
print("whisking rate, variant without the licking and after-touch exclusions (whisking outside touch epochs): medians",
{g: round(base.loc[G[g], "rate_whisk_no_touch_noexcl"].median(), 3) for g in GROUP_ORDER},
"means", {g: round(base.loc[G[g], "rate_whisk_no_touch_noexcl"].mean(), 3) for g in GROUP_ORDER})
L4C2: 256 reading combinations; values within tolerance: min 10, max 15 (of 18 values)
verdict
reproduced 128
partially 128
min max
start
cross 15 15
first 10 14
fraction of combinations in which each value is within tolerance:
ok_spt_first_mean 1.0
ok_spt_first_median 0.0
ok_spt_later_mean 1.0
ok_spt_later_median 0.5
ok_touch_pct_mean 0.75
ok_touch_pct_median 0.75
ok_tp_pct_mean 0.75
ok_tp_pct_median 0.75
ok_depth_mean 1.0
ok_depth_median 0.0
ok_nonwhisk_mean 1.0
ok_nonwhisk_median 1.0
ok_whisk_mean 1.0
ok_whisk_median 0.0
ok_latency_mean 1.0
ok_latency_median 1.0
ok_min_isi_mean 1.0
ok_min_isi_median 1.0
range of each value over the combinations:
min max
spt_first_mean 0.529298 1.282354
ok_spt_first_mean True True
spt_first_median 0.175439 0.732143
ok_spt_first_median False False
spt_later_mean 0.27892 0.684603
ok_spt_later_mean True True
spt_later_median 0.099548 0.331361
ok_spt_later_median False True
touch_pct_mean 38.616734 75.796184
ok_touch_pct_mean False True
touch_pct_median 47.325409 80.34398
ok_touch_pct_median False True
tp_pct_mean 41.935009 81.427983
ok_tp_pct_mean False True
tp_pct_median 48.308271 84.20398
ok_tp_pct_median False True
depth_mean 0.684925 0.684925
ok_depth_mean True True
depth_median 0.858605 0.858605
ok_depth_median False False
nonwhisk_mean 1.303091 1.303091
ok_nonwhisk_mean True True
nonwhisk_median 0.315315 0.315315
ok_nonwhisk_median True True
whisk_mean 1.46234 3.697774
ok_whisk_mean True True
whisk_median 0.375552 2.269603
ok_whisk_median False False
latency_mean 7.290323 8.419355
ok_latency_mean True True
latency_median 7.0 8.0
ok_latency_median True True
min_isi_mean 2.767742 2.767742
ok_min_isi_mean True True
min_isi_median 2.3 2.3
ok_min_isi_median True True
L4out: 256 reading combinations; values within tolerance: min 7, max 11 (of 17 values)
verdict
partially 152
not reproduced 104
min max
start
cross 7 9
first 10 11
fraction of combinations in which each value is within tolerance:
ok_spt_first_mean 0.69
ok_spt_first_median 0.0
ok_spt_later_mean 0.5
ok_spt_later_median 0.0
ok_touch_pct_mean 0.69
ok_touch_pct_median 0.0
ok_tp_pct_mean 0.98
ok_tp_pct_median 0.12
ok_depth_mean 1.0
ok_depth_median 0.0
ok_nonwhisk_mean 1.0
ok_nonwhisk_median 1.0
ok_whisk_mean 1.0
ok_whisk_median 0.0
ok_latency_mean 0.0
ok_latency_median NaN
ok_min_isi_mean 1.0
ok_min_isi_median 1.0
range of each value over the combinations:
min max
spt_first_mean 0.015231 0.372658
ok_spt_first_mean False True
spt_first_median 0.012518 0.138736
ok_spt_first_median False False
spt_later_mean 0.002922 0.228281
ok_spt_later_mean False True
spt_later_median 0.002369 0.086298
ok_spt_later_median False False
touch_pct_mean 1.526468 31.712281
ok_touch_pct_mean False True
touch_pct_median 1.147735 30.715696
ok_touch_pct_median False False
tp_pct_mean 18.9711 65.053496
ok_tp_pct_mean False True
tp_pct_median 3.575994 68.55976
ok_tp_pct_median False True
depth_mean 0.603831 0.603831
ok_depth_mean True True
depth_median 0.548106 0.548106
ok_depth_median False False
nonwhisk_mean 2.182879 2.182879
ok_nonwhisk_mean True True
nonwhisk_median 1.093322 1.093322
ok_nonwhisk_median True True
whisk_mean 2.424642 2.501074
ok_whisk_mean True True
whisk_median 1.463681 1.544407
ok_whisk_median False False
latency_mean 8.142857 10.0
ok_latency_mean False False
latency_median NaN NaN
ok_latency_median None None
min_isi_mean 3.11 3.11
ok_min_isi_mean True True
min_isi_median 2.6 2.6
ok_min_isi_median True True
medians of spikes per first and later touch both 0: {False: 256}
L5: 256 reading combinations; values within tolerance: min 8, max 13 (of 18 values)
verdict
partially 216
not reproduced 40
min max
start
cross 10 13
first 8 11
fraction of combinations in which each value is within tolerance:
ok_spt_first_mean 1.0
ok_spt_first_median 0.06
ok_spt_later_mean 1.0
ok_spt_later_median 0.12
ok_touch_pct_mean 0.94
ok_touch_pct_median 0.12
ok_tp_pct_mean 0.72
ok_tp_pct_median 0.22
ok_depth_mean 1.0
ok_depth_median 0.0
ok_nonwhisk_mean 1.0
ok_nonwhisk_median 1.0
ok_whisk_mean 1.0
ok_whisk_median 0.5
ok_latency_mean 1.0
ok_latency_median 0.0
ok_min_isi_mean 1.0
ok_min_isi_median 0.0
range of each value over the combinations:
min max
spt_first_mean 0.181116 2.303749
ok_spt_first_mean True True
spt_first_median 0.06383 0.980392
ok_spt_first_median False True
spt_later_mean 0.10613 1.48964
ok_spt_later_mean True True
spt_later_median 0.010753 0.613248
ok_spt_later_median False True
touch_pct_mean 5.098896 30.569613
ok_touch_pct_mean False True
touch_pct_median 1.258881 29.053647
ok_touch_pct_median False True
tp_pct_mean 10.202406 48.937939
ok_tp_pct_mean False True
tp_pct_median 2.125531 50.625869
ok_tp_pct_median False True
depth_mean 0.181606 0.181606
ok_depth_mean True True
depth_median 0.10563 0.10563
ok_depth_median False False
nonwhisk_mean 12.038534 12.038534
ok_nonwhisk_mean True True
nonwhisk_median 12.903281 12.903281
ok_nonwhisk_median True True
whisk_mean 16.660252 19.693169
ok_whisk_mean True True
whisk_median 11.521873 12.818982
ok_whisk_median False True
latency_mean 4.9 10.888889
ok_latency_mean True True
latency_median 4.0 10.0
ok_latency_median False False
min_isi_mean 2.627273 2.627273
ok_min_isi_mean True True
min_isi_median 1.5 1.5
ok_min_isi_median False False
modulation depth, variants (the criterion's values use every phase-analysis cycle of the trial):
L4C2: all cells with a fit mean 0.685 median 0.859 (n 27); cells with >= 20 phase spikes mean 0.506 median 0.444 (n 11); exploration time only mean 0.718 median 1.000 (n 25)
L4out: all cells with a fit mean 0.604 median 0.548 (n 10); cells with >= 20 phase spikes mean 0.505 median 0.496 (n 8); exploration time only mean 0.729 median 0.980 (n 10)
L5: all cells with a fit mean 0.182 median 0.106 (n 11); cells with >= 20 phase spikes mean 0.182 median 0.106 (n 11); exploration time only mean 0.212 median 0.141 (n 11)
whisking rate, variant without the licking and after-touch exclusions (whisking outside touch epochs): medians {'L4C2': np.float64(0.86), 'L4out': np.float64(1.748), 'L5': np.float64(10.962)} means {'L4C2': np.float64(2.402), 'L4out': np.float64(2.61), 'L5': np.float64(18.491)}
t = T1[T1.group == "L4C2"]
print("L4 C2: verdicts over the reading combinations:", t.verdict.value_counts().to_dict())
print("medians within tolerance, per combination: min %d, max %d of 9" % (t[[c for c in t.columns if c.startswith("ok_") and c.endswith("median")]].sum(axis=1).min(),
t[[c for c in t.columns if c.startswith("ok_") and c.endswith("median")]].sum(axis=1).max()))
print(t.groupby(["start", "tmax"]).n_ok.agg(["min", "max"]))
print("verdicts by start of the cutoff search (cross is the criterion's, first a variant) and whisking reading:")
print(t.groupby(["start", "whisk_with_touch", "latency_from"]).n_ok.agg(["min", "max"]).to_string())
print(t.groupby("start").verdict.value_counts().unstack(fill_value=0).to_string())
c = t[t.criterion]
print("under the criterion's readings (cross, whisking outside touch epochs): values within tolerance %d to %d over %d combinations;" % (c.n_ok.min(), c.n_ok.max(), len(c)),
"verdicts", c.verdict.value_counts().to_dict())
print(c.groupby(["tmax", "width"]).n_ok.agg(["min", "max"]).to_string())
print("values under the first reading of each open detail (measured, expected, tolerance, within):")
for row in first_reading(t):
print(" ", row)
L4 C2: verdicts over the reading combinations: {'reproduced': 128, 'partially': 128}
medians within tolerance, per combination: min 3, max 6 of 9
min max
start tmax
cross 100 15 15
200 15 15
first 100 10 10
200 14 14
verdicts by start of the cutoff search (cross is the criterion's, first a variant) and whisking reading:
min max
start whisk_with_touch latency_from
cross False lat_all 15 15
lat_first 15 15
True lat_all 15 15
lat_first 15 15
first False lat_all 10 14
lat_first 10 14
True lat_all 10 14
lat_first 10 14
verdict reproduced partially
start
cross 128 0
first 0 128
under the criterion's readings (cross, whisking outside touch epochs): values within tolerance 15 to 15 over 64 combinations; verdicts {'reproduced': 64}
min max
tmax width
100 pm 15 15
total 15 15
200 pm 15 15
total 15 15
values under the first reading of each open detail (measured, expected, tolerance, within):
('spt_first_mean', 1.24, 1.36, 1.32, True)
('spt_first_median', 0.704, 0.87, 0.087, False)
('spt_later_mean', 0.66, 0.74, 0.87, True)
('spt_later_median', 0.314, 0.33, 0.033, True)
('touch_pct_mean', 73.807, 70.6, 20.9, True)
('touch_pct_median', 79.979, 74.7, 7.47, True)
('tp_pct_mean', 81.059, 75.2, 19.0, True)
('tp_pct_median', 83.333, 80.3, 8.03, True)
('depth_mean', 0.685, 0.67, 0.32, True)
('depth_median', 0.859, 0.7, 0.07, False)
('nonwhisk_mean', 1.303, 1.31, 2.88, True)
('nonwhisk_median', 0.315, 0.32, 0.032, True)
('whisk_mean', 1.462, 1.39, 2.33, True)
('whisk_median', 0.376, 0.33, 0.033, False)
('latency_mean', 7.29, 7.8, 3.0, True)
('latency_median', 7.0, 8, 1, True)
('min_isi_mean', 2.768, 2.9, 1.9, True)
('min_isi_median', 2.3, 2.3, 0.23, True)
Result: reproduced. Under the criterion's readings, the transition cutoff with whisking periods outside touch epochs, 15 of the 18 values are within tolerance in all 64 combinations of the other readings (the 100 and 200 ms curves, the shuffle, the closest onset, the phase window width and the latency PCTH), at the threshold for reproduced. The three outside tolerance are medians: spikes per first touch (0.704 under the first reading of each detail, against 0.87 +- 0.087), modulation depth (0.859 against 0.70 +- 0.07) and whisking rate (0.376 against 0.33 +- 0.033). Spikes evoked by touch are 73.8% (median 80.0%) and by touch or phase 81.1% (83.3%) under the first reading. With touch epochs in the whisking rate (a variant) the count is also 15 of 18 in every combination; with the literal cutoff (a variant) it is 10 to 14. Twelve of the 27 cells with a phase fit have a depth of 1; with only the 11 cells that have at least 20 phase spikes, a variant, the depth median is 0.444.
t = T1[T1.group == "L4out"]
print("L4 outside C2: verdicts over the reading combinations:", t.verdict.value_counts().to_dict())
print(t.groupby(["start", "lat_reading"]).n_ok.agg(["min", "max"]))
print("verdicts by start of the cutoff search (cross is the criterion's, first a variant) and whisking reading:")
print(t.groupby(["start", "whisk_with_touch", "latency_from"]).n_ok.agg(["min", "max"]).to_string())
print(t.groupby("start").verdict.value_counts().unstack(fill_value=0).to_string())
c = t[t.criterion]
print("under the criterion's readings (cross, whisking outside touch epochs): values within tolerance %d to %d over %d combinations;" % (c.n_ok.min(), c.n_ok.max(), len(c)),
"verdicts", c.verdict.value_counts().to_dict())
print(c.groupby(["tmax", "width"]).n_ok.agg(["min", "max"]).to_string())
print("values under the first reading of each open detail (measured, expected, tolerance, within):")
for row in first_reading(t):
print(" ", row)
L4 outside C2: verdicts over the reading combinations: {'partially': 152, 'not reproduced': 104}
min max
start lat_reading
cross a 7 9
b 7 9
first a 10 11
b 10 11
verdicts by start of the cutoff search (cross is the criterion's, first a variant) and whisking reading:
min max
start whisk_with_touch latency_from
cross False lat_all 7 9
lat_first 7 9
True lat_all 7 9
lat_first 7 9
first False lat_all 10 11
lat_first 10 11
True lat_all 10 11
lat_first 10 11
verdict not reproduced partially
start
cross 104 24
first 0 128
under the criterion's readings (cross, whisking outside touch epochs): values within tolerance 7 to 9 over 64 combinations; verdicts {'not reproduced': 52, 'partially': 12}
min max
tmax width
100 pm 7 9
total 7 9
200 pm 7 9
total 7 9
values under the first reading of each open detail (measured, expected, tolerance, within):
('spt_first_mean', 0.31, 0.08, 0.19, False)
('spt_later_mean', 0.152, 0.02, 0.04, False)
('touch_pct_mean', 26.167, 9.2, 13.7, False)
('touch_pct_median', 23.759, 3.6, 0.36, False)
('tp_pct_mean', 56.667, 33.0, 30.2, True)
('tp_pct_median', 64.565, 23.8, 2.38, False)
('depth_mean', 0.604, 0.58, 0.36, True)
('depth_median', 0.548, 0.61, 0.061, False)
('nonwhisk_mean', 2.183, 2.2, 2.91, True)
('nonwhisk_median', 1.093, 1.1, 0.11, True)
('whisk_mean', 2.501, 2.41, 3.48, True)
('whisk_median', 1.464, 1.31, 0.131, False)
('latency_mean', 8.143, 18.5, 6.0, False)
('min_isi_mean', 3.11, 3.3, 1.3, True)
('min_isi_median', 2.6, 2.8, 0.28, True)
Result: blocked (methods_underspecified). Under the criterion's readings the count of values within tolerance is 7 to 9 of 17, which falls on both sides of the threshold for partially reproduced (9): 52 combinations give not reproduced and 12 partially, the difference coming from the touch-window readings. The medians of spikes per touch, now counted among the 17, are not 0 under any reading; the means of spikes per touch (0.31 and 0.15 under the first reading, against 0.08 and 0.02), the touch-evoked fraction (26.2%, median 23.8%, against 9.2% and 3.6%), the median modulation depth (0.548 against 0.61 +- 0.061), the median whisking rate (1.464 against 1.31 +- 0.131) and the latency mean (8.1 ms over the cells with a latency, against 18.5) are outside tolerance. Rates, mean depth and the minimum ISI are within it. The literal cutoff, a variant, gives 10 to 11.
t = T1[T1.group == "L5"]
print("L5: verdicts over the reading combinations:", t.verdict.value_counts().to_dict())
print(t.groupby(["start", "lat_reading"]).n_ok.agg(["min", "max"]))
print("verdicts by start of the cutoff search (cross is the criterion's, first a variant) and whisking reading:")
print(t.groupby(["start", "whisk_with_touch", "latency_from"]).n_ok.agg(["min", "max"]).to_string())
print(t.groupby("start").verdict.value_counts().unstack(fill_value=0).to_string())
c = t[t.criterion]
print("under the criterion's readings (cross, whisking outside touch epochs): values within tolerance %d to %d over %d combinations;" % (c.n_ok.min(), c.n_ok.max(), len(c)),
"verdicts", c.verdict.value_counts().to_dict())
print(c.groupby(["tmax", "width"]).n_ok.agg(["min", "max"]).to_string())
print("values under the first reading of each open detail (measured, expected, tolerance, within):")
for row in first_reading(t):
print(" ", row)
L5: verdicts over the reading combinations: {'partially': 216, 'not reproduced': 40}
min max
start lat_reading
cross a 10 13
b 10 13
first a 8 11
b 8 11
verdicts by start of the cutoff search (cross is the criterion's, first a variant) and whisking reading:
min max
start whisk_with_touch latency_from
cross False lat_all 11 13
lat_first 11 13
True lat_all 10 12
lat_first 10 12
first False lat_all 9 11
lat_first 9 11
True lat_all 8 10
lat_first 8 10
verdict not reproduced partially
start
cross 0 128
first 40 88
under the criterion's readings (cross, whisking outside touch epochs): values within tolerance 11 to 13 over 64 combinations; verdicts {'partially': 64}
min max
tmax width
100 pm 12 13
total 11 12
200 pm 12 13
total 11 13
values under the first reading of each open detail (measured, expected, tolerance, within):
('spt_first_mean', 2.046, 1.95, 2.96, True)
('spt_first_median', 0.227, 0.98, 0.098, False)
('spt_later_mean', 1.313, 1.28, 1.88, True)
('spt_later_median', 0.188, 0.42, 0.042, False)
('touch_pct_mean', 27.202, 24.0, 18.6, True)
('touch_pct_median', 22.115, 24.4, 2.44, True)
('tp_pct_mean', 45.029, 31.5, 16.4, True)
('tp_pct_median', 38.078, 36.0, 3.6, True)
('depth_mean', 0.182, 0.18, 0.12, True)
('depth_median', 0.106, 0.13, 0.013, False)
('nonwhisk_mean', 12.039, 12.0, 9.87, True)
('nonwhisk_median', 12.903, 12.8, 1.28, True)
('whisk_mean', 16.66, 16.2, 15.3, True)
('whisk_median', 12.819, 13.5, 1.35, True)
('latency_mean', 4.9, 9.7, 5.2, True)
('latency_median', 4.0, 8, 1, False)
('min_isi_mean', 2.627, 3.7, 2.5, True)
('min_isi_median', 1.5, 2.7, 0.27, False)
Result: partially reproduced. Under the criterion's readings 11 to 13 of the 18 values are within tolerance in every combination, partially reproduced throughout. Rates, mean depth (0.182), the median whisking rate (12.82 against 13.5 +- 1.35), the means of spikes per touch and the touch-evoked fraction (27.2%, median 22.1% under the first reading) are within tolerance. The medians of spikes per first and later touch (0.227 and 0.188 against 0.98 and 0.42), the median depth (0.106 against 0.13 +- 0.013), the median latency (4 against 8 ms) and the median minimum ISI (1.5 against 2.7 ms) are not. The literal cutoff, a variant, gives 8 to 11.
latency-determinable¶
A touch-response latency could be determined for all 31 L4 C2 cells but for only four of the 10 L4 cells outside C2, for lack of a touch response.
Readings of the PCTH: first touches or all touches. The search runs 0 to 50 ms.
for src in ("lat_first", "lat_all"):
has = base[src].notna().groupby(base.group).sum()
print(f"{src}: cells with a latency {has.to_dict()}; within (31, 4): {has['L4C2'] == 31 and has['L4out'] == 4}")
print(base.loc[G["L4out"] + G["L5"], ["group", "lat_first", "lat_all"]].to_string())
lat_first: cells with a latency {'L4C2': 31, 'L4out': 7, 'L5': 10}; within (31, 4): False
lat_all: cells with a latency {'L4C2': 31, 'L4out': 10, 'L5': 9}; within (31, 4): False
group lat_first lat_all
cell
32 L4out NaN 23.0
33 L4out 16.0 16.0
34 L4out NaN 16.0
35 L4out 20.0 3.0
36 L4out 5.0 0.0
37 L4out 3.0 0.0
38 L4out 4.0 0.0
39 L4out 8.0 7.0
40 L4out 1.0 33.0
41 L4out NaN 2.0
42 L5 3.0 12.0
43 L5 11.0 15.0
44 L5 10.0 10.0
45 L5 4.0 9.0
46 L5 0.0 9.0
47 L5 4.0 24.0
48 L5 4.0 7.0
49 L5 3.0 NaN
50 L5 0.0 0.0
51 L5 NaN NaN
52 L5 10.0 12.0
Result: not reproduced. A latency is found for all 31 L4 C2 cells under both readings, but for 7 of the 10 L4 cells outside C2 from first touches and 10 of 10 from all touches, not 4. The two readings agree that the count is not the paper's.
l4-outside-weaker-touch¶
L4 neurons outside C2 respond much more weakly to touch than L4 neurons in C2, for first and for later touches.
Readings of the measure: (a) spikes per touch in the touch window (every window reading), or (b) the peak touch rate, the highest 1 ms bin of the first-touch PCTH and of the later-touch PCTH, 0 to 50 ms after onset.
rows = []
for key, g in W[W.group.isin(["L4C2", "L4out"])].groupby(WKEY):
a, b = g[g.group == "L4C2"], g[g.group == "L4out"]
pf, pl = ranksum(a.spt_first, b.spt_first), ranksum(a.spt_later, b.spt_later)
rows.append(dict(zip(WKEY, key), med_c2_first=a.spt_first.median(), med_out_first=b.spt_first.median(), p_first=pf,
med_c2_later=a.spt_later.median(), med_out_later=b.spt_later.median(), p_later=pl,
both=(pf < 0.05) and (pl < 0.05) and a.spt_first.median() > b.spt_first.median() and a.spt_later.median() > b.spt_later.median()))
LW = pd.DataFrame(rows)
print("(a) spikes per touch: combinations with both tests p < 0.05 and C2 higher: %d of %d" % (LW.both.sum(), len(LW)))
print(LW.groupby("start").both.agg(["sum", "size"]).to_string())
print(LW.groupby("start")[["p_first", "p_later"]].agg(["min", "max"]).to_string(float_format="{:.3g}".format))
peak = {}
for no in G["L4C2"] + G["L4out"]:
c = by_no[no]
T = D[no]["T"]
s = D[no]["s"]
pk = []
for sel in (T.order == 1, T.order > 1):
t_, h, n = hires.pcth(c, s, T.on[sel])
pk.append(h[(t_ >= 0)].max() * 1000)
peak[no] = pk
peak = pd.DataFrame(peak, index=["peak_first", "peak_later"]).T
pa, pb = peak.loc[G["L4C2"]], peak.loc[G["L4out"]]
print("(b) peak touch rate (spikes/s): first touches C2 median %.1f, outside %.1f, p = %.3g; later touches C2 %.1f, outside %.1f, p = %.3g" %
(pa.peak_first.median(), pb.peak_first.median(), ranksum(pa.peak_first, pb.peak_first),
pa.peak_later.median(), pb.peak_later.median(), ranksum(pa.peak_later, pb.peak_later)))
(a) spikes per touch: combinations with both tests p < 0.05 and C2 higher: 32 of 32
sum size
start
cross 16 16
first 16 16
p_first p_later
min max min max
start
cross 0.000102 0.00432 0.0002 0.0152
first 1.63e-05 0.00129 1.74e-05 0.00987
(b) peak touch rate (spikes/s): first touches C2 median 148.1, outside 25.0, p = 1.63e-05; later touches C2 51.3, outside 5.6, p = 1.07e-05
Result: reproduced. Under reading (a) both rank-sum tests give p < 0.05 with the C2 median higher in all 16 window readings of the transition cutoff (first touches p = 1.0e-4 to 0.0043, later touches 2.0e-4 to 0.015; the paper gives 1.4e-5 and 6.32e-6), and in the 16 of the literal variant. Under reading (b) the peak touch rate is 148.1 against 25.0 spikes/s for first touches (p = 1.63e-5) and 51.3 against 5.6 for later touches (p = 1.07e-5).
l5-higher-rates¶
L5 neurons near C2 have much higher spike rates than L4 neurons, with more diverse touch responses.
for wt, col in ((False, "rate_whisk_no_touch"), (True, "rate_whisk_with_touch")):
print(f"whisking read {'with' if wt else 'without'} touch epochs")
allok, c2ok = True, True
for rate in ("rate_nonwhisk", col):
m = {g: base.loc[G[g], rate].median() for g in GROUP_ORDER}
for other in ("L4C2", "L4out"):
p = ranksum(base.loc[G["L5"], rate], base.loc[G[other], rate])
ok = (m["L5"] > m[other]) and p < 0.05
allok &= ok
if other == "L4C2":
c2ok &= ok
print(f" {rate}: L5 median {m['L5']:.2f} vs {other} {m[other]:.2f} spikes/s, rank sum p = {p:.3g}")
print(" all four comparisons hold:", allok, "| those against L4 C2 hold:", c2ok)
fig, axs = plt.subplots(1, 3, figsize=(10, 2.6), sharey=False)
for ax, g in zip(axs, GROUP_ORDER):
for no in G[g]:
T = D[no]["T"]
t_, h, n = hires.pcth(by_no[no], D[no]["s"], T.on[T.order == 1])
ax.plot(t_, h * 1000, lw=0.7)
ax.set_title(g); ax.set_xlabel("ms from first touch")
axs[0].set_ylabel("spikes/s")
plt.tight_layout(); plt.show()
print("plotted: first-touch PCTHs, 1 ms bins, -50 to 50 ms; all values inside the axes")
whisking read without touch epochs rate_nonwhisk: L5 median 12.90 vs L4C2 0.32 spikes/s, rank sum p = 0.000119 rate_nonwhisk: L5 median 12.90 vs L4out 1.09 spikes/s, rank sum p = 0.00918 rate_whisk_no_touch: L5 median 12.82 vs L4C2 0.38 spikes/s, rank sum p = 2.15e-05 rate_whisk_no_touch: L5 median 12.82 vs L4out 1.46 spikes/s, rank sum p = 0.00485 all four comparisons hold: True | those against L4 C2 hold: True whisking read with touch epochs rate_nonwhisk: L5 median 12.90 vs L4C2 0.32 spikes/s, rank sum p = 0.000119 rate_nonwhisk: L5 median 12.90 vs L4out 1.09 spikes/s, rank sum p = 0.00918 rate_whisk_with_touch: L5 median 11.52 vs L4C2 2.27 spikes/s, rank sum p = 0.00157 rate_whisk_with_touch: L5 median 11.52 vs L4out 1.54 spikes/s, rank sum p = 0.00485 all four comparisons hold: True | those against L4 C2 hold: True
plotted: first-touch PCTHs, 1 ms bins, -50 to 50 ms; all values inside the axes
Result: reproduced. L5 medians are 12.90 spikes/s outside whisking and 12.82 (11.52 with touch epochs) during it, against 0.32 and 0.38 (2.27) for L4 C2 and 1.09 and 1.46 (1.54) for L4 outside C2; all four rank-sum tests give p < 0.05 under both readings (largest p = 0.0092). The diversity of touch responses is not judged; the first-touch PCTHs are plotted.
phase-mod-l4-inside-outside and phase-mod-l4-vs-l5¶
Modulation by whisking phase does not differ between L4 neurons inside and outside C2.
L4 neurons, inside and outside C2, are more modulated by whisking phase than L5 neurons.
Rank-sum tests of the per-cell modulation depth, from every phase-analysis cycle of the trial; cells without spikes in those cycles have no depth and are left out (their numbers are printed). Variants: only cells with at least 20 phase spikes, and the cycles counted only in exploration time (the earlier run's).
dm = {g: base.loc[G[g], "depth_mod"] for g in GROUP_ORDER}
print({g: f"n {v.notna().sum()} of {len(v)}, median {v.median():.3f}, mean {v.mean():.3f}" for g, v in dm.items()})
p = ranksum(dm["L4C2"], dm["L4out"])
print(f"L4 C2 vs L4 outside: rank sum p = {p:.3g}; absence holds at 0.05: {p >= 0.05}")
print("cells with depth 1 (no fitted baseline):", {g: int((v > 0.999).sum()) for g, v in dm.items()})
print("phase-analysis spikes per cell, median:", {g: base.loc[G[g], "phase_spikes"].median() for g in GROUP_ORDER})
for name, col, keep in (("cells with >= 20 phase spikes (variant)", "depth_mod", base.phase_spikes >= 20),
("exploration time only (variant)", "depth_mod_x", base.phase_spikes_x > 0)):
a_, b_ = base.loc[G["L4C2"]], base.loc[G["L4out"]]
a_, b_ = a_.loc[keep[a_.index], col], b_.loc[keep[b_.index], col]
print(f"{name}: L4 C2 n {a_.notna().sum()} median {a_.median():.3f}, L4 outside n {b_.notna().sum()} median {b_.median():.3f}, rank sum p = {ranksum(a_, b_):.3g}")
{'L4C2': 'n 27 of 31, median 0.859, mean 0.685', 'L4out': 'n 10 of 10, median 0.548, mean 0.604', 'L5': 'n 11 of 11, median 0.106, mean 0.182'}
L4 C2 vs L4 outside: rank sum p = 0.584; absence holds at 0.05: True
cells with depth 1 (no fitted baseline): {'L4C2': 12, 'L4out': 4, 'L5': 0}
phase-analysis spikes per cell, median: {'L4C2': np.float64(14.0), 'L4out': np.float64(83.0), 'L5': np.float64(923.0)}
cells with >= 20 phase spikes (variant): L4 C2 n 11 median 0.444, L4 outside n 8 median 0.496, rank sum p = 0.934
exploration time only (variant): L4 C2 n 25 median 1.000, L4 outside n 10 median 0.980, rank sum p = 0.971
Result: reproduced. The rank-sum test of modulation depth between L4 C2 (27 cells with a fit, median 0.859) and L4 outside C2 (10 cells, median 0.548) gives p = 0.584, so the paper's test finds no difference (paper p = 0.68). With only cells that have at least 20 phase spikes (11 and 8 cells, medians 0.444 and 0.496) p = 0.934, and with the cycles counted only in exploration time p = 0.798.
for g in ("L4C2", "L4out"):
p = ranksum(dm[g], dm["L5"])
print(f"{g} vs L5: medians {dm[g].median():.3f} vs {dm['L5'].median():.3f}, rank sum p = {p:.3g}, holds: {p < 0.05 and dm[g].median() > dm['L5'].median()}")
# variant: only cells with at least 20 phase-analysis spikes
sel = base.phase_spikes >= 20
print("variant, cells with >= 20 phase spikes:", {g: int(sel[G[g]].sum()) for g in GROUP_ORDER})
for a_, b_ in (("L4C2", "L4out"), ("L4C2", "L5"), ("L4out", "L5")):
x, y = dm[a_][sel[G[a_]]], dm[b_][sel[G[b_]]]
print(f" {a_} {x.median():.3f} vs {b_} {y.median():.3f}: p = {ranksum(x, y):.3g}")
for name, col, keep in (("exploration time only (variant)", "depth_mod_x", base.phase_spikes_x > 0),):
l5 = base.loc[G["L5"]]
l5 = l5.loc[keep[l5.index], col]
for g in ("L4C2", "L4out"):
a_ = base.loc[G[g]]
a_ = a_.loc[keep[a_.index], col]
print(f"{name}: {g} n {a_.notna().sum()} median {a_.median():.3f} vs L5 n {l5.notna().sum()} median {l5.median():.3f}, rank sum p = {ranksum(a_, l5):.3g}")
L4C2 vs L5: medians 0.859 vs 0.106, rank sum p = 8.05e-05, holds: True
L4out vs L5: medians 0.548 vs 0.106, rank sum p = 0.00918, holds: True
variant, cells with >= 20 phase spikes: {'L4C2': 11, 'L4out': 8, 'L5': 11}
L4C2 0.444 vs L4out 0.496: p = 0.934
L4C2 0.444 vs L5 0.106: p = 0.00428
L4out 0.496 vs L5 0.106: p = 0.0318
exploration time only (variant): L4C2 n 25 median 1.000 vs L5 n 11 median 0.141, rank sum p = 0.000291
exploration time only (variant): L4out n 10 median 0.980 vs L5 n 11 median 0.141, rank sum p = 0.0031
Result: reproduced. L4 C2 (median 0.859) and L4 outside C2 (0.548) are both more phase modulated than L5 (0.106): p = 8.05e-5 and 0.00918 (paper 4.3e-5 and 0.012). With only cells that have at least 20 phase spikes (11, 8 and 11 cells), p = 0.00428 and 0.0318; with the cycles counted only in exploration time, p = 2.9e-4 and 0.0031.
whisking-rate-l4-no-change and whisking-rate-l5-increase¶
Across L4 neurons, inside and outside C2, spike rate does not change between whisking and non-whisking periods.
L5 neurons fire modestly but significantly more during whisking than during non-whisking periods.
Wilcoxon signed-rank test of per-cell whisking against non-whisking rate, each group apart,
whisking read with and without touch epochs. Whisking leaves out licking cycles and the 70 ms
after each touch (exclusions True); the earlier run's mask without them is shown as a
variant.
WR = []
for g in GROUP_ORDER:
for ex in (True, False):
for wt in (False, True):
col = ("rate_whisk_with_touch" if wt else "rate_whisk_no_touch") + ("" if ex else "_noexcl")
b = base.loc[G[g]]
diff = b[col] - b.rate_nonwhisk
p = stats.wilcoxon(b[col], b.rate_nonwhisk).pvalue
WR.append(dict(group=g, exclusions=ex, whisking_with_touch=wt, median_whisk=b[col].median(), median_nonwhisk=b.rate_nonwhisk.median(),
median_difference=diff.median(), n_higher=int((diff > 0).sum()), n=len(b), p=p))
WR = pd.DataFrame(WR)
print("exclusions True: licking cycles and the 70 ms after touch left out of whisking (the criterion's); False: the earlier run's mask, a variant")
print(WR.to_string(index=False, float_format="{:.3g}".format))
print("per-cell rates, spikes/s (whisking outside touch epochs, with the exclusions):")
print(base[["group", "rate_nonwhisk", "rate_whisk_no_touch", "rate_whisk_with_touch", "rate_whisk_no_touch_noexcl"]].round(3).to_string())
exclusions True: licking cycles and the 70 ms after touch left out of whisking (the criterion's); False: the earlier run's mask, a variant
group exclusions whisking_with_touch median_whisk median_nonwhisk median_difference n_higher n p
L4C2 True False 0.376 0.315 0.0741 18 31 0.595
L4C2 True True 2.27 0.315 1.18 29 31 1.28e-07
L4C2 False False 0.86 0.315 0.633 26 31 5.02e-05
L4C2 False True 2.28 0.315 1.21 30 31 2.79e-09
L4out True False 1.46 1.09 0.138 6 10 0.625
L4out True True 1.54 1.09 0.235 6 10 0.695
L4out False False 1.75 1.09 0.419 6 10 0.432
L4out False True 1.78 1.09 0.365 6 10 0.432
L5 True False 12.8 12.9 2.1 9 11 0.0186
L5 True True 11.5 12.9 2.43 9 11 0.0322
L5 False False 11 12.9 1.88 9 11 0.0244
L5 False True 11.3 12.9 1.79 10 11 0.0137
per-cell rates, spikes/s (whisking outside touch epochs, with the exclusions):
group rate_nonwhisk rate_whisk_no_touch rate_whisk_with_touch rate_whisk_no_touch_noexcl
cell
1 L4C2 0.076 0.339 2.270 0.778
2 L4C2 1.379 3.123 4.836 3.219
3 L4C2 0.546 2.656 12.626 3.444
4 L4C2 0.385 0.033 1.490 0.148
5 L4C2 0.171 0.024 1.601 0.133
6 L4C2 0.275 0.349 0.646 0.527
7 L4C2 0.022 0.297 0.701 0.262
8 L4C2 1.130 0.376 2.307 0.886
9 L4C2 13.677 10.625 17.120 15.096
10 L4C2 0.719 0.157 2.345 2.067
11 L4C2 2.144 2.294 3.193 2.358
12 L4C2 0.149 0.241 0.583 0.437
13 L4C2 0.337 4.102 11.878 4.767
14 L4C2 0.148 0.201 3.778 0.781
15 L4C2 1.333 0.617 0.949 0.823
16 L4C2 1.065 6.121 6.507 5.353
17 L4C2 0.611 1.596 2.800 3.896
18 L4C2 8.706 4.723 8.016 6.579
19 L4C2 0.075 0.037 0.523 0.228
20 L4C2 0.276 0.175 0.896 0.561
21 L4C2 0.632 0.639 5.067 3.718
22 L4C2 0.043 0.161 2.954 0.291
23 L4C2 0.227 0.777 2.195 1.257
24 L4C2 0.123 0.388 0.569 0.757
25 L4C2 0.192 0.514 3.470 1.897
26 L4C2 0.164 0.016 0.820 0.359
27 L4C2 0.487 0.089 1.372 0.660
28 L4C2 0.103 0.638 0.688 0.860
29 L4C2 0.315 0.138 0.897 0.554
30 L4C2 0.050 0.247 0.400 1.935
31 L4C2 4.838 3.642 11.136 9.841
32 L4out 0.025 0.641 0.470 0.545
33 L4out 0.040 0.187 0.218 0.375
34 L4out 8.922 12.138 11.435 12.300
35 L4out 0.533 0.268 0.286 0.285
36 L4out 0.853 1.465 1.464 2.019
37 L4out 0.266 1.827 1.777 1.658
38 L4out 1.413 0.859 0.915 0.851
39 L4out 5.584 4.067 4.033 3.758
40 L4out 1.334 1.462 1.625 1.837
41 L4out 2.859 2.095 2.024 2.476
42 L5 29.846 59.310 69.098 69.686
43 L5 1.869 11.215 10.901 10.507
44 L5 23.101 26.976 28.320 26.972
45 L5 1.029 1.204 1.540 1.329
46 L5 12.971 15.067 15.401 14.323
47 L5 15.858 17.009 16.332 17.743
48 L5 21.836 24.623 45.087 34.934
49 L5 5.689 9.497 11.522 9.549
50 L5 0.240 0.978 0.741 0.862
51 L5 12.903 12.819 11.357 10.962
52 L5 7.082 4.565 6.326 6.532
Result: blocked (methods_underspecified). The two readings of whisking periods give different outcomes. With touch epochs left out, neither group changes: L4 C2 p = 0.595 (18 of 31 cells higher when whisking, median difference +0.074 spikes/s) and L4 outside C2 p = 0.625, so the claim is reproduced (paper 0.75 and 0.92). With touch epochs in, L4 C2 changes (p = 1.28e-7, 29 of 31 higher, +1.18 spikes/s) and L4 outside C2 does not (p = 0.695), which is partially reproduced. The earlier mask without the licking and after-touch exclusions, a variant, gives a change in L4 C2 under both readings (p = 5.0e-5 and 2.8e-9).
print(WR[(WR.group == "L5")].to_string(index=False, float_format="{:.3g}".format))
for r in WR[WR.group == "L5"].itertuples():
print(f"exclusions {r.exclusions}, whisking with touch epochs {r.whisking_with_touch}: p < 0.05 {r.p < 0.05}, median of the per-cell differences positive {r.median_difference > 0}")
group exclusions whisking_with_touch median_whisk median_nonwhisk median_difference n_higher n p L5 True False 12.8 12.9 2.1 9 11 0.0186 L5 True True 11.5 12.9 2.43 9 11 0.0322 L5 False False 11 12.9 1.88 9 11 0.0244 L5 False True 11.3 12.9 1.79 10 11 0.0137 exclusions True, whisking with touch epochs False: p < 0.05 True, median of the per-cell differences positive True exclusions True, whisking with touch epochs True: p < 0.05 True, median of the per-cell differences positive True exclusions False, whisking with touch epochs False: p < 0.05 True, median of the per-cell differences positive True exclusions False, whisking with touch epochs True: p < 0.05 True, median of the per-cell differences positive True
Result: reproduced. The signed-rank test is significant under both readings of whisking periods (p = 0.0186 without touch epochs and 0.0322 with them, paper 0.019), and the median of the per-cell differences, the measure the answered criterion names, is positive under both (+2.10 and +2.43 spikes/s), with 9 of the 11 cells firing more when whisking under each. The group medians, shown only for comparison, are 12.8 and 11.5 spikes/s when whisking against 12.9 outside whisking (Table 1 gives 13.5 and 12.8). The earlier mask, a variant, also gives p < 0.05 under both readings (0.0244 and 0.0137).
Figure 3C: fig3c-touch-spikes-l4c2, fig3c-touch-spikes-outside, fig3c-touch-spikes-l5¶
In L4 C2 neurons 70.6 +- 20.9% (mean +- s.d.) of spikes during exploration are coupled to touch, within an average of 14.9% of exploration time.
L4 neurons outside C2 have far fewer touch-coupled spikes, 8.4 +- 13.3%, in 3.5% of exploration time.
L5 neurons have fewer touch-coupled spikes than L4 C2 neurons, 24.0 +- 18.6%, in 19.0% of exploration time.
Every touch-window reading is shown. For the outside and L5 groups, the rank-sum test against L4 C2 uses the same reading for both groups; the variant sets cells without a first-touch latency to 0%.
F3 = []
c2 = W[W.group == "L4C2"]
for g, exp in (("L4C2", (70.6, 20.9, 14.9, 1.49)), ("L4out", (8.4, 13.3, 3.5, 0.35)), ("L5", (24.0, 18.6, 19.0, 1.9))):
for key, kk in W[W.group == g].groupby(WKEY):
ref = c2[(c2.lat_reading == "a") & (c2.shuffle == key[1]) & (c2.nearest == key[2]) & (c2.tmax == key[3]) & (c2.start == key[4])]
fs = 100 * kk.frac_spikes
nolat = base.loc[kk.cell, "lat_first"].isna().to_numpy()
fs0 = np.where(nolat, 0.0, fs)
p = ranksum(ref.frac_spikes, kk.frac_spikes) if g != "L4C2" else np.nan
F3.append(dict(group=g, **dict(zip(WKEY, key)), touch_mean=fs.mean(), touch_sd=fs.std(ddof=1), touch_median=fs.median(),
time_mean=100 * kk.frac_time.mean(), p_vs_c2=p, lower_than_c2=fs.median() < 100 * ref.frac_spikes.median(),
ok_values=within(fs.mean(), exp[0], exp[1]) and within(100 * kk.frac_time.mean(), exp[2], exp[3]),
variant_zero_mean=fs0.mean()))
F3 = pd.DataFrame(F3)
F3["ok"] = F3.ok_values & ((F3.group == "L4C2") | ((F3.p_vs_c2 < 0.05) & F3.lower_than_c2))
t = F3[(F3.group == "L4C2") & (F3.lat_reading == "a")]
print(t.drop(columns=["group", "lat_reading", "p_vs_c2", "lower_than_c2", "variant_zero_mean"]).round(2).to_string(index=False))
print("readings within both tolerances: %d of %d" % (t.ok.sum(), len(t)))
print("by start and longest window (cross is the criterion's start, first a variant):")
print(t.groupby(["start", "tmax"]).agg(ok=("ok", "sum"), n=("ok", "size"), touch_min=("touch_mean", "min"), touch_max=("touch_mean", "max"),
time_min=("time_mean", "min"), time_max=("time_mean", "max")).round(2).to_string())
fig, ax = plt.subplots(figsize=(4, 3))
k = W[(W.group == "L4C2") & (W.lat_reading == "a") & (W.shuffle == "uniform") & (W.nearest == "before") & (W.tmax == 100) & (W.start == "cross")]
for g, col in (("L4C2", "red"), ("L4out", "grey"), ("L5", "orange")):
kk = W[(W.group == g) & (W.lat_reading == "a") & (W.shuffle == "uniform") & (W.nearest == "before") & (W.tmax == 100) & (W.start == "cross")]
ax.scatter(100 * kk.frac_time, 100 * kk.frac_spikes, s=12, c=col, label=g)
ax.set_xlabel("% of exploration time"); ax.set_ylabel("% of spikes in touch window"); ax.legend(fontsize=7)
ax.set_title("uniform, before, 100 ms, cross", fontsize=8)
plt.tight_layout(); plt.show()
kk = W[(W.lat_reading == "a") & (W.shuffle == "uniform") & (W.nearest == "before") & (W.tmax == 100) & (W.start == "cross")]
print("plotted ranges: time %.1f to %.1f %%, spikes %.1f to %.1f %%" % (100 * kk.frac_time.min(), 100 * kk.frac_time.max(), 100 * kk.frac_spikes.min(), 100 * kk.frac_spikes.max()))
shuffle nearest tmax start touch_mean touch_sd touch_median time_mean ok_values ok
circular before 100 cross 70.86 17.36 77.48 14.60 True True
circular before 100 first 40.77 37.06 55.10 8.53 False False
circular before 200 cross 73.89 12.99 78.09 16.01 True True
circular before 200 first 61.37 32.34 76.90 11.89 False False
circular either 100 cross 66.88 15.46 68.00 14.92 True True
circular either 100 first 38.62 34.63 47.33 8.89 False False
circular either 200 cross 68.86 14.16 71.85 16.50 False False
circular either 200 first 59.24 29.21 71.85 12.60 False False
uniform before 100 cross 73.81 14.89 79.98 16.27 True True
uniform before 100 first 42.84 37.69 57.63 9.51 False False
uniform before 200 cross 75.80 13.20 80.34 17.48 False False
uniform before 200 first 64.47 31.00 80.34 13.34 False False
uniform either 100 cross 68.56 15.64 71.31 16.00 True True
uniform either 100 first 40.60 34.75 47.92 9.60 False False
uniform either 200 cross 70.28 14.30 71.85 17.58 False False
uniform either 200 first 60.37 29.72 71.85 13.36 False False
readings within both tolerances: 5 of 16
by start and longest window (cross is the criterion's start, first a variant):
ok n touch_min touch_max time_min time_max
start tmax
cross 100 4 4 66.88 73.81 14.60 16.27
200 1 4 68.86 75.80 16.01 17.58
first 100 0 4 38.62 42.84 8.53 9.60
200 0 4 59.24 64.47 11.89 13.36
plotted ranges: time 0.9 to 49.0 %, spikes 0.0 to 94.3 %
Result: blocked (methods_underspecified). Under the transition cutoff, the criterion's, 66.9 to 75.8% of exploration spikes fall in touch windows covering 14.6 to 17.6% of exploration time. All 4 readings with the 100 ms curve put both values in tolerance (70.6 +- 20.9% and 14.9 +- 1.49%), but only 1 of the 4 with the 200 ms curve does: the other three cover 16.5 to 17.6% of exploration time. The readings of the curve's range, which the paper does not give, disagree. The literal cutoff, a variant, gives 38.6 to 64.5% in 8.5 to 13.4% of time.
t = F3[F3.group == "L4out"]
print(t.drop(columns=["group"]).round(3).to_string(index=False, formatters={"p_vs_c2": "{:.3g}".format}))
print("readings meeting the criterion: %d of %d" % (t.ok.sum(), len(t)))
print(t.groupby(["start", "lat_reading", "tmax"]).agg(ok=("ok", "sum"), n=("ok", "size"), touch_min=("touch_mean", "min"), touch_max=("touch_mean", "max"),
time_min=("time_mean", "min"), time_max=("time_mean", "max"), p_max=("p_vs_c2", "max")).round(3).to_string())
print("per cell, first reading of each detail (uniform, before, 100 ms, cross, latency reading a): latency, t_on, fractions, and whether the increment rose above chance:")
print(W[(W.group == "L4out") & (W.lat_reading == "a") & (W.shuffle == "uniform") & (W.nearest == "before") & (W.tmax == 100) & (W.start == "cross")]
[["cell", "latency", "t_on", "frac_spikes", "frac_time", "crossed", "n_explore_spikes"]].assign(has_first_latency=lambda x: base.loc[x.cell, "lat_first"].notna().to_numpy()).round(3).to_string(index=False))
lat_reading shuffle nearest tmax start touch_mean touch_sd touch_median time_mean p_vs_c2 lower_than_c2 ok_values variant_zero_mean ok
a circular before 100 cross 28.186 22.421 29.245 21.209 0 True False 27.014 False
a circular before 100 first 1.526 1.261 1.148 1.250 0.069 True False 1.382 False
a circular before 200 cross 28.827 26.705 25.921 21.673 0 True False 22.300 False
a circular before 200 first 1.758 1.275 1.764 1.564 0.001 True False 1.613 False
a circular either 100 cross 20.540 16.630 17.785 19.713 0 True False 19.521 False
a circular either 100 first 1.526 1.261 1.148 1.250 0.069 True False 1.382 False
a circular either 200 cross 18.772 16.689 15.122 18.949 0 True False 15.544 False
a circular either 200 first 6.894 15.365 1.569 2.616 0 True False 6.749 False
a uniform before 100 cross 26.167 23.694 23.759 18.978 0 True False 23.642 False
a uniform before 100 first 1.526 1.261 1.148 1.250 0.069 True False 1.382 False
a uniform before 200 cross 31.712 25.259 25.517 26.286 0 True False 23.224 False
a uniform before 200 first 1.763 1.475 1.148 1.403 0 True False 1.618 False
a uniform either 100 cross 19.269 17.186 16.824 17.547 0 True False 17.265 False
a uniform either 100 first 1.526 1.261 1.148 1.250 0.052 True False 1.382 False
a uniform either 200 cross 24.746 15.052 21.743 26.001 0 True False 17.831 False
a uniform either 200 first 6.702 15.435 1.148 2.424 0 True False 6.557 False
b circular before 100 cross 29.665 21.069 29.245 22.783 0 True False 27.014 False
b circular before 100 first 2.010 1.230 1.443 1.338 0.249 True False 1.382 False
b circular before 200 cross 30.302 25.157 30.716 23.231 0 True False 22.300 False
b circular before 200 first 2.242 1.140 1.959 1.538 0.002 True False 1.613 False
b circular either 100 cross 21.142 15.957 17.785 20.943 0 True False 19.521 False
b circular either 100 first 2.181 1.271 1.770 1.465 0.249 True False 1.382 False
b circular either 200 cross 19.550 16.506 18.103 19.435 0 True False 15.544 False
b circular either 200 first 7.361 15.181 1.770 2.590 0.001 True False 6.749 False
b uniform before 100 cross 26.318 23.360 22.873 18.983 0 True False 23.642 False
b uniform before 100 first 2.010 1.230 1.443 1.224 0.249 True False 1.382 False
b uniform before 200 cross 31.386 24.825 23.203 25.959 0 True False 23.224 False
b uniform before 200 first 2.247 1.358 1.638 1.376 0.001 True False 1.618 False
b uniform either 100 cross 19.232 16.840 14.529 17.529 0 True False 17.265 False
b uniform either 100 first 2.181 1.271 1.770 1.465 0.182 True False 1.382 False
b uniform either 200 cross 23.159 15.672 18.040 24.664 0 True False 17.831 False
b uniform either 200 first 7.169 15.259 1.443 2.398 0.001 True False 6.557 False
readings meeting the criterion: 0 of 32
ok n touch_min touch_max time_min time_max p_max
start lat_reading tmax
cross a 100 0 4 19.269 28.186 17.547 21.209 0.000
200 0 4 18.772 31.712 18.949 26.286 0.000
b 100 0 4 19.232 29.665 17.529 22.783 0.000
200 0 4 19.550 31.386 19.435 25.959 0.000
first a 100 0 4 1.526 1.526 1.250 1.250 0.069
200 0 4 1.758 6.894 1.403 2.616 0.001
b 100 0 4 2.010 2.181 1.224 1.465 0.249
200 0 4 2.242 7.361 1.376 2.590 0.002
per cell, first reading of each detail (uniform, before, 100 ms, cross, latency reading a): latency, t_on, fractions, and whether the increment rose above chance:
cell latency t_on frac_spikes frac_time crossed n_explore_spikes has_first_latency
32 0.0 1 0.000 0.014 False 52 False
33 16.0 24 0.593 0.095 True 54 True
34 0.0 72 0.215 0.234 True 1242 False
35 20.0 5 0.091 0.023 True 55 True
36 5.0 86 0.609 0.490 True 256 True
37 3.0 51 0.260 0.317 True 196 True
38 4.0 1 0.007 0.011 False 134 True
39 8.0 53 0.300 0.354 True 407 True
40 1.0 56 0.504 0.328 True 127 True
41 0.0 4 0.037 0.032 True 534 False
Result: blocked (methods_underspecified). The claim depends on fig3c-touch-spikes-l4c2, which is blocked. Under every one of its 16 readings with the transition cutoff the time fraction is outside tolerance: 18.8 to 31.7% of spikes in 17.5 to 26.3% of exploration time (paper 8.4 +- 13.3% in 3.5%), with the rank-sum test against L4 C2 at p < 0.001. Per cell, under the first reading, five cells (34, 36, 37, 39, 40) have windows of 51 to 86 ms covering 23 to 49% of exploration time; three of them (36, 37, 40) have latencies of 1 to 5 ms, and cell 34 has no first-touch latency. The literal cutoff, a variant, gives 1.5 to 7.4% in 1.2 to 2.6%.
t = F3[F3.group == "L5"]
print(t.drop(columns=["group"]).round(3).to_string(index=False, formatters={"p_vs_c2": "{:.3g}".format}))
print("readings meeting the criterion: %d of %d" % (t.ok.sum(), len(t)))
print(t.groupby(["start", "lat_reading", "tmax"]).agg(ok=("ok", "sum"), n=("ok", "size"), touch_min=("touch_mean", "min"), touch_max=("touch_mean", "max"),
time_min=("time_mean", "min"), time_max=("time_mean", "max"), p_max=("p_vs_c2", "max")).round(3).to_string())
print("per cell, first reading of each detail (uniform, before, 100 ms, cross, latency reading a): latency, t_on, fractions, and whether the increment rose above chance:")
print(W[(W.group == "L5") & (W.lat_reading == "a") & (W.shuffle == "uniform") & (W.nearest == "before") & (W.tmax == 100) & (W.start == "cross")]
[["cell", "latency", "t_on", "frac_spikes", "frac_time", "crossed", "n_explore_spikes"]].assign(has_first_latency=lambda x: base.loc[x.cell, "lat_first"].notna().to_numpy()).round(3).to_string(index=False))
lat_reading shuffle nearest tmax start touch_mean touch_sd touch_median time_mean p_vs_c2 lower_than_c2 ok_values variant_zero_mean ok
a circular before 100 cross 25.651 21.797 20.673 17.784 0 True True 25.597 True
a circular before 100 first 6.143 15.006 1.259 3.542 0.141 True False 6.089 False
a circular before 200 cross 27.514 22.468 26.437 19.637 0 True True 27.359 True
a circular before 200 first 13.674 23.388 1.706 7.696 0.003 True False 13.519 False
a circular either 100 cross 20.770 16.408 18.391 16.128 0 True False 20.716 False
a circular either 100 first 5.099 11.557 1.259 3.154 0.126 True False 5.045 False
a circular either 200 cross 24.278 18.690 20.305 20.888 0 True True 24.123 True
a circular either 200 first 11.932 19.565 1.706 8.163 0.001 True False 11.776 False
a uniform before 100 cross 27.202 23.027 22.115 19.163 0 True True 27.147 True
a uniform before 100 first 6.318 15.222 1.259 3.696 0.133 True False 6.264 False
a uniform before 200 cross 30.570 21.880 29.054 22.466 0 True False 27.925 False
a uniform before 200 first 14.081 24.569 1.525 8.067 0.001 True False 14.027 False
a uniform either 100 cross 21.358 18.285 19.595 15.400 0 True False 21.303 False
a uniform either 100 first 5.559 12.721 1.259 3.696 0.106 True False 5.505 False
a uniform either 200 cross 21.735 19.227 20.305 16.664 0 True False 21.681 False
a uniform either 200 first 11.932 19.811 1.525 8.231 0.001 True False 11.878 False
b circular before 100 cross 25.651 21.797 20.673 17.784 0 True True 25.597 True
b circular before 100 first 6.143 15.006 1.259 3.542 0.141 True False 6.089 False
b circular before 200 cross 27.514 22.468 26.437 19.637 0 True True 27.359 True
b circular before 200 first 13.674 23.388 1.706 7.696 0.003 True False 13.519 False
b circular either 100 cross 20.770 16.408 18.391 16.128 0 True False 20.716 False
b circular either 100 first 5.099 11.557 1.259 3.154 0.126 True False 5.045 False
b circular either 200 cross 24.278 18.690 20.305 20.888 0 True True 24.123 True
b circular either 200 first 11.932 19.565 1.706 8.163 0.001 True False 11.776 False
b uniform before 100 cross 27.202 23.027 22.115 19.163 0 True True 27.147 True
b uniform before 100 first 6.318 15.222 1.259 3.696 0.133 True False 6.264 False
b uniform before 200 cross 30.570 21.880 29.054 22.466 0 True False 27.925 False
b uniform before 200 first 14.081 24.569 1.525 8.067 0.001 True False 14.027 False
b uniform either 100 cross 21.358 18.285 19.595 15.400 0 True False 21.303 False
b uniform either 100 first 5.559 12.721 1.259 3.696 0.106 True False 5.505 False
b uniform either 200 cross 21.735 19.227 20.305 16.664 0 True False 21.681 False
b uniform either 200 first 11.932 19.811 1.525 8.231 0.001 True False 11.878 False
readings meeting the criterion: 8 of 32
ok n touch_min touch_max time_min time_max p_max
start lat_reading tmax
cross a 100 2 4 20.770 27.202 15.400 19.163 0.000
200 2 4 21.735 30.570 16.664 22.466 0.000
b 100 2 4 20.770 27.202 15.400 19.163 0.000
200 2 4 21.735 30.570 16.664 22.466 0.000
first a 100 0 4 5.099 6.318 3.154 3.696 0.141
200 0 4 11.932 14.081 7.696 8.231 0.003
b 100 0 4 5.099 6.318 3.154 3.696 0.141
200 0 4 11.932 14.081 7.696 8.231 0.003
per cell, first reading of each detail (uniform, before, 100 ms, cross, latency reading a): latency, t_on, fractions, and whether the increment rose above chance:
cell latency t_on frac_spikes frac_time crossed n_explore_spikes has_first_latency
42 3.0 69 0.574 0.486 True 8023 True
43 11.0 1 0.007 0.009 False 894 True
44 10.0 12 0.197 0.122 True 3999 True
45 4.0 15 0.221 0.121 True 208 True
46 0.0 20 0.286 0.173 True 1659 True
47 4.0 4 0.056 0.051 True 1967 True
48 4.0 50 0.688 0.398 True 6355 True
49 3.0 14 0.184 0.165 True 592 True
50 0.0 100 0.253 0.368 True 87 True
51 0.0 1 0.006 0.013 False 1172 False
52 10.0 59 0.520 0.203 True 3595 True
Result: blocked (methods_underspecified). The claim depends on fig3c-touch-spikes-l4c2, which is blocked. Under the transition cutoff 4 of the 8 window readings meet the criterion under both latency readings (20.8 to 30.6% of spikes in 15.4 to 22.5% of time across the 8, paper 24.0 +- 18.6% in 19.0%; p against L4 C2 below 0.001), 2 with the 100 ms curve and 2 with the 200 ms curve. The literal cutoff, a variant, gives 5.1 to 14.1% in 3.2 to 8.2% of time.
Figure 3D: fig3d-touch-phase-l4c2, fig3d-touch-phase-outside, fig3d-touch-phase-l5, fig3d-phase-only¶
In L4 C2 neurons 75.2 +- 19.0% of exploration spikes are coupled to touch or to whisking phase, within an average of 22.0% of exploration time.
L4 neurons outside C2 have fewer touch or phase coupled spikes than L4 C2, 33.0 +- 30.2% in 16.9% of exploration time.
L5 neurons have fewer touch or phase coupled spikes than L4 C2, 31.5 +- 16.4% in 26.1% of exploration time.
Spikes coupled to whisking phase alone are 4.6 +- 8.7% in L4 C2, 23.8 +- 27.9% in L4 outside C2 and 7.6 +- 14.9% in L5.
F3D = []
c2p = PH[PH.group == "L4C2"]
for g, exp in (("L4C2", (75.2, 19.0, 22.0, 2.2)), ("L4out", (33.0, 30.2, 16.9, 1.69)), ("L5", (31.5, 16.4, 26.1, 2.61))):
for key, kk in PH[PH.group == g].groupby(PKEY):
kd = dict(zip(PKEY, key))
ref = c2p[(c2p.lat_reading == "a") & (c2p.shuffle == kd["shuffle"]) & (c2p.nearest == kd["nearest"]) & (c2p.tmax == kd["tmax"]) &
(c2p.start == kd["start"]) & (c2p.width == kd["width"])]
tp = 100 * kk.frac_tp
p = ranksum(ref.frac_tp, kk.frac_tp) if g != "L4C2" else np.nan
F3D.append(dict(group=g, **kd, tp_mean=tp.mean(), tp_sd=tp.std(ddof=1), time_mean=100 * kk.time_tp.mean(),
phase_only_mean=100 * kk.frac_phase_only.mean(), phase_only_sd=100 * kk.frac_phase_only.std(ddof=1),
p_vs_c2=p, lower_than_c2=tp.median() < 100 * ref.frac_tp.median(),
ok_values=within(tp.mean(), exp[0], exp[1]) and within(100 * kk.time_tp.mean(), exp[2], exp[3])))
F3D = pd.DataFrame(F3D)
F3D["ok"] = F3D.ok_values & ((F3D.group == "L4C2") | ((F3D.p_vs_c2 < 0.05) & F3D.lower_than_c2))
t = F3D[(F3D.group == "L4C2") & (F3D.lat_reading == "a")]
print(t.drop(columns=["group", "lat_reading", "p_vs_c2", "lower_than_c2"]).round(2).to_string(index=False))
print("readings within both tolerances: %d of %d" % (t.ok.sum(), len(t)))
print(t.groupby(["start", "lat_reading", "tmax", "width"]).agg(ok=("ok", "sum"), n=("ok", "size"), tp_min=("tp_mean", "min"), tp_max=("tp_mean", "max"),
time_min=("time_mean", "min"), time_max=("time_mean", "max")).round(2).to_string())
shuffle nearest tmax start width tp_mean tp_sd time_mean phase_only_mean phase_only_sd ok_values ok
circular before 100 cross pm 77.82 14.40 25.32 6.96 7.99 False False
circular before 100 cross total 76.37 14.32 22.52 5.51 6.88 True True
circular before 100 first pm 43.72 36.81 10.98 2.94 6.99 False False
circular before 100 first total 43.06 36.49 10.32 2.28 6.52 False False
circular before 200 cross pm 79.99 13.14 26.18 6.10 7.54 False False
circular before 200 cross total 79.70 11.49 25.62 5.82 6.93 False False
circular before 200 first pm 65.01 31.95 15.61 3.65 6.41 False False
circular before 200 first total 63.18 32.29 13.87 1.81 4.21 False False
circular either 100 cross pm 76.83 12.39 26.62 9.95 10.63 False False
circular either 100 cross total 74.89 12.59 24.47 8.01 9.27 False False
circular either 100 first pm 42.96 35.38 12.20 4.35 9.82 False False
circular either 100 first total 41.94 34.95 11.55 3.32 8.83 False False
circular either 200 cross pm 77.73 13.40 28.10 8.88 9.96 False False
circular either 200 cross total 77.62 12.23 26.88 8.76 9.65 False False
circular either 200 first pm 64.83 29.78 17.53 5.59 9.69 False False
circular either 200 first total 62.66 29.49 15.50 3.42 8.47 False False
uniform before 100 cross pm 81.06 10.94 27.11 7.25 9.75 False False
uniform before 100 cross total 78.90 11.76 24.21 5.09 6.72 False False
uniform before 100 first pm 45.35 37.63 11.48 2.51 6.72 False False
uniform before 100 first total 44.65 37.35 10.86 1.82 5.91 False False
uniform before 200 cross pm 81.43 12.75 27.30 5.63 6.87 False False
uniform before 200 cross total 80.66 11.57 25.54 4.86 6.15 False False
uniform before 200 first pm 67.45 30.85 16.45 2.98 6.16 False False
uniform before 200 first total 65.56 31.23 14.42 1.09 3.52 False False
uniform either 100 cross pm 78.78 10.80 28.51 10.22 11.24 False False
uniform either 100 cross total 75.82 11.68 24.59 7.26 8.13 False False
uniform either 100 first pm 44.21 35.30 12.24 3.61 8.78 False False
uniform either 100 first total 43.21 34.83 11.46 2.61 7.41 False False
uniform either 200 cross pm 77.97 13.15 26.44 7.69 9.29 False False
uniform either 200 cross total 78.53 11.63 28.37 8.25 8.67 False False
uniform either 200 first pm 65.33 30.22 17.69 4.96 8.44 False False
uniform either 200 first total 63.20 30.45 15.61 2.83 7.09 False False
readings within both tolerances: 1 of 32
ok n tp_min tp_max time_min time_max
start lat_reading tmax width
cross a 100 pm 0 4 76.83 81.06 25.32 28.51
total 1 4 74.89 78.90 22.52 24.59
200 pm 0 4 77.73 81.43 26.18 28.10
total 0 4 77.62 80.66 25.54 28.37
first a 100 pm 0 4 42.96 45.35 10.98 12.24
total 0 4 41.94 44.65 10.32 11.55
200 pm 0 4 64.83 67.45 15.61 17.69
total 0 4 62.66 65.56 13.87 15.61
Result: blocked (methods_underspecified). The claim depends on fig3c-touch-spikes-l4c2. Under the transition cutoff 74.9 to 81.4% of spikes are touch or phase coupled in 22.5 to 28.5% of exploration time, and 1 of the 16 readings puts both values in tolerance (75.2 +- 19.0% in 22.0 +- 2.2%): circular shuffle, onset before, 100 ms curve, window t_w wide in total, 76.4% in 22.5%. The literal cutoff, a variant, gives 41.9 to 67.5% in 10.3 to 17.7%. With the preferred phase now fitted over the whole trial, the count of readings that meet the criterion fell from 4 to 1.
t = F3D[F3D.group == "L4out"]
print(t.groupby(["start", "lat_reading", "width"])[["tp_mean", "time_mean", "p_vs_c2"]].agg(["min", "max"]).to_string(float_format="{:.3g}".format))
print("readings meeting the criterion: %d of %d" % (t.ok.sum(), len(t)))
print(t.groupby(["start", "lat_reading", "tmax", "width"]).agg(ok=("ok", "sum"), n=("ok", "size"), tp_min=("tp_mean", "min"), tp_max=("tp_mean", "max"),
time_min=("time_mean", "min"), time_max=("time_mean", "max")).round(2).to_string())
tp_mean time_mean p_vs_c2
min max min max min max
start lat_reading width
cross a pm 42.9 57.6 29.1 39.3 8.99e-05 0.00293
total 42.4 65.1 28 42.8 0.000167 0.00293
b pm 44.8 58 30.2 38.6 0.000167 0.00129
total 43 59.1 29 39.1 0.000148 0.000839
first a pm 28 30 11.1 12.8 0.00177 0.485
total 19 26.9 5.89 10.7 0.00217 0.395
b pm 28.5 30.5 11.1 12.6 0.0024 0.627
total 19.5 27.2 5.9 10.7 0.00356 0.495
readings meeting the criterion: 0 of 64
ok n tp_min tp_max time_min time_max
start lat_reading tmax width
cross a 100 pm 0 4 49.27 57.62 32.69 39.30
total 0 4 42.39 50.73 28.03 30.28
200 pm 0 4 42.92 53.55 29.13 36.30
total 0 4 49.00 65.05 35.91 42.83
b 100 pm 0 4 50.88 56.79 32.40 35.70
total 0 4 42.98 50.97 28.98 30.54
200 pm 0 4 44.83 58.03 30.18 38.62
total 0 4 49.82 59.12 34.23 39.12
first a 100 pm 0 4 28.04 28.35 11.10 11.49
total 0 4 18.97 19.35 5.89 6.00
200 pm 0 4 28.25 30.02 11.71 12.80
total 0 4 23.63 26.92 8.88 10.73
b 100 pm 0 4 28.52 29.00 11.15 11.66
total 0 4 19.45 20.00 5.90 6.18
200 pm 0 4 28.73 30.50 11.54 12.63
total 0 4 24.10 27.21 8.82 10.69
Result: blocked (methods_underspecified). The claim depends on fig3d-touch-phase-l4c2. None of the 32 readings with the transition cutoff meets the criterion: the touch or phase fraction is 42.4 to 65.1% (within 33.0 +- 30.2%), but the time fraction is 28.0 to 42.8% against 16.9 +- 1.69%. The rank-sum test against L4 C2 gives p = 9.0e-5 to 0.0029. None of the literal variant meets it either (5.9 to 12.8% of time).
t = F3D[F3D.group == "L5"]
print(t.groupby(["start", "lat_reading", "width"])[["tp_mean", "time_mean", "p_vs_c2"]].agg(["min", "max"]).to_string(float_format="{:.3g}".format))
print("readings meeting the criterion: %d of %d" % (t.ok.sum(), len(t)))
print(t.groupby(["start", "lat_reading", "tmax", "width"]).agg(ok=("ok", "sum"), n=("ok", "size"), tp_min=("tp_mean", "min"), tp_max=("tp_mean", "max"),
time_min=("time_mean", "min"), time_max=("time_mean", "max")).round(2).to_string())
tp_mean time_mean p_vs_c2
min max min max min max
start lat_reading width
cross a pm 38.9 48.9 27.5 36.4 2.91e-06 0.000237
total 32.9 47 23.6 35.6 3.84e-06 9.43e-05
b pm 38.9 48.9 27.5 36.4 2.91e-06 0.000237
total 32.9 47 23.6 35.6 3.84e-06 9.43e-05
first a pm 11.2 24.9 7.83 16.1 0.000629 0.247
total 10.2 22.7 6.71 15.4 0.000458 0.235
b pm 11.2 24.9 7.83 16.1 0.000629 0.247
total 10.2 22.7 6.71 15.4 0.000458 0.235
readings meeting the criterion: 12 of 64
ok n tp_min tp_max time_min time_max
start lat_reading tmax width
cross a 100 pm 1 4 38.88 45.03 27.53 33.61
total 4 4 32.85 38.06 23.59 28.00
200 pm 0 4 40.88 48.94 30.19 36.37
total 1 4 36.29 46.99 26.44 35.64
b 100 pm 1 4 38.88 45.03 27.53 33.61
total 4 4 32.85 38.06 23.59 28.00
200 pm 0 4 40.88 48.94 30.19 36.37
total 1 4 36.29 46.99 26.44 35.64
first a 100 pm 0 4 11.19 14.89 7.83 10.19
total 0 4 10.20 11.31 6.71 7.34
200 pm 0 4 21.90 24.94 14.69 16.14
total 0 4 16.95 22.74 11.86 15.38
b 100 pm 0 4 11.19 14.89 7.83 10.19
total 0 4 10.20 11.31 6.71 7.34
200 pm 0 4 21.90 24.94 14.69 16.14
total 0 4 16.95 22.74 11.86 15.38
Result: blocked (methods_underspecified). The claim depends on fig3d-touch-phase-l4c2. Under the transition cutoff 12 of the 32 readings meet the criterion, 6 under each latency reading: all 4 with the 100 ms curve and the phase window t_w wide in total, 1 with the 100 ms curve and +-t_w, and 1 with the 200 ms curve and the total width. Over the 32 the fraction is 32.9 to 48.9% in 23.6 to 36.4% of time (paper 31.5 +- 16.4% in 26.1 +- 2.61%), p against L4 C2 below 2.4e-4. None of the literal variant does (6.7 to 16.1% of time).
po = F3D.pivot_table(index=PKEY, columns="group", values="phase_only_mean")
po = po.reset_index()
po = po[po.lat_reading == "a"].drop(columns="lat_reading").merge(
F3D[(F3D.lat_reading == "b")].pivot_table(index=PKEY[1:], columns="group", values="phase_only_mean")[["L4out", "L5"]].add_suffix("_latb").reset_index(),
on=PKEY[1:])
po["ok"] = [within(r.L4C2, 4.6, 8.7) and within(r.L4out, 23.8, 27.9) and within(r.L5, 7.6, 14.9) and r.L4out > max(r.L4C2, r.L5) for r in po.itertuples()]
po["ok_latb"] = [within(r.L4C2, 4.6, 8.7) and within(r.L4out_latb, 23.8, 27.9) and within(r.L5_latb, 7.6, 14.9) and r.L4out_latb > max(r.L4C2, r.L5_latb) for r in po.itertuples()]
print(po.round(2).to_string(index=False))
print("readings meeting the criterion: latency reading (a) %d, (b) %d, of %d" % (po.ok.sum(), po.ok_latb.sum(), len(po)))
print(po.groupby("start")[["ok", "ok_latb"]].agg(["sum", "size"]).to_string())
shuffle nearest tmax start width L4C2 L4out L5 L4out_latb L5_latb ok ok_latb
circular before 100 cross pm 6.96 28.39 17.71 21.21 17.71 True True
circular before 100 cross total 5.51 18.57 12.41 18.55 12.41 True True
circular before 100 first pm 2.94 26.51 8.75 26.51 8.75 True True
circular before 100 first total 2.28 17.44 4.32 17.44 4.32 True True
circular before 200 cross pm 6.10 23.07 21.42 23.02 21.42 True True
circular before 200 cross total 5.82 36.23 12.98 28.82 12.98 True True
circular before 200 first pm 3.65 28.26 11.27 28.26 11.27 True True
circular before 200 first total 1.81 25.17 9.07 24.97 9.07 True True
circular either 100 cross pm 9.95 37.08 18.96 29.80 18.96 True True
circular either 100 cross total 8.01 21.85 14.97 21.84 14.97 True True
circular either 100 first pm 4.35 26.66 8.77 26.66 8.77 True True
circular either 100 first total 3.32 17.44 5.10 17.44 5.10 True True
circular either 200 cross pm 8.88 24.15 19.72 25.28 19.72 True True
circular either 200 cross total 8.76 36.72 17.91 30.27 17.91 True True
circular either 200 first pm 5.59 22.19 10.09 22.19 10.09 True True
circular either 200 first total 3.42 19.46 5.21 19.46 5.21 True True
uniform before 100 cross pm 7.25 30.50 17.83 30.47 17.83 True True
uniform before 100 cross total 5.09 24.56 8.34 24.66 8.34 True True
uniform before 100 first pm 2.51 26.82 8.48 26.82 8.48 True True
uniform before 100 first total 1.82 17.64 4.99 17.64 4.99 True True
uniform before 200 cross pm 5.63 21.84 14.74 26.65 14.74 True True
uniform before 200 cross total 4.86 27.46 16.42 26.82 16.42 True True
uniform before 200 first pm 2.98 26.48 10.46 26.48 10.46 True True
uniform before 200 first total 1.09 23.12 5.20 23.12 5.20 True True
uniform either 100 cross pm 10.22 30.00 17.53 34.62 17.53 True True
uniform either 100 cross total 7.26 27.72 11.49 28.07 11.49 True True
uniform either 100 first pm 3.61 26.82 5.63 26.82 5.63 True True
uniform either 100 first total 2.61 17.82 5.06 17.82 5.06 True True
uniform either 200 cross pm 7.69 24.36 19.15 30.03 19.15 True True
uniform either 200 cross total 8.25 24.26 14.56 29.80 14.56 True True
uniform either 200 first pm 4.96 22.78 9.97 22.78 9.97 True True
uniform either 200 first total 2.83 16.93 5.02 16.93 5.02 True True
readings meeting the criterion: latency reading (a) 32, (b) 32, of 32
group ok ok_latb
sum size sum size
start
cross 16 16 16 16
first 16 16 16 16
Result: blocked (methods_underspecified). The claim depends on fig3d-touch-phase-l4c2, which is blocked. Under the transition cutoff the phase-only fractions meet the criterion in all 16 readings under both latency readings (L4 C2 4.9 to 10.2%, L4 outside 18.6 to 37.1%, L5 8.3 to 21.4%, the outside mean the largest in each), and so do all 16 of the literal variant.
Figure 4A: fig4a-touch-detection-c2 and fig4a-touch-detection-outside¶
Pooling the activity of only 15 L4 C2 neurons, integrated over 10 ms, is enough to detect 95% of touches.
Touch detection by neurons from surrounding barrels, outside C2, is poor.
Touch-aligned responses from -30 to +50 ms, leaving out touches with another onset within 50
ms; N neurons drawn with replacement, 1000 responses each with replacement, 10 ms causal
integration, pooled; one threshold per post-touch time point, chosen on half the samples to
maximise percent correct, tested on the other half. Open details from the claims file: (1)
the non-touch epochs, (a) the 30 ms before touch in the same responses, or (b) 10 ms windows of
exploration time away from touch; (2) the measure at the best post-touch time point, (a)
fraction of touches detected or (b) percent correct, both with the threshold that maximises
percent correct on the training half. Within (1a), after the answer to the question, the
criterion lists three samples: (a1) the window ending at onset (pre_end, -10 to 0 ms), (a2)
the earliest (pre_start, -30 to -20 ms) and (a3) every pre-touch window pooled (pre_all, the
21 windows from -30 to -20 ms through -10 to 0 ms). Reading (1b) uses windows with no touch
onset from 60 ms before to 100 ms after them (a standard choice).
The fraction of touches detected is the hit rate at the time point of best percent correct (the
same time point as measure (2b)); thresholds that tie for the best training percent correct are
broken at the middle one, a standard choice, and the lowest one (the earlier run's rule) is shown
as a variant. 50 repetitions per population size give the mean and the 95% bounds.
t0 = time.time()
DEC = {}
for g in ("L4C2", "L4out"):
I_list, Q_list = [], []
for no in G[g]:
c = by_no[no]
d = D[no]
R = hires.touch_responses(c, d["s"], d["T"])
Ii = hires.causal_integrate(R.astype(np.int32))
I_list.append(Ii)
Q_list.append({"pre_end": Ii[:, 29], "pre_start": Ii[:, 9], "pre_all": Ii[:, 9:30].ravel(), "explore": hires.nontouch_windows(c, d["s"], d["T"], d["E"])})
DEC[g] = (I_list, Q_list)
print(g, "touch responses per cell:", [len(x) for x in I_list], "| non-touch windows per cell (1b):", [len(q["explore"]) for q in Q_list])
N_GRID = list(range(1, 21)) + [25, 30, 40, 60, 100, 200]
REPS4A = 50
rng4a = np.random.default_rng(4100)
det_rows = []
for g in ("L4C2", "L4out"):
I_list, Q_list = DEC[g]
for N in N_GRID:
res = hires.touch_detection(I_list, Q_list, N, rng4a, reps=REPS4A)
for rd, arr in res.items():
for tie, off in (("middle", 0), ("lowest", 2)):
hit, pc = arr[:, :, off], arr[:, :, off + 1]
tb = pc.argmax(axis=1) # the best post-touch time point, by percent correct
for measure, best in (("detected", hit[np.arange(len(tb)), tb]), ("percent_correct", pc.max(axis=1))):
det_rows.append(dict(group=g, N=N, nontouch=rd, tie=tie, measure=measure, mean=best.mean(),
lo=np.percentile(best, 2.5), hi=np.percentile(best, 97.5)))
DET_ALL = pd.DataFrame(det_rows)
DET = DET_ALL[DET_ALL.tie == "middle"].drop(columns="tie")
print(f"{time.time() - t0:.0f} s")
L4C2 touch responses per cell: [207, 102, 157, 263, 130, 114, 61, 164, 574, 1146, 532, 866, 343, 1274, 256, 561, 1214, 245, 404, 1013, 444, 69, 797, 408, 512, 1085, 1585, 431, 157, 410, 804] | non-touch windows per cell (1b): [2457, 1376, 1283, 1917, 1466, 1195, 1913, 1051, 3710, 7304, 3131, 4180, 2721, 4603, 3772, 4507, 11740, 2215, 2320, 7478, 2854, 1453, 5018, 1686, 2764, 5741, 11938, 2960, 1202, 5273, 3453] L4out touch responses per cell: [1244, 569, 411, 659, 564, 777, 650, 636, 320, 1245] | non-touch windows per cell (1b): [9135, 7650, 6259, 10925, 2757, 5082, 6426, 3009, 1983, 6563]
89 s
# check of the detector on synthetic neurons whose touch and non-touch counts come from the same
# distribution (Poisson, mean 0.3), N = 1: the fraction detected should not reach 1
rng_syn = np.random.default_rng(7)
I_syn = [rng_syn.poisson(0.3, (300, 80)) for _ in range(5)]
Q_syn = [{"same": rng_syn.poisson(0.3, 2000)} for _ in range(5)]
r_syn = hires.touch_detection(I_syn, Q_syn, 1, np.random.default_rng(8), reps=50)["same"]
for tie, off in (("middle", 0), ("lowest", 2)):
hit, pc = r_syn[:, :, off], r_syn[:, :, off + 1]
tb = pc.argmax(axis=1)
print(f"ties at the {tie} threshold: fraction detected {hit[np.arange(len(tb)), tb].mean():.3f}, percent correct {pc.max(axis=1).mean():.3f}")
print(f"earlier run's rule (lowest threshold, largest hit rate over time points): fraction detected {r_syn[:, :, 2].max(axis=1).mean():.3f}")
ties at the middle threshold: fraction detected 0.329, percent correct 0.541 ties at the lowest threshold: fraction detected 0.329, percent correct 0.541 earlier run's rule (lowest threshold, largest hit rate over time points): fraction detected 0.987
d4 = DET[DET.group == "L4C2"]
n95 = d4[d4["mean"] >= 0.95].groupby(["nontouch", "measure"]).N.min().reindex(pd.MultiIndex.from_product([["pre_end", "pre_start", "pre_all", "explore"], ["detected", "percent_correct"]]))
print(f"{REPS4A} repetitions per population size; the fraction detected is the hit rate at the time point of best percent correct;")
print("thresholds tied for the best training percent correct are broken at the middle one")
print("smallest N whose mean performance reaches 95%:")
print(n95.to_string())
print("within 15 +- 1.5:", {k: within(v, 15, 1.5) for k, v in n95.items()})
print(d4[d4.N.isin([1, 5, 10, 13, 14, 15, 16, 17, 20, 30])].pivot_table(index="N", columns=["nontouch", "measure"], values="mean").round(3).to_string())
print("95% bounds at N = 13 to 17 (lo, hi):")
print(d4[d4.N.isin([13, 14, 15, 16, 17])].pivot_table(index="N", columns=["nontouch", "measure"], values=["lo", "hi"]).round(3).to_string())
dv = DET_ALL[(DET_ALL.group == "L4C2") & (DET_ALL.tie == "lowest")]
print("variant, ties broken at the lowest threshold (the earlier run's rule), smallest N reaching 95%:")
print(dv[dv["mean"] >= 0.95].groupby(["nontouch", "measure"]).N.min().to_string())
50 repetitions per population size; the fraction detected is the hit rate at the time point of best percent correct;
thresholds tied for the best training percent correct are broken at the middle one
smallest N whose mean performance reaches 95%:
pre_end detected 19
percent_correct 19
pre_start detected 30
percent_correct 25
pre_all detected 25
percent_correct 25
explore detected 14
percent_correct 12
within 15 +- 1.5: {('pre_end', 'detected'): False, ('pre_end', 'percent_correct'): False, ('pre_start', 'detected'): False, ('pre_start', 'percent_correct'): False, ('pre_all', 'detected'): False, ('pre_all', 'percent_correct'): False, ('explore', 'detected'): True, ('explore', 'percent_correct'): False}
nontouch explore pre_all pre_end pre_start
measure detected percent_correct detected percent_correct detected percent_correct detected percent_correct
N
1 0.296 0.642 0.296 0.621 0.296 0.628 0.296 0.618
5 0.784 0.854 0.746 0.793 0.763 0.808 0.734 0.776
10 0.927 0.935 0.854 0.864 0.871 0.886 0.846 0.848
13 0.947 0.956 0.888 0.895 0.909 0.916 0.868 0.878
14 0.954 0.963 0.912 0.911 0.922 0.929 0.895 0.893
15 0.966 0.969 0.909 0.914 0.927 0.935 0.897 0.898
16 0.961 0.972 0.918 0.925 0.929 0.944 0.910 0.909
17 0.968 0.974 0.932 0.931 0.945 0.949 0.918 0.920
20 0.983 0.986 0.943 0.949 0.952 0.962 0.933 0.937
30 0.994 0.994 0.969 0.971 0.981 0.981 0.961 0.963
95% bounds at N = 13 to 17 (lo, hi):
hi lo
nontouch explore pre_all pre_end pre_start explore pre_all pre_end pre_start
measure detected percent_correct detected percent_correct detected percent_correct detected percent_correct detected percent_correct detected percent_correct detected percent_correct detected percent_correct
N
13 0.995 0.990 0.969 0.967 0.980 0.977 0.963 0.952 0.876 0.894 0.724 0.794 0.741 0.824 0.721 0.786
14 0.990 0.992 0.974 0.974 0.988 0.979 0.966 0.963 0.889 0.928 0.772 0.842 0.813 0.861 0.748 0.826
15 0.996 0.995 0.985 0.980 0.992 0.986 0.966 0.961 0.892 0.930 0.767 0.786 0.798 0.810 0.767 0.785
16 0.994 0.992 0.976 0.975 0.982 0.982 0.976 0.962 0.866 0.917 0.824 0.843 0.833 0.864 0.824 0.824
17 0.998 0.995 0.984 0.980 0.988 0.990 0.977 0.967 0.912 0.943 0.851 0.880 0.863 0.888 0.829 0.868
variant, ties broken at the lowest threshold (the earlier run's rule), smallest N reaching 95%:
nontouch measure
explore detected 14
percent_correct 12
pre_all detected 25
percent_correct 25
pre_end detected 19
percent_correct 19
pre_start detected 25
percent_correct 25
Result: blocked (methods_underspecified). The smallest population whose mean performance reaches 95% depends on the readings: 14 (fraction detected) or 12 (percent correct) neurons with non-touch windows from exploration time; 19, 25 and 25 with percent correct and the pre-touch windows (a1), (a2) and (a3); 19, 30 and 25 with the fraction detected and the same windows. Only the exploration windows with the fraction detected (14) is within 15 +- 1.5. Breaking threshold ties at the lowest threshold, a variant, changes only the (a2) fraction-detected value (25). With one neuron, percent correct is 0.618 to 0.642. On synthetic neurons whose touch and non-touch counts come from one distribution, one neuron gives a fraction detected of 0.329 at percent correct 0.541 under either tie rule, against 0.987 under the earlier run's rule (the largest hit rate over time points).
print('performance at N = 1:'); print(DET[DET.N == 1].pivot_table(index=['nontouch', 'measure'], columns='group', values='mean').round(3).to_string())
d15 = DET[DET.N == 15].pivot_table(index=["nontouch", "measure"], columns="group", values=["mean", "hi"])
print(d15.round(3).to_string())
for (rd, ms), g in DET.groupby(["nontouch", "measure"]):
c2_15 = g[(g.group == "L4C2") & (g.N == 15)]["mean"].iloc[0]
o = g[g.group == "L4out"]
o15 = o[o.N == 15].iloc[0]
reach = o[o["mean"] >= c2_15].N.min()
print(f"{rd:9s} {ms:15s}: N=15 outside {o15['mean']:.3f} (upper bound {o15.hi:.3f}) vs C2 {c2_15:.3f}; below C2 and bound < 0.95: {o15['mean'] < c2_15 and o15.hi < 0.95};"
f" outside N for C2's N=15 performance: {reach if pd.notna(reach) else '> 200'} (ratio {reach / 15 if pd.notna(reach) else np.inf:.1f})")
fig, ax = plt.subplots(1, 2, figsize=(8, 3), sharey=True)
for a_, ms in zip(ax, ("detected", "percent_correct")):
for g, col in (("L4C2", "red"), ("L4out", "grey")):
x = DET[(DET.group == g) & (DET.nontouch == "explore") & (DET.measure == ms)]
a_.plot(x.N, x["mean"], color=col, label=g)
a_.fill_between(x.N, x.lo, x.hi, color=col, alpha=0.2)
a_.set_xscale("log"); a_.axhline(0.95, ls="--", c="k", lw=0.5); a_.set_xlabel("neurons"); a_.set_title(f"explore, {ms}", fontsize=8)
ax[0].set_ylabel("performance"); ax[0].legend(fontsize=7)
plt.tight_layout(); plt.show()
print("plotted range %.3f to %.3f" % (DET[DET.nontouch == "explore"].lo.min(), DET[DET.nontouch == "explore"].hi.max()))
performance at N = 1:
group L4C2 L4out
nontouch measure
explore detected 0.296 0.042
percent_correct 0.642 0.511
pre_all detected 0.296 0.046
percent_correct 0.621 0.510
pre_end detected 0.296 0.055
percent_correct 0.628 0.514
pre_start detected 0.296 0.034
percent_correct 0.618 0.510
hi mean
group L4C2 L4out L4C2 L4out
nontouch measure
explore detected 0.996 0.473 0.966 0.252
percent_correct 0.995 0.583 0.969 0.534
pre_all detected 0.985 0.458 0.909 0.196
percent_correct 0.980 0.557 0.914 0.522
pre_end detected 0.992 0.481 0.927 0.351
percent_correct 0.986 0.570 0.935 0.541
pre_start detected 0.966 0.455 0.897 0.197
percent_correct 0.961 0.564 0.898 0.529
explore detected : N=15 outside 0.252 (upper bound 0.473) vs C2 0.966; below C2 and bound < 0.95: True; outside N for C2's N=15 performance: > 200 (ratio inf)
explore percent_correct: N=15 outside 0.534 (upper bound 0.583) vs C2 0.969; below C2 and bound < 0.95: True; outside N for C2's N=15 performance: > 200 (ratio inf)
pre_all detected : N=15 outside 0.196 (upper bound 0.458) vs C2 0.909; below C2 and bound < 0.95: True; outside N for C2's N=15 performance: > 200 (ratio inf)
pre_all percent_correct: N=15 outside 0.522 (upper bound 0.557) vs C2 0.914; below C2 and bound < 0.95: True; outside N for C2's N=15 performance: > 200 (ratio inf)
pre_end detected : N=15 outside 0.351 (upper bound 0.481) vs C2 0.927; below C2 and bound < 0.95: True; outside N for C2's N=15 performance: > 200 (ratio inf)
pre_end percent_correct: N=15 outside 0.541 (upper bound 0.570) vs C2 0.935; below C2 and bound < 0.95: True; outside N for C2's N=15 performance: > 200 (ratio inf)
pre_start detected : N=15 outside 0.197 (upper bound 0.455) vs C2 0.897; below C2 and bound < 0.95: True; outside N for C2's N=15 performance: > 200 (ratio inf)
pre_start percent_correct: N=15 outside 0.529 (upper bound 0.564) vs C2 0.898; below C2 and bound < 0.95: True; outside N for C2's N=15 performance: > 200 (ratio inf)
plotted range 0.000 to 1.000
Result: blocked (methods_underspecified). The claim depends on fig4a-touch-detection-c2, which is blocked. Its own criterion holds under every reading: with 15 outside neurons percent correct is 0.522 to 0.541 (upper bound 0.557 to 0.583) and the fraction of touches detected 0.196 to 0.351 (upper bound 0.455 to 0.481), against 0.897 to 0.969 in C2, and no outside population up to 200 reaches C2's level at 15. One outside neuron performs at 0.510 to 0.514 percent correct.
fig4b-touch-time-precision¶
A naive Bayes decoder reading 200 L4 C2 neurons integrated over 10 ms decodes the time elapsed since touch with high precision, a minimum uncertainty of 0.55 ms (95% confidence) at 10 ms after touch onset.
Multinomial naive Bayes over the 10 ms integrated responses of 200 neurons, 100 bootstrap runs. Open details from the claims file: (1) the uncertainty, (a) RMS error, (b) half width of the 95% bound of decoded time, (c) 1.96 x RMS; (2) training and test sets, (a) halves of 1000 draws or (b) independent sets of 1000 draws; (3) the times the decoder chooses among, added after the answer to the question, (c-all) every 1 ms where a full 10 ms window fits, from -21 ms, or (c-post) every 1 ms from onset on; (4) the smoothing of zero counts, add-one (alpha = 1), with a smaller pseudocount (alpha = 0.01) as a variant.
t0 = time.time()
I_c2 = DEC["L4C2"][0]
tb_rows = []
for si, split in enumerate(("half", "independent")):
for ci, t0b in enumerate((9, 30)):
for ai, alpha in enumerate((1.0, 0.01)): # add-one smoothing (criterion), a smaller pseudocount (variant)
rng4b = np.random.default_rng([4200, si, ci] + ([ai] if ai else []))
dec, tru, rep = hires.time_decoding(I_c2, rng4b, N=200, reps=100, split=split, t0_bin=t0b, alpha=alpha)
for t in range(0, 21):
sel = tru == t
err = dec[sel] - t
rms = np.sqrt(np.mean(err ** 2))
lo, hi = np.percentile(dec[sel], [2.5, 97.5])
tb_rows.append(dict(split=split, classes_from=t0b - 30, alpha=alpha, t=t, rms=rms, half95=(hi - lo) / 2,
rms196=1.96 * rms, exact=np.mean(err == 0)))
TB = pd.DataFrame(tb_rows)
print(f"{time.time() - t0:.0f} s")
118 s
print(TB[TB.t.isin([6, 8, 9, 10, 11, 12, 14]) & (TB.alpha == 1.0)].round(3).to_string(index=False))
near = TB[(TB.t >= 5) & (TB.t <= 15)]
mins = near.groupby(["alpha", "split", "classes_from"])[["rms", "half95", "rms196"]].min()
print("minimum over 5 to 15 ms after onset:")
print(mins.round(3).to_string())
print("minimum over all times 0 to 20 ms, and where:")
for (al, sp, cf), g in TB.groupby(["alpha", "split", "classes_from"]):
print(al, sp, cf, {m: (round(g[m].min(), 3), int(g.t[g[m].idxmin()])) for m in ("rms", "half95", "rms196")})
print("within 0.55 +- 0.055:", {k: [m for m in ("rms", "half95", "rms196") if within(v[m], 0.55, 0.055)] for k, v in mins.iterrows()})
rng4bi = np.random.default_rng(4210)
inset = []
for N in (10, 25, 50, 100, 200):
dec, tru, rep = hires.time_decoding(I_c2, rng4bi, N=N, reps=10, split="half", t0_bin=9)
sel = tru == 10
inset.append((N, np.sqrt(np.mean((dec[sel] - 10) ** 2))))
print("RMS error at 10 ms by population size (10 runs each):", [(n, round(e, 2)) for n, e in inset])
split classes_from alpha t rms half95 rms196 exact
half -21 1.0 6 9.009 22.0 17.657 0.363
half -21 1.0 8 1.396 2.0 2.737 0.544
half -21 1.0 9 1.643 3.0 3.220 0.440
half -21 1.0 10 1.620 3.0 3.176 0.409
half -21 1.0 11 1.717 3.0 3.366 0.284
half -21 1.0 12 1.819 3.5 3.565 0.289
half -21 1.0 14 1.792 3.5 3.512 0.218
half 0 1.0 6 6.888 15.0 13.501 0.401
half 0 1.0 8 1.124 2.0 2.202 0.546
half 0 1.0 9 1.444 3.0 2.830 0.442
half 0 1.0 10 1.531 3.0 3.000 0.405
half 0 1.0 11 1.620 3.0 3.175 0.280
half 0 1.0 12 1.729 3.5 3.389 0.285
half 0 1.0 14 1.748 3.5 3.426 0.218
independent -21 1.0 6 8.484 21.0 16.628 0.383
independent -21 1.0 8 1.307 2.0 2.562 0.547
independent -21 1.0 9 1.583 3.0 3.103 0.450
independent -21 1.0 10 1.585 3.0 3.107 0.414
independent -21 1.0 11 1.643 3.0 3.221 0.288
independent -21 1.0 12 1.730 3.5 3.392 0.294
independent -21 1.0 14 1.768 3.0 3.465 0.218
independent 0 1.0 6 6.600 12.5 12.935 0.416
independent 0 1.0 8 1.093 2.0 2.142 0.546
independent 0 1.0 9 1.409 3.0 2.761 0.451
independent 0 1.0 10 1.485 3.0 2.910 0.415
independent 0 1.0 11 1.582 3.0 3.101 0.290
independent 0 1.0 12 1.699 3.5 3.331 0.294
independent 0 1.0 14 1.724 3.0 3.380 0.221
minimum over 5 to 15 ms after onset:
rms half95 rms196
alpha split classes_from
0.01 half -21 1.481 2.0 2.903
0 1.314 2.0 2.575
independent -21 1.283 2.0 2.514
0 1.175 2.0 2.303
1.00 half -21 1.396 2.0 2.737
0 1.124 2.0 2.202
independent -21 1.307 2.0 2.562
0 1.093 2.0 2.142
minimum over all times 0 to 20 ms, and where:
0.01 half -21 {'rms': (np.float64(1.481), 8), 'half95': (np.float64(2.0), 8), 'rms196': (np.float64(2.903), 8)}
0.01 half 0 {'rms': (np.float64(1.124), 18), 'half95': (np.float64(2.0), 8), 'rms196': (np.float64(2.203), 18)}
0.01 independent -21 {'rms': (np.float64(1.283), 8), 'half95': (np.float64(2.0), 8), 'rms196': (np.float64(2.514), 8)}
0.01 independent 0 {'rms': (np.float64(1.098), 18), 'half95': (np.float64(2.0), 8), 'rms196': (np.float64(2.152), 18)}
1.0 half -21 {'rms': (np.float64(1.396), 8), 'half95': (np.float64(2.0), 8), 'rms196': (np.float64(2.737), 8)}
1.0 half 0 {'rms': (np.float64(1.113), 18), 'half95': (np.float64(2.0), 8), 'rms196': (np.float64(2.182), 18)}
1.0 independent -21 {'rms': (np.float64(1.307), 8), 'half95': (np.float64(2.0), 8), 'rms196': (np.float64(2.562), 8)}
1.0 independent 0 {'rms': (np.float64(1.093), 8), 'half95': (np.float64(2.0), 8), 'rms196': (np.float64(2.142), 8)}
within 0.55 +- 0.055: {(0.01, 'half', -21): [], (0.01, 'half', 0): [], (0.01, 'independent', -21): [], (0.01, 'independent', 0): [], (1.0, 'half', -21): [], (1.0, 'half', 0): [], (1.0, 'independent', -21): [], (1.0, 'independent', 0): []}
RMS error at 10 ms by population size (10 runs each): [(10, np.float64(19.88)), (25, np.float64(11.58)), (50, np.float64(6.42)), (100, np.float64(3.21)), (200, np.float64(1.55))]
Result: not reproduced. With add-one smoothing the smallest uncertainty 5 to 15 ms after onset, at 8 ms after onset, is an RMS error of 1.09 to 1.40 ms, a half width of the 95% bound of 2.0 ms, or 2.14 to 2.74 ms as 1.96 x RMS, under both splits and both sets of time classes; none is within 0.55 +- 0.055 ms. The smaller pseudocount, a variant, gives 1.18 to 1.48 ms RMS and the same 2.0 ms half width. The error falls with population size (19.9 ms at 10 neurons to 1.55 ms at 200, at 10 ms after onset), as the Figure 4B inset shows.
fig4c-phase-decoding¶
Even 1000 L4 neurons give only a rough estimate of whisking phase, a resolution of 24 deg in C2 and 33 deg outside C2 at d' = 1 (76% correct discrimination).
Per cell, spike counts in 120 phase bins for every whisk cycle used for phase tuning, as the answered criterion says: complete cycles with amplitude above 2.5 deg throughout, monotonic phase, no lick and not within 70 ms after a touch, over the whole trial. 1000 neurons drawn with replacement and circularly shifted to uniformly distributed preferred phases (from the cosine fit; 0 for the cells without phase spikes); a Poisson naive Bayes decision between two phases, 50 test samples per phase and phase difference, phase differences of 6 to 180 deg in 6 deg steps, 100 runs, median. Open details from the claims file: (1) responses, (a) counts of held-out cycles or (b) Poisson counts with the cycle-averaged rates; (2) shifts, (a) random or (b) even; (3) empty bins, add-one smoothing of the counts, or under (1b) a floor rate, here one spike over all of a cell's cycles and bins (the paper gives no value). Two variants: the same cycles inside the sample epoch only (the earlier run's), and every cycle peaking above 2.5 deg with one spike spread over the 120 bins (the first run's).
t0 = time.time()
CB, CB_all, CB_x = {}, {}, {}
for g in ("L4C2", "L4out"):
Cl, Ca, Cx, pref = [], [], [], []
for no in G[g]:
cyc = D[no]["cyc"]
cp = cyc[cyc.phase_ok] # the cycles of phase tuning (criterion)
Cl.append(hires.cycle_bin_counts(by_no[no], D[no]["s"], cp))
ca = cyc[cyc.amp_max > hires.WHISK_AMP] # every cycle peaking above 2.5 deg (earlier run, variant)
Ca.append(hires.cycle_bin_counts(by_no[no], D[no]["s"], ca))
cx = cyc[cyc.phase_ok & cyc.in_sample] # phase-tuning cycles inside the sample epoch (earlier run, variant)
Cx.append(hires.cycle_bin_counts(by_no[no], D[no]["s"], cx))
phi = base.loc[no, "phi_pref"]
pref.append(0.0 if pd.isna(phi) else phi)
CB[g] = (Cl, np.array(pref))
CB_all[g] = (Ca, np.array(pref))
CB_x[g] = (Cx, np.array(pref))
for name, L in (("phase-tuning cycles", Cl), ("phase-tuning cycles in the sample epoch (variant)", Cx), ("all cycles > 2.5 deg (variant)", Ca)):
print(g, name, "| cycles per cell: median %d, range %d to %d; spikes per cycle, mean over cells %.3f; cells with no spike %d of %d" % (
np.median([len(x) for x in L]), min(len(x) for x in L), max(len(x) for x in L),
np.mean([x.sum() / max(len(x), 1) for x in L]), sum(x.sum() == 0 for x in L), len(L)))
READ4C = [("cycles", "add_one"), ("poisson", "add_one"), ("poisson", "floor")]
pd_rows = []
for g in ("L4C2", "L4out"):
for mi, (model, smooth) in enumerate(READ4C):
for shi, shifts in enumerate(("random", "even")):
rng4c = np.random.default_rng([4300, mi, shi, GROUP_ORDER.index(g), 1])
res, cur = hires.phase_discrimination(CB[g][0], CB[g][1], rng4c, N=1000, reps=100, model=model, shifts=shifts, smooth=smooth)
pd_rows.append(dict(cycles="phase tuning", group=g, model=model, smooth=smooth, shifts=shifts,
median_deg=np.nanmedian(res) if np.isfinite(res).any() else np.nan,
runs_reaching=int(np.isfinite(res).sum()), lo=np.nanpercentile(res, 2.5) if np.isfinite(res).any() else np.nan,
hi=np.nanpercentile(res, 97.5) if np.isfinite(res).any() else np.nan, pc_at_180=cur[:, -1].mean()))
# variant: the phase-tuning cycles inside the sample epoch, as the earlier run, Poisson rates, add-one, random shifts
rng4c = np.random.default_rng([4300, 8, 0, GROUP_ORDER.index(g)])
res, cur = hires.phase_discrimination(CB_x[g][0], CB_x[g][1], rng4c, N=1000, reps=100, model="poisson", shifts="random", smooth="add_one")
pd_rows.append(dict(cycles="phase tuning, sample epoch (variant)", group=g, model="poisson", smooth="add_one", shifts="random",
median_deg=np.nanmedian(res) if np.isfinite(res).any() else np.nan, runs_reaching=int(np.isfinite(res).sum()),
lo=np.nanpercentile(res, 2.5) if np.isfinite(res).any() else np.nan,
hi=np.nanpercentile(res, 97.5) if np.isfinite(res).any() else np.nan, pc_at_180=cur[:, -1].mean()))
# variant: the earlier run's cycles and smoothing, Poisson rates, random shifts
rng4c = np.random.default_rng([4300, 9, 0, GROUP_ORDER.index(g)])
res, cur = hires.phase_discrimination(CB_all[g][0], CB_all[g][1], rng4c, N=1000, reps=100, model="poisson", shifts="random", smooth="spread")
pd_rows.append(dict(cycles="all > 2.5 deg (variant)", group=g, model="poisson", smooth="spread", shifts="random",
median_deg=np.nanmedian(res) if np.isfinite(res).any() else np.nan, runs_reaching=int(np.isfinite(res).sum()),
lo=np.nanpercentile(res, 2.5) if np.isfinite(res).any() else np.nan,
hi=np.nanpercentile(res, 97.5) if np.isfinite(res).any() else np.nan, pc_at_180=cur[:, -1].mean()))
PDC = pd.DataFrame(pd_rows)
print(f"{time.time() - t0:.0f} s")
L4C2 phase-tuning cycles | cycles per cell: median 494, range 92 to 1778; spikes per cycle, mean over cells 0.107; cells with no spike 4 of 31 L4C2 phase-tuning cycles in the sample epoch (variant) | cycles per cell: median 137, range 23 to 804; spikes per cycle, mean over cells 0.100; cells with no spike 6 of 31 L4C2 all cycles > 2.5 deg (variant) | cycles per cell: median 2869, range 680 to 8697; spikes per cycle, mean over cells 0.245; cells with no spike 0 of 31 L4out phase-tuning cycles | cycles per cell: median 966, range 461 to 1481; spikes per cycle, mean over cells 0.163; cells with no spike 0 of 10 L4out phase-tuning cycles in the sample epoch (variant) | cycles per cell: median 283, range 84 to 710; spikes per cycle, mean over cells 0.158; cells with no spike 0 of 10 L4out all cycles > 2.5 deg (variant) | cycles per cell: median 4170, range 2089 to 6651; spikes per cycle, mean over cells 0.179; cells with no spike 0 of 10
130 s
print(PDC.round(3).to_string(index=False))
for (cy, model, smooth, shifts), g in PDC.groupby(["cycles", "model", "smooth", "shifts"]):
a_, b_ = g.set_index("group").loc["L4C2"], g.set_index("group").loc["L4out"]
print(f"{cy} | {model}/{smooth}/{shifts}: C2 {a_.median_deg:.1f} deg (24 +- 2.4: {within(a_.median_deg, 24, 2.4)}), "
f"outside {b_.median_deg:.1f} deg (33 +- 3.3: {within(b_.median_deg, 33, 3.3)}), C2 better: {a_.median_deg < b_.median_deg}")
cycles group model smooth shifts median_deg runs_reaching lo hi pc_at_180
phase tuning L4C2 cycles add_one random NaN 0 NaN NaN 0.532
phase tuning L4C2 cycles add_one even NaN 0 NaN NaN 0.531
phase tuning L4C2 poisson add_one random 102.750 42 12.265 171.900 0.660
phase tuning L4C2 poisson add_one even 98.250 38 6.000 179.207 0.667
phase tuning L4C2 poisson floor random 72.000 91 7.235 161.769 0.697
phase tuning L4C2 poisson floor even 66.000 91 6.000 154.125 0.708
phase tuning, sample epoch (variant) L4C2 poisson add_one random 60.000 73 10.982 168.000 0.669
all > 2.5 deg (variant) L4C2 poisson spread random 58.450 100 17.096 111.059 0.739
phase tuning L4out cycles add_one random NaN 0 NaN NaN 0.582
phase tuning L4out cycles add_one even NaN 0 NaN NaN 0.579
phase tuning L4out poisson add_one random 90.000 91 11.700 176.450 0.699
phase tuning L4out poisson add_one even 90.000 89 11.435 173.520 0.708
phase tuning L4out poisson floor random 41.250 100 6.000 106.181 0.747
phase tuning L4out poisson floor even 36.000 100 6.000 122.100 0.743
phase tuning, sample epoch (variant) L4out poisson add_one random 47.143 100 6.000 154.462 0.716
all > 2.5 deg (variant) L4out poisson spread random 138.000 9 77.867 167.880 0.654
all > 2.5 deg (variant) | poisson/spread/random: C2 58.4 deg (24 +- 2.4: False), outside 138.0 deg (33 +- 3.3: False), C2 better: True
phase tuning | cycles/add_one/even: C2 nan deg (24 +- 2.4: False), outside nan deg (33 +- 3.3: False), C2 better: False
phase tuning | cycles/add_one/random: C2 nan deg (24 +- 2.4: False), outside nan deg (33 +- 3.3: False), C2 better: False
phase tuning | poisson/add_one/even: C2 98.2 deg (24 +- 2.4: False), outside 90.0 deg (33 +- 3.3: False), C2 better: False
phase tuning | poisson/add_one/random: C2 102.8 deg (24 +- 2.4: False), outside 90.0 deg (33 +- 3.3: False), C2 better: False
phase tuning | poisson/floor/even: C2 66.0 deg (24 +- 2.4: False), outside 36.0 deg (33 +- 3.3: True), C2 better: False
phase tuning | poisson/floor/random: C2 72.0 deg (24 +- 2.4: False), outside 41.2 deg (33 +- 3.3: False), C2 better: False
phase tuning, sample epoch (variant) | poisson/add_one/random: C2 60.0 deg (24 +- 2.4: False), outside 47.1 deg (33 +- 3.3: False), C2 better: False
Result: not reproduced. Under every reading C2 is worse than outside C2 and far from 24 deg. With held-out cycle counts and add-one smoothing, neither group reaches 76% correct in any run (0.53 and 0.58 at 180 deg). With Poisson rates and add-one smoothing, C2 reaches it in only 42 (random shifts) and 38 (even shifts) of the 100 runs, at a median of 102.8 and 98.2 deg over those runs, so in most runs C2 never reaches 76% correct; outside C2 reaches it at 90.0 deg (in 91 and 89 runs); with a floor rate, C2 at 72.0 or 66.0 deg (91 runs each) and outside C2 at 41.2 or 36.0 deg. Only the outside value with the floor rate and even shifts (36.0) is within 33 +- 3.3. The readings no longer straddle the paper's value, as they did when the cycles were restricted to the sample epoch; that restriction, now a variant, gives 60.0 deg in C2 and 47.1 outside. The phase-tuning cycles hold 0.107 spikes per cycle in C2 (median 494 cycles per cell, 4 of 31 cells with no spike).
fig5c-modulation-index¶
Touch response magnitude of L4 C2 neurons is strongly modulated by touch order, whisker velocity before touch and maximum curvature during touch, with mean modulation indices of 0.71 +- 0.25, 0.71 +- 0.27 and 0.72 +- 0.23.
Per cell, touches sorted by the property and split into 10 bins of equal numbers; mean spikes per touch in the cell's touch window per bin; index (max - min) / (max + min). Open details from the claims file: velocity at the onset frame or in the frame before; curvature over the whole touch or 0 to 20 ms; touch order by ranks with ties broken at random or by bins of whole order values. Every touch-window reading is run.
def mi_bins(values, counts, n_bins=10, rng=None, whole=False):
"""Mean count per bin of equal numbers of touches sorted by value; returns the bin means."""
values = np.asarray(values, float)
if whole:
u, inv, cnt = np.unique(values, return_inverse=True, return_counts=True)
mid = (np.cumsum(cnt) - cnt / 2) / len(values)
b = np.minimum((mid * n_bins).astype(int), n_bins - 1)[inv]
else:
key = values + (rng.random(len(values)) * 1e-6 if rng is not None else 0)
order = np.argsort(key, kind="stable")
b = np.empty(len(values), int)
b[order] = np.arange(len(values)) * n_bins // len(values)
return np.array([counts[b == k].mean() for k in range(n_bins) if np.any(b == k)])
def mod_index(m):
return (m.max() - m.min()) / (m.max() + m.min()) if m.max() + m.min() > 0 else np.nan
t0 = time.time()
mi_rows, curves = [], []
for no in G["L4C2"]:
c = by_no[no]
d = D[no]
T = d["T"]
props = {("order", "ranks"): T.order.astype(float), ("order", "whole"): T.order.astype(float),
("velocity", "onset"): hires.touch_velocity(c, T, "onset"), ("velocity", "before"): hires.touch_velocity(c, T, "before"),
("curvature", "touch"): hires.touch_max_dkappa(c, T, "touch"), ("curvature", "0-20"): hires.touch_max_dkappa(c, T, "0-20")}
for key, r in W[(W.cell == no) & (W.lat_reading == "a")].set_index(WKEY[1:]).iterrows():
n = hires.spikes_per_touch(c, d["s"], T, r.latency, r.t_on)
for (prop, rd), v in props.items():
rng5 = np.random.default_rng([5100, no])
m = mi_bins(v, n, rng=rng5, whole=(rd == "whole"))
mi_rows.append(dict(cell=no, **dict(zip(WKEY[1:], key)), prop=prop, reading=rd, mi=mod_index(m), n_bins=len(m)))
curves.append(dict(cell=no, **dict(zip(WKEY[1:], key)), prop=prop, reading=rd, curve=m / m.max() if m.max() > 0 else m))
MI = pd.DataFrame(mi_rows)
print(f"{time.time() - t0:.0f} s; bins per curve: {MI.n_bins.value_counts().to_dict()}")
1 s; bins per curve: {10: 2864, 9: 64, 8: 48}
mis = MI.groupby(WKEY[1:] + ["prop", "reading"]).mi.agg(["mean", "std", "count"]).reset_index()
exp5 = {"order": (0.71, 0.25), "velocity": (0.71, 0.27), "curvature": (0.72, 0.23)}
mis["ok"] = [within(r["mean"], *exp5[r.prop]) for _, r in mis.iterrows()]
print(mis.groupby(["prop", "reading"])["mean"].agg(["min", "max"]).round(3).to_string())
print(mis.groupby(["prop", "reading", "start"]).ok.agg(["sum", "size"]).to_string())
print("range of the means by start of the cutoff search:")
print(mis.groupby(["prop", "start"])["mean"].agg(["min", "max"]).round(3).to_string())
sel = (mis.shuffle == "uniform") & (mis.nearest == "before") & (mis.tmax == 100) & (mis.start == "cross")
print("first reading of each open detail (uniform, before, 100 ms, transition cutoff; ranks, onset, whole touch):",
{p: round(float(mis[sel & (mis.prop == p) & (mis.reading == r)]["mean"].iloc[0]), 3) for p, r in (("order", "ranks"), ("velocity", "onset"), ("curvature", "touch"))})
print("cells with an undefined index (no spikes in any bin):", MI[MI.mi.isna()].groupby(["prop"]).cell.unique().apply(list).to_dict())
min max
prop reading
curvature 0-20 0.686 0.848
touch 0.664 0.832
order ranks 0.666 0.854
whole 0.676 0.871
velocity before 0.712 0.894
onset 0.703 0.861
sum size
prop reading start
curvature 0-20 cross 8 8
first 8 8
touch cross 8 8
first 8 8
order ranks cross 8 8
first 8 8
whole cross 8 8
first 8 8
velocity before cross 8 8
first 8 8
onset cross 8 8
first 8 8
range of the means by start of the cutoff search:
min max
prop start
curvature cross 0.664 0.698
first 0.759 0.848
order cross 0.666 0.692
first 0.744 0.871
velocity cross 0.703 0.735
first 0.776 0.894
first reading of each open detail (uniform, before, 100 ms, transition cutoff; ranks, onset, whole touch): {'order': 0.677, 'velocity': 0.713, 'curvature': 0.672}
cells with an undefined index (no spikes in any bin): {}
Result: reproduced. Under the transition cutoff every reading of the window, velocity, curvature and touch order gives mean modulation indices within tolerance: touch order 0.666 to 0.692, velocity 0.703 to 0.735 and curvature 0.664 to 0.698 (paper 0.71 +- 0.25, 0.71 +- 0.27, 0.72 +- 0.23); under the first reading of each detail 0.677, 0.713 and 0.672. The literal cutoff, a variant, gives 0.744 to 0.894, also within tolerance.
fig5-modulation-correlations¶
Across cells, responsiveness to pre-touch velocity and to maximum curvature covary strongly (r = 0.80), and responsiveness to touch order is somewhat less correlated with both (r = 0.70 and 0.65).
Readings: (a) Pearson correlation across the 31 cells of the modulation indices; (b) Pearson correlation of the peak-normalised 10-bin curves pooled over cells and bins. Under every reading of fig5c.
CV = pd.DataFrame(curves)
cor_rows = []
for key, g in MI.groupby(WKEY[1:]):
for ro in ("ranks", "whole"):
for rv in ("onset", "before"):
for rc in ("touch", "0-20"):
sel = {"order": ro, "velocity": rv, "curvature": rc}
x = {p: g[(g.prop == p) & (g.reading == sel[p])].set_index("cell").mi for p in sel}
df = pd.DataFrame(x).dropna()
cv = CV[(CV.shuffle == key[0]) & (CV.nearest == key[1]) & (CV.tmax == key[2]) & (CV.start == key[3])]
cc = {p: cv[(cv.prop == p) & (cv.reading == sel[p])].set_index("cell").curve for p in sel}
common = [k for k in cc["order"].index if all(len(cc[p][k]) == 10 for p in sel)]
pooled = {p: np.concatenate([cc[p][k] for k in common]) for p in sel}
for rd, data in (("(a) indices", df), ("(b) curves", pd.DataFrame(pooled).dropna())):
r_vc, p_vc = stats.pearsonr(data.velocity, data.curvature)
r_ov, p_ov = stats.pearsonr(data.order, data.velocity)
r_oc, p_oc = stats.pearsonr(data.order, data.curvature)
cor_rows.append(dict(**dict(zip(WKEY[1:], key)), order=ro, velocity=rv, curvature=rc, reading=rd, n=len(data),
r_vc=r_vc, r_ov=r_ov, r_oc=r_oc, p_max=max(p_vc, p_ov, p_oc),
ok=within(r_vc, 0.80, 0.08) and within(r_ov, 0.70, 0.07) and within(r_oc, 0.65, 0.065)
and max(p_vc, p_ov, p_oc) < 0.05 and r_vc >= max(r_ov, r_oc)))
COR = pd.DataFrame(cor_rows)
print(COR.groupby(["reading", "start"])[["r_vc", "r_ov", "r_oc"]].agg(["min", "max"]).round(3).to_string())
print(COR.groupby(["reading", "start"]).ok.agg(["sum", "size"]).to_string())
r_vc r_ov r_oc
min max min max min max
reading start
(a) indices cross 0.618 0.788 0.706 0.788 0.588 0.641
first 0.549 0.783 0.621 0.830 0.534 0.681
(b) curves cross 0.272 0.508 0.012 0.159 -0.015 0.162
first 0.250 0.441 0.114 0.253 0.071 0.250
sum size
reading start
(a) indices cross 16 64
first 8 64
(b) curves cross 0 64
first 0 64
Result: blocked (methods_underspecified). Under reading (a), correlations of the indices across cells, with the transition cutoff r is 0.618 to 0.788 for velocity and curvature, 0.706 to 0.788 for order and velocity, and 0.588 to 0.641 for order and curvature, and the criterion holds in 16 of 64 combinations of the fig5c readings; where it fails, at least one r is outside its tolerance or velocity and curvature is not the largest. Under reading (b), pooled curves, r is 0.27 to 0.51 for velocity and curvature and the criterion never holds. The fig5c readings give different outcomes. The literal cutoff, a variant, holds in 8 of 64 under (a).
fig6s1b-clusters-less-variable¶
The clustering of touches into five groups finds subsets with significantly less variability in touch properties than random samples of touches.
Per L4 C2 cell, touches with inter-contact interval > 250 ms; z-scored pre-touch velocity (the frame before onset, fig5c's "pretouch" reading) and maximum curvature change 0 to 20 ms after onset; OPTICS ordering; five groups of N touches, N = 20% of the selected touches and at least 20. Open details from the claims file: OPTICS minimum points 5 or N; interval from the previous offset or onset. The first group's s.d. of z-scored velocity is compared with 1000 random samples of N of the cell's selected touches.
t0 = time.time()
GRP = {}
cl_rows = []
for oi, ici_from in enumerate(("offset", "onset")):
for mi_, minpts in enumerate(("5", "N")):
for no in G["L4C2"]:
c = by_no[no]
T = D[no]["T"]
sel = np.flatnonzero(hires.ici(T, c, ici_from) > 250)
v = hires.touch_velocity(c, T, "before")[sel]
k = hires.touch_max_dkappa(c, T, "0-20")[sel]
X = np.c_[stats.zscore(v), stats.zscore(k)]
X = np.nan_to_num(X)
N = max(20, int(round(0.2 * len(sel))))
groups = hires.optics_groups(X, N, 5, 5 if minpts == "5" else N) if len(sel) >= N else []
GRP[(ici_from, minpts, no)] = [sel[g] for g in groups]
rng6 = np.random.default_rng([6100, oi, mi_, no])
if groups:
g1 = X[groups[0], 0].std(ddof=1)
rnd = np.array([X[rng6.choice(len(X), N, replace=False), 0].std(ddof=1) for _ in range(1000)])
p5 = np.percentile(rnd, 5)
else:
g1, p5 = np.nan, np.nan
cl_rows.append(dict(ici_from=ici_from, minpts=minpts, cell=no, n_selected=len(sel), N=N, n_groups=len(groups),
sd_group1=g1, random_p5=p5, less=bool(g1 < p5) if groups else False))
CL = pd.DataFrame(cl_rows)
print(f"{time.time() - t0:.0f} s")
40 s
print(CL.groupby(["ici_from", "minpts"]).agg(cells_with_groups=("n_groups", lambda x: int((x > 0).sum())), five_groups=("n_groups", lambda x: int((x == 5).sum())),
less=("less", "sum"), n=("less", "size"), median_selected=("n_selected", "median")).to_string())
for key, g in CL.groupby(["ici_from", "minpts"]):
print(key, "fraction of the 31 cells with the first group below the 5th percentile: %.2f" % g.less.mean())
cells_with_groups five_groups less n median_selected
ici_from minpts
offset 5 28 8 28 31 84.0
N 28 8 28 31 84.0
onset 5 29 8 28 31 96.0
N 29 8 28 31 96.0
('offset', '5') fraction of the 31 cells with the first group below the 5th percentile: 0.90
('offset', 'N') fraction of the 31 cells with the first group below the 5th percentile: 0.90
('onset', '5') fraction of the 31 cells with the first group below the 5th percentile: 0.90
('onset', 'N') fraction of the 31 cells with the first group below the 5th percentile: 0.90
Result: reproduced. Under all four readings, 28 of the 31 L4 C2 cells (0.90) have a first group whose s.d. of z-scored velocity is below the 5th percentile of 1000 random samples of the same size; 28 or 29 cells have enough touches for at least one group, and 8 for all five.
example-cell-ff-around-touch¶
For the example neuron, the 10 ms sliding-window Fano factor is about 1 before and after touch, where the mean spike count is about 0.03, and falls near the peak of the touch response, where the mean count is about 0.7.
Windows of 10 ms slid in 1 ms steps from -50 to +40 ms around onset, across touches with inter-contact interval > 250 ms (from the previous offset or onset). "Before touch" is read as the windows that end at or before onset (-50 to -10 ms starts), a standard choice; the peak is the window of highest mean count. Readings of the touches: (a) pooled, (b) each of the five groups (every clustering reading). Readings of the cell: the Figure 1 cell, or any L4 C2 cell.
starts = np.arange(-50, 41)
pre = starts <= -10
def ff_profile(c, s, onsets):
X = hires.window_counts(c, s, onsets, starts, 10)
m = X.mean(0)
v = X.var(0, ddof=1)
with np.errstate(invalid="ignore", divide="ignore"):
ff = v / m
pk = int(np.argmax(m))
return dict(n=len(onsets), m_pre=m[pre].mean(), ff_pre=np.nanmean(ff[pre]), m_peak=m[pk], ff_peak=ff[pk], peak_start=int(starts[pk]))
rows = []
for no in G["L4C2"]:
c = by_no[no]
T = D[no]["T"]
for ici_from in ("offset", "onset"):
sel = np.flatnonzero(hires.ici(T, c, ici_from) > 250)
rows.append(dict(cell=no, touches="pooled", ici_from=ici_from, minpts="-", group=0, **ff_profile(c, D[no]["s"], T.on[sel])))
for minpts in ("5", "N"):
for gi, gsel in enumerate(GRP[(ici_from, minpts, no)]):
rows.append(dict(cell=no, touches="group", ici_from=ici_from, minpts=minpts, group=gi + 1, **ff_profile(c, D[no]["s"], T.on[gsel])))
FFP = pd.DataFrame(rows)
FFP["ok"] = [within(r.m_pre, 0.03, 0.003) and within(r.m_peak, 0.7, 0.07) and abs(r.ff_pre - 1) <= 0.1 and r.ff_peak < 1 for r in FFP.itertuples()]
print(FFP[FFP.cell == EXC].round(3).to_string(index=False))
gr = FFP[FFP.touches == "group"]
print("units tested: pooled, %d (31 cells x 2 interval readings); grouped, %d cell x group units over the 4 clustering readings, from %d cells" %
((FFP.touches == "pooled").sum(), len(gr), gr.cell.nunique()))
print("rows that match:")
print(FFP[FFP.ok].round(4).to_string(index=False))
print("L4 C2 cells matching under pooled touches:", sorted(FFP[(FFP.touches == "pooled") & FFP.ok].cell.unique().tolist()),
"| under a group:", sorted(FFP[(FFP.touches == "group") & FFP.ok].cell.unique().tolist()))
cell touches ici_from minpts group n m_pre ff_pre m_peak ff_peak peak_start ok 23 pooled offset - 0 156 0.003 1.000 0.641 1.428 12 False 23 group offset 5 1 31 0.000 NaN 0.548 1.075 10 False 23 group offset 5 2 31 0.008 1.000 0.806 1.109 12 False 23 group offset 5 3 31 0.008 1.000 0.484 1.498 13 False 23 group offset 5 4 31 0.000 NaN 0.839 1.597 12 False 23 group offset 5 5 31 0.000 NaN 0.710 1.803 12 False 23 group offset N 1 31 0.000 NaN 0.516 1.146 10 False 23 group offset N 2 31 0.008 1.000 0.806 0.613 14 False 23 group offset N 3 31 0.008 1.000 0.419 1.395 13 False 23 group offset N 4 31 0.000 NaN 0.290 1.652 12 False 23 group offset N 5 31 0.000 NaN 1.355 1.257 12 False 23 pooled onset - 0 193 0.004 0.999 0.539 1.507 12 False 23 group onset 5 1 39 0.000 NaN 0.436 1.183 10 False 23 group onset 5 2 39 0.000 NaN 0.154 0.868 13 False 23 group onset 5 3 39 0.006 1.000 0.718 1.096 12 False 23 group onset 5 4 39 0.006 1.000 1.205 1.274 12 False 23 group onset N 1 39 0.000 NaN 0.256 1.379 13 False 23 group onset N 2 39 0.000 NaN 0.769 1.058 10 False 23 group onset N 3 39 0.006 1.000 0.462 1.123 13 False 23 group onset N 4 39 0.006 1.000 0.872 1.822 12 False units tested: pooled, 62 (31 cells x 2 interval readings); grouped, 394 cell x group units over the 4 clustering readings, from 29 cells rows that match: cell touches ici_from minpts group n m_pre ff_pre m_peak ff_peak peak_start ok 11 group offset N 4 20 0.0305 0.99 0.65 0.3684 5 True 11 group onset 5 3 20 0.0280 1.00 0.70 0.3158 5 True L4 C2 cells matching under pooled touches: [] | under a group: [11]
Result: blocked (methods_underspecified). The Figure 1 cell (23) matches under no reading: its mean 10 ms count before touch is 0.000 to 0.008 (paper about 0.03), its peak mean count 0.15 to 1.36, and its FF at the peak 0.61 to 1.82. With pooled touches no L4 C2 cell matches (62 units); with grouped touches cell 11 does, in 2 of 394 cell and group units (groups of 20 touches: mean counts 0.0305 and 0.0280 before touch and 0.65 and 0.70 at the peak, FF 0.99 and 1.00 before and 0.37 and 0.32 at the peak). The two readings of the touches give different outcomes under the criterion's rule for other cells.
fig6cd-ff-near-binomial¶
Across L4 C2 neurons, the Fano factor following touch drops below one, lying close to the binomial limit (the theoretical minimum) for most cells.
Per cell and touch group, the FF in 10 ms windows starting 0 to 20 ms after onset; its minimum per group (1a) or across groups (1b). "Close": (2a) the binomial limit lies inside the 95% interval of 1000 resamples of the touches, or (2b) the FF is nearer the binomial limit than 1. The binomial limit for a mean count m is f(1 - f)/m with f the fractional part of m (1 - m for m < 1). For (1a) each cell and group counts as one unit.
starts20 = np.arange(0, 21)
t0 = time.time()
gff_rows = []
for oi, ici_from in enumerate(("offset", "onset")):
for mi_, minpts in enumerate(("5", "N")):
for no in G["L4C2"]:
c = by_no[no]
T = D[no]["T"]
for gi, gsel in enumerate(GRP[(ici_from, minpts, no)]):
X = hires.window_counts(c, D[no]["s"], T.on[gsel], starts20, 10)
m = X.mean(0)
with np.errstate(invalid="ignore", divide="ignore"):
ff = X.var(0, ddof=1) / m
j = int(np.nanargmin(ff)) if np.isfinite(ff).any() else 0
rngb = np.random.default_rng([6300, oi, mi_, no, gi])
idx = rngb.integers(0, len(X), (1000, len(X)))
xb = X[:, j][idx]
with np.errstate(invalid="ignore", divide="ignore"):
ffb = xb.var(1, ddof=1) / xb.mean(1)
lo, hi = np.nanpercentile(ffb, [2.5, 97.5]) if np.isfinite(ffb).any() else (np.nan, np.nan)
bl = float(hires.binomial_limit(m[j]))
gff_rows.append(dict(ici_from=ici_from, minpts=minpts, cell=no, group=gi + 1, n=len(gsel), start=int(starts20[j]),
m=m[j], ff=ff[j], binom=bl, lo=lo, hi=hi, ff_by_start=ff, m_by_start=m))
GFF = pd.DataFrame(gff_rows)
GFF["close_a"] = (GFF.binom >= GFF.lo) & (GFF.binom <= GFF.hi)
GFF["close_b"] = (GFF.ff - GFF.binom).abs() < (GFF.ff - 1).abs()
print(f"{time.time() - t0:.0f} s; groups with no spikes in any window (FF undefined): {int(GFF.ff.isna().sum())} of {len(GFF)}")
0 s; groups with no spikes in any window (FF undefined): 6 of 394
rows = []
for (ici_from, minpts), g in GFF.groupby(["ici_from", "minpts"]):
g = g.dropna(subset=["ff"])
best = g.loc[g.groupby("cell").ff.idxmin()]
for which, u in (("(1a) per group", g), ("(1b) across groups", best)):
for close in ("close_a", "close_b"):
rows.append(dict(ici_from=ici_from, minpts=minpts, minimum=which, close=close, n_units=len(u), mean_min_ff=u.ff.mean(),
fraction_close=u[close].mean(), ok=(u.ff.mean() < 1) and (u[close].mean() > 0.5)))
F6 = pd.DataFrame(rows)
print(F6.round(3).to_string(index=False))
print("combinations meeting the criterion: %d of %d" % (F6.ok.sum(), len(F6)))
ici_from minpts minimum close n_units mean_min_ff fraction_close ok offset 5 (1a) per group close_a 92 0.643 0.880 True offset 5 (1a) per group close_b 92 0.643 0.891 True offset 5 (1b) across groups close_a 28 0.504 0.821 True offset 5 (1b) across groups close_b 28 0.504 0.964 True offset N (1a) per group close_a 93 0.652 0.871 True offset N (1a) per group close_b 93 0.652 0.903 True offset N (1b) across groups close_a 28 0.509 0.821 True offset N (1b) across groups close_b 28 0.509 0.964 True onset 5 (1a) per group close_a 101 0.657 0.851 True onset 5 (1a) per group close_b 101 0.657 0.901 True onset 5 (1b) across groups close_a 29 0.522 0.828 True onset 5 (1b) across groups close_b 29 0.522 1.000 True onset N (1a) per group close_a 102 0.676 0.863 True onset N (1a) per group close_b 102 0.676 0.882 True onset N (1b) across groups close_a 29 0.548 0.793 True onset N (1b) across groups close_b 29 0.548 0.966 True combinations meeting the criterion: 16 of 16
Result: reproduced. Under all 16 combinations of readings the population mean of the minimum FF after touch is below 1 (0.504 to 0.676) and the FF is close to the binomial limit in more than half of the units (0.79 to 1.00).
fig6d-not-refractory¶
The low Fano factor after touch cannot be explained by a Poisson process with a 2.3 ms refractory period and rates matched to each neuron.
Per cell, spike trains simulated at 0.1 ms resolution as a Poisson process with an absolute refractory period of 2.3 ms, one train per real touch; the FF in the same 10 ms windows. After the answer to the question, the model's average touch response is matched to the data: its rate is the free rate of Berry and Meister (1998), lambda(t) = r(t) / (1 - integral of r over [t - 2.3 ms, t)), with r the cell's touch PSTH in 1 ms bins from the same touches, starting 10 ms before onset so that the integral is defined in the first window. The PSTH used directly as the rate, the earlier run's model, is shown as a variant. Open details from the claims file: (1) touches pooled (> 250 ms touches) or per group; (2) the comparison, (a) the population mean data FF below the model's at the window of the data minimum by more than the sum of the two s.e.m., or (b) the per-cell data FF below the model's in more than half the cells, sign test p < 0.05.
def free_rate(rate_ms, dead=2.3, dt=0.1, cap=0.95):
"""Berry and Meister (1998): the rate of the process before the dead time, lambda(t) = r(t) /
(1 - integral of r over [t - dead, t)), so that the output rate is r. r is given per 1 ms bin
and spread over steps of dt; the integral is capped at `cap` where it would reach 1."""
r = np.repeat(rate_ms, int(round(1 / dt)))
k = int(round(dead / dt))
c = np.r_[0, np.cumsum(r * dt)]
i = np.arange(len(r))
integ = c[i] - c[np.maximum(i - k, 0)]
return r / (1 - np.minimum(integ, cap)), int((integ > cap).sum())
def simulate_refractory(rate_fine, n_trains, rng, dead=2.3, dt=0.1):
"""Spike times (ms) of n_trains inhomogeneous Poisson trains with an absolute dead time;
rate_fine is in spikes/ms per step of dt."""
p = rate_fine * dt
last = np.full(n_trains, -np.inf)
out = [[] for _ in range(n_trains)]
for i, pi in enumerate(p):
t = i * dt
fire = (rng.random(n_trains) < pi) & (t - last >= dead)
last[fire] = t
for k in np.flatnonzero(fire):
out[k].append(t)
return out
def ff_windows(trains, starts, w=10):
X = np.array([[np.sum((np.asarray(tr) >= s) & (np.asarray(tr) < s + w)) for s in starts] for tr in trains])
m = X.mean(0)
with np.errstate(invalid="ignore", divide="ignore"):
return X.var(0, ddof=1) / m
PRE6D = 10 # ms of PSTH before onset, so that the integral of the free rate starts before the first window
t0 = time.time()
ref_rows = []
for oi, ici_from in enumerate(("offset", "onset")):
sets = {("pooled", "-"): None}
for minpts in ("5", "N"):
sets[("group", minpts)] = minpts
for (touches, minpts), _ in sets.items():
for no in G["L4C2"]:
c = by_no[no]
T = D[no]["T"]
s = D[no]["s"]
units = [np.flatnonzero(hires.ici(T, c, ici_from) > 250)] if touches == "pooled" else GRP[(ici_from, minpts, no)]
for gi, sel in enumerate(units):
on = T.on[sel]
X = hires.window_counts(c, s, on, np.arange(-PRE6D, 31), 1) # 1 ms PSTH, -10 to 30 ms
rate = X.mean(0)
Xd = hires.window_counts(c, s, on, starts20, 10)
with np.errstate(invalid="ignore", divide="ignore"):
ff_data = Xd.var(0, ddof=1) / Xd.mean(0)
lam, n_cap = free_rate(rate)
for mi, (model, rf) in enumerate((("compensated", lam), ("uncompensated", np.repeat(rate, 10)))):
rngr = np.random.default_rng([6400, oi, 0 if touches == "pooled" else (1 if minpts == "5" else 2), no, gi, mi])
sim = [list(np.asarray(tr) - PRE6D) for tr in simulate_refractory(rf, len(on), rngr)]
Xm = np.array([[np.sum((np.asarray(tr) >= a) & (np.asarray(tr) < a + 10)) for a in starts20] for tr in sim])
ref_rows.append(dict(ici_from=ici_from, touches=touches, minpts=minpts, cell=no, group=gi + 1, model=model, ff_data=ff_data,
ff_model=ff_windows(sim, starts20), mean_data=Xd.mean(0), mean_model=Xm.mean(0), n_capped=n_cap))
REF = pd.DataFrame(ref_rows)
print(f"{time.time() - t0:.0f} s")
17 s
rows = []
for (model, ici_from, touches, minpts), g in REF.groupby(["model", "ici_from", "touches", "minpts"]):
Dd = np.vstack(g.ff_data.to_numpy())
Mm = np.vstack(g.ff_model.to_numpy())
md, mm = np.nanmean(Dd, 0), np.nanmean(Mm, 0)
sd_ = np.nanstd(Dd, 0, ddof=1) / np.sqrt(np.sum(np.isfinite(Dd), 0))
sm_ = np.nanstd(Mm, 0, ddof=1) / np.sqrt(np.sum(np.isfinite(Mm), 0))
j = int(np.nanargmin(md))
a_ok = (mm[j] - md[j]) > (sd_[j] + sm_[j])
jc = np.array([np.nanargmin(r) if np.isfinite(r).any() else -1 for r in Dd])
ok = jc >= 0
below = np.array([Dd[i, jc[i]] < Mm[i, jc[i]] for i in np.flatnonzero(ok)])
below = below[np.isfinite(Mm[np.flatnonzero(ok), jc[ok]])]
k, n = int(below.sum()), len(below)
p = stats.binomtest(k, n, 0.5, alternative="greater").pvalue if n else np.nan
rows.append(dict(model=model, ici_from=ici_from, touches=touches, minpts=minpts, units=len(g), start_of_min=int(starts20[j]), data_ff=md[j], model_ff=mm[j],
sem_sum=sd_[j] + sm_[j], a_ok=a_ok, cells_below=k, n=n, sign_p=p, b_ok=(k > n / 2) and p < 0.05))
F6D = pd.DataFrame(rows)
print(F6D.round(4).to_string(index=False))
# does the model's mean count match the data's? ratio of the population mean counts in the 10 ms windows
for model, g in REF.groupby("model"):
md = np.vstack(g.mean_data.to_numpy()).mean(0)
mm = np.vstack(g.mean_model.to_numpy()).mean(0)
print(f"{model}: model / data mean count over the 21 windows, {np.min(mm / md):.3f} to {np.max(mm / md):.3f} (at the 0 to 10 ms window {mm[0] / md[0]:.3f});"
f" steps where the free-rate integral was capped: {int(g.n_capped.sum())}")
model ici_from touches minpts units start_of_min data_ff model_ff sem_sum a_ok cells_below n sign_p b_ok compensated offset group 5 94 3 0.7466 0.7445 0.0601 False 61 89 0.0003 True compensated offset group N 94 4 0.7530 0.7133 0.0588 False 52 91 0.1041 False compensated offset pooled - 31 3 0.7912 0.7633 0.0855 False 17 31 0.3601 False compensated onset group 5 103 3 0.7657 0.7775 0.0625 False 61 100 0.0176 True compensated onset group N 103 3 0.7755 0.7155 0.0546 False 59 99 0.0350 True compensated onset pooled - 31 3 0.8190 0.7862 0.0870 False 22 31 0.0147 True uncompensated offset group 5 94 3 0.7466 0.7600 0.0547 False 66 86 0.0000 True uncompensated offset group N 94 4 0.7530 0.7782 0.0521 False 67 89 0.0000 True uncompensated offset pooled - 31 3 0.7912 0.7766 0.0777 False 22 31 0.0147 True uncompensated onset group 5 103 3 0.7657 0.7687 0.0528 False 77 97 0.0000 True uncompensated onset group N 103 3 0.7755 0.8135 0.0467 False 80 102 0.0000 True uncompensated onset pooled - 31 3 0.8190 0.8031 0.0751 False 21 31 0.0354 True compensated: model / data mean count over the 21 windows, 1.001 to 1.034 (at the 0 to 10 ms window 1.001); steps where the free-rate integral was capped: 0 uncompensated: model / data mean count over the 21 windows, 0.796 to 0.891 (at the 0 to 10 ms window 0.807); steps where the free-rate integral was capped: 0
Result: blocked (methods_underspecified). The compensated model matches the data's mean count (model over data 1.001 to 1.034 across the windows; 0.80 to 0.89 without compensation) and its FF at the data minimum is 0.71 to 0.79. Under reading (2a) the population mean data FF (0.75 to 0.82) is not below the model's by more than the sum of the two s.e.m. in any of the 6 touch readings. Under reading (2b) the per-cell data FF is below the model's in 4 of the 6 (61 of 89 units, 61 of 100, 59 of 99, and 22 of 31 pooled cells counted from onsets, sign test p = 0.0003 to 0.035), but not with groups of N points counted from offsets (52 of 91, p = 0.10) or pooled touches counted from offsets (17 of 31, p = 0.36). The readings disagree. The uncompensated variant, whose FF is higher, gives (2b) in all 6.
Figure 6E: fig6e-sample-epoch-ff, fig6e-random-38ms-ff, fig6e-touch-38ms-ff, fig6e-touch-10ms-ff, fig6e-similar-touch-ff¶
Counting spikes over the whole sample epoch, ignoring the millisecond structure of behavior, gives a very large Fano factor of 6.57 +- 5.41 (mean +- s.d. over L4 C2 cells).
Counting spikes in random 38 ms windows, as many per trial as there are touches, gives a Fano factor of 1.59 +- 0.49.
Counting spikes in 38 ms windows aligned to touch onset (after a 6 ms latency) gives a Fano factor of 1.37 +- 0.74.
A narrower 10 ms window after touch reduces the Fano factor to 0.73 +- 0.19.
Computed only across similar touches, the 10 ms Fano factor falls further to 0.54 +- 0.27, so that L4 responses to touch are almost noiseless.
x = ffs[G["L4C2"]]
print(f"sample-epoch FF over the 31 L4 C2 cells: mean {x.mean():.2f}, s.d. {x.std():.2f}, median {x.median():.2f}, range {x.min():.2f} to {x.max():.2f}")
print("within 6.57 +- 5.41 and above 1:", within(x.mean(), 6.57, 5.41) and x.mean() > 1)
sample-epoch FF over the 31 L4 C2 cells: mean 6.58, s.d. 5.40, median 4.50, range 1.41 to 21.78 within 6.57 +- 5.41 and above 1: True
Result: reproduced. The sample-epoch FF over the 31 L4 C2 cells is 6.58 +- 5.40 (paper 6.57 +- 5.41).
rng6e = np.random.default_rng(6500)
rr = {"(a) sample epoch": [], "(b) whole trial": []}
for no in G["L4C2"]:
c = by_no[no]
T = D[no]["T"]
s = np.sort(D[no]["s"])
ntouch = np.bincount(T.trial, minlength=len(c.trials))
e = EP_by[no].set_index("trial")
nseg = np.diff(np.r_[c.seg_start, len(c.vt)])
for rd in rr:
ffd = []
for _ in range(100):
st = []
for k in np.flatnonzero(ntouch):
if rd.startswith("(a)"):
if k not in e.index or e.loc[k, "off_ms"] - e.loc[k, "on_ms"] < 38:
continue
lo, hi = e.loc[k, "on_ms"], e.loc[k, "off_ms"] - 38
else:
lo, hi = 0, nseg[k] - 38
st.append(k * hires.TRIAL_STRIDE_MS + rng6e.uniform(lo, hi, ntouch[k]))
st = np.concatenate(st)
ffd.append(hires.fano(hires.counts_in(s, st, 38)))
rr[rd].append(np.nanmean(ffd))
for rd, v in rr.items():
v = np.array(v)
print(f"{rd}: mean {np.nanmean(v):.2f}, s.d. {np.nanstd(v, ddof=1):.2f}; within 1.59 +- 0.49: {within(np.nanmean(v), 1.59, 0.49)}")
FF2 = {rd: np.nanmean(v) for rd, v in rr.items()}
(a) sample epoch: mean 1.56, s.d. 0.47; within 1.59 +- 0.49: True (b) whole trial: mean 1.59, s.d. 0.47; within 1.59 +- 0.49: True
Result: reproduced. Random 38 ms windows give 1.56 +- 0.47 within the sample epoch and 1.59 +- 0.47 anywhere in the trial (paper 1.59 +- 0.49).
FF3, FF4 = {}, {}
for no in G["L4C2"]:
c = by_no[no]
T = D[no]["T"]
s = D[no]["s"]
for tr in ("all", "ici>250 offset", "ici>250 onset"):
sel = np.arange(len(T.on)) if tr == "all" else np.flatnonzero(hires.ici(T, c, tr.split()[1]) > 250)
X38 = hires.window_counts(c, s, T.on[sel], np.array([6.0]), 38)[:, 0]
FF3.setdefault(tr, []).append(hires.fano(X38))
X10 = hires.window_counts(c, s, T.on[sel], starts20, 10)
with np.errstate(invalid="ignore", divide="ignore"):
ff = X10.var(0, ddof=1) / X10.mean(0)
for rd, lim in (("0-20", 20), ("0-10", 10)):
v = ff[starts20 <= lim]
FF4.setdefault((tr, rd), []).append(np.nanmin(v) if np.isfinite(v).any() else np.nan)
for tr, v in FF3.items():
v = np.array(v)
print(f"touches {tr}: FF 6-44 ms mean {np.nanmean(v):.2f}, s.d. {np.nanstd(v, ddof=1):.2f}, n {np.isfinite(v).sum()}; within 1.37 +- 0.74: {within(np.nanmean(v), 1.37, 0.74)}")
touches all: FF 6-44 ms mean 1.39, s.d. 0.60, n 31; within 1.37 +- 0.74: True touches ici>250 offset: FF 6-44 ms mean 1.34, s.d. 0.74, n 31; within 1.37 +- 0.74: True touches ici>250 onset: FF 6-44 ms mean 1.36, s.d. 0.74, n 31; within 1.37 +- 0.74: True
Result: reproduced. Windows 6 to 44 ms after touch give 1.39 +- 0.60 over all touches and 1.34 or 1.36 +- 0.74 over touches with an interval above 250 ms (paper 1.37 +- 0.74).
for key, v in FF4.items():
v = np.array(v)
print(f"touches {key[0]}, window starts {key[1]} ms: minimum FF mean {np.nanmean(v):.2f}, s.d. {np.nanstd(v, ddof=1):.2f}, n {np.isfinite(v).sum()}; within 0.73 +- 0.19: {within(np.nanmean(v), 0.73, 0.19)}")
touches all, window starts 0-20 ms: minimum FF mean 0.85, s.d. 0.13, n 31; within 0.73 +- 0.19: True touches all, window starts 0-10 ms: minimum FF mean 0.86, s.d. 0.13, n 31; within 0.73 +- 0.19: True touches ici>250 offset, window starts 0-20 ms: minimum FF mean 0.70, s.d. 0.20, n 31; within 0.73 +- 0.19: True touches ici>250 offset, window starts 0-10 ms: minimum FF mean 0.71, s.d. 0.21, n 31; within 0.73 +- 0.19: True touches ici>250 onset, window starts 0-20 ms: minimum FF mean 0.72, s.d. 0.19, n 31; within 0.73 +- 0.19: True touches ici>250 onset, window starts 0-10 ms: minimum FF mean 0.74, s.d. 0.20, n 31; within 0.73 +- 0.19: True
Result: reproduced. The minimum 10 ms FF after touch is 0.85 or 0.86 over all touches and 0.70 to 0.74 over touches with an interval above 250 ms, within 0.73 +- 0.19 under all six readings.
rows = []
for (ici_from, minpts), g in GFF.groupby(["ici_from", "minpts"]):
for rd, lim in (("0-20", 20), ("0-10", 10)):
per = []
for no, gg in g.groupby("cell"):
mins = [np.nanmin(f[starts20 <= lim]) if np.isfinite(f[starts20 <= lim]).any() else np.nan for f in gg.ff_by_start]
per.append((np.nanmin(mins) if np.isfinite(mins).any() else np.nan, np.nanmean(mins)))
per = np.array(per)
for vi, val in enumerate(("(a) min across groups", "(b) mean of group minima")):
v = per[:, vi]
m5 = np.nanmean(v)
ff1, ff2 = ffs[G["L4C2"]].mean(), FF2["(a) sample epoch"]
ff3 = np.nanmean(FF3[f"ici>250 {ici_from}"])
ff4 = np.nanmean(FF4[(f"ici>250 {ici_from}", rd)])
rows.append(dict(ici_from=ici_from, minpts=minpts, starts=rd, value=val, mean=m5, sd=np.nanstd(v, ddof=1), n=int(np.isfinite(v).sum()),
ok_value=within(m5, 0.54, 0.27), ordered=bool(ff1 > ff2 > ff3 > ff4 > m5), ff1=ff1, ff2=ff2, ff3=ff3, ff4=ff4))
F6E = pd.DataFrame(rows)
F6E["ok"] = F6E.ok_value & F6E.ordered
print(F6E.round(3).to_string(index=False))
print("combinations meeting the criterion: %d of %d" % (F6E.ok.sum(), len(F6E)))
print("order with the whole-trial random windows (2b) and all touches (3a, 4a):",
[round(v, 2) for v in (ffs[G["L4C2"]].mean(), FF2["(b) whole trial"], np.nanmean(FF3["all"]), np.nanmean(FF4[("all", "0-20")]))])
ici_from minpts starts value mean sd n ok_value ordered ff1 ff2 ff3 ff4 ok offset 5 0-20 (a) min across groups 0.504 0.233 28 True True 6.58 1.557 1.337 0.697 True offset 5 0-20 (b) mean of group minima 0.628 0.233 28 True True 6.58 1.557 1.337 0.697 True offset 5 0-10 (a) min across groups 0.548 0.253 28 True True 6.58 1.557 1.337 0.709 True offset 5 0-10 (b) mean of group minima 0.655 0.236 28 True True 6.58 1.557 1.337 0.709 True offset N 0-20 (a) min across groups 0.509 0.233 28 True True 6.58 1.557 1.337 0.697 True offset N 0-20 (b) mean of group minima 0.636 0.222 28 True True 6.58 1.557 1.337 0.697 True offset N 0-10 (a) min across groups 0.547 0.252 28 True True 6.58 1.557 1.337 0.709 True offset N 0-10 (b) mean of group minima 0.654 0.231 28 True True 6.58 1.557 1.337 0.709 True onset 5 0-20 (a) min across groups 0.522 0.224 29 True True 6.58 1.557 1.363 0.722 True onset 5 0-20 (b) mean of group minima 0.652 0.209 29 True True 6.58 1.557 1.363 0.722 True onset 5 0-10 (a) min across groups 0.559 0.246 29 True True 6.58 1.557 1.363 0.737 True onset 5 0-10 (b) mean of group minima 0.679 0.217 29 True True 6.58 1.557 1.363 0.737 True onset N 0-20 (a) min across groups 0.548 0.202 29 True True 6.58 1.557 1.363 0.722 True onset N 0-20 (b) mean of group minima 0.670 0.203 29 True True 6.58 1.557 1.363 0.722 True onset N 0-10 (a) min across groups 0.581 0.227 29 True True 6.58 1.557 1.363 0.737 True onset N 0-10 (b) mean of group minima 0.691 0.211 29 True True 6.58 1.557 1.363 0.737 True combinations meeting the criterion: 16 of 16 order with the whole-trial random windows (2b) and all touches (3a, 4a): [np.float64(6.58), np.float64(1.59), np.float64(1.39), np.float64(0.85)]
Result: reproduced. Within the touch groups the minimum 10 ms FF is 0.504 to 0.691 under all 16 readings (paper 0.54 +- 0.27), and the five conditions fall in the paper's order (6.58 > 1.56 > 1.34 to 1.36 > 0.70 to 0.74 > 0.50 to 0.69) under each.
whisk-cycle-spikes¶
L4 neurons fire few spikes during whisker movements, on average fewer than 0.15 spikes per whisk cycle, which hinders decoding of whisking phase.
Spikes per whisk cycle over the cycles used for phase analysis (amplitude > 2.5 deg, monotonic phase, no lick, not within 70 ms after touch), over the whole trial. Readings: L4 C2 or all 41 L4 cells; mean over cells or over all cycles pooled. As a variant, the same cycles inside the sample epoch only (the earlier run's).
rows = []
for no in G["L4C2"] + G["L4out"]:
c = by_no[no]
cyc = D[no]["cyc"]
si = np.sort(hires.spike_samples(c, D[no]["s"]))
for rd, cc in (("phase-analysis cycles", cyc[cyc.phase_ok]), ("in the sample epoch (variant)", cyc[cyc.phase_ok & cyc.in_sample])):
n = np.searchsorted(si, cc.stop.to_numpy()) - np.searchsorted(si, cc.start.to_numpy())
rows.append(dict(cell=no, group=base.loc[no, "group"], cycles=rd, n_cycles=len(cc), spikes=int(n.sum())))
WCS = pd.DataFrame(rows)
WCS["per_cycle"] = WCS.spikes / WCS.n_cycles
for rd, g in WCS.groupby("cycles"):
for pop, gg in (("(a) L4 C2", g[g.group == "L4C2"]), ("(b) all L4", g)):
print(f"{rd}, {pop}: mean over cells {gg.per_cycle.mean():.3f}, pooled {gg.spikes.sum() / gg.n_cycles.sum():.3f} spikes/cycle "
f"({gg.n_cycles.sum()} cycles); below 0.15: {gg.per_cycle.mean() < 0.15}, {gg.spikes.sum() / gg.n_cycles.sum() < 0.15}")
print("cells above 0.15 spikes/cycle (phase-analysis cycles):", WCS[(WCS.cycles == "phase-analysis cycles") & (WCS.per_cycle > 0.15)].cell.tolist())
in the sample epoch (variant), (a) L4 C2: mean over cells 0.100, pooled 0.092 spikes/cycle (5377 cycles); below 0.15: True, True in the sample epoch (variant), (b) all L4: mean over cells 0.114, pooled 0.114 spikes/cycle (8627 cycles); below 0.15: True, True phase-analysis cycles, (a) L4 C2: mean over cells 0.107, pooled 0.116 spikes/cycle (17669 cycles); below 0.15: True, True phase-analysis cycles, (b) all L4: mean over cells 0.121, pooled 0.124 spikes/cycle (27769 cycles); below 0.15: True, True cells above 0.15 spikes/cycle (phase-analysis cycles): [2, 3, 9, 11, 13, 16, 18, 31, 34, 37, 39]
Result: reproduced. L4 cells fire 0.107 (L4 C2, mean over cells) to 0.121 (all L4) spikes per phase-analysis whisk cycle, 0.116 and 0.124 pooled, all below 0.15; eleven cells are above 0.15. Inside the sample epoch only, a variant, the values are 0.092 to 0.114.
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 |
|---|---|
| recordings-counts | reproduced |
| recordings-mice | not reproduced |
| behavior-performance | reproduced |
| behavior-sample-epoch | reproduced |
| behavior-whisking-bouts | not reproduced |
| behavior-touches-reaction | not reproduced |
| example-cell-rates | reproduced |
| example-cell-sample-ff | reproduced |
| example-cell-touch-window | not reproduced |
| example-cell-ff-around-touch | blocked: methods underspecified |
| text-spikes-per-touch | reproduced |
| text-whisking-rates | not reproduced |
| example-cell-phase-modulation | reproduced |
| table1-l4-c2 | reproduced |
| table1-l4-outside | blocked: methods underspecified |
| table1-l5 | partially reproduced |
| latency-determinable | not reproduced |
| l4-outside-weaker-touch | reproduced |
| l5-higher-rates | reproduced |
| phase-mod-l4-inside-outside | reproduced |
| phase-mod-l4-vs-l5 | reproduced |
| whisking-rate-l4-no-change | blocked: methods underspecified |
| whisking-rate-l5-increase | reproduced |
| fig3c-touch-spikes-l4c2, fig3c-touch-spikes-outside, fig3c-touch-spikes-l5 | blocked: methods underspecified |
| fig3d-touch-phase-l4c2, fig3d-touch-phase-outside, fig3d-touch-phase-l5, fig3d-phase-only | blocked: methods underspecified |
| fig4a-touch-detection-c2, fig4a-touch-detection-outside | blocked: methods underspecified |
| fig4b-touch-time-precision | not reproduced |
| fig4c-phase-decoding | not reproduced |
| fig5c-modulation-index | reproduced |
| fig5-modulation-correlations | blocked: methods underspecified |
| fig6s1b-clusters-less-variable | reproduced |
| fig6cd-ff-near-binomial | reproduced |
| fig6d-not-refractory | blocked: methods underspecified |
| fig6e-sample-epoch-ff, fig6e-random-38ms-ff, fig6e-touch-38ms-ff, fig6e-touch-10ms-ff, fig6e-similar-touch-ff | reproduced |
| whisk-cycle-spikes | reproduced |
Of the 45 claims, 22 are reproduced, 1 partially reproduced, 8 not reproduced and 14 blocked as methods underspecified. After the first verification round five outcomes changed. behavior-whisking-bouts, behavior-touches-reaction and example-cell-touch-window go from partially reproduced to not reproduced: their criteria require every value within tolerance and give no partial rule, and the bout duration, the touch count and the time fraction fail under every reading. whisking-rate-l4-no-change goes from partially reproduced to blocked: with whisking leaving out licking and the 70 ms after touch, as the claims file defines it, L4 C2 shows no change when touch epochs are left out (p = 0.595) and a change when they are in (p = 1.28e-7). fig4c-phase-decoding goes from blocked to not reproduced: with the phase-tuning cycles of the whole trial every reading puts C2 at 66 deg or more (or never at 76% correct) and worse than outside C2. The second verification round agreed with all 45 outcomes and changed none.
How the analysis differs from the paper's.
- The paper's analysis code is not public (the claims file records that the repository it links holds only the contact detection). Every procedure is rebuilt from the claims file.
- Touch onsets and offsets, the Hilbert amplitude and phase, and the lick and pole series are
the dandiset's own. The stored amplitude and phase match a 6 to 60 Hz fourth-order
Butterworth band-pass and Hilbert transform of
theta_at_base. - Spikes are placed on the
SpikeTraintimestamps, which put the L4 C2 touch-response peak at a median of 8 ms in the files without the offset and 12 ms in those with it; the voltage timestamps would put it at 8 and 22 ms. The residual 4 ms is not explained. - The details asked about in
questions/000013.yaml(where the cutoff search starts, the longest window of the cutoff curve, the pre-touch window of the touch detector, the decoder's time classes and smoothing, the phase decoder's cycles and empty bins, the rate of the refractory model, the measure of the L5 whisking increase) are now set or listed as readings by the claims file and run as it says. The values chosen here where the answers list a reading without a number: the floor rate of the phase decoder (one spike over all of a cell's cycles and bins) and the smaller pseudocount of the time decoder (0.01). Whisking periods for Table 1 are read as the claims file defines them, outside touch epochs, leaving out licking cycles and the 70 ms after each touch. Standard choices that are not expected to change a result: the polynomial's conditioning, the guard around touches for non-touch windows, the random seeds. Choices that are shown both ways: how ties between thresholds of the touch detector are broken (the middle one, the lowest as a variant), and the cycle holding a lick onset standing for "during licking", as the files hold lick onsets only. - Fewer repetitions than a reader might expect where the paper does not give a number: 50 per population size for touch detection.
What was found about the dataset.
- It holds all 52 recordings with the paper's group labels (31, 10, 11) and 115 +- 60 trials per neuron; the performance, the sample epoch and three of the five Figure 6E Fano factors come out as the paper's to two or three digits (the 10 ms touch and similar-touch Fano factors, 0.70 to 0.86 and 0.50 to 0.69 against 0.73 and 0.54, are within tolerance but vary more across readings), and the example cell of Figure 1 can be identified (cell 23).
- The files fall in batches converted differently. In 36 files the
SpikeTraintimestamps run 10 ms before the voltage timestamps of the same samples. In 38 files (those 36 and cells 20 and 40)theta_filtis not the band-passed angle the stored amplitude and phase come from, and the voltage, with its stated conversion of 0.001 to volts, spans only microvolts. In the other 14 files none of these holds. - The subjects (23) and genotypes do not match the paper's 21 mice and four lines: three subjects (six cells) have no genotype, and the C57BL6 and Ai32 x PV-Cre counts differ by one cell each.
pole_positionis described in mm, but its values range over several scales from file to file: 0.2 to 1.6 in some, 992 to 13888 (motor steps) in others, about 1e8 in others, 0 and 5952 in one, and up to 1.4e12 in another.set_pointis described as the Hilbert amplitude but follows the angle minus its band-passed part (r 0.957 to 0.999 where it has no missing values; 23 files hold missing values).- File names say
ogen, but no file holds an optogenetic stimulus andstim_presentis 0 in every trial. first_lick_timedisagrees with the lick series in some trials (printed in the data check).
print(f"notebook run time: {(time.time() - T_START) / 60:.1f} min")
notebook run time: 10.9 min