Reproducing Severson et al. (2017) from dandiset 000226¶
Paper. Severson KS, Xu D, Van de Loo M, Bai L, Ginty DD, O'Connor DH (2017). Active touch and self-motion encoding by Merkel cell-associated afferents. Neuron 94:666-676. doi:10.1016/j.neuron.2017.03.045. This attempt used the full text of the accepted author manuscript in PubMed Central (PMC5528144) and the publisher's Supplemental Information (Table S1 and the supplemental figure legends). The version of record is not open access. The paper's analysis code is not public; the Methods say it is available on request, and no repository for it was found.
Dandiset. DANDI:000226, version 0.230607.1747. The paper cites no DANDI version (it predates the deposit). This is the only published version; the draft has the same 60 assets and size, so the most recent published version is pinned.
What the dandiset holds. One NWB file per recording: 56 trigeminal afferents recorded with
high-speed whisker video, and 4 EMG sessions, which match the 60 recordings of Table S1. Each file
holds the spike times of one unit (none in three EMG files) and, for every 5 s trial, the whisker
kinematics, contact, distance to the pole, curvature change, and for most recordings with a pole,
forces and bending moment, at 500 Hz. The files do not say which experiment a recording belongs to,
the afferent type, or which afferents were identified as Merkel afferents; these come from Table S1
(see table-s1-link).
Streaming. Each file is read whole into memory from the archive and never written to disk
(see merkel.py). The run takes about 21 minutes, most of it reading 13.7 GB and fitting the spline models.
| Claim | What | Status |
|---|---|---|
| table-s1-link | Recordings match Table S1 | tested |
| video-frames | Frame counts (54.6 M, 65.3 M) | tested |
| merkel-minutes | Minutes of Merkel recordings | tested (outside labels) |
| fig2c-baseline-rate | Non-whisking rate 0.0 Hz | tested |
| fig2c-touch-rate | Contact rate 48.0 Hz, all touch-responsive | tested |
| fig2c-wt-star | 13 WT* afferents, whisking rate 19.1 Hz | tested |
| fig3a-direction-selective | Merkel and SA contact responses direction selective | tested (outside labels) |
| fig3d-gam-full | Full model r = 0.71 | tested |
| fig3d-gam-single | M0 alone 0.45, M0' alone 0.42 | tested |
| fig3d-gam-best-pair | M0 and M0' the best pair (0.63) | tested |
| fig4c-viscoelastic | Viscoelastic model comparable to the M0, M0' model | tested (outside labels) |
| fig5-phase-modulation | WT* phase modulation 2.1 against 51.2 Hz | tested |
| fig5e-phase-span | Preferred phases span the cycle | tested |
| figs3-phase-invariance | Preferred phase invariant to amplitude, frequency, setpoint | tested |
| figs6-kinematic-gam | Kinematic models predict phase tuning (r = 0.97) | tested on 15 of 28 |
| fig7g-phase-span-28 | 28 WT* phases, more in retraction | blocked |
| fig6-cut-eliminated | Cutting eliminated responses in 6 of 13 | blocked |
| fig6-cut-remained | Responses remained in 7 of 13 | blocked |
| fig6a-final-length | Final whisker length 0.37 mm | blocked |
| fig7d-cut-emg | Cut-whisker tuning matches EMG (r = 0.66) | blocked |
| fig7e-lost-rate-kinematics | Lost rate matches acceleration or jerk (r = 0.71) | blocked |
| fig7f-precut-emg | Pre-cut afferents match EMG (r = 0.76) | blocked |
| figs4-nasolabialis | Nasolabialis motor units follow negative velocity | tested |
| figs4-intrinsic-protractor | Intrinsic protractor EMG follows positive velocity | blocked |
| sigma-kappa | Curvature noise 7.8% | tested |
| figs1-optotag-latency | Light-evoked spike latency 4.9 ms | blocked |
| figs1-labelling-specificity | 99.2% of labelled endings are Merkel | blocked |
| contact-classifier | Contact classifier 99.5% correct | blocked |
In [1]:
import time
import warnings
from concurrent.futures import ThreadPoolExecutor
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from scipy import stats
import merkel
pd.set_option("display.width", 220)
pd.set_option("display.max_columns", 40)
pd.set_option("display.max_rows", 200)
warnings.filterwarnings("ignore", category=RuntimeWarning)
started = time.time()
table = pd.DataFrame(merkel.table_s1()).set_index("rec")
assets = merkel.assets()
print(f"dandiset {merkel.DANDISET_ID} version {merkel.DANDISET_VERSION}: {len(assets)} assets, "
f"{sum(a[1] for a in assets) / 1e9:.2f} GB")
dandiset 000226 version 0.230607.1747: 60 assets, 13.75 GB
The recordings to keep frames for are chosen from Table S1 before loading: the touch recordings
(for the touch models) and the 15 WT* recordings of Figure S3 (for the whisking-in-air models).
Their identifiers are matched to the files in table-s1-link; here only the file name is needed,
through the identifier read while loading.
In [2]:
WT_STAR_15 = sorted(list(table.index[(table.group == "touch") & (table.cls == "WT*")]) +
list(table.index[table.group == "whiskonly"]))
TOUCH_33 = list(table.index[table.group == "touch"])
import re
def rec_of(identifier):
# KS0240A -> 240A; KS0007A -> 7A; KS0125_ -> 125_ (matched by number in table-s1-link)
m = re.fullmatch(r"KS0*(\d+)([A-Z_])", identifier)
return f"{int(m.group(1))}{m.group(2)}"
def number_of(identifier):
return int(re.fullmatch(r"KS0*(\d+)[A-Z_]", identifier).group(1))
# The touch and WT* recordings by mouse and date, to decide what to keep before the identifier is read.
def wanted(path):
sub = path.split("/")[0].replace("sub-", "")
date = path.split("ses-")[1][:8]
rows = table[(table.mouse == sub) & (table.date.map(lambda d: pd.to_datetime(d, format="%m/%d/%y").strftime("%Y%m%d")) == date)]
return set(rows.index)
def reduce_one(item):
i, (path, size, url) = item
cands = wanted(path)
return merkel.reduce(path, url, seed=1000 + i,
keep_touch=bool(cands & set(TOUCH_33)), keep_wia=bool(cands & set(WT_STAR_15)))
with ThreadPoolExecutor(3) as ex:
recs = list(ex.map(reduce_one, enumerate(assets)))
print(f"loaded and reduced {len(recs)} files in {(time.time() - started) / 60:.1f} min")
loaded and reduced 60 files in 11.6 min
Data check¶
In [3]:
info = pd.DataFrame([{
"path": r["path"].split("/")[1][:48], "identifier": r["identifier"], "subject": r["subject_id"], "sex": r["sex"],
"strain": r["strain"], "date": r["session_start"][:10], "trials": r["n_trials"],
"spikes": r.get("n_spikes"), "rate_hz": r["rates"], "frames_per_series": r["lengths"],
"trial_dur_s": r["trial_duration"], "MOI": r["MOI_values"][:2], "MOI_nan": r["MOI_nan"],
"cutting": r["cutting_groups"], "has_M0": "M0" in r["series_present"],
"has_dpole": "distanceToPoleCenter" in r["series_present"]} for r in recs])
print("Files:", len(info), "; subjects:", info.subject.nunique(), "; files with spike times:", int(info.spikes.notna().sum()))
print("Archive summary: 60 files, 43 subjects (assetsSummary).")
print("Series per file at", sorted({x for r in recs for x in r["rates"]}), "Hz with", sorted({x for r in recs for x in r["lengths"]}), "frames per trial")
print("Trial table columns:", sorted({tuple(r["trial_columns"]) for r in recs}))
print("Acquisition groups present in any file:", sorted({x for r in recs for x in r["acquisition"]}))
print("Cutting groups in any file:", sorted({x for r in recs for x in r["cutting_groups"]}))
print("Trial durations in the trial table (s):", sorted({tuple(np.round(r["trial_duration"], 3)) for r in recs}),
"; series start equals trial start (max offset s):", max(r["start_offset"] for r in recs))
print("related_publications:", sorted({p for r in recs for p in r["related_publications"]}))
print("Clock times of session_start_time:", sorted({r["session_start"][11:19] for r in recs}))
print()
print(info.sort_values("path").to_string(index=False))
Files: 60 ; subjects: 43 ; files with spike times: 57
Archive summary: 60 files, 43 subjects (assetsSummary).
Series per file at [500.0] Hz with [2500] frames per trial
Trial table columns: [('start_time', 'stop_time', 'MOI', 'cuttingGroup', 'barCoor', 'barRadii', 'barAzim')]
Acquisition groups present in any file: []
Cutting groups in any file: ['intact']
Trial durations in the trial table (s): [(np.float64(6.998), np.float64(6.998))] ; series start equals trial start (max offset s): 0.0
related_publications: ['DOI: 10.1016/j.neuron.2017.03.045']
Clock times of session_start_time: ['12:00:00']
path identifier subject sex strain date trials spikes rate_hz frames_per_series trial_dur_s MOI MOI_nan cutting has_M0 has_dpole
sub-KSm11_ses-20160616T120000_behavior.nwb KS0011A KSm11 F Rosa Ai9 2016-06-16 378 13851.0 [500.0] [2500] (6.997999999999934, 6.998000000000047) [9.627381365394453e-13] 0 [intact] False False
sub-KSm30_ses-20160831T120000_behavior.nwb KS0278B KSm30 M Rosa Ai32 2016-08-31 150 5897.0 [500.0] [2500] (6.997999999999934, 6.998000000000047) [9.627381365394453e-13] 0 [intact] False False
sub-KSm31_ses-20160824T120000_behavior.nwb KS0277A KSm31 M Rosa Ai32 2016-08-24 128 7676.0 [500.0] [2500] (6.9979999999999905, 6.998000000000047) [9.627381365394453e-13] 0 [intact] False False
sub-KSm6_ses-20160602T120000_behavior.nwb KS0007A KSm6 M Rosa Ai32 2016-06-02 305 NaN [500.0] [2500] (6.997999999999934, 6.998000000000047) [9.627381365394453e-13] 0 [intact] False False
sub-KSm8_ses-20160604T120000_behavior.nwb KS0008A KSm8 M Rosa Ai32 2016-06-04 213 NaN [500.0] [2500] (6.997999999999934, 6.998000000000047) [9.627381365394453e-13] 0 [intact] False False
sub-KSm9_ses-20160610T120000_behavior.nwb KS0009A KSm9 M CCK CreER ;Rosa Ai32 2016-06-10 202 NaN [500.0] [2500] (6.997999999999934, 6.998000000000047) [9.627381365394453e-13] 0 [intact] False False
sub-KSt101_ses-20150511T120000_behavior.nwb KS0223B KSt101 F TrKC CreER ;Rosa Ai32 2015-05-11 1059 104155.0 [500.0] [2500] (6.9979999999995925, 6.998000000000047) [6.47569466309281e-12] 0 [intact] True True
sub-KSt106_ses-20150514T120000_obj-1pi120h_behav KS0221B KSt106 F TrKC CreER ;Rosa Ai32 2015-05-14 467 257945.0 [500.0] [2500] (6.997999999999934, 6.998000000000047) [9.376322888922841e-12] 0 [intact] True True
sub-KSt106_ses-20150514T120000_obj-g7vu6s_behavi KS0221C KSt106 F TrKC CreER ;Rosa Ai32 2015-05-14 726 22234.0 [500.0] [2500] (6.9979999999995925, 6.998000000000047) [9.376322888922841e-12] 0 [intact] True True
sub-KSt107_ses-20150514T120000_behavior.nwb KS0225D KSt107 F TrKC CreER ;Rosa Ai32 2015-05-14 843 241793.0 [500.0] [2500] (6.9979999999995925, 6.998000000000047) [1.4120907205937285e-12] 0 [intact] True True
sub-KSt107_ses-20150515T120000_obj-1lzhy6l_behav KS0227C KSt107 F TrKC CreER ;Rosa Ai32 2015-05-15 345 636.0 [500.0] [2500] (6.997999999999934, 6.998000000000047) [6.368891979190124e-12] 0 [intact] True True
sub-KSt107_ses-20150515T120000_obj-1yadv1x_behav KS0227A KSt107 F TrKC CreER ;Rosa Ai32 2015-05-15 266 28389.0 [500.0] [2500] (6.997999999999934, 6.998000000000047) [2.5570905640246913e-12] 0 [intact] True True
sub-KSt109_ses-20150521T120000_behavior.nwb KS0231A KSt109 F TrKC CreER ;Rosa Ai32 2015-05-21 903 38388.0 [500.0] [2500] (6.9979999999995925, 6.998000000000047) [2.236398921329255e-12] 0 [intact] True True
sub-KSt111_ses-20150529T120000_behavior.nwb KS0234A KSt111 M TrKC CreER ;Rosa Ai32 2015-05-29 942 47990.0 [500.0] [2500] (6.9979999999995925, 6.998000000000047) [8.823708484362019e-14] 0 [intact] True True
sub-KSt112_ses-20150606T120000_behavior.nwb KS0240A KSt112 M TrKC CreER ;Rosa Ai32 2015-06-06 58 19824.0 [500.0] [2500] (6.9979999999999905, 6.998000000000019) [2.6808257310632973e-12] 0 [intact] True True
sub-KSt112_ses-20150610T120000_behavior.nwb KS0243B KSt112 M TrKC CreER ;Rosa Ai32 2015-06-10 465 11392.0 [500.0] [2500] (6.997999999999934, 6.998000000000047) [2.46877044004752e-12] 0 [intact] True True
sub-KSt115_ses-20150622T120000_behavior.nwb KS0250A KSt115 F TrKC CreER ;Rosa Ai32 2015-06-22 701 237376.0 [500.0] [2500] (6.9979999999995925, 6.998000000000047) [2.8568692113804697e-12] 0 [intact] True True
sub-KSt117_ses-20150629T120000_behavior.nwb KS0251B KSt117 M CCK Cre ; Rosa Ai32 2015-06-29 1169 131041.0 [500.0] [2500] (6.9979999999995925, 6.998000000000047) [1.516523587976499e-12] 0 [intact] True True
sub-KSt119_ses-20150702T120000_behavior.nwb KS0253A KSt119 M TrKC CreER ;Rosa Ai32 2015-07-02 63 16532.0 [500.0] [2500] (6.9979999999999905, 6.998000000000019) [] 63 [intact] False True
sub-KSt119_ses-20150703T120000_behavior.nwb KS0254C KSt119 M TrKC CreER ;Rosa Ai32 2015-07-03 210 19187.0 [500.0] [2500] (6.997999999999934, 6.998000000000047) [4.000592989152262e-12] 0 [intact] True True
sub-KSt119_ses-20150706T120000_behavior.nwb KS0256B KSt119 M TrKC CreER ;Rosa Ai32 2015-07-06 430 118563.0 [500.0] [2500] (6.997999999999934, 6.998000000000047) [5.168252697511919e-12] 0 [intact] True True
sub-KSt122_ses-20150719T120000_behavior.nwb KS0265A KSt122 M TrKC CreER ;Rosa Ai32 2015-07-19 1058 35605.0 [500.0] [2500] (6.9979999999995925, 6.998000000000047) [2.1203026311995754e-13] 0 [intact] True True
sub-KSt123_ses-20150729T120000_behavior.nwb KS0271A KSt123 F TrKC CreER ;Rosa Ai32 2015-07-29 392 8230.0 [500.0] [2500] (6.997999999999934, 6.998000000000047) [3.964449366158483e-12] 0 [intact] True True
sub-KSt124_ses-20150721T120000_behavior.nwb KS0266B KSt124 F TrKC CreER ;Rosa Ai32 2015-07-21 879 10899.0 [500.0] [2500] (6.9979999999995925, 6.998000000000047) [5.883765665694027e-12] 0 [intact] True True
sub-KSt126_ses-20160831T120000_behavior.nwb KS0282A KSt126 M TrKC CreER ;Rosa Ai32 2016-08-31 118 5025.0 [500.0] [2500] (6.9979999999999905, 6.998000000000047) [5.5331587397828323e-14] 0 [intact] True True
sub-KSt127_ses-20160909T120000_behavior.nwb KS0286A KSt127 F TrKC CreER ;Rosa Ai32 2016-09-09 937 179963.0 [500.0] [2500] (6.9979999999995925, 6.998000000000047) [9.696257458150666e-13] 0 [intact] True True
sub-KSt128_ses-20160913T120000_behavior.nwb KS0288C KSt128 M TrKC CreER ;Rosa Ai32 2016-09-13 189 1087.0 [500.0] [2500] (6.997999999999934, 6.998000000000047) [9.627381365394453e-13] 0 [intact] False False
sub-KSt130_ses-20160923T120000_behavior.nwb KS0295A KSt130 M TrKC CreER ;Rosa Ai32 2016-09-23 213 3380.0 [500.0] [2500] (6.997999999999934, 6.998000000000047) [9.627381365394453e-13] 0 [intact] False False
sub-KSt131_ses-20160923T120000_behavior.nwb KS0295C KSt131 M TrKC CreER ;Rosa Ai32 2016-09-23 327 6537.0 [500.0] [2500] (6.997999999999934, 6.998000000000047) [9.627381365394453e-13] 0 [intact] False False
sub-KSt135_ses-20161024T120000_behavior.nwb KS0309B KSt135 M TrKC CreER ;Rosa Ai32 2016-10-24 277 6350.0 [500.0] [2500] (6.997999999999934, 6.998000000000047) [9.627381365394453e-13] 0 [intact] False False
sub-KSt135_ses-20161027T120000_behavior.nwb KS0311C KSt135 M TrKC CreER ;Rosa Ai32 2016-10-27 402 166850.0 [500.0] [2500] (6.997999999999934, 6.998000000000047) [9.627381365394453e-13] 0 [intact] False False
sub-KSt137_ses-20161015T120000_behavior.nwb KS0303A KSt137 M Rosa Ai32 2016-10-15 391 725.0 [500.0] [2500] (6.997999999999934, 6.998000000000047) [9.627381365394453e-13] 0 [intact] False False
sub-KSt137_ses-20161017T120000_behavior.nwb KS0304B KSt137 M Rosa Ai32 2016-10-17 183 820.0 [500.0] [2500] (6.997999999999934, 6.998000000000047) [9.627381365394453e-13] 0 [intact] False False
sub-KSt137_ses-20161018T120000_obj-1cimbwp_behav KS0305A KSt137 M Rosa Ai32 2016-10-18 82 3590.0 [500.0] [2500] (6.9979999999999905, 6.998000000000047) [9.627381365394453e-13] 0 [intact] False False
sub-KSt137_ses-20161018T120000_obj-1jbvpxr_behav KS0305B KSt137 M Rosa Ai32 2016-10-18 288 65413.0 [500.0] [2500] (6.997999999999934, 6.998000000000047) [9.627381365394453e-13] 0 [intact] False False
sub-KSt138_ses-20161120T120000_behavior.nwb KS0313A KSt138 M TrKC CreER ;Rosa Ai32 2016-11-20 332 2534.0 [500.0] [2500] (6.997999999999934, 6.998000000000047) [9.627381365394453e-13] 0 [intact] False False
sub-KSt140_ses-20161201T120000_behavior.nwb KS0318A KSt140 M TrKC CreER ;Rosa Ai32 2016-12-01 366 5459.0 [500.0] [2500] (6.997999999999934, 6.998000000000047) [9.627381365394453e-13] 0 [intact] False False
sub-KSt140_ses-20161205T120000_obj-13pr4ey_behav KS0321B KSt140 M TrKC CreER ;Rosa Ai32 2016-12-05 250 2403.0 [500.0] [2500] (6.997999999999934, 6.998000000000047) [2.015732973379472e-12] 0 [intact] False False
sub-KSt140_ses-20161205T120000_obj-ofsdfk_behavi KS0321A KSt140 M TrKC CreER ;Rosa Ai32 2016-12-05 450 58164.0 [500.0] [2500] (6.997999999999934, 6.998000000000047) [2.015732973379472e-12] 0 [intact] False False
sub-KSt66_ses-20140711T120000_behavior.nwb KS0125_ KSt66 M TrKC CreER ;Rosa Ai32 2014-07-11 591 2016.0 [500.0] [2500] (6.9979999999995925, 6.998000000000047) [5.485906282110278e-12] 0 [intact] True True
sub-KSt66_ses-20140718T120000_behavior.nwb KS0128_ KSt66 M TrKC CreER ;Rosa Ai32 2014-07-18 984 58165.0 [500.0] [2500] (6.9979999999995925, 6.998000000000047) [8.131375988068089e-12] 0 [intact] True True
sub-KSt69_ses-20140725T120000_behavior.nwb KS0132A KSt69 M TrKC CreER ;Rosa Ai32 2014-07-25 167 40281.0 [500.0] [2500] (6.997999999999934, 6.998000000000047) [9.573041163130305e-12] 0 [intact] True True
sub-KSt76_ses-20140927T120000_behavior.nwb KS0145B KSt76 M TrKC CreER ;Rosa Ai32 2014-09-27 216 4100.0 [500.0] [2500] (6.997999999999934, 6.998000000000047) [3.124442570484683e-13] 0 [intact] True True
sub-KSt77_ses-20141110T120000_obj-k3liaq_behavio KS0163B KSt77 M TrKC CreER ;Rosa Ai32 2014-11-10 366 6542.0 [500.0] [2500] (6.997999999999934, 6.998000000000047) [7.114635290530249e-12] 0 [intact] True True
sub-KSt77_ses-20141110T120000_obj-maz7dk_behavio KS0161A KSt77 M TrKC CreER ;Rosa Ai32 2014-11-10 499 48973.0 [500.0] [2500] (6.997999999999934, 6.998000000000047) [5.373143465195956e-12] 0 [intact] True True
sub-KSt80_ses-20141025T120000_behavior.nwb KS0154B KSt80 F TrKC CreER ;Rosa Ai32 2014-10-25 605 35801.0 [500.0] [2500] (6.9979999999995925, 6.998000000000047) [7.365107458580644e-12] 0 [intact] True True
sub-KSt82_ses-20141208T120000_behavior.nwb KS0167A KSt82 M TrKC CreER ;Rosa Ai32 2014-12-08 433 32194.0 [500.0] [2500] (6.997999999999934, 6.998000000000047) [5.753993221585802e-12] 0 [intact] True True
sub-KSt88_ses-20150320T120000_behavior.nwb KS0194A KSt88 M TrKC CreER ;Rosa Ai32 2015-03-20 425 6647.0 [500.0] [2500] (6.997999999999934, 6.998000000000047) [1.0100276055640553e-11] 0 [intact] True True
sub-KSt89_ses-20150331T120000_behavior.nwb KS0198B KSt89 F TrKC CreER ;Rosa Ai32 2015-03-31 651 5860.0 [500.0] [2500] (6.9979999999995925, 6.998000000000047) [5.368259023735076e-12] 0 [intact] True True
sub-KSt89_ses-20150401T120000_behavior.nwb KS0199B KSt89 F TrKC CreER ;Rosa Ai32 2015-04-01 173 7669.0 [500.0] [2500] (6.997999999999934, 6.998000000000047) [3.3510321638291216e-13] 0 [intact] True True
sub-KSt91_ses-20150413T120000_obj-2btolv_behavio KS0206A KSt91 M TrKC CreER ;Rosa Ai32 2015-04-13 62 13013.0 [500.0] [2500] (6.9979999999999905, 6.998000000000019) [9.627381365394453e-13] 0 [intact] False False
sub-KSt91_ses-20150413T120000_obj-u72o7l_behavio KS0207A KSt91 M TrKC CreER ;Rosa Ai32 2015-04-13 739 91811.0 [500.0] [2500] (6.9979999999995925, 6.998000000000047) [7.269901863698197e-14] 0 [intact] True True
sub-KSt92_ses-20150405T120000_behavior.nwb KS0202A KSt92 M TrKC CreER ;Rosa Ai32 2015-04-05 851 10241.0 [500.0] [2500] (6.9979999999995925, 6.998000000000047) [1.5835269847415939e-12] 0 [intact] True True
sub-KSt93_ses-20150416T120000_behavior.nwb KS0209A KSt93 M TrKC CreER ;Rosa Ai32 2015-04-16 36 1460.0 [500.0] [2500] (6.9979999999999905, 6.998000000000005) [3.3510321638291216e-13] 0 [intact] True True
sub-KSt94_ses-20150424T120000_behavior.nwb KS0214A KSt94 F TrKC CreER ;Rosa Ai32 2015-04-24 181 214.0 [500.0] [2500] (6.997999999999934, 6.998000000000047) [2.597704214658853e-12] 0 [intact] True True
sub-KSt95_ses-20150429T120000_behavior.nwb KS0217F KSt95 M TrKC CreER ;Rosa Ai32 2015-04-29 986 178935.0 [500.0] [2500] (6.9979999999995925, 6.998000000000047) [6.063252415836189e-12] 0 [intact] True True
sub-KSt96_ses-20150420T120000_behavior.nwb KS0211A KSt96 F TrKC CreER ;Rosa Ai32 2015-04-20 1384 47444.0 [500.0] [2500] (6.9979999999995925, 6.998000000000047) [4.78094818607619e-13] 0 [intact] True True
sub-KSt97_ses-20150427T120000_behavior.nwb KS0216D KSt97 M CCK Cre ; Rosa Ai32 2015-04-27 634 86220.0 [500.0] [2500] (6.9979999999995925, 6.998000000000047) [2.219660763133158e-13] 0 [intact] True True
sub-KSt99_ses-20150506T120000_obj-1wybaw_behavio KS0220C KSt99 M TrKC CreER ;Rosa Ai32 2015-05-06 950 5426.0 [500.0] [2500] (6.9979999999995925, 6.998000000000047) [2.1768206065022113e-13] 0 [intact] True True
sub-KSt99_ses-20150506T120000_obj-r9qsp2_behavio KS0220B KSt99 M TrKC CreER ;Rosa Ai32 2015-05-06 210 44145.0 [500.0] [2500] (6.997999999999934, 6.998000000000047) [4.756721480804156e-12] 0 [intact] True True
Each 5 s trial holds 2,500 frames at 500 Hz, as the Methods state, but stop_time in every trial
table is 7 s after start_time. Every trial's series start at the trial's start_time. All trials
of all files are labelled intact in cuttingGroup, including the 13 progressive cutting
recordings, and MOI (moment of inertia) has one value per file; one file (KS0253A) has MOI
missing in every trial. Session start times are all at 12:00, so only their dates carry
information. The archive's count of 43 subjects matches the files.
In [4]:
units = {}
for r in recs:
for v, u in r["units"].items():
units.setdefault(v, set()).add(u)
print("Stated units:", {v: sorted(u) for v, u in units.items()})
rng_tab = pd.DataFrame([{"variable": v,
"min": min(r["ranges"][v][0] for r in recs if v in r["ranges"]),
"max": max(r["ranges"][v][1] for r in recs if v in r["ranges"]),
"files": sum(v in r["ranges"] for r in recs),
"finite_fraction_min": min(r["ranges"][v][2] for r in recs if v in r["ranges"])}
for v in merkel.VARIABLES])
print(rng_tab.to_string(index=False))
missing = {r["identifier"]: [v for v in ("M0", "Faxial", "distanceToPoleCenter") if v not in r["series_present"]] for r in recs}
print("\nFiles without forces and moment:", sorted(k for k, v in missing.items() if "M0" in v))
print("Files without distance to pole:", sorted(k for k, v in missing.items() if "distanceToPoleCenter" in v))
Stated units: {'position': ['degree'], 'velocity': ['degree/s'], 'acceleration': ['degree/s^2'], 'jerk': ['degree/s^3'], 'phase': ['radians'], 'amplitude': ['degree'], 'setpoint': ['degree'], 'frequency': ['Hz'], 'contact': ['boolean'], 'deltaKappa': ['1/mm'], 'Mi': ['N*m'], 'dMi': ['N*m/s'], 'M0': ['N*m'], 'dM0': ['N*m/s'], 'Faxial': ['N'], 'Flateral': ['N'], 'dFaxial': ['N/s'], 'dFlateral': ['N/s'], 'distanceToPoleCenter': ['unknown']}
variable min max files finite_fraction_min
position -984.727295 596.020752 60 0.940523
velocity -965.753601 832.439575 60 0.949328
acceleration -1849.780762 453.450958 60 0.955824
jerk -1695.111694 1511.652832 60 0.960827
phase -3.141593 3.141593 60 0.923546
amplitude 0.000009 75.579056 60 0.923546
setpoint -69.105934 78.830322 60 0.923546
frequency -249.858047 249.981522 60 0.918912
M0 -0.522661 0.186033 39 0.999448
dM0 -0.163885 0.155027 39 0.999341
Faxial -4.362401 12.256233 39 0.999448
Flateral -121.551895 43.264404 39 0.999448
dFaxial -3.635336 3.843059 39 0.999341
dFlateral -38.113747 36.053642 39 0.999341
contact 0.000000 1.000000 60 1.000000
distanceToPoleCenter -0.257661 20.450836 40 0.793651
deltaKappa -0.410638 0.500954 60 0.923519
Files without forces and moment: ['KS0007A', 'KS0008A', 'KS0009A', 'KS0011A', 'KS0206A', 'KS0253A', 'KS0277A', 'KS0278B', 'KS0288C', 'KS0295A', 'KS0295C', 'KS0303A', 'KS0304B', 'KS0305A', 'KS0305B', 'KS0309B', 'KS0311C', 'KS0313A', 'KS0318A', 'KS0321A', 'KS0321B']
Files without distance to pole: ['KS0007A', 'KS0008A', 'KS0009A', 'KS0011A', 'KS0206A', 'KS0277A', 'KS0278B', 'KS0288C', 'KS0295A', 'KS0295C', 'KS0303A', 'KS0304B', 'KS0305A', 'KS0305B', 'KS0309B', 'KS0311C', 'KS0313A', 'KS0318A', 'KS0321A', 'KS0321B']
Phase spans -pi to pi radians, amplitude is in degrees, and contact is 0 or 1. The ranges include
extreme values (position from -985 to 596 degrees, frequency from -250 to 250 Hz, M0 down to
-0.52 N m) that the Glossary's outlier removal excludes. The unit of distanceToPoleCenter is
"unknown"; the next cells check it against contact, and the units of the derivatives.
In [5]:
spk = [r for r in recs if r["has_spikes"]]
place = pd.DataFrame([{"identifier": r["identifier"], **r["spike_placement"], "sorted": r["spikes_sorted"]} for r in spk])
print("Spikes in the 2,500 frames of a trial:", place.in_frames.sum(), "; inside a trial interval but after its frames:",
place.in_trial_after_frames.sum(), "; outside every trial:", place.outside_trials.sum())
print("Files whose spike times are not sorted:", list(place.identifier[~place.sorted]))
sig = pd.DataFrame([{"identifier": r["identifier"], "pos": r["signature"][0], "spk": r["signature"][1]} for r in recs])
dup_pos = sig[sig.duplicated("pos", keep=False)]
dup_spk = sig[sig.duplicated("spk", keep=False) & sig.identifier.isin([r["identifier"] for r in spk])]
print("Files with identical whisker position in their first 5 trials:", list(dup_pos.identifier) or "none")
print("Files with identical spike trains:", list(dup_spk.identifier) or "none")
same_session = info.groupby(["subject", "date"]).identifier.apply(list)
print("Sessions with more than one file (different units; checked above for duplicates):",
[v for v in same_session if len(v) > 1])
Spikes in the 2,500 frames of a trial: 2613060 ; inside a trial interval but after its frames: 0 ; outside every trial: 0 Files whose spike times are not sorted: ['KS0277A'] Files with identical whisker position in their first 5 trials: none Files with identical spike trains: none Sessions with more than one file (different units; checked above for duplicates): [['KS0221B', 'KS0221C'], ['KS0227C', 'KS0227A'], ['KS0305A', 'KS0305B'], ['KS0321B', 'KS0321A'], ['KS0163B', 'KS0161A'], ['KS0206A', 'KS0207A'], ['KS0220C', 'KS0220B']]
No two files share their whisker traces or spike trains, so no recording is deposited twice; the seven sessions with two files hold two different units each, as Table S1 lists them. Spike times in one file are not sorted (named above). No spike falls outside the 5 s of video of its trial.
In [6]:
# distanceToPoleCenter against contact, and the sign of M0 against the direction of motion at contact onset.
rows = []
for r in recs:
if "episodes" not in r:
continue
ep = r["episodes"]
ok = np.isfinite(ep[:, 0]) & (ep[:, 3] != 0)
rows.append({"identifier": r["identifier"], "contacts": len(ep),
"M0_sign_matches_velocity_at_onset": np.mean(np.sign(ep[ok, 0]) == np.sign(ep[ok, 3]))})
sgn = pd.DataFrame(rows)
print("Per touch recording, fraction of contacts whose mean M0 has the sign of the whisker velocity at contact onset:")
print(sgn.M0_sign_matches_velocity_at_onset.describe().round(3).to_string())
secs = pd.DataFrame([{"identifier": r["identifier"], **{k: round(v / 60, 1) for k, v in r["category_seconds"].items()}} for r in recs])
print("\nMinutes per category over all files:", secs.drop(columns="identifier").sum().round(0).to_dict())
dp = pd.DataFrame([{"identifier": r["identifier"], "contact_p1": r["dpole_contact_pct"][0], "contact_p50": r["dpole_contact_pct"][1],
"contact_p99": r["dpole_contact_pct"][2], "free_p1": (r["dpole_free_pct"] or [np.nan])[0],
"free_p50": (r["dpole_free_pct"] or [np.nan, np.nan])[1]} for r in recs if "dpole_contact_pct" in r])
print("\ndistanceToPoleCenter in frames with contact = 1 (percentiles 1, 50, 99) and without contact (1, 50), per file:")
print(dp.round(2).to_string(index=False))
kp = pd.DataFrame({"phase": recs[0]["kinematics_by_phase"]["centres"].round(2),
**{f"{q}_{k}": np.mean([r["kinematics_by_phase"][k] for r in recs if "kinematics_by_phase" in r], axis=0)
for k in ("position",) for q in ("mean",)},
"mean_velocity": np.mean([r["kinematics_by_phase"]["velocity"] for r in recs if "kinematics_by_phase" in r], axis=0)})
print("\nMean position (minus its mean, degrees) and velocity (as stored; per ms, see the ratio check below) by phase, whisking frames, averaged over files:")
print(kp.round(2).to_string(index=False))
neg_vel = [np.mean(r["kinematics_by_phase"]["velocity"][r["kinematics_by_phase"]["centres"] < 0]) for r in recs if "kinematics_by_phase" in r]
print("Files where the mean velocity over phases below 0 is negative:", int(np.sum(np.array(neg_vel) < 0)), "of", len(neg_vel))
pdir = pd.DataFrame([r["phase_direction"] for r in recs])
for k in ("sin", "cos"):
print(f"Correlation of {k}(stored phase) with {k}(Hilbert phase of band-passed position), per file: median",
round(float(pdir[k].median()), 4), ", range", round(float(pdir[k].min()), 4), "to", round(float(pdir[k].max()), 4))
print("(sin -1 and cos +1: the stored phase is the negative of the Hilbert phase; both -1: shifted by pi)")
print("Files where mean position is highest in the bin nearest phase 0:",
int(sum(np.argmin(np.abs(r["kinematics_by_phase"]["centres"])) in np.argsort(r["kinematics_by_phase"]["position"])[-2:]
for r in recs if "kinematics_by_phase" in r)), "of", sum("kinematics_by_phase" in r for r in recs),
"(bins of width pi/6; the two bins either side of 0 count)")
ratios = pd.DataFrame([{"identifier": r["identifier"], **r["derivative_ratios"]} for r in recs])
print("\nStored derivative over the change of its variable per 2 ms frame (0.5 means per ms, 500 per s):")
print(ratios.drop(columns="identifier").describe().loc[["count", "min", "50%", "max"]].round(3).to_string())
Per touch recording, fraction of contacts whose mean M0 has the sign of the whisker velocity at contact onset:
count 34.000
mean 0.864
std 0.041
min 0.769
25% 0.841
50% 0.863
75% 0.889
max 0.925
Minutes per category over all files: {'contact': 358.0, 'wia': 708.0, 'nonwhisking': 387.0, 'whisking': 1322.0}
distanceToPoleCenter in frames with contact = 1 (percentiles 1, 50, 99) and without contact (1, 50), per file:
identifier contact_p1 contact_p50 contact_p99 free_p1 free_p50
KS0223B -0.16 -0.08 0.17 -0.00 2.11
KS0221B -0.09 0.06 0.14 0.05 3.35
KS0221C -0.08 0.06 0.11 0.12 2.90
KS0225D -0.16 -0.01 0.10 0.09 3.08
KS0227C -0.12 0.08 0.17 0.12 3.82
KS0227A -0.15 0.07 0.14 0.13 3.30
KS0231A -0.13 -0.01 0.09 0.11 3.84
KS0234A -0.18 -0.06 0.17 0.06 1.58
KS0243B -0.17 -0.01 0.14 0.04 3.50
KS0250A -0.07 -0.00 0.10 0.11 3.85
KS0251B -0.18 -0.03 0.06 0.09 2.84
KS0254C -0.06 -0.00 0.14 0.10 3.01
KS0256B -0.10 0.00 0.08 0.07 2.93
KS0265A -0.17 -0.06 0.14 0.09 1.92
KS0271A -0.14 0.03 0.21 0.08 3.13
KS0266B -0.12 0.02 0.09 0.14 3.97
KS0282A -0.19 -0.01 0.10 0.10 1.77
KS0286A -0.11 0.00 0.08 0.10 2.68
KS0125_ -0.05 0.02 0.06 0.15 3.39
KS0128_ -0.17 0.09 0.56 0.20 3.08
KS0132A -0.14 0.08 0.29 0.34 5.52
KS0145B -0.13 0.02 0.34 0.15 2.68
KS0163B -0.04 0.05 0.33 0.26 4.29
KS0161A -0.11 0.06 0.37 0.23 4.77
KS0154B -0.03 0.05 0.33 0.11 2.17
KS0167A -0.09 0.05 0.27 0.20 4.84
KS0194A -0.11 0.05 0.13 0.15 3.99
KS0198B -0.08 0.04 0.15 0.10 2.71
KS0199B -0.08 -0.04 0.03 1.07 9.97
KS0207A -0.21 -0.02 0.21 0.21 2.72
KS0202A -0.10 -0.02 0.07 0.14 4.12
KS0209A 5.57 6.02 6.49 6.39 11.36
KS0214A -0.12 0.01 0.08 0.12 3.12
KS0217F -0.12 0.06 0.16 0.07 3.30
KS0211A -0.16 -0.01 0.06 0.07 2.34
KS0216D -0.17 -0.05 0.08 0.07 2.15
KS0220C -0.17 0.02 0.12 0.12 2.47
KS0220B -0.13 -0.05 0.04 0.07 4.42
Mean position (minus its mean, degrees) and velocity (as stored; per ms, see the ratio check below) by phase, whisking frames, averaged over files:
phase mean_position mean_velocity
-2.88 -11.80 -0.11
-2.36 -9.66 -0.63
-1.83 -4.39 -1.18
-1.31 3.33 -1.36
-0.79 10.31 -0.85
-0.26 13.06 -0.06
0.26 11.97 0.43
0.79 8.54 0.66
1.31 3.39 0.79
1.83 -2.49 0.80
2.36 -7.74 0.64
2.88 -11.10 0.31
Files where the mean velocity over phases below 0 is negative: 60 of 60
Correlation of sin(stored phase) with sin(Hilbert phase of band-passed position), per file: median -0.9999 , range -1.0 to -0.9962
Correlation of cos(stored phase) with cos(Hilbert phase of band-passed position), per file: median 0.9999 , range 0.9976 to 1.0
(sin -1 and cos +1: the stored phase is the negative of the Hilbert phase; both -1: shifted by pi)
Files where mean position is highest in the bin nearest phase 0: 60 of 60 (bins of width pi/6; the two bins either side of 0 count)
Stored derivative over the change of its variable per 2 ms frame (0.5 means per ms, 500 per s):
velocity acceleration jerk dM0 dFaxial dFlateral
count 60.000 60.000 60.000 37.000 37.000 37.000
min 0.481 0.451 0.431 0.491 0.484 0.490
50% 0.497 0.485 0.467 0.502 0.494 0.502
max 0.499 0.493 0.478 0.507 0.501 0.507
Frames with contact = 1 have distances near 0 and frames without contact have larger ones (table above), consistent with the Methods' d_pole in mm, measured to the pole's edge, though the files do not state the unit. Protraction (positive velocity) contacts mostly have M0 > 0, as Figure 3F has it, so the sign of M0 gives the contact direction.
Mean position is highest in the bins next to phase 0 in all 60 files, so phase 0 is full protraction, as in the Hilbert convention the paper cites (Hill et al. 2011). But the stored phase runs backwards in time: the sine of the stored phase is anticorrelated (median -0.9999) and its cosine correlated (median 0.9999) with those of the Hilbert phase computed here from the band-passed position, so the stored phase is the negative of the Hilbert phase, not a shift by pi. Velocity is negative at phases below 0 in every file. In the files, phases from -pi to 0 are retraction and 0 to pi protraction, the reverse of the usual convention. Neither the files nor the paper state which convention the stored phase follows.
The stored velocity, acceleration, jerk and the rates of change of moment and forces are per millisecond (ratio near 0.5 to the change per 2 ms frame), while the files state per second ("degree/s", "N*m/s"). Phase curves and the spline models do not depend on these scales; the viscoelastic model is given dM0 in N m per second.
In [7]:
# Clock: responses with the true spike times against spikes shifted by 0.5 s within each trial.
clk = pd.DataFrame([{"identifier": r["identifier"], "rec": rec_of(r["identifier"]),
"contact_minus_wia_true": r["clock"]["true"]["contact_minus_wia"],
"contact_minus_wia_shifted": r["clock"]["shifted"]["contact_minus_wia"],
"phase_depth_true": r["clock"]["true"]["phase_depth"],
"phase_depth_shifted": r["clock"]["shifted"]["phase_depth"]} for r in spk])
clk["touch_contrast_holds"] = clk.contact_minus_wia_true > 2 * np.abs(clk.contact_minus_wia_shifted)
clk["phase_contrast_holds"] = clk.phase_depth_true > 1.5 * clk.phase_depth_shifted
clk["established"] = clk.touch_contrast_holds.fillna(False) | clk.phase_contrast_holds
print(clk.round(2).to_string(index=False))
print("\nClock shown by one of the two contrasts:", int(clk.established.sum()), "of", len(clk))
print("Not established:", list(clk.rec[~clk.established]))
identifier rec contact_minus_wia_true contact_minus_wia_shifted phase_depth_true phase_depth_shifted touch_contrast_holds phase_contrast_holds established KS0011A 11A NaN NaN 36.90 2.12 False True True KS0278B 278B NaN NaN 17.29 2.64 False True True KS0277A 277A NaN NaN 69.11 4.53 False True True KS0223B 223B 88.180000 35.230000 1.19 14.09 True False True KS0221B 221B 156.580002 56.209999 145.45 60.53 True True True KS0221C 221C 17.350000 3.320000 2.17 6.64 True False True KS0225D 225D 195.860001 27.500000 48.28 42.50 True False True KS0227C 227C -0.180000 0.210000 7.93 0.97 False True True KS0227A 227A 83.510002 15.770000 1.85 18.08 True False True KS0231A 231A 24.430000 16.230000 33.88 3.62 False True True KS0234A 234A 33.720001 20.110001 0.18 2.02 False False False KS0240A 240A NaN NaN 314.59 24.39 False True True KS0243B 243B 29.049999 2.160000 0.72 3.43 True False True KS0250A 250A 101.459999 17.420000 117.50 17.30 True True True KS0251B 251B 112.519997 43.139999 5.83 8.97 True False True KS0253A 253A NaN NaN 97.40 14.59 False True True KS0254C 254C 95.370003 35.349998 2.03 16.37 True False True KS0256B 256B 59.810001 12.890000 115.62 22.01 True True True KS0265A 265A 37.900002 5.190000 3.09 2.74 True False True KS0271A 271A 16.240000 3.030000 0.76 4.16 True False True KS0266B 266B 11.550000 2.320000 0.66 0.89 True False True KS0282A 282A 41.549999 19.020000 0.40 11.48 True False True KS0286A 286A 161.080002 41.450001 15.13 26.50 True False True KS0288C 288C NaN NaN 5.55 1.15 False True True KS0295A 295A NaN NaN 29.38 2.12 False True True KS0295C 295C NaN NaN 66.54 1.36 False True True KS0309B 309B NaN NaN 24.29 1.60 False True True KS0311C 311C NaN NaN 142.27 4.64 False True True KS0303A 303A NaN NaN 2.48 0.96 False True True KS0304B 304B NaN NaN 7.56 0.77 False True True KS0305A 305A NaN NaN 58.55 3.80 False True True KS0305B 305B NaN NaN 106.02 6.18 False True True KS0313A 313A NaN NaN 20.16 1.09 False True True KS0318A 318A NaN NaN 10.67 1.19 False True True KS0321B 321B NaN NaN 35.86 1.81 False True True KS0321A 321A NaN NaN 49.29 3.52 False True True KS0125_ 125_ 3.590000 2.020000 0.00 1.18 False False False KS0128_ 128_ 23.459999 6.400000 48.18 6.33 True True True KS0132A 132A 91.970001 38.430000 328.45 22.66 True True True KS0145B 145B 18.570000 2.840000 4.05 6.48 True False True KS0163B 163B 11.720000 4.600000 0.00 2.31 True False True KS0161A 161A 58.299999 15.800000 10.62 12.68 True False True KS0154B 154B 47.970001 20.889999 4.31 5.71 True False True KS0167A 167A 46.450001 11.920000 7.66 8.71 True False True KS0194A 194A 11.720000 1.970000 0.06 2.22 True False True KS0198B 198B 4.370000 1.940000 0.96 2.00 True False True KS0199B 199B 61.139999 13.820000 53.81 4.85 True True True KS0206A 206A NaN NaN 187.46 14.80 False True True KS0207A 207A 96.489998 52.470001 1.22 8.16 False False False KS0202A 202A 10.700000 3.450000 1.39 2.34 True False True KS0209A 209A NaN NaN 41.92 8.79 False True True KS0214A 214A 2.030000 1.080000 0.42 0.55 False False False KS0217F 217F 112.080002 48.270000 36.37 11.53 True True True KS0211A 211A 21.139999 8.480000 0.04 3.60 True False True KS0216D 216D 101.379997 38.459999 0.42 11.88 True False True KS0220C 220C 3.730000 1.430000 0.53 0.81 True False True KS0220B 220B 202.979996 53.959999 54.54 17.13 True True True Clock shown by one of the two contrasts: 53 of 57 Not established: ['234A', '125_', '207A', '214A']
The spikes and the video series share a clock when the touch response (contact minus whisking in air) or the phase modulation depth is clearly larger with the true spike times than with spikes shifted by 0.5 s. Recordings with no touch response and little phase tuning cannot show it either way; they are named above.
In [8]:
# Progressive cutting recordings: whisking-in-air rate over the session.
prog = [r for r in spk if rec_of(r["identifier"]) in set(table.index[table.group == "progressive"])]
print("Progressive cutting recordings:", len(prog), "; cutting groups in their trial tables:",
sorted({g for r in prog for g in r["cutting_groups"]}), "; MOI values per file:", sorted({len(r["MOI_values"]) for r in prog}))
print("(The rate by fifth of the session is shown under fig6-cut-eliminated.)")
Progressive cutting recordings: 13 ; cutting groups in their trial tables: ['intact'] ; MOI values per file: [1] (The rate by fifth of the session is shown under fig6-cut-eliminated.)
The clock is shown for 53 of 57 recordings with spikes. For four touch recordings with little whisking-in-air activity (named above), neither contrast is clear, because contacts last long enough that spikes shifted by 0.5 s still fall in contact; their timing is taken as it is.
table-s1-link¶
The dandiset's recordings can be matched one to one to the recordings of Table S1, and the matched records agree on mouse, sex, mouse line and date.
In [9]:
rows = []
for r in recs:
rec = rec_of(r["identifier"])
if rec in table.index:
match, how = rec, "id"
else:
cands = [t for t in table.index if re.fullmatch(r"\d+", t[:-1]) and int(t[:-1]) == number_of(r["identifier"])
and table.loc[t, "mouse"] == r["subject_id"]]
match, how = (cands[0], "number and mouse") if len(cands) == 1 else (None, "none")
row = {"identifier": r["identifier"], "matched": match, "how": how}
if match:
t = table.loc[match]
norm = lambda x: x.replace(" ", "").lower()
row.update({"mouse_ok": t.mouse == r["subject_id"], "sex_ok": t.sex == r["sex"],
"line_ok": norm(t.line) == norm(r["strain"]),
"date_ok": pd.to_datetime(t.date, format="%m/%d/%y").strftime("%Y-%m-%d") == r["session_start"][:10],
"group": t.group})
rows.append(row)
link = pd.DataFrame(rows)
link["all_ok"] = link[["mouse_ok", "sex_ok", "line_ok", "date_ok"]].all(axis=1)
REC = dict(zip(link.identifier, link.matched))
print(link.to_string(index=False))
n = len(link)
print(f"\nmatched: {link.matched.notna().sum()} of {n}; distinct: {link.matched.nunique()}; by number and mouse only: {(link.how == 'number and mouse').sum()}")
print(f"agree on mouse, sex, line and date: {link.all_ok.sum()} of {n} ({100 * link.all_ok.mean():.1f}%)")
print("Table S1 recordings without a file:", sorted(set(table.index) - set(link.matched.dropna())))
print("Files per Table S1 group:", link.group.value_counts().to_dict())
identifier matched how mouse_ok sex_ok line_ok date_ok group all_ok
KS0011A 11A id True True True True emg True
KS0278B 278B id True True True True progressive True
KS0277A 277A id True True True True cutonly True
KS0007A 7A id True True True True emg True
KS0008A 8A id True True True True emg True
KS0009A 9A id True True True True emg True
KS0223B 223B id True True True True touch True
KS0221B 221B id True True True True touch True
KS0221C 221C id True True True True touch True
KS0225D 225D id True True True True touch True
KS0227C 227C id True True True True excluded True
KS0227A 227A id True True True True touch True
KS0231A 231A id True True True True touch True
KS0234A 234A id True True True True touch True
KS0240A 240A id True True True True whiskonly True
KS0243B 243B id True True True True touch True
KS0250A 250A id True True True True touch True
KS0251B 251B id True True True True touch True
KS0253A 253A id True True True True whiskonly True
KS0254C 254C id True True True True touch True
KS0256B 256B id True True True True touch True
KS0265A 265A id True True True True touch True
KS0271A 271A id True True True True excluded True
KS0266B 266B id True True True True touch True
KS0282A 282A id True True True True touch True
KS0286A 286A id True True True True touch True
KS0288C 288C id True True True True progressive True
KS0295A 295A id True True True True progressive True
KS0295C 295C id True True True True progressive True
KS0309B 309B id True True True True progressive True
KS0311C 311C id True True True True progressive True
KS0303A 303A id True True True True progressive True
KS0304B 304B id True True True True progressive True
KS0305A 305A id True True True True cutonly True
KS0305B 305B id True True True True progressive True
KS0313A 313A id True True True True progressive True
KS0318A 318A id True True True True progressive True
KS0321B 321B id True True True True progressive True
KS0321A 321A id True True True True progressive True
KS0125_ 125B number and mouse True True True True touch True
KS0128_ 128D number and mouse True True True True touch True
KS0132A 132A id True True True True touch True
KS0145B 145B id True True True True touch True
KS0163B 163B id True True True True touch True
KS0161A 161A id True True True True touch True
KS0154B 154B id True True True True touch True
KS0167A 167A id True True True True touch True
KS0194A 194A id True True True True touch True
KS0198B 198B id True True True True touch True
KS0199B 199B id True True True True cutonly True
KS0206A 206A id True True True True cutonly True
KS0207A 207A id True True True True touch True
KS0202A 202A id True True True True touch True
KS0209A 209A id True True True True cutonly True
KS0214A 214A id True True True True excluded True
KS0217F 217F id True True True True touch True
KS0211A 211A id True True True True touch True
KS0216D 216D id True True True True touch True
KS0220C 220C id True True True True touch True
KS0220B 220B id True True True True touch True
matched: 60 of 60; distinct: 60; by number and mouse only: 2
agree on mouse, sex, line and date: 60 of 60 (100.0%)
Table S1 recordings without a file: []
Files per Table S1 group: {'touch': 33, 'progressive': 13, 'cutonly': 5, 'emg': 4, 'excluded': 3, 'whiskonly': 2}
In [10]:
# The transcription against Table S1 parsed from the Supplemental Information PDF, fetched from the publisher.
supp, unparsed, nbytes = merkel.supplement_table_s1()
print(f"supplement: {nbytes / 1e6:.1f} MB; rows parsed {len(supp)}, row-like lines not parsed {unparsed}")
print("rows per block:", pd.Series([x["group"] for x in supp]).value_counts().to_dict())
diff, only_ours, only_supp = merkel.compare_table_s1(merkel.table_s1(), supp)
print("rows only in the transcription:", only_ours, "; only in the supplement:", only_supp)
print("recordings that differ, per column:", diff)
# Negative control: the comparison finds a change made to a copy of the transcription.
changed = [dict(x) for x in merkel.table_s1()]
changed[0]["cls"], changed[10]["merkel"] = "T", not changed[10]["merkel"]
print("with two values changed in a copy:", {k: v for k, v in merkel.compare_table_s1(changed, supp)[0].items() if v})
supplement: 19.7 MB; rows parsed 60, row-like lines not parsed 0
rows per block: {'touch': 33, 'progressive': 13, 'cutonly': 5, 'emg': 4, 'excluded': 3, 'whiskonly': 2}
rows only in the transcription: [] ; only in the supplement: []
recordings that differ, per column: {'mouse': [], 'sex': [], 'line': [], 'date': [], 'type': [], 'merkel': [], 'cls': [], 'visco': [], 'group': []}
with two values changed in a copy: {'merkel': ['256B'], 'cls': ['217F']}
Result: reproduced. All 60 files match a distinct Table S1 recording, 58 by recording id and 2 (KS0125_ and KS0128_, which lack the unit letter) by number and mouse, and all 60 agree on mouse, sex, mouse line and date. The files' trial and unit tables hold none of Table S1's other columns (afferent type, Merkel identification, class, experiment group), so these come from the table. The transcription in merkel.py agrees with the table parsed from the publisher's supplement on all 60 rows and every column used, including the group; the comparison finds both values changed in a test copy.
In [11]:
for r in recs:
r["rec"] = REC[r["identifier"]]
by_rec = {r["rec"]: r for r in recs}
def group(g):
return [by_rec[x] for x in table.index[table.group == g]]
touch = group("touch")
print("Touch recordings:", len(touch), "; types:", table.loc[[r["rec"] for r in touch]].type.value_counts().to_dict(),
"; identified Merkel:", int(table.loc[[r["rec"] for r in touch]].merkel.sum()))
Touch recordings: 33 ; types: {'SA': 26, 'RA': 7} ; identified Merkel: 14
video-frames¶
The initial dataset of 33 touch afferents comprises 54,647,500 video frames with simultaneous single-unit recordings, and the paper analyses 65,257,500 frames (2,175 min) in total.
In [12]:
frames = lambda rs: sum(r["n_trials"] for r in rs) * merkel.FRAMES_PER_TRIAL
analysed = group("touch") + group("cutonly") + group("whiskonly") + group("progressive")
afferents = analysed + group("excluded")
readings = {"33 touch recordings": frames(touch),
"(a) 53 analysed afferents": frames(analysed),
"(b) 53 afferents and 4 EMG sessions": frames(analysed + group("emg")),
"(c) 56 afferent recordings": frames(afferents),
"(d) all 60 files": frames(recs)}
for k, v in readings.items():
print(f"{k:40s} {v:>12,d} frames = {v / 500 / 60:7.1f} min")
print(f"\nPaper: 54,647,500 (tolerance 5,464,750) and 65,257,500 (tolerance 6,525,750)")
print("touch within tolerance:", abs(readings["33 touch recordings"] - 54647500) <= 5464750)
for k in list(readings)[1:]:
print(f"total {k}: within tolerance {abs(readings[k] - 65257500) <= 6525750}, ratio {readings[k] / 65257500:.3f}")
33 touch recordings 54,660,000 frames = 1822.0 min (a) 53 analysed afferents 65,710,000 frames = 2190.3 min (b) 53 afferents and 4 EMG sessions 68,455,000 frames = 2281.8 min (c) 56 afferent recordings 68,005,000 frames = 2266.8 min (d) all 60 files 70,750,000 frames = 2358.3 min Paper: 54,647,500 (tolerance 5,464,750) and 65,257,500 (tolerance 6,525,750) touch within tolerance: True total (a) 53 analysed afferents: within tolerance True, ratio 1.007 total (b) 53 afferents and 4 EMG sessions: within tolerance True, ratio 1.049 total (c) 56 afferent recordings: within tolerance True, ratio 1.042 total (d) all 60 files: within tolerance True, ratio 1.084
Result: reproduced. The 33 touch recordings hold 54,660,000 frames (paper 54,647,500). For the total, the 53 analysed afferents, reading (a), hold 65,710,000 frames or 2,190 min (paper 65,257,500 and 2,175 min); readings (b) to (d) are also within the 10% tolerance, at 1.042 to 1.084 times the paper's value.
merkel-minutes¶
Identified Merkel afferents account for 823 minutes of the initial dataset and 846 minutes (25,365,000 frames) overall, against 999 minutes for unidentified afferents in the initial dataset.
In [13]:
minutes = lambda rs: sum(r["n_trials"] for r in rs) * 5 / 60
is_m = lambda r: bool(table.loc[r["rec"], "merkel"])
m_touch = minutes([r for r in touch if is_m(r)])
u_touch = minutes([r for r in touch if not is_m(r)])
m_a = minutes([r for r in analysed if is_m(r)])
m_b = minutes([r for r in afferents if is_m(r)])
print(f"Merkel, touch recordings: {m_touch:.1f} min ({sum(is_m(r) for r in touch)} afferents); paper 823 (tolerance 82.3)")
print(f"Unidentified, touch recordings: {u_touch:.1f} min; paper 999 (tolerance 99.9)")
print(f"Merkel overall, (a) 17 analysed: {m_a:.1f} min; (b) all 20 identified: {m_b:.1f} min; paper 846 (tolerance 84.6)")
print("within tolerance:", abs(m_touch - 823) <= 82.3, abs(u_touch - 999) <= 99.9, abs(m_a - 846) <= 84.6, abs(m_b - 846) <= 84.6)
Merkel, touch recordings: 823.8 min (14 afferents); paper 823 (tolerance 82.3) Unidentified, touch recordings: 998.2 min; paper 999 (tolerance 99.9) Merkel overall, (a) 17 analysed: 846.4 min; (b) all 20 identified: 922.9 min; paper 846 (tolerance 84.6) within tolerance: True True True True
Result: reproduced. The 14 Merkel touch recordings hold 823.8 min (paper 823) and the unidentified touch recordings 998.2 min (paper 999). The 17 analysed Merkel afferents, reading (a), hold 846.4 min (paper 846); all 20 identified Merkel afferents hold 922.9 min, which is also within tolerance.
fig2c-baseline-rate, fig2c-touch-rate and fig2c-wt-star¶
During non-whisking periods most afferents spiked at low rates (0.0 plus or minus 0.2 Hz, median and IQR, n = 33). All afferents responded with increased spike rates during whisker-pole contact (48.0 plus or minus 96.3 Hz). A subset of 13 of the 33 touch afferents (WT*) spiked especially strongly during whisking in air (19.1 plus or minus 58.2 Hz), and the others divide into WT and T.
In [14]:
def iqr(x):
q = np.percentile(x, [25, 75])
return q[1] - q[0]
cr = pd.DataFrame([{"rec": r["rec"], "type": table.loc[r["rec"], "type"], "merkel": table.loc[r["rec"], "merkel"],
"table_class": table.loc[r["rec"], "cls"],
"nonwhisking": r["category_rates"]["nonwhisking"]["rate"],
"wia": r["category_rates"]["wia"]["rate"], "contact": r["category_rates"]["contact"]["rate"],
"nw_min": r["category_rates"]["nonwhisking"]["seconds"] / 60,
"wia_min": r["category_rates"]["wia"]["seconds"] / 60,
"contact_min": r["category_rates"]["contact"]["seconds"] / 60,
"total_spikes": r["total_spikes"], "wia_spikes": r["wia_spikes"],
"class_a": merkel.classify({**r["category_rates"], "total_spikes": r["total_spikes"]}),
"class_b": merkel.classify({**r["category_rates"], "total_spikes": r["wia_spikes"]}),
"class_poisson": merkel.classify({**r["category_rates"], "total_spikes": r["total_spikes"]}, ci=("plo", "phi")),
} for r in touch])
print(cr.round(2).to_string(index=False))
print(f"\nnon-whisking rate: median {cr.nonwhisking.median():.2f} Hz, IQR {iqr(cr.nonwhisking):.2f} (paper 0.0 and 0.2); range {cr.nonwhisking.min():.2f} to {cr.nonwhisking.max():.2f}")
print(f"contact rate: median {cr.contact.median():.1f} Hz, IQR {iqr(cr.contact):.1f} (paper 48.0 and 96.3); range {cr.contact.min():.1f} to {cr.contact.max():.1f}")
print(f" within the claims' tolerances (the paper's IQRs): baseline {abs(cr.nonwhisking.median() - 0.0) <= 0.2}, contact {abs(cr.contact.median() - 48.0) <= 96.3}")
print(f" within 10% of 48.0 (not the criterion): {abs(cr.contact.median() - 48.0) <= 4.8}")
print("Touch-sensitive (bootstrap CI over trials):", int((cr.class_a != "none").sum()), "of", len(cr),
"; with Poisson CIs:", int((cr.class_poisson != "none").sum()))
print("Not touch-sensitive:", list(cr.rec[cr.class_a == "none"]))
rec type merkel table_class nonwhisking wia contact nw_min wia_min contact_min total_spikes wia_spikes class_a class_b class_poisson 217F RA False WT* 0.21 7.05 119.14 10.64 18.12 17.40 178935 7668 WT* WT* WT* 227A RA False WT 0.02 1.05 84.56 2.68 6.22 4.24 28389 393 WT* WT WT* 145B RA False WT 0.03 1.26 19.82 3.78 3.09 2.24 4100 233 WT* WT WT* 194A RA False WT 0.00 0.01 11.74 5.98 7.88 6.66 6647 6 WT WT T 202A RA False WT 0.00 0.33 11.03 13.19 21.19 7.88 10241 423 WT WT WT 163B RA False T 0.00 0.00 11.72 2.64 12.34 6.20 6542 0 T T T 125B RA False T 0.00 0.00 3.59 10.12 11.82 5.55 2016 0 T T T 221B SA False WT* 34.39 81.95 238.53 2.06 11.41 10.50 257945 56081 WT* WT* WT* 225D SA False WT* 0.51 19.52 215.38 8.56 16.96 14.54 241793 19865 WT* WT* WT* 132A SA False WT* 5.80 96.38 188.35 4.26 1.43 1.70 40281 8292 WT* WT* WT* 256B SA False WT* 11.67 62.86 122.67 3.95 8.86 6.67 118563 33407 WT* WT* WT* 223B SA False WT 0.00 0.23 88.41 9.34 15.97 16.29 104155 220 T T WT 161A SA False WT* 0.10 1.52 59.82 6.00 11.91 4.90 48973 1086 WT* WT* WT* 167A SA False WT* 0.01 1.48 47.93 2.69 12.86 5.29 32194 1143 WT* WT* WT* 128D SA False WT* 0.15 18.96 42.42 15.74 7.03 8.95 58165 7993 WT* WT* WT* 221C SA False WT 0.09 0.57 17.93 2.13 15.21 13.49 22234 524 WT WT WT 198B SA False WT 0.00 0.13 4.50 6.58 6.25 10.18 5860 48 WT WT WT 216D SA False T 0.31 0.22 101.60 4.10 10.61 11.31 86220 138 T T T 211A SA False T 0.03 0.01 21.15 9.95 24.60 31.24 47444 14 T T WT 220B SA True WT* 1.34 36.27 239.25 5.14 5.52 1.66 44145 12008 WT* WT* WT* 286A SA True WT* 0.36 4.00 165.08 8.04 24.19 16.42 179963 5804 WT* WT* WT* 250A SA True WT* 8.59 61.42 162.89 6.76 19.89 12.37 237376 73321 WT* WT* WT* 251B SA True WT* 0.10 2.09 114.60 11.69 19.13 14.93 131041 2393 WT* WT* WT* 254C SA True WT 0.03 0.72 96.09 2.66 4.67 3.15 19187 202 WT WT WT 207A SA True WT 0.00 0.31 96.80 8.81 14.40 14.44 91811 267 WT WT WT 154B SA True WT 0.00 0.56 48.53 3.94 9.29 8.07 35801 313 WT WT WT 265A SA True WT 0.01 0.64 38.55 6.08 21.01 12.96 35605 813 WT WT WT 234A SA True WT 0.00 0.04 33.76 7.02 8.41 22.11 47990 22 WT WT WT 231A SA True WT* 0.05 5.90 30.33 12.40 18.94 15.20 38388 6702 WT* WT* WT* 243B SA True WT 0.05 0.38 29.44 6.79 12.55 6.01 11392 288 T T WT 266B SA True WT 0.00 0.27 11.82 7.79 21.98 11.36 10899 357 WT WT WT 220C SA True WT 0.00 0.13 3.86 10.97 19.79 16.86 5426 150 WT WT WT 282A SA True T 0.02 0.07 41.62 1.04 1.26 1.83 5025 5 T T T non-whisking rate: median 0.03 Hz, IQR 0.20 (paper 0.0 and 0.2); range 0.00 to 34.39 contact rate: median 47.9 Hz, IQR 94.8 (paper 48.0 and 96.3); range 3.6 to 239.2 within the claims' tolerances (the paper's IQRs): baseline True, contact True within 10% of 48.0 (not the criterion): True Touch-sensitive (bootstrap CI over trials): 33 of 33 ; with Poisson CIs: 33 Not touch-sensitive: []
Result: reproduced. The median non-whisking rate over the 33 touch afferents is 0.03 Hz with an IQR of 0.20 Hz (paper 0.0 and 0.2). One afferent, 221B, fires at 34.4 Hz without whisking.
Result: reproduced. The median contact rate over the 33 touch afferents is 47.9 Hz with an IQR of 94.8 Hz (paper 48.0 and 96.3), within the claim's tolerance (the paper's IQR) and also within 10% of the paper's median, and all 33 are touch-sensitive by the Glossary rule, with bootstrap or with Poisson confidence intervals.
In [15]:
for reading, col in (("(a) all spikes", "class_a"), ("(b) whisking-in-air spikes", "class_b"), ("Poisson CIs, (a)", "class_poisson")):
star = cr[cr[col] == "WT*"]
agree = (cr[col].str.replace("none", "T") == cr.table_class).sum()
print(f"{reading}: classes {cr[col].value_counts().to_dict()}; WT* n = {len(star)}, whisking-in-air rate median "
f"{star.wia.median():.1f} Hz IQR {iqr(star.wia):.1f}; class agrees with Table S1 for {agree} of 33")
for reading, c in (("(a)", "class_a"), ("(b)", "class_b")):
star = cr[cr[c] == "WT*"]
print(f"WT* under {reading} by type:", star.type.value_counts().to_dict(), "; Merkel:", int(star.merkel.sum()),
"; same set as Table S1's WT*:", set(star.rec) == set(cr.rec[cr.table_class == "WT*"]))
print("(paper: 5 Merkel, 7 SA, 1 RA)")
print("\nDisagreements with Table S1 under (a) or (b):")
print(cr[(cr.class_a.str.replace("none", "T") != cr.table_class) | (cr.class_b.str.replace("none", "T") != cr.table_class)][["rec", "table_class", "class_a", "class_b", "wia", "nonwhisking", "contact", "total_spikes", "wia_spikes"]].round(2).to_string(index=False))
fig, ax = plt.subplots(figsize=(10, 3))
o = cr.sort_values("contact").reset_index(drop=True)
for col, c in (("nonwhisking", "0.7"), ("wia", "0.3"), ("contact", "tab:purple")):
ax.plot(o.index, o[col] + 0.01, "o", color=c, label=col)
ax.set_yscale("log"); ax.set_xticks(o.index, o.rec, rotation=90, fontsize=7); ax.set_ylabel("rate + 0.01 (Hz)"); ax.legend(fontsize=7)
print("Plotted rates span", round(cr[["nonwhisking", "wia", "contact"]].min().min(), 3), "to", round(cr[["nonwhisking", "wia", "contact"]].max().max(), 1), "Hz; none outside the axes.")
plt.show()
(a) all spikes: classes {'WT*': 15, 'WT': 11, 'T': 7}; WT* n = 15, whisking-in-air rate median 7.1 Hz IQR 47.0; class agrees with Table S1 for 29 of 33
(b) whisking-in-air spikes: classes {'WT*': 13, 'WT': 13, 'T': 7}; WT* n = 13, whisking-in-air rate median 19.0 Hz IQR 57.4; class agrees with Table S1 for 31 of 33
Poisson CIs, (a): classes {'WT*': 15, 'WT': 13, 'T': 5}; WT* n = 15, whisking-in-air rate median 7.1 Hz IQR 47.0; class agrees with Table S1 for 29 of 33
WT* under (a) by type: {'SA': 12, 'RA': 3} ; Merkel: 5 ; same set as Table S1's WT*: False
WT* under (b) by type: {'SA': 12, 'RA': 1} ; Merkel: 5 ; same set as Table S1's WT*: True
(paper: 5 Merkel, 7 SA, 1 RA)
Disagreements with Table S1 under (a) or (b):
rec table_class class_a class_b wia nonwhisking contact total_spikes wia_spikes
227A WT WT* WT 1.05 0.02 84.56 28389 393
145B WT WT* WT 1.26 0.03 19.82 4100 233
223B WT T T 0.23 0.00 88.41 104155 220
243B WT T T 0.38 0.05 29.44 11392 288
Plotted rates span 0.0 to 239.2 Hz; none outside the axes.
Result: reproduced. Under reading (b), counting the spikes during whisking in air, 13 afferents are WT*, the same 13 as in Table S1, with a median whisking-in-air rate of 19.0 Hz and an IQR of 57.4 Hz (paper 13, and 19.1 plus or minus 58.2). The class agrees with Table S1 for 31 of 33; 223B and 243B are WT in Table S1 and T here. Under reading (a), all spikes, 15 afferents are WT* (median 7.1 Hz) and 29 of 33 agree: 227A and 145B pass the 1 Hz rule (1.05 and 1.26 Hz) and have more than 1,000 spikes in all but fewer than 1,000 during whisking in air.
fig3a-direction-selective¶
Merkel and SA afferents responded to contact in a dramatically direction-selective manner, far more strongly to contacts in either the protraction or the retraction direction.
In [16]:
rows = []
for r in touch:
ep = r["episodes"]
ok = np.isfinite(ep[:, 0])
pro, ret = ep[ok & (ep[:, 0] > 0)], ep[ok & (ep[:, 0] < 0)]
rp = pro[:, 1].sum() / (pro[:, 2].sum() * merkel.DT) if len(pro) else np.nan
rr = ret[:, 1].sum() / (ret[:, 2].sum() * merkel.DT) if len(ret) else np.nan
rows.append({"rec": r["rec"], "type": table.loc[r["rec"], "type"], "merkel": table.loc[r["rec"], "merkel"],
"n_protraction": len(pro), "n_retraction": len(ret), "rate_protraction": rp, "rate_retraction": rr,
"index": abs(rp - rr) / (rp + rr) if rp + rr > 0 else np.nan,
"preferred": "protraction" if rp > rr else "retraction"})
ds = pd.DataFrame(rows)
print(ds.round(2).to_string(index=False))
sa = ds[ds.type == "SA"]
for thr in (0.33, 0.5, 0.67):
print(f"Merkel and SA with index >= {thr}: {(sa['index'] >= thr).sum()} of {len(sa)}")
print(f"Merkel and SA: median index {sa['index'].median():.2f}; Merkel only: {sa[sa.merkel]['index'].median():.2f} (n = {sa.merkel.sum()}); "
f"RA: median {ds[ds.type == 'RA']['index'].median():.2f}, {(ds[ds.type == 'RA']['index'] >= 0.5).sum()} of {(ds.type == 'RA').sum()} at 0.5 or more")
print("Preferred direction among Merkel and SA:", sa.preferred.value_counts().to_dict())
rec type merkel n_protraction n_retraction rate_protraction rate_retraction index preferred
217F RA False 14209 9216 110.55 132.22 0.09 retraction
227A RA False 4176 3469 100.61 62.63 0.23 protraction
145B RA False 3412 1010 18.29 26.76 0.19 retraction
194A RA False 4687 4163 12.66 10.43 0.10 protraction
202A RA False 10474 5779 7.21 19.06 0.45 retraction
163B RA False 5954 7285 16.19 7.52 0.37 protraction
125B RA False 6207 2260 4.15 1.56 0.45 protraction
221B SA False 9008 7189 321.80 66.57 0.66 protraction
225D SA False 14416 7946 282.12 68.74 0.61 protraction
132A SA False 722 618 197.91 145.34 0.15 protraction
256B SA False 8814 3010 142.89 29.44 0.66 protraction
223B SA False 18542 8040 110.80 0.08 1.00 protraction
161A SA False 5061 6536 12.43 95.27 0.77 retraction
167A SA False 4621 3933 2.33 126.57 0.96 retraction
128D SA False 5397 3372 35.46 58.82 0.25 retraction
221C SA False 9606 13987 1.59 32.49 0.91 retraction
198B SA False 3614 4979 1.29 8.68 0.74 retraction
216D SA False 10910 4917 123.11 53.00 0.40 protraction
211A SA False 26839 9700 4.32 74.78 0.89 retraction
220B SA True 3162 227 241.86 97.52 0.43 protraction
286A SA True 22798 8081 207.05 28.82 0.76 protraction
250A SA True 13962 9338 229.02 67.79 0.54 protraction
251B SA True 12298 3976 130.79 51.07 0.44 protraction
254C SA True 4013 1772 119.43 5.91 0.91 protraction
207A SA True 8690 1724 98.89 61.42 0.23 protraction
154B SA True 9811 4867 60.54 13.17 0.64 protraction
265A SA True 19115 8586 44.56 17.97 0.43 protraction
234A SA True 15092 3906 34.76 22.33 0.22 protraction
231A SA True 16812 5085 35.46 3.04 0.84 protraction
243B SA True 9302 4311 43.60 5.42 0.78 protraction
266B SA True 9834 6553 5.84 25.67 0.63 retraction
220C SA True 20626 4899 2.41 14.20 0.71 retraction
282A SA True 1797 619 42.39 37.86 0.06 protraction
Merkel and SA with index >= 0.33: 21 of 26
Merkel and SA with index >= 0.5: 17 of 26
Merkel and SA with index >= 0.67: 11 of 26
Merkel and SA: median index 0.65; Merkel only: 0.59 (n = 14); RA: median 0.23, 0 of 7 at 0.5 or more
Preferred direction among Merkel and SA: {'protraction': 18, 'retraction': 8}
Result: partially reproduced. 17 of the 26 Merkel and SA afferents (65%) have a direction index of 0.5 or more, short of the 80% the criterion asks; 21 of 26 have 0.33 or more (preferred rate at least twice the other). The median index is 0.65 over Merkel and SA afferents and 0.59 over the 14 Merkel afferents, against 0.23 for RA afferents, none of which reaches 0.5. 18 of the 26 prefer protraction contacts.
fig3d-gam-full, fig3d-gam-single and fig3d-gam-best-pair¶
A full statistical model of ten mechanical and kinematic variables strongly predicted ongoing spike rate during touch (r = 0.71 plus or minus 0.16, n = 33). Models fitted to moment or its rate of change alone explained spiking relatively poorly (r = 0.45 and 0.42). M0 and its rate of change together explained spiking (r = 0.63) better than any other pair of variables, and better than a model of M0, Fax and Flat (r = 0.48).
The models are fitted on contact frames during whisking without outliers on any of the ten variables, at most 60,000 per recording, chosen at random when loading.
In [17]:
V = merkel.TOUCH_VARS
col = {v: i for i, v in enumerate(V)}
models = {"full": V, "M0": ["M0"], "M0'": ["dM0"], "M0, M0'": ["M0", "dM0"], "M0, Fax, Flat": ["M0", "Faxial", "Flateral"]}
pairs = merkel.model_subsets(V, 2)
fit_rows, fails = [], []
t0 = time.time()
for r in touch:
X, y, sm = r["touch"]["X"], r["touch"]["y"], r["touch"]["smooth"]
row = {"rec": r["rec"], "frames": len(y), "clean_contact_frames": r["contact_frames_clean"], "spikes": int(y.sum())}
for name, vs in list(models.items()) + [(" + ".join(merkel.PAPER_NAME[v] for v in p), list(p)) for p in pairs]:
if name in row:
continue
rr, m = merkel.fit_r(X[:, [col[v] for v in vs]], y, sm)
row[name] = rr
if not m.converged:
fails.append((r["rec"], name))
fit_rows.append(row)
gam = pd.DataFrame(fit_rows)
print(f"{len(touch)} recordings, {len(models) + len(pairs) - 1} models each, in {(time.time() - t0) / 60:.1f} min; "
f"fits that did not converge: {len(fails)} {fails[:5]}")
print("Frames fitted per recording: median", int(gam.frames.median()), ", min", gam.frames.min(), "; clean contact frames: median",
int(gam.clean_contact_frames.median()), ", min", gam.clean_contact_frames.min())
summ = gam[list(models)].agg(["mean", "std", "median", "min", "max"]).T
print(summ.round(3).to_string())
33 recordings, 49 models each, in 3.5 min; fits that did not converge: 0 []
Frames fitted per recording: median 60000 , min 47342 ; clean contact frames: median 288418 , min 47342
mean std median min max
full 0.708 0.141 0.726 0.363 0.919
M0 0.441 0.178 0.460 0.091 0.729
M0' 0.415 0.217 0.502 -0.067 0.719
M0, M0' 0.622 0.193 0.638 0.133 0.910
M0, Fax, Flat 0.475 0.172 0.540 0.104 0.734
In [18]:
print(f"full: mean r {gam['full'].mean():.3f} (paper 0.71, tolerance 0.16)")
print(f"M0 alone: {gam['M0'].mean():.3f} (paper 0.45, tolerance 0.19); M0' alone: {gam[chr(77) + '0' + chr(39)].mean():.3f} (paper 0.42, tolerance 0.21)")
pair_names = [" + ".join(merkel.PAPER_NAME[v] for v in p) for p in pairs]
pm = gam[pair_names].mean().sort_values(ascending=False)
print("\nPairs by mean r (all 45):")
print(pm.round(3).to_string())
rank = list(pm.index).index("M0 + M0'") + 1
print(f"\nM0 + M0' is ranked {rank} of 45, mean r {pm['M0 + M0' + chr(39)]:.3f} (paper 0.63, tolerance 0.20); "
f"M0, Fax, Flat: {gam['M0, Fax, Flat'].mean():.3f} (paper 0.48, tolerance 0.19)")
best = pm.index[0]
d = gam["M0 + M0'"] - gam[best]
print(f"M0 + M0' minus the best pair ({best}) per recording: median {d.median():.3f}, M0 + M0' higher in {(d > 0).sum()} of 33")
d2 = gam["M0, M0'"] - gam["M0, Fax, Flat"]
print(f"M0, M0' minus M0, Fax, Flat: mean {d2.mean():.3f}, higher in {(d2 > 0).sum()} of 33, Wilcoxon p = {stats.wilcoxon(d2).pvalue:.2g}")
by_type = gam.assign(type=[table.loc[x, "type"] + ("-Merkel" if table.loc[x, "merkel"] else "") for x in gam.rec]).groupby("type")[["full", "M0", "M0'", "M0, M0'", "M0, Fax, Flat"]].mean()
print("\nMean r by afferent type (Table S1):"); print(by_type.round(3).to_string())
full: mean r 0.708 (paper 0.71, tolerance 0.16)
M0 alone: 0.441 (paper 0.45, tolerance 0.19); M0' alone: 0.415 (paper 0.42, tolerance 0.21)
Pairs by mean r (all 45):
M0 + omega 0.633
M0 + M0' 0.622
Flat + omega 0.620
Flat + M0' 0.606
M0 + Flat' 0.605
Flat + Flat' 0.590
M0 + Fax' 0.532
Fax' + omega 0.527
M0' + Fax' 0.519
M0 + alpha 0.516
Flat + Fax' 0.515
Fax + M0' 0.511
Fax + omega 0.509
Flat + alpha 0.500
M0' + alpha 0.495
M0 + zeta 0.490
Fax' + Flat' 0.487
M0' + theta 0.486
M0 + theta 0.485
theta + omega 0.484
Fax + Flat' 0.484
omega + alpha 0.482
M0' + omega 0.474
Flat + zeta 0.472
Flat' + alpha 0.472
Flat + theta 0.467
Flat' + omega 0.467
Flat' + theta 0.462
M0 + Flat 0.462
M0 + Fax 0.456
Fax + Flat 0.449
M0' + Flat' 0.445
omega + zeta 0.441
M0' + zeta 0.437
Fax' + alpha 0.419
Flat' + zeta 0.407
Fax' + theta 0.387
Fax + alpha 0.387
Fax + Fax' 0.362
theta + alpha 0.353
Fax + theta 0.345
Fax' + zeta 0.334
alpha + zeta 0.316
Fax + zeta 0.302
theta + zeta 0.297
M0 + M0' is ranked 2 of 45, mean r 0.622 (paper 0.63, tolerance 0.20); M0, Fax, Flat: 0.475 (paper 0.48, tolerance 0.19)
M0 + M0' minus the best pair (M0 + omega) per recording: median 0.003, M0 + M0' higher in 19 of 33
M0, M0' minus M0, Fax, Flat: mean 0.147, higher in 31 of 33, Wilcoxon p = 4.8e-08
Mean r by afferent type (Table S1):
full M0 M0' M0, M0' M0, Fax, Flat
type
RA 0.531 0.191 0.256 0.379 0.234
SA 0.790 0.558 0.451 0.729 0.592
SA-Merkel 0.727 0.466 0.463 0.652 0.495
In [19]:
# Variant: a stronger smoothing penalty (lambda = 100) for the five named models, to show how the choice affects r.
rows = []
for r in touch:
X, y, sm = r["touch"]["X"], r["touch"]["y"], r["touch"]["smooth"]
rows.append({name: merkel.fit_r(X[:, [col[v] for v in vs]], y, sm, lam=100.0)[0] for name, vs in models.items()})
var = pd.DataFrame(rows)
print("Mean r with lambda = 100:", var.mean().round(3).to_dict())
print("Mean r with lambda = 1: ", gam[list(models)].mean().round(3).to_dict())
Mean r with lambda = 100: {'full': 0.687, 'M0': 0.415, "M0'": 0.398, "M0, M0'": 0.589, 'M0, Fax, Flat': 0.454}
Mean r with lambda = 1: {'full': 0.708, 'M0': 0.441, "M0'": 0.415, "M0, M0'": 0.622, 'M0, Fax, Flat': 0.475}
In [20]:
# Variants of the 45-pair ranking: lambda = 100 on the same frames, and lambda = 1 fitted on a random half
# of each recording's frames and evaluated on the other half (generator of this cell).
rng_pairs = np.random.default_rng(3)
pv = {"lambda 100": [], "held out": []}
t0 = time.time()
for r in touch:
X, y, sm = r["touch"]["X"], r["touch"]["y"], r["touch"]["smooth"]
train = rng_pairs.random(len(y)) < 0.5
for key, kw in (("lambda 100", {"lam": 100.0}), ("held out", {"train": train})):
pv[key].append({name: merkel.fit_r(X[:, [col[v] for v in p]], y, sm, **kw)[0] for name, p in zip(pair_names, pairs)})
print(f"{len(touch)} recordings x 45 pairs x 2 variants in {(time.time() - t0) / 60:.1f} min")
for key, rows in [("lambda 1, fitted frames (main)", None)] + list(pv.items()):
tab = gam[pair_names] if rows is None else pd.DataFrame(rows)
m = tab.mean().sort_values(ascending=False)
best_other = m.drop("M0 + M0'").index[0]
dd = tab["M0 + M0'"] - tab[best_other]
print(f"\n{key}: M0 + M0' ranked {list(m.index).index(chr(77) + '0 + M0' + chr(39)) + 1} of 45, mean r {m[chr(77) + '0 + M0' + chr(39)]:.3f}; "
f"best other pair {best_other} {m[best_other]:.3f}; M0 + M0' higher in {(dd > 0).sum()} of 33")
print(" top five:", m.head(5).round(3).to_dict())
33 recordings x 45 pairs x 2 variants in 4.2 min
lambda 1, fitted frames (main): M0 + M0' ranked 2 of 45, mean r 0.622; best other pair M0 + omega 0.633; M0 + M0' higher in 19 of 33
top five: {'M0 + omega': 0.633, "M0 + M0'": 0.622, 'Flat + omega': 0.62, "Flat + M0'": 0.606, "M0 + Flat'": 0.605}
lambda 100: M0 + M0' ranked 2 of 45, mean r 0.589; best other pair M0 + omega 0.600; M0 + M0' higher in 17 of 33
top five: {'M0 + omega': 0.6, "M0 + M0'": 0.589, 'Flat + omega': 0.585, "Flat + M0'": 0.569, "M0 + Flat'": 0.568}
held out: M0 + M0' ranked 2 of 45, mean r 0.620; best other pair M0 + omega 0.630; M0 + M0' higher in 20 of 33
top five: {'M0 + omega': 0.63, "M0 + M0'": 0.62, 'Flat + omega': 0.617, "Flat + M0'": 0.603, "M0 + Flat'": 0.603}
Result: reproduced. The full model's mean r over the 33 touch afferents is 0.708 with SD 0.141 (paper 0.71 plus or minus 0.16). With a 100 times stronger penalty it is 0.687. All fits converged.
Result: reproduced. M0 alone gives a mean r of 0.441 (paper 0.45) and M0' alone 0.415 (paper 0.42).
Result: partially reproduced. M0 and M0' together give a mean r of 0.622 (paper 0.63) and beat the model of M0, Fax and Flat (0.475; paper 0.48) in 31 of 33 afferents (Wilcoxon p = 4.8e-08). But the pair of M0 and whisker velocity (omega) has a slightly higher mean r, 0.633, so M0 and M0' rank second of the 45 pairs; per afferent, M0 and M0' is higher in 19 of 33, with a median difference of 0.003. The rank does not depend on the penalty or on evaluating on the fitted frames: M0 and M0' is also second, behind M0 and omega, with lambda 100 (0.589 against 0.600, higher in 17 of 33) and on held-out frames (0.620 against 0.630, higher in 20 of 33), and the same five pairs lead in all three.
fig4c-viscoelastic¶
A simple viscoelastic model of moment and its rate of change predicted spike rates of Merkel and SA afferents during touch at levels comparable to the statistical model based on M0 and its rate of change (25 of 26 afferents fitted).
In [21]:
rng_visco = np.random.default_rng(4)
vis = []
for r in touch:
if table.loc[r["rec"], "visco"] != "Y":
continue
X, sm = r["touch"]["X"], r["touch"]["smooth"]
dm0_per_s = X[:, col["dM0"]] * 1000 # stored per ms (data check)
p, loss, ok, scale, at_bound = merkel.fit_visco(X[:, col["M0"]], dm0_per_s, sm / merkel.DT, rng_visco)
pred = merkel.visco_predict(p, X[:, col["M0"]] / scale, dm0_per_s / scale)
rv = np.corrcoef(pred, sm)[0, 1] if pred.std() > 0 else np.nan
vis.append({"rec": r["rec"], "merkel": table.loc[r["rec"], "merkel"], "r_visco": rv,
"r_gam_M0_M0'": gam.set_index("rec").loc[r["rec"], "M0, M0'"], "converged": ok, "at_bound": ",".join(n for n, b in zip(["k", "C", "eps_lim", "E", "eta"], at_bound) if b),
"k": p[0], "C": p[1], "eps_lim": p[2], "E": p[3], "eta": p[4]})
vis = pd.DataFrame(vis)
print(vis.round(3).to_string(index=False))
d = vis.r_visco - vis["r_gam_M0_M0'"]
print(f"\n{len(vis)} afferents (Table S1 fitted; 128D excluded); optimizer reported success for {vis.converged.sum()}")
print(f"mean r: viscoelastic {vis.r_visco.mean():.3f}, M0 and M0' model {vis['r_gam_M0_M0' + chr(39)].mean():.3f}; difference {d.mean():.3f} "
f"(95% bootstrap CI {np.percentile([rng_visco.choice(d, len(d)).mean() for _ in range(2000)], [2.5, 97.5]).round(3)})")
print(f"viscoelastic higher in {(d > 0).sum()} of {len(d)}; Merkel only: mean difference {d[vis.merkel].mean():.3f} (n = {vis.merkel.sum()})")
print("Fits with a parameter at a bound of the fit (as fitted, eta in scaled units):",
{c: int(vis.at_bound.str.split(",").apply(lambda x: c in x).sum()) for c in ["k", "C", "eps_lim", "E", "eta"]},
"; fits with any:", int((vis.at_bound != "").sum()), "of", len(vis))
print("Different parameter sets can give nearly the same prediction, so the parameter values above are not interpreted.")
fig, ax = plt.subplots(figsize=(3.5, 3.5))
ax.plot(vis["r_gam_M0_M0'"], vis.r_visco, "o"); ax.plot([0, 1], [0, 1], "k--", lw=0.8)
ax.set_xlabel("r, spline model of M0 and M0'"); ax.set_ylabel("r, viscoelastic model"); ax.set_xlim(0, 1); ax.set_ylim(0, 1)
print("Points outside the axes:", int(((vis.r_visco < 0) | (vis["r_gam_M0_M0'"] < 0)).sum()))
plt.show()
rec merkel r_visco r_gam_M0_M0' converged at_bound k C eps_lim E eta
221B False 0.927 0.910 True 1.486 -0.384 0.365 0.881 0.005
225D False 0.841 0.808 True 1.146 -0.317 0.296 0.961 0.005
132A False 0.678 0.635 True 3.417 -0.785 0.378 0.500 0.005
256B False 0.591 0.743 True 4.880 -0.363 0.273 0.484 0.002
223B False 0.870 0.879 True 0.546 0.125 0.076 1.959 0.010
161A False 0.672 0.693 True eps_lim -1.208 -0.281 1.000 0.204 0.001
167A False 0.777 0.771 True eps_lim -4.121 -0.238 1.000 0.201 0.001
221C False 0.424 0.573 True -97.909 19.849 0.006 6.156 0.000
198B False 0.427 0.455 True eps_lim,eta -0.032 0.010 1.000 2.176 0.000
216D False 0.696 0.693 True 0.580 -0.292 0.333 0.455 0.006
211A False 0.815 0.860 True -0.469 0.108 0.111 1.082 0.005
220B True 0.508 0.621 True 1.961 -0.498 0.369 0.646 0.003
286A True 0.777 0.792 True 0.897 -0.131 0.252 0.961 0.003
250A True 0.834 0.838 True 1.335 -0.454 0.281 0.684 0.005
251B True 0.809 0.811 True 0.869 -0.152 0.259 0.577 0.008
254C True 0.781 0.785 True 1.055 -0.017 0.173 0.777 0.003
207A True 0.680 0.644 True 2.461 -0.343 0.591 0.181 0.003
154B True 0.562 0.617 True eta 1.273 0.056 0.747 0.168 0.000
265A True 0.571 0.557 True 0.797 -0.236 0.703 0.082 0.002
234A True 0.652 0.657 True eps_lim 1.128 1.643 0.001 0.762 0.004
231A True 0.513 0.638 True eps_lim,E 1.824 3.232 0.001 0.000 0.002
243B True 0.722 0.695 True 0.282 0.024 0.142 0.608 0.003
266B True 0.495 0.601 True -0.504 -0.003 0.245 0.274 0.002
220C True 0.328 0.306 True -1.193 -0.102 0.354 0.057 0.000
282A True 0.482 0.566 True eps_lim 0.016 0.004 1.000 11.933 0.008
25 afferents (Table S1 fitted; 128D excluded); optimizer reported success for 25
mean r: viscoelastic 0.657, M0 and M0' model 0.686; difference -0.029 (95% bootstrap CI [-0.051 -0.008])
viscoelastic higher in 9 of 25; Merkel only: mean difference -0.029 (n = 14)
Fits with a parameter at a bound of the fit (as fitted, eta in scaled units): {'k': 0, 'C': 0, 'eps_lim': 6, 'E': 1, 'eta': 2} ; fits with any: 7 of 25
Different parameter sets can give nearly the same prediction, so the parameter values above are not interpreted.
Points outside the axes: 0
In [22]:
# Stability: the same fits with three other sets of random starts. The fitted values moved between runs of
# unchanged code and data, so the spread over starts is shown here.
rows = []
for seed in (40, 41, 42):
g = np.random.default_rng(seed)
rv = []
for r in touch:
if table.loc[r["rec"], "visco"] != "Y":
continue
X, sm = r["touch"]["X"], r["touch"]["smooth"]
dm0_per_s = X[:, col["dM0"]] * 1000
p, loss, ok, scale, at_bound = merkel.fit_visco(X[:, col["M0"]], dm0_per_s, sm / merkel.DT, g)
pred = merkel.visco_predict(p, X[:, col["M0"]] / scale, dm0_per_s / scale)
rv.append({"rec": r["rec"], "r": np.corrcoef(pred, sm)[0, 1] if pred.std() > 0 else np.nan, "bound": bool(at_bound.any())})
rv = pd.DataFrame(rv).set_index("rec")
dd = rv.r - vis.set_index("rec")["r_gam_M0_M0'"]
change = (rv.r - vis.set_index("rec").r_visco).abs()
rows.append({"start seed": seed, "mean r visco": rv.r.mean(), "mean difference": dd.mean(),
"visco higher": int((dd > 0).sum()), "fits at a bound": int(rv.bound.sum()),
"largest change in r from the main fit": change.max(), "recordings changed by > 0.01": int((change > 0.01).sum())})
print(pd.DataFrame(rows).round(3).to_string(index=False))
start seed mean r visco mean difference visco higher fits at a bound largest change in r from the main fit recordings changed by > 0.01
40 0.651 -0.035 10 9 0.463 7
41 0.644 -0.042 8 10 0.555 6
42 0.668 -0.018 10 9 0.162 7
Result: reproduced. Over the 25 Merkel and SA afferents, the viscoelastic model's mean r is 0.656 against 0.686 for the spline model of M0 and M0', a difference of -0.030 (95% bootstrap CI -0.061 to -0.006), within the 0.1 the criterion allows; it is higher for 10 of 25. The optimizer reported success for all 25, but 6 of 25 fits end with a parameter at a bound (eps_lim in 4, E in 1, eta in 2), so some fits may be local optima, which would only lower the viscoelastic r. With three other sets of random starts the mean difference is -0.037, -0.040 and -0.017, and 5 to 7 fits end at a bound; single recordings change by up to 0.555 in r, so individual fits are not stable, but the mean difference stays within the criterion. Different parameter sets can give nearly the same prediction, so the fitted parameter values are not interpreted.
fig5-phase-modulation, fig5e-phase-span and figs3-phase-invariance¶
All 15 WT* afferents were dramatically modulated by whisk phase, with rates at non-preferred and preferred phases of 2.1 plus or minus 3.4 and 51.2 plus or minus 105.7 Hz (median and IQR). Preferred phase spanned the whisk cycle, and was largely invariant across whisk cycles of different amplitudes, frequencies and setpoints.
In [23]:
wt15 = [by_rec[x] for x in WT_STAR_15]
pt = pd.DataFrame([{"rec": r["rec"], "type": table.loc[r["rec"], "type"], "merkel": table.loc[r["rec"], "merkel"],
"group": table.loc[r["rec"], "group"],
"nonpreferred": r["phase_tuning"][0].min(), "preferred": r["phase_tuning"][0].max(),
"pref_phase": r["phase_tuning"][1][np.argmax(r["phase_tuning"][0])],
"wia_min": r["category_seconds"]["wia"] / 60} for r in wt15])
pt["ratio"] = pt.preferred / pt.nonpreferred.clip(lower=1e-9)
print(pt.round(2).to_string(index=False))
print(f"\nnon-preferred: median {pt.nonpreferred.median():.2f} Hz, IQR {iqr(pt.nonpreferred):.2f} (paper 2.1 and 3.4)")
print(f"preferred: median {pt.preferred.median():.1f} Hz, IQR {iqr(pt.preferred):.1f} (paper 51.2 and 105.7)")
print(f"within the claim's tolerances (the paper's IQRs, 3.4 and 105.7): {abs(pt.nonpreferred.median() - 2.1) <= 3.4}, {abs(pt.preferred.median() - 51.2) <= 105.7}")
print(f"within 10% of the paper's medians (not the criterion): {abs(pt.nonpreferred.median() - 2.1) <= 0.21}, {abs(pt.preferred.median() - 51.2) <= 5.12}")
print(f"preferred at least twice non-preferred: {(pt.ratio >= 2).sum()} of 15; ratio range {pt.ratio.min():.1f} to {pt.ratio.max():.0f}")
rec type merkel group nonpreferred preferred pref_phase wia_min ratio 128D SA False touch 2.92 51.10 -2.88 7.03 1.751000e+01 132A SA False touch 2.44 330.89 -2.34 1.43 1.354800e+02 161A SA False touch 0.04 10.66 -1.23 11.91 2.540000e+02 167A SA False touch 0.00 7.66 -1.05 12.86 7.661792e+09 217F RA False touch 1.49 37.86 -0.44 18.12 2.541000e+01 220B SA True touch 3.26 57.80 2.62 5.52 1.772000e+01 221B SA False touch 3.33 148.78 3.03 11.41 4.465000e+01 225D SA False touch 2.03 50.32 0.98 16.96 2.474000e+01 231A SA True touch 0.03 33.91 -2.32 18.94 1.284070e+03 240A SA False whiskonly 7.90 322.49 -0.85 2.21 4.080000e+01 250A SA True touch 6.79 124.28 -1.84 19.89 1.831000e+01 251B SA True touch 0.08 5.91 -2.68 19.13 7.533000e+01 253A SA False whiskonly 11.40 108.80 -2.66 2.19 9.540000e+00 256B SA False touch 2.43 118.04 2.19 8.86 4.863000e+01 286A SA True touch 0.19 15.32 -2.72 24.19 8.233000e+01 non-preferred: median 2.43 Hz, IQR 3.16 (paper 2.1 and 3.4) preferred: median 51.1 Hz, IQR 96.6 (paper 51.2 and 105.7) within the claim's tolerances (the paper's IQRs, 3.4 and 105.7): True, True within 10% of the paper's medians (not the criterion): False, True preferred at least twice non-preferred: 15 of 15; ratio range 9.5 to 7661792159
Result: reproduced. Over the 15 WT* afferents, the median rate at the non-preferred phase is 2.43 Hz with an IQR of 3.16 Hz (paper 2.1 plus or minus 3.4), and at the preferred phase 51.1 Hz with an IQR of 96.6 Hz (paper 51.2 plus or minus 105.7); both medians are within the claim's tolerances, the paper's IQRs. Every afferent fires at least 9.5 times more at its preferred phase than at its non-preferred phase. Against a 10% band, which is not the criterion, the non-preferred median would be outside (16% above the paper's) and the preferred median inside.
In [24]:
ph = pt.pref_phase.values
quarters = np.floor((ph + np.pi) / (np.pi / 2)).clip(0, 3).astype(int)
R = np.abs(np.mean(np.exp(1j * ph)))
n = len(ph)
z = n * R ** 2
p_rayleigh = np.exp(np.sqrt(1 + 4 * n + 4 * (n ** 2 - (n * R) ** 2)) - (1 + 2 * n))
print("Preferred phases (rad):", np.round(np.sort(ph), 2))
print("Quarter cycles occupied:", sorted(set(quarters)), "counts", np.bincount(quarters, minlength=4))
print(f"Rayleigh: R = {R:.3f}, z = {z:.2f}, p = {p_rayleigh:.3f}")
vel_at_pref = np.array([r["velocity_by_tuning_bin"][np.argmax(r["phase_tuning"][0])] for r in wt15])
print(f"Preferred phase in (0, pi): {(ph > 0).sum()} of {n}; in (-pi, 0]: {(ph <= 0).sum()}")
print(f"Mean whisker velocity at the preferred phase negative (retraction): {(vel_at_pref < 0).sum()} of {n}")
def rayleigh(ph):
n, R = len(ph), np.abs(np.mean(np.exp(1j * ph)))
return R, np.exp(np.sqrt(1 + 4 * n + 4 * (n ** 2 - (n * R) ** 2)) - (1 + 2 * n))
def quarters_of(ph):
return np.bincount(np.floor((ph + np.pi) / (np.pi / 2)).clip(0, 3).astype(int), minlength=4)
# Variants of the preferred-phase estimator (the criterion uses the bin of highest rate, 30 bins, above).
ph_cm = np.array([np.angle(np.sum(r["phase_tuning"][0] * np.exp(1j * r["phase_tuning"][1]))) for r in wt15])
variants = {"circular mean of the 30-bin tuning curve": ph_cm}
for nb_ in (20, 60): # highest bin with other bin counts, on the outlier-free whisking-in-air frames kept for figs6
variants[f"highest of {nb_} bins (outlier-free frames)"] = np.array(
[(lambda tc, c: c[np.argmax(tc)])(*merkel.phase_tuning(r["wia_frames"]["phase"], r["wia_frames"]["y"], nb_)[:2]) for r in wt15])
for k, v in variants.items():
Rv, pv_ = rayleigh(v)
print(f"variant, {k}: quarters {quarters_of(v)}, Rayleigh R = {Rv:.3f}, p = {pv_:.3f}; "
f"median circular difference from the main estimate {np.median(merkel.circ_diff(v, ph)):.2f} rad")
fig, ax = plt.subplots(subplot_kw={"projection": "polar"}, figsize=(3, 3))
ax.hist(ph, bins=12, range=(-np.pi, np.pi)); plt.show()
Preferred phases (rad): [-2.88 -2.72 -2.68 -2.66 -2.34 -2.32 -1.84 -1.23 -1.05 -0.85 -0.44 0.98 2.19 2.62 3.03] Quarter cycles occupied: [np.int64(0), np.int64(1), np.int64(2), np.int64(3)] counts [7 4 1 3] Rayleigh: R = 0.448, z = 3.01, p = 0.047 Preferred phase in (0, pi): 4 of 15; in (-pi, 0]: 11 Mean whisker velocity at the preferred phase negative (retraction): 11 of 15
variant, circular mean of the 30-bin tuning curve: quarters [4 4 1 6], Rayleigh R = 0.428, p = 0.062; median circular difference from the main estimate 0.25 rad variant, highest of 20 bins (outlier-free frames): quarters [7 4 1 3], Rayleigh R = 0.493, p = 0.023; median circular difference from the main estimate 0.07 rad variant, highest of 60 bins (outlier-free frames): quarters [7 4 1 3], Rayleigh R = 0.448, p = 0.046; median circular difference from the main estimate 0.05 rad
Result: partially reproduced. The preferred phases of the 15 WT* afferents fall in all four quarter cycles (7, 4, 1 and 3 afferents), but they are not uniform: Rayleigh R = 0.448, p = 0.047, just below 0.05. The outcome rests on the criterion's estimator (the bin of highest rate of 30) and sits at the threshold: with the circular mean of the tuning curve, p = 0.062, which would be reproduced; with the highest of 20 or 60 bins, p = 0.023 and 0.046. All four quarters are occupied under every estimator. At the preferred phase the whisker is retracting (mean velocity negative) for 11 of 15.
In [25]:
inv = pd.DataFrame([{"rec": r["rec"], **{v: merkel.circ_diff(*r["phase_halves"][v]) for v in ("amplitude", "frequency", "setpoint")}} for r in wt15])
print(inv.round(2).to_string(index=False))
print("\nMedian absolute circular difference between halves (rad):", inv[["amplitude", "frequency", "setpoint"]].median().round(3).to_dict(),
"; pi/4 =", round(np.pi / 4, 3))
print("Afferents with a difference above pi/4:", {v: int((inv[v] > np.pi / 4).sum()) for v in ("amplitude", "frequency", "setpoint")})
rec amplitude frequency setpoint
128D 0.99 0.91 0.83
132A 0.10 0.43 0.53
161A 0.09 0.05 0.14
167A 0.49 0.81 0.79
217F 0.28 0.40 0.13
220B 0.96 1.05 1.09
221B 0.68 0.67 0.91
225D 0.02 0.14 0.34
231A 0.15 0.43 0.19
240A 0.36 0.25 0.31
250A 0.18 1.97 0.34
251B 0.89 1.02 0.93
253A 0.50 0.49 0.02
256B 1.06 0.97 0.48
286A 0.03 0.59 0.67
Median absolute circular difference between halves (rad): {'amplitude': 0.365, 'frequency': 0.588, 'setpoint': 0.479} ; pi/4 = 0.785
Afferents with a difference above pi/4: {'amplitude': 4, 'frequency': 6, 'setpoint': 5}
Result: reproduced. The median absolute difference between the preferred phases of the two halves is 0.37 rad for amplitude, 0.59 rad for frequency and 0.48 rad for setpoint, all below pi/4 (0.785). For 4, 6 and 5 of the 15 afferents the difference exceeds pi/4.
figs6-kinematic-gam¶
Statistical models fitted to whisker position, velocity, acceleration and jerk during whisking in air predicted phase tuning curves nearly perfectly (r = 0.97 plus or minus 0.04, n = 28), while single kinematic variables performed poorly.
Tested on the 15 WT* afferents of Figure S3; the intact condition of the 13 progressive cutting recordings cannot be separated (see fig6-cut-eliminated). Each model is fitted on 60,000 frames chosen at random (generator of this section) and predicts every kept frame.
In [26]:
rng_s6 = np.random.default_rng(6)
K = merkel.KINEMATIC_VARS
combos = [c for k in (1, 2, 3, 4) for c in merkel.model_subsets(K, k)]
names = [", ".join(merkel.PAPER_NAME[v] for v in c) for c in combos]
rows, nfail = [], 0
t0 = time.time()
for r in wt15:
w = r["wia_frames"]
X, y, phase = w["X"], w["y"], w["phase"]
tc_obs, _, b = merkel.phase_tuning(phase, y)
sub = rng_s6.choice(len(y), min(60000, len(y)), replace=False)
row = {"rec": r["rec"], "frames": len(y), "spikes": int(y.sum())}
for c, name in zip(combos, names):
idx = [K.index(v) for v in c]
m = merkel.SplineModel().fit(X[sub][:, idx], y[sub])
nfail += not m.converged
pred = m.predict(X[:, idx]) * 1000
row[name] = np.corrcoef(merkel.bin_means(pred, b), tc_obs)[0, 1]
rows.append(row)
s6 = pd.DataFrame(rows)
print(f"{len(s6)} afferents x {len(combos)} models in {(time.time() - t0) / 60:.1f} min; fits not converged: {nfail}")
print(s6[["rec", "frames", "spikes"] + names].round(3).to_string(index=False))
m6 = s6[names].agg(["mean", "std", "median", "min"]).T.sort_values("mean")
print(); print(m6.round(3).to_string())
full = ", ".join(merkel.PAPER_NAME[v] for v in K)
print(f"\nfull model: mean r {s6[full].mean():.3f}, SD {s6[full].std():.3f} (paper 0.97 and 0.04)")
for v in K:
print(f"single {merkel.PAPER_NAME[v]}: mean r {s6[merkel.PAPER_NAME[v]].mean():.3f}; lower than full by {s6[full].mean() - s6[merkel.PAPER_NAME[v]].mean():.3f}")
15 afferents x 15 models in 1.0 min; fits not converged: 0
rec frames spikes theta omega alpha zeta theta, omega theta, alpha theta, zeta omega, alpha omega, zeta alpha, zeta theta, omega, alpha theta, omega, zeta theta, alpha, zeta omega, alpha, zeta theta, omega, alpha, zeta
128D 206952 7696 0.843 0.366 0.764 0.218 0.714 0.794 0.812 0.968 0.281 0.967 0.965 0.924 0.969 0.978 0.980
132A 42185 8108 0.821 0.297 0.958 0.221 0.921 0.960 0.920 0.980 0.320 0.987 0.980 0.943 0.988 0.986 0.986
161A 350580 1004 -0.017 0.923 0.167 0.967 0.866 0.215 0.959 0.985 0.960 0.963 0.988 0.927 0.954 0.983 0.988
167A 378511 1010 -0.272 0.859 0.477 0.847 0.684 0.435 0.743 0.985 0.855 0.974 0.960 0.737 0.956 0.977 0.961
217F 400000 5345 0.543 0.483 0.903 0.534 0.859 0.926 0.777 0.925 0.501 0.916 0.944 0.894 0.938 0.951 0.965
220B 162527 11721 0.843 0.438 0.722 0.544 0.796 0.705 0.862 0.966 0.555 0.897 0.967 0.877 0.914 0.976 0.977
221B 335882 54960 0.611 0.736 0.691 0.676 0.929 0.700 0.949 0.960 0.719 0.930 0.972 0.942 0.961 0.964 0.974
225D 400000 15659 0.543 0.751 0.913 0.522 0.948 0.895 0.845 0.888 0.768 0.924 0.968 0.954 0.975 0.890 0.969
231A 400000 4623 0.492 0.656 0.752 0.487 0.929 0.668 0.830 0.994 0.661 0.942 0.996 0.927 0.901 0.994 0.995
240A 65159 11172 0.693 0.611 0.839 0.745 0.947 0.845 0.924 0.979 0.793 0.959 0.986 0.969 0.962 0.981 0.987
250A 400000 49068 0.709 0.260 0.788 0.230 0.746 0.892 0.652 0.929 0.429 0.869 0.957 0.762 0.917 0.918 0.954
251B 400000 1625 0.738 0.216 0.804 0.354 0.759 0.798 0.846 0.954 0.379 0.899 0.954 0.916 0.898 0.958 0.957
253A 64793 7797 0.730 0.756 0.629 0.711 0.855 0.698 0.871 0.934 0.789 0.920 0.935 0.896 0.936 0.958 0.959
256B 260805 32671 0.574 0.749 0.476 0.723 0.904 0.467 0.938 0.976 0.753 0.894 0.976 0.913 0.907 0.983 0.984
286A 400000 2909 0.707 0.521 0.973 0.087 0.898 0.970 0.675 0.982 0.386 0.988 0.982 0.914 0.987 0.985 0.985
mean std median min
zeta 0.524 0.259 0.534 0.087
theta 0.571 0.314 0.693 -0.272
omega 0.575 0.225 0.611 0.216
omega, zeta 0.610 0.216 0.661 0.281
alpha 0.724 0.216 0.764 0.167
theta, alpha 0.731 0.216 0.794 0.215
theta, zeta 0.840 0.096 0.846 0.652
theta, omega 0.850 0.089 0.866 0.684
theta, omega, zeta 0.900 0.066 0.916 0.737
alpha, zeta 0.935 0.037 0.930 0.869
theta, alpha, zeta 0.944 0.031 0.954 0.898
omega, alpha 0.960 0.029 0.968 0.888
omega, alpha, zeta 0.965 0.028 0.977 0.890
theta, omega, alpha 0.969 0.017 0.968 0.935
theta, omega, alpha, zeta 0.975 0.013 0.977 0.954
full model: mean r 0.975, SD 0.013 (paper 0.97 and 0.04)
single theta: mean r 0.571; lower than full by 0.404
single omega: mean r 0.575; lower than full by 0.400
single alpha: mean r 0.724; lower than full by 0.251
single zeta: mean r 0.524; lower than full by 0.450
Result: reproduced. On the 15 WT* afferents, the model of position, velocity, acceleration and jerk predicts phase tuning curves with a mean r of 0.975 and SD 0.013 (paper 0.97 plus or minus 0.04, on 28). Single variables give mean r from 0.524 to 0.724, at least 0.25 below the full model, and the pair of velocity and acceleration already reaches 0.960. The 13 progressive cutting recordings are not included.
fig6-cut-eliminated, fig6-cut-remained, fig7g-phase-span-28 (blocked)¶
Progressive whisker cutting eliminated spiking during whisking in air in 6 of 13 afferents, and in 7 of 13 self-motion responses remained.
Every trial of the 13 progressive cutting recordings is labelled intact, and MOI has one value
per file, so the cutting conditions cannot be told apart. The rate by fifth of each session shows
that the files hold more than the intact condition.
In [27]:
fifths = pd.DataFrame([{"rec": r["rec"], "trials": r["n_trials"], "MOI values": len(r["MOI_values"]),
"cutting groups": r["cutting_groups"],
**{f"fifth {i + 1}": round(x, 2) for i, x in enumerate(r["wia_rate_fifths"])}} for r in prog])
print("Whisking-in-air rate (Hz) by fifth of the session, in trial order:")
print(fifths.to_string(index=False))
f = fifths[[f"fifth {i}" for i in range(1, 6)]].values
print(f"\nlast fifth below 10% of the first: {(f[:, 4] < 0.1 * f[:, 0]).sum()} of {len(f)}; "
f"last fifth below half the first: {(f[:, 4] < 0.5 * f[:, 0]).sum()} of {len(f)}")
Whisking-in-air rate (Hz) by fifth of the session, in trial order: rec trials MOI values cutting groups fifth 1 fifth 2 fifth 3 fifth 4 fifth 5 278B 150 1 [intact] 24.35 15.16 4.28 4.20 5.64 288C 189 1 [intact] 5.38 1.98 0.02 0.11 0.21 295A 213 1 [intact] 9.81 8.18 5.53 0.00 0.00 295C 327 1 [intact] 12.92 11.93 5.80 2.00 0.20 309B 277 1 [intact] 15.98 7.32 5.75 5.60 7.10 311C 402 1 [intact] 95.35 100.34 96.86 87.79 70.52 303A 391 1 [intact] 3.24 1.24 0.30 0.07 0.02 304B 183 1 [intact] 2.62 0.87 0.99 0.33 0.37 305B 288 1 [intact] 83.36 82.84 34.67 29.59 28.52 313A 332 1 [intact] 6.40 5.47 0.55 0.18 0.34 318A 366 1 [intact] 2.32 3.95 4.90 2.28 2.09 321B 250 1 [intact] 3.07 4.07 5.11 1.93 2.60 321A 450 1 [intact] 55.28 42.17 56.59 27.44 22.12 last fifth below 10% of the first: 5 of 13; last fifth below half the first: 10 of 13
Result: blocked (metadata insufficient). Every trial of the 13 progressive cutting recordings is labelled intact and each file has one moment of inertia, so the cutting conditions cannot be told apart. The whisking-in-air rate falls over most of these sessions: the last fifth is below half the first in 10 of 13 and below 10% of it in 5 of 13, which fits the files holding the cut conditions under the intact label.
What the archive holds for the other blocked claims¶
The archive's own description of each asset (variableMeasured), and every group and dataset in
each file (per-trial series collapsed to one entry), checked for the data the blocked claims need:
EMG, optogenetic stimulation, whisker length and manual contact labels.
In [28]:
vm = merkel.asset_variables()
print("variableMeasured across assets:", pd.Series([tuple(v) for v in vm.values()]).value_counts().to_dict())
for x in ("7A", "8A", "9A", "11A"):
r = by_rec[x]
print(f"{x} ({r['identifier']}): variableMeasured {vm[r['path']]}; neurodata types {r['neurodata_types']}")
print(" paths outside processing/behavior, general and specifications:",
[p for p in r["paths"] if not p.startswith(("processing/behavior/", "general", "specifications"))])
print(" series in processing/behavior:", r["series_present"])
other = sorted({p for r in recs for p in r["paths"] if not p.startswith(("processing/behavior/", "specifications"))})
print("\nPaths outside processing/behavior and specifications in any file:", other)
print("Members of acquisition, analysis or stimulus in any file:",
sorted({p for r in recs for p in r["paths"] if p.startswith(("acquisition/", "analysis/", "stimulus/presentation/", "stimulus/templates/"))}) or "none")
print("Neurodata types in any file:", sorted({t for r in recs for t in r["neurodata_types"]}))
print("Series in processing/behavior in any file:", sorted({x for r in recs for x in r["series_present"]}))
words = ("emg", "optogen", "stim", "laser", "light", "length", "curat", "manual", "label")
print("Paths or series names containing", words, ":",
sorted({p for r in recs for p in r["paths"] + r["series_present"] if any(w in p.lower() for w in words)}) or "none")
print("Per-trial size (MB): EMG files", {x: round(by_rec[x]["bytes"] / by_rec[x]["n_trials"] / 1e6, 3) for x in ("7A", "8A", "9A", "11A")},
"; progressive files, median", round(float(np.median([r["bytes"] / r["n_trials"] / 1e6 for r in prog])), 3))
variableMeasured across assets: {('BehavioralTimeSeries', 'ProcessingModule', 'Units'): 57, ('BehavioralTimeSeries', 'ProcessingModule'): 3}
7A (KS0007A): variableMeasured ['BehavioralTimeSeries', 'ProcessingModule']; neurodata types ['BehavioralTimeSeries', 'ElementIdentifiers', 'ProcessingModule', 'Subject', 'TimeIntervals', 'TimeSeries', 'VectorData']
paths outside processing/behavior, general and specifications: ['acquisition', 'analysis', 'file_create_date', 'identifier', 'intervals', 'intervals/trials', 'intervals/trials/MOI', 'intervals/trials/barAzim', 'intervals/trials/barCoor', 'intervals/trials/barRadii', 'intervals/trials/cuttingGroup', 'intervals/trials/id', 'intervals/trials/start_time', 'intervals/trials/stop_time', 'processing', 'processing/behavior', 'session_description', 'session_start_time', 'stimulus', 'stimulus/presentation', 'stimulus/templates', 'timestamps_reference_time']
series in processing/behavior: ['Mi', 'acceleration', 'amplitude', 'contact', 'dMi', 'deltaKappa', 'follicleCoordsX', 'follicleCoordsY', 'frequency', 'jerk', 'phase', 'position', 'setpoint', 'velocity']
8A (KS0008A): variableMeasured ['BehavioralTimeSeries', 'ProcessingModule']; neurodata types ['BehavioralTimeSeries', 'ElementIdentifiers', 'ProcessingModule', 'Subject', 'TimeIntervals', 'TimeSeries', 'VectorData']
paths outside processing/behavior, general and specifications: ['acquisition', 'analysis', 'file_create_date', 'identifier', 'intervals', 'intervals/trials', 'intervals/trials/MOI', 'intervals/trials/barAzim', 'intervals/trials/barCoor', 'intervals/trials/barRadii', 'intervals/trials/cuttingGroup', 'intervals/trials/id', 'intervals/trials/start_time', 'intervals/trials/stop_time', 'processing', 'processing/behavior', 'session_description', 'session_start_time', 'stimulus', 'stimulus/presentation', 'stimulus/templates', 'timestamps_reference_time']
series in processing/behavior: ['Mi', 'acceleration', 'amplitude', 'contact', 'dMi', 'deltaKappa', 'follicleCoordsX', 'follicleCoordsY', 'frequency', 'jerk', 'phase', 'position', 'setpoint', 'velocity']
9A (KS0009A): variableMeasured ['BehavioralTimeSeries', 'ProcessingModule']; neurodata types ['BehavioralTimeSeries', 'ElementIdentifiers', 'ProcessingModule', 'Subject', 'TimeIntervals', 'TimeSeries', 'VectorData']
paths outside processing/behavior, general and specifications: ['acquisition', 'analysis', 'file_create_date', 'identifier', 'intervals', 'intervals/trials', 'intervals/trials/MOI', 'intervals/trials/barAzim', 'intervals/trials/barCoor', 'intervals/trials/barRadii', 'intervals/trials/cuttingGroup', 'intervals/trials/id', 'intervals/trials/start_time', 'intervals/trials/stop_time', 'processing', 'processing/behavior', 'session_description', 'session_start_time', 'stimulus', 'stimulus/presentation', 'stimulus/templates', 'timestamps_reference_time']
series in processing/behavior: ['Mi', 'acceleration', 'amplitude', 'contact', 'dMi', 'deltaKappa', 'follicleCoordsX', 'follicleCoordsY', 'frequency', 'jerk', 'phase', 'position', 'setpoint', 'velocity']
11A (KS0011A): variableMeasured ['BehavioralTimeSeries', 'ProcessingModule', 'Units']; neurodata types ['BehavioralTimeSeries', 'ElementIdentifiers', 'ProcessingModule', 'Subject', 'TimeIntervals', 'TimeSeries', 'Units', 'VectorData', 'VectorIndex']
paths outside processing/behavior, general and specifications: ['acquisition', 'analysis', 'file_create_date', 'identifier', 'intervals', 'intervals/trials', 'intervals/trials/MOI', 'intervals/trials/barAzim', 'intervals/trials/barCoor', 'intervals/trials/barRadii', 'intervals/trials/cuttingGroup', 'intervals/trials/id', 'intervals/trials/start_time', 'intervals/trials/stop_time', 'processing', 'processing/behavior', 'session_description', 'session_start_time', 'stimulus', 'stimulus/presentation', 'stimulus/templates', 'timestamps_reference_time', 'units', 'units/id', 'units/obs_intervals', 'units/obs_intervals_index', 'units/spike_times', 'units/spike_times_index']
series in processing/behavior: ['Mi', 'acceleration', 'amplitude', 'contact', 'dMi', 'deltaKappa', 'follicleCoordsX', 'follicleCoordsY', 'frequency', 'jerk', 'phase', 'position', 'setpoint', 'velocity']
Paths outside processing/behavior and specifications in any file: ['acquisition', 'analysis', 'file_create_date', 'general', 'general/experimenter', 'general/institution', 'general/lab', 'general/related_publications', 'general/subject', 'general/subject/date_of_birth', 'general/subject/description', 'general/subject/sex', 'general/subject/species', 'general/subject/strain', 'general/subject/subject_id', 'identifier', 'intervals', 'intervals/trials', 'intervals/trials/MOI', 'intervals/trials/barAzim', 'intervals/trials/barCoor', 'intervals/trials/barRadii', 'intervals/trials/cuttingGroup', 'intervals/trials/id', 'intervals/trials/start_time', 'intervals/trials/stop_time', 'processing', 'processing/behavior', 'session_description', 'session_start_time', 'stimulus', 'stimulus/presentation', 'stimulus/templates', 'timestamps_reference_time', 'units', 'units/id', 'units/obs_intervals', 'units/obs_intervals_index', 'units/spike_times', 'units/spike_times_index']
Members of acquisition, analysis or stimulus in any file: none
Neurodata types in any file: ['BehavioralTimeSeries', 'ElementIdentifiers', 'ProcessingModule', 'Subject', 'TimeIntervals', 'TimeSeries', 'Units', 'VectorData', 'VectorIndex']
Series in processing/behavior in any file: ['Faxial', 'Flateral', 'Fnorm', 'M0', 'Mi', 'acceleration', 'amplitude', 'contact', 'dFaxial', 'dFlateral', 'dM0', 'dMi', 'deltaKappa', 'distanceToPoleCenter', 'follicleCoordsX', 'follicleCoordsY', 'frequency', 'jerk', 'phase', 'position', 'setpoint', 'thetaAtContact', 'velocity']
Paths or series names containing ('emg', 'optogen', 'stim', 'laser', 'light', 'length', 'curat', 'manual', 'label') : ['stimulus', 'stimulus/presentation', 'stimulus/templates']
Per-trial size (MB): EMG files {'7A': 0.309, '8A': 0.31, '9A': 0.311, '11A': 0.311} ; progressive files, median 0.311
In [29]:
# The EMG files' undescribed Mi and dMi series, by phase (30 equal-count bins of whisking frames),
# against positive and negative velocity, to see whether either holds the intrinsic protractor EMG.
rows = []
for x in ("7A", "8A", "9A", "11A"):
e = by_rec[x].get("extra_by_phase")
if e is None:
rows.append({"rec": x, "note": "no Mi/dMi or too few whisking frames"})
continue
row = {"rec": x}
for v in merkel.EXTRA_VARS:
row[f"{v}_finite"] = e[v + "_finite"]
row[f"{v}_min"], row[f"{v}_max"] = e[v + "_range"]
row[f"{v}_curve_sd"] = float(np.std(e[v]))
for d in ("pos", "neg"):
row[f"r_{v}_vel_{d}"] = np.corrcoef(e[v], e["vel_" + d])[0, 1] if np.std(e[v]) > 0 else np.nan
rows.append(row)
print("stored units:", {v: sorted({r["units"].get(v) for r in recs if v in r["units"]}) for v in merkel.EXTRA_VARS})
print(pd.DataFrame(rows).to_string(index=False, float_format=lambda v: f"{v:.3g}"))
allr = pd.DataFrame([{"rec": r["rec"], "group": table.loc[r["rec"], "group"],
**{f"r_{v}_vel_pos": np.corrcoef(r["extra_by_phase"][v], r["extra_by_phase"]["vel_pos"])[0, 1]
for v in merkel.EXTRA_VARS}}
for r in recs if "extra_by_phase" in r and min(np.std(r["extra_by_phase"][v]) for v in merkel.EXTRA_VARS) > 0])
print("\nFor comparison, the same correlations in every file with whisking, median by group:")
print(allr.groupby("group")[["r_Mi_vel_pos", "r_dMi_vel_pos"]].agg(["median", "count"]).round(3).to_string())
stored units: {'Mi': ['N*m'], 'dMi': ['N*m/s']}
rec Mi_finite Mi_min Mi_max Mi_curve_sd r_Mi_vel_pos r_Mi_vel_neg dMi_finite dMi_min dMi_max dMi_curve_sd r_dMi_vel_pos r_dMi_vel_neg
7A 1 -5.73e-09 5.35e-09 7.64e-10 0.00247 0.00234 1 -8.39e-10 9.98e-10 7.85e-11 -0.543 0.948
8A 1 -4.96e-09 3.99e-09 9.59e-10 -0.0257 -0.017 1 -6.11e-10 8.01e-10 1.12e-10 -0.57 0.911
9A 1 -8.76e-09 6.94e-09 1.13e-09 -0.0167 -0.017 1 -9.58e-10 1.34e-09 1.37e-10 -0.534 0.961
11A 1 -2.77e-08 2.32e-08 1.46e-09 0.000219 -0.00375 1 -4.74e-09 4.92e-09 1.62e-10 -0.379 0.866
For comparison, the same correlations in every file with whisking, median by group:
r_Mi_vel_pos r_dMi_vel_pos
median count median count
group
cutonly -0.025 5 -0.518 5
emg -0.008 4 -0.538 4
excluded -0.014 3 -0.485 3
progressive -0.026 13 -0.536 13
touch -0.015 33 -0.535 33
whiskonly -0.016 1 -0.523 1
Result: blocked (data missing) for figs4-intrinsic-protractor, fig7d-cut-emg, fig6a-final-length, figs1-optotag-latency, figs1-labelling-specificity and contact-classifier. The archive lists BehavioralTimeSeries and ProcessingModule for every asset and Units for 57; the intrinsic protractor files 7A, 8A and 9A have no Units, and their per-trial size (0.31 MB) equals that of 11A and of the progressive cutting files. No file has a member in acquisition, analysis or stimulus, or a series or path named for EMG, optogenetic stimulation, light, whisker length or manual labels. The trial tables hold MOI, cuttingGroup and the bar position only. The EMG files hold the tracked whisker series that the afferent files hold, including Mi (unit Nm) and dMi (Nm/s), with values below 3e-8 in absolute size. Their phase curves do not follow positive velocity: r with positive velocity is -0.026 to 0.002 for Mi and -0.570 to -0.534 for dMi in 7A, 8A and 9A, close to the medians in every group of afferent files (-0.026 to -0.014 for Mi, -0.536 to -0.485 for dMi), so neither series is the intrinsic protractor EMG. fig7d-cut-emg also needs the cut condition (see fig6-cut-eliminated).
fig7f-precut-emg (blocked): the nasolabialis curve alone¶
In 5 afferents recorded with the whisker already cut near its base, phase tuning was strongly correlated with intrinsic protractor or nasolabialis EMG (r = 0.76 plus or minus 0.06).
The best match is chosen between the intrinsic protractor and nasolabialis curves, and the intrinsic protractor EMG is not in the dandiset, so the claim stays blocked. The correlation with the nasolabialis curve alone (11A) is a lower bound on each afferent's best match. Both curves are on whisking frames in 30 equal-count phase bins; the 11A curve is interpolated (circularly) to each afferent's bin centres.
In [30]:
w11 = by_rec["11A"]["whisk_phase_curves"]
o11 = np.argsort(w11["centres"])
rows = []
for r in group("cutonly"):
w = r["whisk_phase_curves"]
naso = np.interp(w["centres"], w11["centres"][o11], w11["rate"][o11], period=2 * np.pi)
rows.append({"rec": r["rec"], "whisking_min": r["category_seconds"]["whisking"] / 60, "spikes": r["n_spikes"],
"rate_min": w["rate"].min(), "rate_max": w["rate"].max(),
"r_nasolabialis": np.corrcoef(w["rate"], naso)[0, 1],
"r_own_neg_velocity": np.corrcoef(w["rate"], w["vel_neg"])[0, 1],
"r_own_pos_velocity": np.corrcoef(w["rate"], w["vel_pos"])[0, 1],
"contact_flag_frames": r["contact_flag_frames"], "has_dpole": "distanceToPoleCenter" in r["series_present"],
"contact_min": r["category_seconds"]["contact"] / 60})
f7 = pd.DataFrame(rows)
print(f7.round(3).to_string(index=False))
print(f"\nmean r with the nasolabialis curve {f7.r_nasolabialis.mean():.3f} (SD {f7.r_nasolabialis.std():.3f}); "
f"paper's mean best-match r 0.76 plus or minus 0.06, which is at least this")
f7["r_best_own_velocity"] = f7[["r_own_neg_velocity", "r_own_pos_velocity"]].max(axis=1)
print(f"Not the criterion: the better of each afferent's correlation with its own positive and negative velocity curves "
f"(Figure S4 relates the two muscles to these) has mean {f7.r_best_own_velocity.mean():.3f}, range "
f"{f7.r_best_own_velocity.min():.3f} to {f7.r_best_own_velocity.max():.3f}; positive velocity is the better for "
f"{(f7.r_own_pos_velocity > f7.r_own_neg_velocity).sum()} of {len(f7)}")
rec whisking_min spikes rate_min rate_max r_nasolabialis r_own_neg_velocity r_own_pos_velocity contact_flag_frames has_dpole contact_min 199B 5.758 7669 0.347 61.046 -0.625 -0.593 0.795 5014 True 0.166 206A 3.213 13013 1.245 188.706 -0.615 -0.563 0.753 0 False 0.000 209A 1.480 1460 0.000 42.230 -0.294 -0.268 0.680 32 True 0.000 277A 7.401 7676 0.405 69.518 -0.544 -0.484 0.691 0 False 0.000 305A 4.347 3590 0.115 58.661 0.695 0.848 -0.647 0 False 0.000 mean r with the nasolabialis curve -0.277 (SD 0.559); paper's mean best-match r 0.76 plus or minus 0.06, which is at least this Not the criterion: the better of each afferent's correlation with its own positive and negative velocity curves (Figure S4 relates the two muscles to these) has mean 0.753, range 0.680 to 0.848; positive velocity is the better for 4 of 5
Result: blocked (data missing). The nasolabialis curve does not bound the claim usefully: 4 of the 5 cut-whisker afferents are anticorrelated with it (mean r -0.277), and only 305A correlates (0.695). These four are better matched by their own positive velocity curve, the variable the paper relates to the intrinsic protractor, whose EMG is not in the dandiset. As a stand-in that is not the criterion, the better of each afferent's correlations with its own positive and negative velocity curves has a mean of 0.753 (range 0.680 to 0.848), near the paper's 0.76. 199B and 209A hold distance-to-pole series and some contact = 1 frames (5,014 and 32) although they were recorded with the whisker cut; at 209A's contact = 1 frames the distance to the pole is 5.6 to 6.5 (1st and 99th percentiles, data check).
figs4-nasolabialis¶
Motor unit activity of m. nasolabialis was phase modulated and correlated with negative whisker velocity (r = 0.93) and anti-correlated with positive velocity (r = -0.75) within whisk cycles.
In [31]:
emg = by_rec["11A"]
w = emg["whisk_phase_curves"]
r_neg = np.corrcoef(w["rate"], w["vel_neg"])[0, 1]
r_pos = np.corrcoef(w["rate"], w["vel_pos"])[0, 1]
print(f"recording 11A ({emg['identifier']}, {emg['subject_id']}): {emg['n_spikes']} spike times, whisking {emg['category_seconds']['whisking'] / 60:.1f} min")
print(f"motor unit rate by phase: {w['rate'].min():.1f} to {w['rate'].max():.1f} Hz, peak at phase {w['centres'][np.argmax(w['rate'])]:.2f} rad, "
f"where the mean velocity is {'negative' if w['vel_neg'][np.argmax(w['rate'])] > w['vel_pos'][np.argmax(w['rate'])] else 'positive'}")
print(f"r with negative velocity {r_neg:.3f} (paper 0.93, tolerance 0.093); with positive velocity {r_pos:.3f} (paper -0.75, tolerance 0.075)")
print("within tolerance:", abs(r_neg - 0.93) <= 0.093, abs(r_pos + 0.75) <= 0.075)
print("Other EMG files (7A, 8A, 9A) have spike times:", [by_rec[x]["has_spikes"] for x in ("7A", "8A", "9A")],
"; acquisition groups:", [by_rec[x]["acquisition"] for x in ("7A", "8A", "9A")])
fig, ax = plt.subplots(figsize=(4, 2.5))
ax.plot(w["centres"], w["rate"] / w["rate"].max(), label="motor units")
ax.plot(w["centres"], w["vel_neg"] / w["vel_neg"].max(), label="-velocity")
ax.plot(w["centres"], w["vel_pos"] / w["vel_pos"].max(), label="+velocity"); ax.legend(fontsize=7); ax.set_xlabel("phase (rad)")
plt.show()
recording 11A (KS0011A, KSm11): 13851 spike times, whisking 23.3 min motor unit rate by phase: 0.7 to 37.6 Hz, peak at phase -0.91 rad, where the mean velocity is negative r with negative velocity 0.924 (paper 0.93, tolerance 0.093); with positive velocity -0.716 (paper -0.75, tolerance 0.075) within tolerance: True True Other EMG files (7A, 8A, 9A) have spike times: [False, False, False] ; acquisition groups: [[], [], []]
Result: reproduced. Across 30 phase bins the motor unit spikes of recording 11A follow negative velocity (r = 0.924; paper 0.93) and anti-correlate with positive velocity (r = -0.716; paper -0.75), both within tolerance. The other three EMG files hold no EMG signal and no spike times.
sigma-kappa¶
The noise level of curvature estimates was 7.8 plus or minus 5.0% (mean and SD) across 33 sessions.
In [32]:
sk = pd.DataFrame([{"rec": r["rec"], "a": r["sigma_kappa"]["a"], "b": r["sigma_kappa"]["b"]} for r in touch])
print(sk.round(2).to_string(index=False))
for k, lab in (("a", "(a) contact = 1 and distance < 1 mm"), ("b", "(b) contact = 1")):
print(f"{lab}: mean {sk[k].mean():.2f}%, SD {sk[k].std():.2f}, range {sk[k].min():.2f} to {sk[k].max():.2f} (paper 7.8 and 5.0)")
rec a b 217F 10.320000 10.320000 227A 8.570000 8.570000 145B 13.840000 13.820000 194A 6.230000 6.230000 202A 9.640000 9.640000 163B 11.310000 11.310000 125B 11.520000 11.510000 221B 8.020000 8.020000 225D 6.360000 6.360000 132A 8.630000 8.630000 256B 5.830000 5.830000 223B 6.660000 6.660000 161A 15.000000 15.000000 167A 6.180000 6.170000 128D 9.450000 9.430000 221C 8.670000 8.670000 198B 4.650000 4.650000 216D 9.000000 9.000000 211A 5.060000 5.060000 220B 6.070000 6.060000 286A 6.600000 6.600000 250A 6.490000 6.490000 251B 10.570000 10.570000 254C 15.710000 15.710000 207A 8.280000 8.280000 154B 9.090000 9.090000 265A 28.490000 28.490000 234A 8.330000 8.330000 231A 6.080000 6.080000 243B 39.889999 39.889999 266B 9.200000 9.200000 220C 6.300000 6.300000 282A 4.230000 4.230000 (a) contact = 1 and distance < 1 mm: mean 10.01%, SD 6.97, range 4.23 to 39.89 (paper 7.8 and 5.0) (b) contact = 1: mean 10.01%, SD 6.97, range 4.23 to 39.89 (paper 7.8 and 5.0)
Result: reproduced. The mean curvature noise over the 33 touch recordings is 10.01% with SD 6.97 (paper 7.8 plus or minus 5.0), within the paper's spread but above its mean; two recordings are high, 243B (39.89%) and 265A (28.49%). Both readings of contact frames give the same values.
Summary¶
| Claim | Outcome |
|---|---|
| table-s1-link | reproduced |
| video-frames | reproduced |
| merkel-minutes (outside labels) | reproduced |
| fig2c-baseline-rate | reproduced |
| fig2c-touch-rate | reproduced |
| fig2c-wt-star | reproduced (reading b) |
| fig3a-direction-selective (outside labels) | partially reproduced |
| fig3d-gam-full | reproduced |
| fig3d-gam-single | reproduced |
| fig3d-gam-best-pair | partially reproduced |
| fig4c-viscoelastic (outside labels) | reproduced |
| fig5-phase-modulation | reproduced |
| fig5e-phase-span | partially reproduced |
| figs3-phase-invariance | reproduced |
| figs6-kinematic-gam | reproduced (15 of 28 afferents) |
| figs4-nasolabialis | reproduced |
| sigma-kappa | reproduced |
| fig7g-phase-span-28, fig6-cut-eliminated, fig6-cut-remained, fig7e-lost-rate-kinematics | blocked: metadata insufficient |
| fig6a-final-length, fig7d-cut-emg, fig7f-precut-emg, figs4-intrinsic-protractor, figs1-optotag-latency, figs1-labelling-specificity, contact-classifier | blocked: data missing |
A failure to reproduce here may be a failure of this attempt, not of the paper: see the reasoning and the differences from the paper's method.
How the analysis differs from the paper's. The paper's GAMs (mgcv, 1 ms bins, UBRE smoothing) are replaced by penalized B-spline models fitted as binomial models with two trials per 2 ms frame, on at most 60,000 frames per recording; the viscoelastic model is fitted with L-BFGS-B from eight starts. Confidence intervals for the Glossary classes are bootstrapped over trials. The kinematic models of Figure S6 are tested on 15 of the paper's 28 afferents. Which recording belongs to which experiment, and the afferent types, come from Table S1. Two outcomes sit near their thresholds: the Rayleigh p of fig5e-phase-span moves across 0.05 with the preferred-phase estimator, and M0 and M0' rank second of the 45 pairs, 0.01 behind M0 and velocity, under every variant tried.
What was found about the dandiset. The 60 files match Table S1 one to one, and most values come close to the paper's frame counts, rates and model performance. The exceptions are the median non-preferred phase rate (16% above the paper's), the mean curvature noise (28% above), and direction selectivity (65% of Merkel and SA afferents at an index of 0.5 or more, against the 80% the criterion asks). But the files do not say which experiment each recording belongs to, which three touch recordings the paper excluded, the afferent type or Merkel identity, or the meaning of their series (every series is "no description"). The progressive cutting recordings label every trial intact, though their rates fall over the session. The intrinsic protractor EMG, the optogenetic tagging data and the whisker cut lengths are not in the dandiset. The derivatives are stored per millisecond but labelled per second, the stored whisk phase is the negative of the Hilbert phase, so it runs backwards in time, trial intervals are 7 s for 5 s of data, one file's spike times are unsorted, and two identifiers lack the unit letter. Two recordings made with the whisker cut (199B, 209A) hold distance-to-pole series and some contact frames, which the files do not explain.
In [33]:
print(f"total run time {(time.time() - started) / 60:.1f} min")
total run time 27.4 min