Reproducing DeWeese and Zador (2004) from CRCNS data set ac-2¶
Paper. DeWeese MR, Zador AM (2004). Shared and private variability in the auditory cortex.
J Neurophysiol 92:1840-1855. doi:10.1152/jn.00197.2004.
The claims file (claims/crcns-ac-2/10.1152-jn.00197.2004.yaml) was written from the version of
record (text extracted from the journal PDF, with its Methods and figure legends; the paper has no
supplement and links no code).
Data set. CRCNS ac-2, doi:10.6080/K0G44N6R, read from
the public bucket crcnsarchive (ac-2/crcns-ac-2.tar.gz and
ac-2/crcns-ac-2-docs_and_scripts.tar.gz). CRCNS data sets are not versioned, so the version is
recorded as unversioned; the only published content is this one, and the paper cites none. It
is pinned by the MD5 sums of the data set's checksums.md5: crcns-ac-2.tar.gz
af006ca3530e7f615a37d0548ba86cf2 and crcns-ac-2-docs_and_scripts.tar.gz
6ffd46575e78b3f8a2ad3efb9f95008a. The notebook computes both sums from the bytes it streams and
compares them (below), so a later change of the files is detected. The data tar is a gzipped tar,
which cannot be read by range; it is streamed once into memory (about 30 s) and nothing is
written to disk.
What the data set holds. The 33 anesthetized whole-cell recordings of the variability analysis, 29 of them with a simultaneous LFP; six recordings after TTX (the paper's mEPSPs); one awake recording that belongs to a later paper and is not used.
| Claim | Tested | Outcome |
|---|---|---|
| fig2-stereotyped-shape | yes | partially reproduced |
| fig3a-reliable-neuron | yes | reproduced |
| fig3b-total-variability-distribution | yes | reproduced (assumed) |
| total-variability-population | yes | reproduced (assumed) |
| lfp-culling-reduction | yes | reproduced (assumed) |
| fig5b-residual-lfp-psp-correlation | yes | reproduced |
| fig7a-lfp-culling-estimation | yes | partially reproduced (assumed) |
| psp-culling-reduction | yes | partially reproduced (assumed) |
| fig4c-psp-culling-can-increase | yes | partially reproduced (assumed) |
| fig7b-private-best-estimate | yes | reproduced (assumed) |
| culling-controls | yes | not reproduced |
| fig7d-large-responses | yes | reproduced (assumed) |
| large-responses-controls | yes | not reproduced |
| cv-large-responses | yes | partially reproduced (assumed) |
| lfp-cv-not-different | yes | reproduced |
| anesthesia-comparison | yes | reproduced (assumed) |
| fig8-spontaneous-veto | yes | reproduced |
| overfitting-correction-example | yes | reproduced (assumed) |
| private-shared-ranges | yes | partially reproduced |
| quantal-comparison | yes | reproduced |
| mepsp-size | yes | partially reproduced (assumed) |
| recording-properties | yes | reproduced (assumed) |
| access-resistance-depth-no-effect | yes | reproduced (assumed) |
"Assumed" marks outcomes that rest on a choice the paper does not make and that would differ under another reading: the spike rule for most, a reading of the open culling details or the measuring of heights once before culling for some, the input resistance for recording-properties and the mEPSP detector for mepsp-size.
import time
import warnings
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from scipy import stats
import ac2
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()
D = ac2.load()
N = D["neurons"]
print(f"data set {ac2.DATASET} ({ac2.DATASET_VERSION}, doi:{ac2.DATASET_DOI}) loaded in {time.time() - started:.0f} s")
print("MD5 of the streamed files:", D["md5"])
print("MD5 in checksums.md5: ", ac2.MD5)
print("Same:", D["md5"] == ac2.MD5)
data set ac-2 (unversioned, doi:10.6080/K0G44N6R) loaded in 33 s
MD5 of the streamed files: {'crcns-ac-2-docs_and_scripts.tar.gz': '6ffd46575e78b3f8a2ad3efb9f95008a', 'crcns-ac-2.tar.gz': 'af006ca3530e7f615a37d0548ba86cf2'}
MD5 in checksums.md5: {'crcns-ac-2.tar.gz': 'af006ca3530e7f615a37d0548ba86cf2', 'crcns-ac-2-docs_and_scripts.tar.gz': '6ffd46575e78b3f8a2ad3efb9f95008a'}
Same: True
Data check¶
Each check of the documentation against the values is stated under its output.
info = pd.DataFrame(D["info"]).T
print("Records in neuron_recording_info.mat:", len(info))
print("Anesthetized neurons with traces:", len(N), "; TTX recordings:", len(D["ttx"]),
"; files of the awake recording (not used):", len(D["awake_files"]))
print("Other files in the data tar:", D["other_files"])
print("Tone order read from Anesthetized_cell_info.txt (Hz):", ac2.TONE_ORDER.tolist())
steps = np.diff(np.log(np.sort(ac2.TONE_ORDER)))
print(f"{len(ac2.TONE_ORDER)} tones, {ac2.TONE_ORDER.min()} to {ac2.TONE_ORDER.max()} Hz; log steps {steps.min():.4f} to {steps.max():.4f}")
rows = []
for name, n in N.items():
st = ac2.stable(n)
p = n["params"]
rows.append({
"neuron": name, "anesthetic": n["info"]["which_anesthetic"], "LFP": int(n["info"]["was_LFP_recorded"]),
"lfp_gain": n["info"]["axogaincorrection_LFP"], "sweeps": n["n_sweeps"], "traces": len(n["r"]),
"cycles": int(len(n["r"]) // 32), "first_cycle": int(p["first_cycle"]), "last_cycle": int(p["last_cycle"]),
"reps": int(p["last_cycle"] - p["first_cycle"] + 1), "stable_traces": int(st.sum()),
"spikethresh": p["spikethresh"], "spikes_cell": int((ac2.spikes(n, "cell") & st).sum()),
"spikes_claims_rule": int((ac2.spikes(n, "claims") & st).sum()),
"R_access": n["info"]["R_access"], "depth_um": n["info"]["Sutter_diag_depth"] / np.sqrt(2),
"max_baseline_offset": float(np.abs(n["r"][:, :ac2.PRE].mean(1)).max()),
"trace_vs_sweep_maxdiff": float(n["trace_match_error"].max()),
"onset_shift_samples": int(np.abs(n["onset_offset"]).max()),
})
cells = pd.DataFrame(rows).set_index("neuron")
print(cells.round(3).to_string())
print("\nNeurons by anesthetic:", cells.anesthetic.value_counts().to_dict(), "(paper: 13 P, 3 P with diazepam, 17 urethane)")
print("Neurons with the LFP flag:", int(cells.LFP.sum()), "(paper: 29 with a stable LFP)")
print("Traces = 16 x sweeps in every neuron:", bool((cells.traces == 16 * cells.sweeps).all()))
print("Largest |mean of the first 15 ms| of any trace (should be 0):", cells.max_baseline_offset.max())
print("Largest difference between a trace and its place in the sweep (mV):", cells.trace_vs_sweep_maxdiff.max(),
"; largest onset shift found (samples):", cells.onset_shift_samples.max())
Records in neuron_recording_info.mat: 40
Anesthetized neurons with traces: 33 ; TTX recordings: 6 ; files of the awake recording (not used): 88
Other files in the data tar: ['crcns/ac-2/Readme.txt']
Tone order read from Anesthetized_cell_info.txt (Hz): [4074, 38134, 10172, 7498, 16909, 3681, 31118, 2713, 15275, 3325, 20722, 2451, 28110, 18719, 4993, 42214, 6773, 13799, 46731, 11260, 2000, 9188, 2214, 8300, 12465, 25393, 3003, 4510, 6119, 34448, 22939, 5527]
32 tones, 2000 to 46731 Hz; log steps 0.1014 to 0.1019
anesthetic LFP lfp_gain sweeps traces cycles first_cycle last_cycle reps stable_traces spikethresh spikes_cell spikes_claims_rule R_access depth_um max_baseline_offset trace_vs_sweep_maxdiff onset_shift_samples
neuron
020701md01b P 1 2.0 152 2432 76 8 70 63 2016 50.0 0 0 15.0 767.211 0.0 0.0 0
020901md01 P_with_Diaz 1 10.0 89 1424 44 3 36 34 1088 50.0 0 0 44.0 319.612 0.0 0.0 0
021001md01 P_with_Diaz 1 10.0 197 3152 98 3 15 13 416 40.0 1 1 105.0 NaN 0.0 0.0 0
021401md01b P 1 10.0 116 1856 58 30 55 26 832 45.0 0 0 37.0 445.477 0.0 0.0 0
021401md02 P 1 10.0 24 384 12 1 11 11 352 50.0 0 0 23.0 381.131 0.0 0.0 0
021401md03 P 1 10.0 22 352 11 1 10 10 320 50.0 1 1 200.0 NaN 0.0 0.0 0
021401md05 P 1 10.0 50 800 25 4 17 14 448 44.0 3 3 60.0 833.679 0.0 0.0 0
021901md01 P 1 10.0 54 864 27 1 27 27 864 29.0 9 5 110.0 424.264 0.0 0.0 0
022301md01 P_with_Diaz 1 10.0 61 976 30 3 30 28 896 30.0 32 10 320.0 517.602 0.0 0.0 0
022801md01b P 1 10.0 46 736 23 3 22 20 640 40.0 0 0 20.0 547.301 0.0 0.0 0
022801md02b P 1 10.0 51 816 25 2 12 11 352 50.0 0 0 10.0 607.405 0.0 0.0 0
030201md01b P 1 10.0 54 864 27 5 24 20 640 50.0 34 34 120.0 577.706 0.0 0.0 0
030201md03b P 1 10.0 82 1312 41 1 17 17 544 50.0 0 0 40.0 675.994 0.0 0.0 0
030301md01b P 1 10.0 40 640 20 7 20 14 448 55.0 0 0 32.0 179.605 0.0 0.0 0
030301md02b P 1 10.0 30 480 15 2 13 12 384 50.0 0 0 32.0 465.276 0.0 0.0 0
030301md03b P 1 10.0 44 704 22 4 19 16 512 35.0 4 3 47.0 762.968 0.0 0.0 0
100300md01 U 0 NaN 67 1072 33 10 25 16 512 45.0 0 0 54.0 388.909 0.0 0.0 0
100300md02 U 0 NaN 29 464 14 2 11 10 320 40.0 19 0 17.0 180.312 0.0 0.0 0
100500md01 U 0 NaN 38 608 19 3 13 11 352 40.0 15 5 160.0 474.469 0.0 0.0 0
100500md02b U 0 NaN 30 480 15 1 14 14 448 45.0 5 0 113.0 NaN 0.0 0.0 0
102400md01 U 1 2.0 84 1344 42 8 24 17 544 45.0 5 3 98.0 530.330 0.0 0.0 0
102500md02 U 1 2.0 162 2592 81 5 48 44 1408 45.0 34 10 128.0 NaN 0.0 0.0 0
102600md01 U 1 2.0 44 704 22 8 19 12 384 45.0 0 0 28.0 NaN 0.0 0.0 0
102600md03 U 1 2.0 80 1280 40 16 27 12 384 40.0 2 0 32.0 NaN 0.0 0.0 0
102800md01 U 1 2.0 60 960 30 9 22 14 448 45.0 0 0 223.0 452.548 0.0 0.0 0
102800md02 U 1 2.0 100 1600 50 20 35 16 512 35.0 31 0 20.0 NaN 0.0 0.0 0
102900md02 U 1 2.0 61 976 30 9 25 17 544 45.0 23 13 159.0 516.895 0.0 0.0 0
102900md04 U 1 2.0 29 464 14 2 12 11 352 32.0 29 2 161.0 NaN 0.0 0.0 0
110700md01 U 1 2.0 82 1312 41 4 15 12 384 28.0 13 1 214.0 258.801 0.0 0.0 0
110900md05 U 1 2.0 86 1376 43 10 22 13 416 40.0 5 0 30.0 622.254 0.0 0.0 0
111200md01 U 1 2.0 104 1664 52 8 28 21 672 45.0 19 17 145.0 464.569 0.0 0.0 0
111400md01 U 1 2.0 45 720 22 10 21 12 384 38.0 37 37 239.0 361.332 0.0 0.0 0
111500md01 U 1 2.0 37 592 18 1 10 10 320 45.0 21 17 110.0 418.607 0.0 0.0 0
Neurons by anesthetic: {'U': 17, 'P': 13, 'P_with_Diaz': 3} (paper: 13 P, 3 P with diazepam, 17 urethane)
Neurons with the LFP flag: 29 (paper: 29 with a stable LFP)
Traces = 16 x sweeps in every neuron: True
Largest |mean of the first 15 ms| of any trace (should be 0): 1.3605753338197246e-05
Largest difference between a trace and its place in the sweep (mV): 1.0172526060614473e-05 ; largest onset shift found (samples): 0
The data set's files agree with its documentation on the counts: 33 anesthetized neurons, 13
under pentobarbital, 3 with diazepam and 17 under urethane (the paper's group sizes), and 29 with
the flag was_LFP_recorded (the paper's 29 with a stable LFP). Each trace is 600 samples (150 ms
at the documented 4 kHz) and starts at its nominal place in the sweeps (200 ms + 500 ms per tone,
minus 15 ms): every trace equals its piece of the sweep minus its stored baseline, to floating
point, and no shift is needed. The documentation says the stimulus-locked traces were corrected
for a mismatch of the tone onset times; the traces are not shifted against the stored sweeps, so
either the sweeps carry the same correction or none was applied; this cannot be told from the
files. Every trace is relative to the mean of its first 15 ms, as documented. Each neuron's
cell_specific_parameters.m gives a stable epoch (first_cycle to last_cycle, cycles of 32
tones) and a spike threshold; the files do not describe them for these neurons (the TTX
neurons' copies of the file describe both fields).
reps = cells.reps
print("Repetitions per tone in the stable epoch: mean %.1f, SE %.1f, range %d to %d, n = %d"
% (reps.mean(), reps.std(ddof=1) / np.sqrt(len(reps)), reps.min(), reps.max(), len(reps)))
print("Paper: 18 +- 2 (mean +- SE), range 10 to 63, n = 33")
# Recordings deposited twice: compare the first sweep of every pair of neurons.
first = {name: n["r"][:64].astype(float) for name, n in N.items()}
dup = [(a, b) for i, a in enumerate(first) for b in list(first)[i + 1:]
if first[a].shape == first[b].shape and np.allclose(first[a], first[b])]
print("Pairs of neurons with identical first two sweeps:", dup)
print("Animals by recording date prefix:", len({c[:6] for c in N}), "(paper: 21 animals)")
print("Ages and dates are in the docs table only; the date prefix of each name is the recording day.")
Repetitions per tone in the stable epoch: mean 18.1, SE 1.9, range 10 to 63, n = 33 Paper: 18 +- 2 (mean +- SE), range 10 to 63, n = 33
Pairs of neurons with identical first two sweeps: [] Animals by recording date prefix: 21 (paper: 21 animals) Ages and dates are in the docs table only; the date prefix of each name is the recording day.
Repetitions per tone are taken as the length of each neuron's stable epoch. They match the paper: mean 18, range 10 to 63 (printed above). No two neurons share their traces. The neuron names carry 21 recording days, the paper's 21 animals; names with the same day are listed as the same animal in the docs table ("same").
# The clock of the traces against the tone sequence: the tone explains the response after onset
# (peak in the 50 ms after onset), and not the 100 ms before onset, taken from the sweeps. The
# null is the same measure with the tone labels permuted at random over the trials.
def eta2(h, labels):
df = pd.DataFrame({"h": h, "g": labels})
ss_b = (df.groupby("g").h.transform("mean") - df.h.mean()) ** 2
return ss_b.sum() / ((df.h - df.h.mean()) ** 2).sum()
rng_clock = np.random.default_rng(1)
clock = []
for name, n in N.items():
idx = np.flatnonzero(ac2.stable(n) & ~ac2.spikes(n, "cell"))
labels = n["tone"][idx]
post = n["r"][idx, ac2.PRE:ac2.PRE + 200].max(1)
pre = n["pre_wc"][idx].max(1) - n["pre_wc"][idx].mean(1)
null = [eta2(post, rng_clock.permutation(labels)) for _ in range(20)]
row = {"neuron": name, "post_onset": eta2(post, labels), "pre_onset": eta2(pre, labels),
"null_median": np.median(null), "null_max": np.max(null)}
if ac2.has_lfp(n):
lpost = np.abs(n["lfp"][idx, ac2.PRE:ac2.PRE + 200]).max(1)
lpre = np.abs(n["pre_lfp"][idx] - n["pre_lfp"][idx].mean(1, keepdims=True)).max(1)
row.update(lfp_post=eta2(lpost, labels), lfp_pre=eta2(lpre, labels),
lfp_null_max=max(eta2(lpost, rng_clock.permutation(labels)) for _ in range(20)))
row["evoked_mV"] = float(n["r"][idx, ac2.PRE + 40:ac2.PRE + 200].mean())
# Response onset: first time after onset where the mean over all tones exceeds 1 mV, and the
# largest |mean| in the 15 ms before onset and the first 5 ms after it.
grand = n["r"][idx].mean(0)
above = np.flatnonzero(grand[ac2.PRE:] > 1.0)
row["onset_latency_ms"] = above[0] / ac2.FS * 1000 if len(above) else np.nan
row["max_abs_mean_-15_to_5ms"] = float(np.abs(grand[:ac2.PRE + 20]).max())
clock.append(row)
clock = pd.DataFrame(clock).set_index("neuron")
print("Fraction of variance explained by the tone (eta squared), after and before onset, with the null:")
print(clock.round(3).to_string())
print("\nWhole cell: post-onset eta2 above the largest of 20 permuted values in",
int((clock.post_onset > clock.null_max).sum()), "of", len(clock), "neurons; pre-onset above it in",
int((clock.pre_onset > clock.null_max).sum()))
lc = clock.dropna(subset=["lfp_post"])
print("LFP: post-onset above the largest permuted value in", int((lc.lfp_post > lc.lfp_null_max).sum()), "of", len(lc),
"; pre-onset above it in", int((lc.lfp_pre > lc.lfp_null_max).sum()))
Fraction of variance explained by the tone (eta squared), after and before onset, with the null:
post_onset pre_onset null_median null_max lfp_post lfp_pre lfp_null_max evoked_mV onset_latency_ms max_abs_mean_-15_to_5ms
neuron
020701md01b 0.573 0.028 0.015 0.022 0.523 0.024 0.021 1.581 26.25 0.010
020901md01 0.229 0.118 0.029 0.043 0.253 0.061 0.041 2.608 17.75 0.199
021001md01 0.367 0.153 0.069 0.089 0.357 0.117 0.130 0.526 38.50 0.035
021401md01b 0.292 0.107 0.033 0.053 0.231 0.044 0.060 1.285 26.25 0.029
021401md02 0.275 0.105 0.085 0.137 0.315 0.071 0.141 0.425 NaN 0.026
021401md03 0.386 0.153 0.105 0.141 0.518 0.084 0.162 1.074 27.75 0.064
021401md05 0.563 0.082 0.071 0.110 0.216 0.093 0.101 5.169 16.75 0.042
021901md01 0.113 0.062 0.037 0.068 0.119 0.030 0.060 2.732 18.50 0.089
022301md01 0.322 0.217 0.040 0.069 0.515 0.071 0.050 3.430 18.00 0.172
022801md01b 0.839 0.135 0.044 0.065 0.376 0.041 0.077 3.485 18.25 0.027
022801md02b 0.769 0.105 0.090 0.118 0.547 0.052 0.120 7.721 16.75 0.038
030201md01b 0.777 0.559 0.047 0.085 0.470 0.050 0.100 1.451 16.25 0.027
030201md03b 0.561 0.124 0.055 0.122 0.684 0.071 0.078 4.281 14.25 0.055
030301md01b 0.500 0.094 0.067 0.110 0.408 0.065 0.149 2.664 20.50 0.104
030301md02b 0.837 0.179 0.073 0.112 0.101 0.039 0.108 3.678 18.25 0.096
030301md03b 0.236 0.076 0.060 0.079 0.640 0.062 0.074 1.791 22.75 0.043
100300md01 0.478 0.340 0.063 0.081 NaN NaN NaN 1.283 30.75 0.025
100300md02 0.436 0.244 0.104 0.174 NaN NaN NaN 6.845 23.00 0.113
100500md01 0.743 0.308 0.090 0.139 NaN NaN NaN 4.413 23.25 0.067
100500md02b 0.577 0.557 0.062 0.095 NaN NaN NaN 1.371 29.75 0.016
102400md01 0.339 0.104 0.054 0.085 0.090 0.055 0.083 2.765 26.50 0.158
102500md02 0.257 0.037 0.022 0.032 0.231 0.019 0.029 3.863 19.50 0.153
102600md01 0.373 0.219 0.080 0.143 0.651 0.075 0.136 1.379 29.75 0.068
102600md03 0.373 0.226 0.082 0.130 0.560 0.104 0.121 0.644 34.25 0.035
102800md01 0.372 0.274 0.072 0.109 0.304 0.145 0.110 1.688 22.50 0.053
102800md02 0.416 0.204 0.062 0.080 0.428 0.122 0.085 1.682 29.00 0.116
102900md02 0.402 0.134 0.057 0.091 0.354 0.105 0.089 3.934 22.00 0.129
102900md04 0.263 0.208 0.096 0.158 0.598 0.138 0.124 3.855 20.50 0.158
110700md01 0.388 0.119 0.085 0.135 0.620 0.150 0.148 3.390 23.50 0.093
110900md05 0.267 0.134 0.073 0.098 0.084 0.066 0.112 1.482 22.25 0.038
111200md01 0.365 0.122 0.048 0.075 0.470 0.130 0.071 7.388 18.75 0.036
111400md01 0.424 0.145 0.085 0.119 0.338 0.244 0.128 6.063 18.25 0.248
111500md01 0.710 0.169 0.095 0.142 0.447 0.247 0.180 9.964 17.00 0.086
Whole cell: post-onset eta2 above the largest of 20 permuted values in 33 of 33 neurons; pre-onset above it in 26
LFP: post-onset above the largest permuted value in 27 of 29 ; pre-onset above it in 11
The tone sequence of the documentation is aligned with the traces: in every neuron the tone
explains the response after onset beyond the largest of 20 permuted labellings, and the mean
response rises above 1 mV only after onset (onset_latency_ms, printed above), with the mean
trace flat from 15 ms before to 5 ms after onset. The 100 ms before onset is also tone-dependent
beyond the null in most neurons; this is expected with a fixed tone order, since each tone is
always preceded by the same tone 500 ms earlier, whose slow after-effects fall in that window.
So this contrast cannot separate an offset from a response, and the alignment rests on the
onset latency and the post-onset tuning. The LFP and the whole-cell traces are cut from the same
sweeps at the same samples (the sweeps file holds both channels), so they share the clock.
# Units and ranges of the variables used.
allr = np.concatenate([n["r"][ac2.stable(n)].ravel()[::50] for n in N.values()])
alll = np.concatenate([n["lfp"][ac2.stable(n)].ravel()[::50] for n in N.values() if ac2.has_lfp(n)])
print("r (mV, relative to the 15-ms baseline), 0.1 / 50 / 99.9 percentiles:", np.percentile(allr, [0.1, 50, 99.9]).round(2))
print("lfp (mV after the gain correction), 0.1 / 50 / 99.9 percentiles:", np.percentile(alll, [0.1, 50, 99.9]).round(3))
base = pd.Series({name: np.median(n["rbaseline"][ac2.stable(n)]) for name, n in N.items()})
print("Median recorded baseline per neuron (mV, uncorrected): %.1f to %.1f" % (base.min(), base.max()))
gains = info.loc[list(N), ["axogaincorrection_wholecell", "axogaincorrection_LFP"]]
print("Whole-cell gain corrections:", gains.axogaincorrection_wholecell.unique(),
"; LFP gain corrections:", pd.Series(gains.axogaincorrection_LFP.astype(float)).value_counts(dropna=False).to_dict())
lfp_param = {name: n["params"].get("Axo_gain_correction_of_lfp") for name, n in N.items()}
mism = [name for name in N if not (np.isnan(float(lfp_param[name])) and np.isnan(float(info.loc[name, "axogaincorrection_LFP"])))
and float(lfp_param[name]) != float(info.loc[name, "axogaincorrection_LFP"])]
print("Neurons whose LFP gain in cell_specific_parameters.m differs from neuron_recording_info.mat:", mism)
print("TTX recordings:", {k: (len(v["ep"]), v["params"]["first_trial"], v["params"]["last_trial"], v["gain"]) for k, v in D["ttx"].items()})
r (mV, relative to the 15-ms baseline), 0.1 / 50 / 99.9 percentiles: [-28.61 0.56 48.02]
lfp (mV after the gain correction), 0.1 / 50 / 99.9 percentiles: [-0.098 -0. 0.051]
Median recorded baseline per neuron (mV, uncorrected): -80.0 to -46.4
Whole-cell gain corrections: [1.0] ; LFP gain corrections: {10.0: 15, 2.0: 14, nan: 4}
Neurons whose LFP gain in cell_specific_parameters.m differs from neuron_recording_info.mat: []
TTX recordings: {'030401md02h': (18, 1.0, 18.0, 4.0), '050301md01': (24, 1.0, 8.0, 1.0), '050401md02': (5, 1.0, 1.0, 1.0), '050501md01': (37, 12.0, 21.0, 1.0), '050501md04': (81, 1.0, 36.0, 1.0), '050801md01c': (44, 5.0, 26.0, 1.0)}
The paper's procedure, under every reading¶
All claims on variability use the same per-tone values: PSP heights, q, and q after each stage of Figure 6 (PSP-based culling, LFP-based culling, estimation with its overfitting correction, and both cullings with estimation). They are computed here once per reading of the open details (claims file header, details 1 to 7), and each claim's section reads them. The main reading:
- trials: the stable epoch of each neuron (
first_cycletolast_cycle), without spike trials; - spikes: the data set's per-neuron
spikethreshabove the trace's baseline (an assumption; see below), with the claims file's rule, a fast rise to above -20 mV after the 12 mV junction potential correction, as a reading (6); - peak of the mean PSP searched over the whole trace after onset (detail 1a), the culling window as a reading (1b);
- peak-to-trough height over the culling window (2a), the whole trace as a reading (2b);
- variance over N - 1 (3), N as a reading;
- negative scale factors removed in PSP-based culling too (4a), not removed as a reading (4b);
- the 3 mV and 15 mV criteria on the mean before culling (5a), after it as a reading (5b);
- overfitting correction with the mean over 20 shuffles (7); per neuron and per tone give the same tone average, so one is run.
PSP-based culling comes first and LFP-based culling is then run on the trials it keeps, its mean
LFP recomputed over them, as the Figure 6 legend has it ("we repeat the above procedure after
first performing the PSP-based culling procedure"); deciding each culling on the full ensemble
and keeping the trials neither removes is a reading ("cullings intersected"). Heights are
measured once on the full ensemble (a choice stated in ac2.py); the Methods say only that q is
"recalculated" after culling, so heights taken again against the peak of the mean of the trials
each culling keeps are a reading ("heights re-measured after culling").
SPK = "spikes: fast rise above -20 mV (6)"
SPK_ANY = "spikes: above -20 mV, any rise"
SPK_SLOW = "spikes: spikethresh, slow events kept"
READINGS = {
"main": {},
SPK: {"spike_rule": "claims"},
SPK_ANY: {"spike_rule": "corrected"},
SPK_SLOW: {"spike_rule": "cell_fast"},
"peak in culling window (1b)": {"peak_window": "culling"},
"peak-to-trough over trace (2b)": {"span": "trace"},
"variance over N (3)": {"ddof": 0},
"no negative-f rule in PSP culling (4b)": {"psp_negative": False},
"cullings intersected": {"combine": "intersection"},
"heights re-measured after culling": {"remeasure": True},
}
SUMMARY_READINGS = {"criteria after culling (5b)": {"criteria": "after"}}
t0 = time.time()
PT = {label: ac2.population(N, seed=2004, **opts) for label, opts in READINGS.items()}
S = {label: ac2.summarise(PT[label], N) for label in READINGS}
for label, opts in SUMMARY_READINGS.items():
S[label] = ac2.summarise(PT["main"], N, **opts)
print(f"{len(S)} readings computed in {time.time() - t0:.0f} s")
main = S["main"]
lfp_neurons = main.index[main.lfp]
def mse(x, digits=2):
x = pd.Series(x).dropna()
return f"{x.mean():.{digits}f} +- {x.std(ddof=1) / np.sqrt(len(x)):.{digits}f} (n = {len(x)})"
overview = pd.DataFrame({label: {
"<q> 33": mse(df.q), "<q> 29 LFP": mse(df.q[df.lfp]), "<q PSP-culled>": mse(df.q_psp),
"<q LFP-culled>": mse(df.q_lfp[df.lfp]), "<q LFP-culled&est>": mse(df.q_lfp_est[df.lfp]),
"<q PSP,LFP-culled&est>": mse(df.q_both_est)} for label, df in S.items()}).T
print(overview.to_string())
11 readings computed in 31 s
<q> 33 <q> 29 LFP <q PSP-culled> <q LFP-culled> <q LFP-culled&est> <q PSP,LFP-culled&est>
main 4.58 +- 0.47 (n = 33) 4.80 +- 0.52 (n = 29) 2.19 +- 0.31 (n = 33) 3.19 +- 0.43 (n = 29) 2.69 +- 0.33 (n = 29) 1.48 +- 0.17 (n = 28)
spikes: fast rise above -20 mV (6) 5.35 +- 0.54 (n = 33) 5.41 +- 0.59 (n = 29) 2.72 +- 0.35 (n = 33) 3.75 +- 0.56 (n = 29) 3.22 +- 0.48 (n = 29) 2.01 +- 0.31 (n = 28)
spikes: above -20 mV, any rise 4.43 +- 0.43 (n = 33) 4.63 +- 0.47 (n = 29) 2.42 +- 0.32 (n = 33) 3.13 +- 0.40 (n = 29) 2.68 +- 0.31 (n = 29) 1.56 +- 0.17 (n = 28)
spikes: spikethresh, slow events kept 5.11 +- 0.51 (n = 33) 5.22 +- 0.56 (n = 29) 2.54 +- 0.34 (n = 33) 3.41 +- 0.48 (n = 29) 2.85 +- 0.37 (n = 29) 1.72 +- 0.24 (n = 28)
peak in culling window (1b) 4.22 +- 0.44 (n = 33) 4.43 +- 0.48 (n = 29) 2.09 +- 0.30 (n = 33) 2.99 +- 0.41 (n = 29) 2.53 +- 0.31 (n = 29) 1.42 +- 0.16 (n = 28)
peak-to-trough over trace (2b) 4.58 +- 0.47 (n = 33) 4.80 +- 0.52 (n = 29) 2.22 +- 0.31 (n = 33) 3.12 +- 0.36 (n = 29) 2.66 +- 0.30 (n = 29) 1.57 +- 0.17 (n = 28)
variance over N (3) 4.29 +- 0.45 (n = 33) 4.50 +- 0.50 (n = 29) 2.02 +- 0.28 (n = 33) 2.94 +- 0.40 (n = 29) 2.48 +- 0.30 (n = 29) 1.36 +- 0.15 (n = 28)
no negative-f rule in PSP culling (4b) 4.58 +- 0.47 (n = 33) 4.80 +- 0.52 (n = 29) 2.69 +- 0.26 (n = 33) 3.19 +- 0.43 (n = 29) 2.67 +- 0.32 (n = 29) 1.78 +- 0.22 (n = 28)
cullings intersected 4.58 +- 0.47 (n = 33) 4.80 +- 0.52 (n = 29) 2.19 +- 0.31 (n = 33) 3.19 +- 0.43 (n = 29) 2.68 +- 0.32 (n = 29) 1.51 +- 0.21 (n = 28)
heights re-measured after culling 4.58 +- 0.47 (n = 33) 4.80 +- 0.52 (n = 29) 2.18 +- 0.29 (n = 33) 3.07 +- 0.37 (n = 29) 2.62 +- 0.30 (n = 29) 1.51 +- 0.17 (n = 28)
criteria after culling (5b) 4.58 +- 0.47 (n = 33) 4.80 +- 0.52 (n = 29) 1.92 +- 0.13 (n = 32) 3.04 +- 0.35 (n = 29) 2.59 +- 0.29 (n = 29) 1.47 +- 0.14 (n = 28)
Spike trials. The claims file's header treats spike detection as a standard choice, "a fast rise
to above -20 mV", that "can change only the few trials with spikes" (detail 6). That rule is run
here as written: the trace, corrected by the paper's 12 mV junction potential, exceeds -20 mV,
and its largest rise between samples exceeds 10 mV/ms (a standard choice for the "fast rise").
The data set gives each neuron its own threshold (spikethresh in cell_specific_parameters.m,
28 to 55 mV above the trace's baseline), which the TTX neurons' copies of the same file describe
as the "height of threshold (mV) above BASELINE/REST for spike detection/removal". The table
below gives, per neuron, the trials each rule removes and their largest rise, and the figure
shows trials that only one rule removes.
# Trials each spike rule removes, per neuron (stable epoch), with the largest rise between samples.
rows = []
for name, n in N.items():
st = ac2.stable(n)
rise = ac2.max_rise(n["r"])
recorded = (n["r"] + n["rbaseline"][:, None]).max(1)
masks = {"above -20 mV uncorrected": (recorded > -20) & st,
"above -20 mV corrected": (recorded - ac2.JUNCTION_MV > -20) & st,
"fast rise and above -20 mV corrected (claims file)": ac2.spikes(n, "claims") & st,
"spikethresh (data set)": ac2.spikes(n, "cell") & st}
row = {"neuron": name, "spikethresh": n["params"]["spikethresh"]}
for label, m in masks.items():
row[label] = int(m.sum())
row[label + ", max rise mV/ms"] = f"{rise[m].min():.1f} to {rise[m].max():.1f}" if m.any() else ""
row["spikethresh only"] = int((masks["spikethresh (data set)"] & ~masks["fast rise and above -20 mV corrected (claims file)"]).sum())
row["claims file only"] = int((~masks["spikethresh (data set)"] & masks["fast rise and above -20 mV corrected (claims file)"]).sum())
row["max rise, 99th percentile of all trials"] = float(np.percentile(rise[st], 99))
rows.append(row)
spk = pd.DataFrame(rows).set_index("neuron")
print(spk.round(1).to_string())
print("\nTotals:", spk.select_dtypes("number").drop(columns=["spikethresh", "max rise, 99th percentile of all trials"]).sum().to_dict())
# Example traces of the six neurons where the rules disagree most; potential corrected by -12 mV.
show_n = (spk["spikethresh only"] + spk["claims file only"]).nlargest(6).index
fig, ax = plt.subplots(2, 3, figsize=(15, 7))
tms = (np.arange(ac2.TRACE) - ac2.PRE) / ac2.FS * 1000
plotted = []
for a, name in zip(ax.ravel(), show_n):
n = N[name]
st = ac2.stable(n)
v = n["r"] + n["rbaseline"][:, None] - ac2.JUNCTION_MV
cell, claims = ac2.spikes(n, "cell") & st, ac2.spikes(n, "claims") & st
for mask, colour in ((cell & ~claims, "tab:blue"), (claims & ~cell, "tab:red"), (cell & claims, "k"), (st & ~cell & ~claims, "0.6")):
for i in np.flatnonzero(mask)[:5]:
a.plot(tms, v[i], color=colour, lw=0.7)
plotted.append(v[i])
a.axhline(-20, ls=":", color="k")
a.set(title=f"{name}: spikethresh only {int((cell & ~claims).sum())}, claims only {int((claims & ~cell).sum())}",
xlabel="ms from onset", ylabel="mV (corrected)")
fig.suptitle("blue: removed by spikethresh only; red: by the claims file's rule only; black: by both; grey: kept")
plt.tight_layout()
plt.show()
print("Plotted range %.1f to %.1f mV (all within the axes); up to 5 traces per class and neuron" % (np.min(plotted), np.max(plotted)))
# The trials only spikethresh removes, by kind: the corrected potential crosses -20 mV with a rise
# below the 10 mV/ms of the claims file's rule; it stays below -20 mV with a fast rise (above 5
# mV/ms); it stays below -20 mV with a slow rise (at most 5 mV/ms).
kinds = []
for name, n in N.items():
st = ac2.stable(n)
only = ac2.spikes(n, "cell") & ~ac2.spikes(n, "claims") & st
if not only.any():
continue
peak = (n["r"] + n["rbaseline"][:, None] - ac2.JUNCTION_MV).max(1)[only]
rise = ac2.max_rise(n["r"])[only]
kinds.append({"neuron": name, "spikethresh only": int(only.sum()),
"above -20 mV, rise < 10": int((peak > -20).sum()),
"below -20 mV, rise > 5": int(((peak <= -20) & (rise > ac2.SLOW_DVDT)).sum()),
"below -20 mV, rise <= 5": int(((peak <= -20) & (rise <= ac2.SLOW_DVDT)).sum()),
"peak mV (corrected)": f"{peak.min():.1f} to {peak.max():.1f}"})
kinds = pd.DataFrame(kinds).set_index("neuron")
print(kinds.to_string())
print("Totals:", kinds.select_dtypes("number").sum().to_dict())
# Trials removed by the two further spike rules shown as readings: the corrected potential above
# -20 mV at any rise, and spikethresh with the slow events kept (largest rise above 5 mV/ms).
for label, rule in ((SPK_ANY, "corrected"), (SPK_SLOW, "cell_fast")):
print(f"[{label}] trials removed: {sum(int((ac2.spikes(n, rule) & ac2.stable(n)).sum()) for n in N.values())}")
spikethresh above -20 mV uncorrected above -20 mV uncorrected, max rise mV/ms above -20 mV corrected above -20 mV corrected, max rise mV/ms fast rise and above -20 mV corrected (claims file) fast rise and above -20 mV corrected (claims file), max rise mV/ms spikethresh (data set) spikethresh (data set), max rise mV/ms spikethresh only claims file only max rise, 99th percentile of all trials
neuron
020701md01b 50.0 0 0 0 0 0 0 2.0
020901md01 50.0 25 4.4 to 16.1 0 0 0 0 0 14.7
021001md01 40.0 1 96.2 to 96.2 1 96.2 to 96.2 1 96.2 to 96.2 1 96.2 to 96.2 0 0 2.9
021401md01b 45.0 2 2.4 to 2.9 0 0 0 0 0 3.4
021401md02 50.0 0 0 0 0 0 0 2.4
021401md03 50.0 1 45.4 to 45.4 1 45.4 to 45.4 1 45.4 to 45.4 1 45.4 to 45.4 0 0 2.9
021401md05 44.0 15 2.4 to 38.1 3 10.7 to 38.1 3 10.7 to 38.1 3 10.7 to 38.1 0 0 7.8
021901md01 29.0 5 23.4 to 46.4 5 23.4 to 46.4 5 23.4 to 46.4 9 2.0 to 46.4 4 0 9.5
022301md01 30.0 35 8.3 to 26.4 11 8.3 to 26.4 10 16.1 to 26.4 32 2.4 to 26.4 24 2 23.0
022801md01b 40.0 0 0 0 0 0 0 4.4
022801md02b 50.0 0 0 0 0 0 0 7.8
030201md01b 50.0 34 30.5 to 61.1 34 30.5 to 61.1 34 30.5 to 61.1 34 30.5 to 61.1 0 0 56.6
030201md03b 50.0 0 0 0 0 0 0 7.3
030301md01b 55.0 1 3.9 to 3.9 1 3.9 to 3.9 0 0 0 0 3.9
030301md02b 50.0 2 6.3 to 7.3 0 0 0 0 0 7.3
030301md03b 35.0 12 1.0 to 43.0 3 27.8 to 43.0 3 27.8 to 43.0 4 9.3 to 43.0 1 0 6.3
100300md01 45.0 1 1.0 to 1.0 0 0 0 0 0 2.4
100300md02 40.0 28 1.0 to 4.9 21 1.0 to 4.9 0 19 2.0 to 4.9 19 0 6.3
100500md01 40.0 14 20.0 to 56.6 5 39.6 to 56.6 5 39.6 to 56.6 15 13.2 to 56.6 10 0 40.7
100500md02b 45.0 5 3.9 to 7.3 5 3.9 to 7.3 0 5 3.9 to 7.3 5 0 3.9
102400md01 45.0 73 0.5 to 51.8 14 0.5 to 51.8 3 39.6 to 51.8 5 2.9 to 43.9 3 1 8.4
102500md02 45.0 51 1.0 to 38.1 46 2.4 to 38.1 10 12.2 to 38.1 34 4.4 to 38.1 25 1 8.8
102600md01 45.0 1 0.5 to 0.5 1 0.5 to 0.5 0 0 0 0 2.9
102600md03 40.0 3 1.0 to 2.4 1 2.0 to 2.0 0 2 2.0 to 2.0 2 0 2.4
102800md01 45.0 0 0 0 0 0 0 2.4
102800md02 35.0 36 1.5 to 2.9 26 1.5 to 2.9 0 31 1.5 to 2.9 31 0 2.9
102900md02 45.0 28 1.5 to 41.0 26 3.9 to 41.0 13 10.3 to 41.0 23 2.4 to 37.1 12 2 26.6
102900md04 32.0 42 1.5 to 35.6 29 2.9 to 35.6 2 26.4 to 35.6 29 2.9 to 35.6 27 0 8.1
110700md01 28.0 9 5.9 to 27.3 1 27.3 to 27.3 1 27.3 to 27.3 13 2.9 to 27.3 12 0 21.8
110900md05 40.0 6 2.0 to 3.9 3 3.4 to 3.9 0 5 2.0 to 3.9 5 0 2.4
111200md01 45.0 18 6.8 to 58.6 18 6.8 to 58.6 17 28.8 to 58.6 19 2.9 to 58.6 2 0 46.0
111400md01 38.0 43 1.0 to 38.6 38 4.4 to 38.6 37 17.1 to 38.6 37 17.1 to 38.6 1 1 37.2
111500md01 45.0 22 2.9 to 91.3 18 5.4 to 91.3 17 25.4 to 91.3 21 4.9 to 91.3 4 0 77.0
Totals: {'above -20 mV uncorrected': 513, 'above -20 mV corrected': 311, 'fast rise and above -20 mV corrected (claims file)': 162, 'spikethresh (data set)': 342, 'spikethresh only': 187, 'claims file only': 7}
Plotted range -79.4 to 3.4 mV (all within the axes); up to 5 traces per class and neuron
spikethresh only above -20 mV, rise < 10 below -20 mV, rise > 5 below -20 mV, rise <= 5 peak mV (corrected)
neuron
021901md01 4 0 0 4 -47.3 to -43.7
022301md01 24 1 22 1 -42.4 to -14.9
030301md03b 1 0 1 0 -20.1 to -20.1
100300md02 19 14 0 5 -32.9 to -2.7
100500md01 10 0 10 0 -34.8 to -21.2
100500md02b 5 5 0 0 -4.7 to 3.1
102400md01 3 3 0 0 -16.6 to 21.9
102500md02 25 25 0 0 -14.0 to 7.2
102600md03 2 1 0 1 -43.7 to -10.9
102800md02 31 24 0 7 -29.0 to -9.3
102900md02 12 11 0 1 -29.5 to -3.7
102900md04 27 26 0 1 -35.6 to -1.1
110700md01 12 0 8 4 -43.6 to -21.0
110900md05 5 3 0 2 -30.4 to -3.3
111200md01 2 0 0 2 -41.1 to -37.4
111400md01 1 0 1 0 -21.0 to -21.0
111500md01 4 1 2 1 -44.1 to -7.8
Totals: {'spikethresh only': 187, 'above -20 mV, rise < 10': 114, 'below -20 mV, rise > 5': 44, 'below -20 mV, rise <= 5': 29}
[spikes: above -20 mV, any rise] trials removed: 311
[spikes: spikethresh, slow events kept] trials removed: 247
Of 342 trials that spikethresh removes over the 33 neurons, the claims file's rule (fast rise
and above -20 mV after correction) also removes 155. It removes 162 in all, 7 that spikethresh
keeps (table above). The 187 trials that only spikethresh removes are of three kinds (second
table). 114 cross -20 mV after correction with a largest rise below the 10 mV/ms chosen here for
the "fast rise": broad spikes that peak between -20 and +7 mV, as in 102500md02 (25 trials),
102900md04 (26) and 102900md02 (11), and the slow depolarizations of 102800md02 (24) and
100300md02 (14). In 102800md02 these rise by 1.5 to 2.9 mV/ms, no faster than that neuron's 99th
percentile over all trials (2.9 mV/ms), and go from about -70 to -10 mV over about 50 ms; they
may be regenerative events that QX-314 leaves in place of fast spikes, or large PSPs, which the
data cannot tell apart. 44 stay below -20 mV with a rise above 5 mV/ms, narrow spikes as in
022301md01 (22 trials) and 100500md01 (10). 29 stay below -20 mV with a rise of at most 5 mV/ms.
The spike rule changes the results (overview above). Two further rules bound the choice: the
corrected potential above -20 mV at any rise (311 trials) gives 4.43 mV, close to
spikethresh's 4.58, and spikethresh with the slow events kept (largest rise above 5 mV/ms,
247 trials) gives 5.11 mV, about two thirds of the way to the claims file's rule with the 10 mV/ms
rise (5.35 mV, 162 trials). So the difference between the claims file's rule and spikethresh
comes from the 10 mV/ms rise criterion, and most of it from the slow events. Applied without the
correction and without the fast rise, as in the earlier version of this notebook, -20 mV marks
513 trials, many of them subthreshold PSPs (in 102400md01, 73 trials with largest rises from 0.5
mV/ms). Under step 1.6 of the protocol all four rules are shown for every claim. The per-neuron
threshold, the authors' own parameter, is adopted as an assumption, and outcomes that differ
between the rules are marked as resting on it. If the slow events are large PSPs, removing them
selects on the response and lowers q.
fig3a-reliable-neuron¶
For the most reliable neuron, the responses to 20 presentations of a 32-kHz tone cluster tightly about the mean (q = 0.23 mV), and its variability averaged over the 26 tones with a mean PSP above 3 mV is 0.88 mV.
This section comes first because fig2-stereotyped-shape uses the neuron it identifies.
def tone_table(label, name):
return pd.DataFrame([{k: v for k, v in r.items() if np.isscalar(v)} for r in PT[label][name]]).set_index("tone")
for label in ("main", SPK, SPK_ANY, SPK_SLOW, "peak in culling window (1b)", "variance over N (3)"):
df = S[label]
best = df.q.idxmin()
near = (df.q - 0.88).abs().idxmin()
tt = tone_table(label, best)
print(f"[{label}] lowest <q>tones: {best} <q> = {df.q[best]:.3f} mV over {df.n_tones[best]} tones;"
f" neuron nearest 0.88 mV: {near} ({df.q[near]:.3f}); repetitions {int(cells.reps[best])};"
f" q at 31,118 Hz {tt.q[31118]:.3f}, at 34,448 Hz {tt.q[34448]:.3f}")
EXAMPLE = S["main"].q.idxmin()
print("\nLowest five <q>tones (main):", S["main"].q.nsmallest(5).round(3).to_dict())
tt = tone_table("main", EXAMPLE)
print(tt[["n", "mean", "q", "peak_ms"]].round(3).to_string())
# As Figure 3A: the 20 responses to 31,118 Hz with their mean, and the histogram of q over the
# included tones (inset of Figure 3A); Figure 3B, in the main section below.
r31 = next(r for r in PT["main"][EXAMPLE] if r["tone"] == 31118)
tr = N[EXAMPLE]["r"][r31["trials"]]
tms = (np.arange(ac2.TRACE) - ac2.PRE) / ac2.FS * 1000
fig, ax = plt.subplots(1, 2, figsize=(10, 3.5))
ax[0].plot(tms, tr.T, color="k", lw=0.6)
ax[0].plot(tms, tr.mean(0), color="0.6", lw=3)
ax[0].set(xlabel="ms from onset", ylabel="mV", title=f"{EXAMPLE}, 31,118 Hz, {len(tr)} trials, q = {r31['q']:.2f} mV")
q_inc = [r["q"] for r in PT["main"][EXAMPLE] if ac2.included(r)]
ax[1].hist(q_inc, bins=np.arange(0, 2.51, 0.25), color="0.5", edgecolor="k")
ax[1].axvline(np.mean(q_inc), ls="--", color="k")
ax[1].set(xlabel="q (mV)", ylabel="tones", title=f"N = {len(q_inc)} tones")
plt.tight_layout()
plt.show()
print(f"Plotted: traces {tr.min():.2f} to {tr.max():.2f} mV; q {min(q_inc):.2f} to {max(q_inc):.2f} mV (all within the axes);"
f" counts per 0.25-mV bin {np.histogram(q_inc, bins=np.arange(0, 2.51, 0.25))[0].tolist()}")
[main] lowest <q>tones: 022801md01b <q> = 0.876 mV over 26 tones; neuron nearest 0.88 mV: 022801md01b (0.876); repetitions 20; q at 31,118 Hz 0.229, at 34,448 Hz 0.648
[spikes: fast rise above -20 mV (6)] lowest <q>tones: 022801md01b <q> = 0.876 mV over 26 tones; neuron nearest 0.88 mV: 022801md01b (0.876); repetitions 20; q at 31,118 Hz 0.229, at 34,448 Hz 0.648
[spikes: above -20 mV, any rise] lowest <q>tones: 022801md01b <q> = 0.876 mV over 26 tones; neuron nearest 0.88 mV: 022801md01b (0.876); repetitions 20; q at 31,118 Hz 0.229, at 34,448 Hz 0.648
[spikes: spikethresh, slow events kept] lowest <q>tones: 022801md01b <q> = 0.876 mV over 26 tones; neuron nearest 0.88 mV: 022801md01b (0.876); repetitions 20; q at 31,118 Hz 0.229, at 34,448 Hz 0.648
[peak in culling window (1b)] lowest <q>tones: 022801md01b <q> = 0.876 mV over 26 tones; neuron nearest 0.88 mV: 022801md01b (0.876); repetitions 20; q at 31,118 Hz 0.229, at 34,448 Hz 0.648
[variance over N (3)] lowest <q>tones: 022801md01b <q> = 0.832 mV over 26 tones; neuron nearest 0.88 mV: 022801md01b (0.832); repetitions 20; q at 31,118 Hz 0.217, at 34,448 Hz 0.615
Lowest five <q>tones (main): {'022801md01b': 0.876, '021401md03': 1.2, '021001md01': 1.21, '030201md01b': 1.231, '030201md03b': 1.575}
n mean q peak_ms
tone
4074 20 7.712 0.507 33.00
38134 20 18.536 0.479 28.00
10172 20 8.682 1.015 31.50
7498 20 5.401 0.878 36.75
16909 20 2.120 1.059 37.25
3681 20 4.495 0.629 37.00
31118 20 15.507 0.229 26.50
2713 20 3.001 0.701 38.00
15275 20 3.122 1.568 34.50
3325 20 2.940 0.591 38.75
20722 20 7.747 0.835 30.50
2451 20 2.099 0.323 37.25
28110 20 11.309 0.494 26.75
18719 20 3.574 1.282 32.25
4993 20 3.678 1.412 34.25
42214 20 17.485 0.310 28.50
6773 20 5.023 0.642 33.25
13799 20 2.680 1.152 34.00
46731 20 9.171 0.378 35.75
11260 20 12.004 0.475 30.25
2000 20 3.421 1.073 40.50
9188 20 11.458 0.860 30.00
2214 20 3.087 1.004 39.50
8300 20 5.378 0.903 36.00
12465 20 2.100 1.429 34.75
25393 20 4.162 0.819 30.75
3003 20 3.168 2.042 41.50
4510 20 3.250 1.249 36.25
6119 20 1.417 0.516 38.00
34448 20 12.372 0.648 28.00
22939 20 4.852 0.606 31.00
5527 20 3.235 1.742 35.00
Plotted: traces -2.69 to 18.82 mV; q 0.23 to 2.04 mV (all within the axes); counts per 0.25-mV bin [1, 5, 6, 5, 4, 2, 2, 0, 1, 0]
Result: reproduced. The neuron with the lowest tones, 022801md01b, has
tones = 0.88 mV
over 26 tones with a mean PSP above 3 mV and 20 repetitions per tone, as in the paper, under the
main reading and with the claims file's spike rule (no spike trials in this neuron). Its q for the
31,118-Hz tone is 0.229 mV (reading a of "32 kHz", the nearest tone); for 34,448 Hz (reading b)
it is 0.648 mV. The neuron and the 0.88 mV do not depend on the spike rule or the peak window;
with the variance over N,
tones is 0.832 mV and q at 31,118 Hz 0.217 mV, both within
tolerance. The two further spike rules give the same values. The panels match Figure 3A, which is
labelled 022801md01b: the 20 responses to 31 kHz rise about 15 ms after onset to about 15 mV and
decay over about 60 ms, two of them falling below baseline late in the trace, and the q of the
26 tones spread from 0.23 to 2.04 mV around their mean of 0.88, as in the inset (0 to 2 mV).
fig2-stereotyped-shape¶
Tones evoke brief, short-latency PSPs in all neurons. In the example neuron the mean responses to the 32 tones peak 27 to 42 ms after onset with heights of 1.2 to 18 mV, and are similar in shape once normalized and aligned.
def mean_traces(n):
idx_ok = ac2.stable(n) & ~ac2.spikes(n, "cell")
return {int(t): n["r"][idx_ok & (n["tone"] == t)].astype(float).mean(0) for t in ac2.TONE_ORDER}
for window, (lo, hi) in {"after onset (1a)": (ac2.PRE, ac2.TRACE), "culling window (1b)": (ac2.PRE, 400)}.items():
mt = mean_traces(N[EXAMPLE])
lat = {t: (lo + np.argmax(m[lo:hi]) - ac2.PRE) / ac2.FS * 1000 for t, m in mt.items()}
hgt = {t: m[lo:hi].max() for t, m in mt.items()}
print(f"[{window}] {EXAMPLE}: latency to peak {min(lat.values()):.2f} to {max(lat.values()):.2f} ms;"
f" peak height {min(hgt.values()):.2f} to {max(hgt.values()):.2f} mV (32 tones)")
lat = pd.Series({t: (ac2.PRE + np.argmax(m[ac2.PRE:]) - ac2.PRE) / ac2.FS * 1000 for t, m in mt.items()})
hgt = pd.Series({t: m[ac2.PRE:].max() for t, m in mt.items()})
print("Latencies (ms) by tone:", lat.round(2).sort_index().to_dict())
print("Tones with latency outside 24.3-46.2 ms:", lat[(lat < 24.3) | (lat > 46.2)].round(2).to_dict())
# Shape: normalize each mean trace to its peak and align it on the peak; window = the culling
# window placed for a tone at the median latency.
med = ac2.PRE + int(round(np.median(lat.values) * ac2.FS / 1000))
before, after = med, 400 - med
aligned = []
for t, m in mt.items():
pk = ac2.PRE + int(np.argmax(m[ac2.PRE:]))
padded = np.pad(m, (before, after), mode="edge") # edge values where the window runs off the trace
aligned.append(padded[pk:pk + before + after] / m[pk])
aligned = np.array(aligned)
avg = aligned.mean(0)
cors = np.array([np.corrcoef(x, avg)[0, 1] for x in aligned])
print(f"Aligned traces: {len(aligned)} of 32, window {-before / 4:.2f} to {after / 4:.2f} ms around the peak; correlation with their average: median {np.median(cors):.3f},"
f" min {cors.min():.3f}, 10th percentile {np.percentile(cors, 10):.3f}")
# Unaligned (no shift), for reference.
raw = np.array([m[:400] / m[ac2.PRE:].max() for m in mt.values()])
cr = np.array([np.corrcoef(x, raw.mean(0))[0, 1] for x in raw])
print(f"Without alignment: median correlation {np.median(cr):.3f}")
# All neurons: latency of the mean response to each neuron's largest tone, with the peak (and the
# largest tone) searched after onset (1a) or in the culling window, up to 85 ms after onset (1b),
# and the response onset of the data check (mean over all tones first above 1 mV).
big = []
for name, n in N.items():
m = mean_traces(n)
row = {"neuron": name}
for label, hi in (("1a", ac2.TRACE), ("1b", 400)):
t_big = max(m, key=lambda t: m[t][ac2.PRE:hi].max())
row[f"tone ({label})"] = t_big
row[f"peak mV ({label})"] = m[t_big][ac2.PRE:hi].max()
row[f"latency ms ({label})"] = np.argmax(m[t_big][ac2.PRE:hi]) / ac2.FS * 1000
row["response onset ms"] = clock.onset_latency_ms[name]
big.append(row)
big = pd.DataFrame(big).set_index("neuron")
for label in ("1a", "1b"):
x = big[f"latency ms ({label})"]
print(f"\nLargest tone per neuron ({label}): latency {x.min():.2f} to {x.max():.2f} ms (median {x.median():.2f});"
f" within 50 ms in {int((x <= 50).sum())} of {len(x)}")
o = big["response onset ms"]
print(f"Response onset: {o.min():.2f} to {o.max():.2f} ms (median {o.median():.2f}); within 50 ms in {int((o <= 50).sum())} of {len(o)}")
print(big[big["latency ms (1a)"] > 50].round(2).to_string())
fig, ax = plt.subplots(1, 2, figsize=(10, 3.5))
tms = (np.arange(ac2.TRACE) - ac2.PRE) / ac2.FS * 1000
for m in mt.values():
ax[0].plot(tms, m, lw=0.7)
for x in aligned:
ax[1].plot((np.arange(before + after) - before) / ac2.FS * 1000, x, lw=0.7)
ax[0].set(xlabel="ms from onset", ylabel="mV", title=f"{EXAMPLE}: mean per tone (all within axes)")
ax[1].set(xlabel="ms from peak", ylabel="normalized", title="normalized and aligned")
plt.tight_layout()
plt.show()
print("Plotted ranges: mean traces %.2f to %.2f mV; normalized %.2f to %.2f" % (
min(m.min() for m in mt.values()), max(m.max() for m in mt.values()), aligned.min(), aligned.max()))
[after onset (1a)] 022801md01b: latency to peak 26.50 to 41.50 ms; peak height 1.19 to 18.12 mV (32 tones)
[culling window (1b)] 022801md01b: latency to peak 26.50 to 41.50 ms; peak height 1.19 to 18.12 mV (32 tones)
Latencies (ms) by tone: {2000: 40.5, 2214: 39.5, 2451: 37.25, 2713: 38.0, 3003: 41.5, 3325: 38.75, 3681: 37.0, 4074: 33.0, 4510: 36.25, 4993: 34.25, 5527: 35.0, 6119: 38.0, 6773: 33.25, 7498: 36.75, 8300: 36.0, 9188: 30.0, 10172: 31.5, 11260: 30.25, 12465: 34.75, 13799: 34.0, 15275: 34.5, 16909: 37.25, 18719: 32.25, 20722: 30.5, 22939: 31.0, 25393: 30.75, 28110: 26.75, 31118: 26.5, 34448: 28.0, 38134: 28.0, 42214: 28.5, 46731: 35.75}
Tones with latency outside 24.3-46.2 ms: {}
Aligned traces: 32 of 32, window -49.50 to 50.50 ms around the peak; correlation with their average: median 0.991, min 0.933, 10th percentile 0.974
Without alignment: median correlation 0.973
Largest tone per neuron (1a): latency 20.75 to 115.50 ms (median 45.25); within 50 ms in 19 of 33
Largest tone per neuron (1b): latency 20.75 to 84.75 ms (median 45.25); within 50 ms in 19 of 33
Response onset: 14.25 to 38.50 ms (median 22.12); within 50 ms in 32 of 33
tone (1a) peak mV (1a) latency ms (1a) tone (1b) peak mV (1b) latency ms (1b) response onset ms
neuron
021401md01b 12465 12.38 66.00 12465 12.38 66.00 26.25
100300md01 8300 10.89 91.50 8300 9.81 84.75 30.75
100300md02 7498 26.09 91.25 7498 25.43 73.00 23.00
100500md02b 46731 17.97 96.75 46731 17.27 84.75 29.75
102400md01 4074 24.97 115.50 4074 22.51 84.75 26.50
102600md01 10172 9.20 68.75 10172 9.20 68.75 29.75
102600md03 6773 7.55 53.00 6773 7.55 53.00 34.25
102800md02 9188 16.12 95.50 9188 15.27 84.75 29.00
102900md02 42214 17.82 82.50 42214 17.82 82.50 22.00
102900md04 31118 16.78 88.25 31118 16.01 84.75 20.50
110700md01 18719 20.16 68.75 18719 20.16 68.75 23.50
110900md05 6773 10.82 58.25 6773 10.82 58.25 22.25
111200md01 6119 20.16 65.00 6119 20.16 65.00 18.75
111500md01 46731 36.72 77.50 46731 36.72 77.50 17.00
Plotted ranges: mean traces -1.07 to 18.12 mV; normalized -0.27 to 1.00
Result: partially reproduced. In 022801md01b (20 trials per tone, the Figure 2 legend's 20 presentations) the mean responses to the 32 tones peak 26.5 to 41.5 ms after onset with heights of 1.19 to 18.12 mV, within the widened ranges (24.3 to 46.2 ms; 1.08 to 19.8 mV), under both peak windows. Normalized and aligned on their peaks, the 32 mean traces correlate with their average with a median r of 0.991 (minimum 0.933), above the 0.9 the criterion sets. The "all neurons" part does not hold under the claims file's 50-ms bound: the mean response to each neuron's largest tone peaks within 50 ms of onset in 19 of 33 neurons, with the peak searched after onset (20.75 to 115.5 ms, median 45.25) or in the culling window (20.75 to 84.75 ms). The responses begin early in every neuron: the mean over all tones first exceeds 1 mV 14.25 to 38.5 ms after onset in 32 neurons (021401md02 never exceeds it). The late peaks of the other 14 come from responses that keep rising after onset, which the 50-ms bound on the peak counts against. The panels match Figure 2 (labelled 022801md01b, the neuron identified in fig3a): mean responses of 1 to 18 mV that rise about 15 ms after onset, and normalized, aligned traces that overlap closely, with a few decaying below zero.
fig3b-total-variability-distribution¶
Across the 33 neurons the tone-averaged variability index ranges over more than an order of magnitude, from below 1 to above 10 mV; about a fifth (7/33) have low variability (1 to 2 mV) and about half (18/33) have high variability (4 to 12 mV).
rows = []
for label, df in S.items():
if label in SUMMARY_READINGS:
continue
q = df.q
rows.append({"reading": label, "min": q.min(), "max": q.max(), "1<q<2": int(((q > 1) & (q < 2)).sum()),
"q<2": int((q < 2).sum()), "4<q<12": int(((q > 4) & (q < 12)).sum())})
print(pd.DataFrame(rows).set_index("reading").round(2).to_string())
print("\n<q>tones per neuron (main), sorted:", main.q.sort_values().round(2).tolist())
print("Neurons by anesthetic within 4<q<12 (main):", main[(main.q > 4) & (main.q < 12)].anesthetic.value_counts().to_dict())
# As Figure 3B: <q>tones of the 33 neurons in 2-mV bins.
counts, edges = np.histogram(main.q, bins=np.arange(0, 12.1, 2))
fig, ax = plt.subplots(figsize=(4.5, 3))
ax.hist(main.q, bins=edges, color="0.5", edgecolor="k")
ax.set(xlabel="<q>tones (mV)", ylabel="neurons", title="N = 33 neurons (main)")
plt.tight_layout()
plt.show()
print("Neurons per 2-mV bin, 0 to 12 mV (main):", counts.tolist(), "; all 33 within the axes:", int(counts.sum()) == 33)
min max 1<q<2 q<2 4<q<12
reading
main 0.88 11.21 6 7 18
spikes: fast rise above -20 mV (6) 0.88 11.13 5 6 20
spikes: above -20 mV, any rise 0.88 10.94 5 6 18
spikes: spikethresh, slow events kept 0.88 11.33 6 7 20
peak in culling window (1b) 0.88 10.15 7 8 15
peak-to-trough over trace (2b) 0.88 11.21 6 7 18
variance over N (3) 0.83 10.61 6 7 16
no negative-f rule in PSP culling (4b) 0.88 11.21 6 7 18
cullings intersected 0.88 11.21 6 7 18
heights re-measured after culling 0.88 11.21 6 7 18
<q>tones per neuron (main), sorted: [0.88, 1.2, 1.21, 1.23, 1.57, 1.6, 1.7, 2.2, 2.48, 2.51, 2.71, 3.76, 3.79, 3.87, 3.92, 4.19, 4.33, 4.79, 4.87, 5.1, 5.14, 5.19, 5.4, 5.45, 5.91, 6.19, 6.3, 7.15, 7.36, 8.15, 8.9, 10.94, 11.21]
Neurons by anesthetic within 4<q<12 (main): {'U': 10, 'P': 6, 'P_with_Diaz': 2}
Neurons per 2-mV bin, 0 to 12 mV (main): [7, 8, 10, 4, 2, 2] ; all 33 within the axes: True
Result: reproduced (resting on the spike-rule assumption and on detail 1). Under the main
reading tones ranges from 0.88 to 11.21 mV (below 1 and above 10 mV). 18 neurons lie between 4
and 12 mV. 6 lie strictly between 1 and 2 mV, outside 7 +- 0.7. Under the criterion's second
reading of the low count (q < 2 mV, which includes the 0.88 mV example neuron) it is 7. Both
counts hold under that reading. With the claims file's spike rule the counts are 5 (6 below 2 mV)
and 20, both outside, which would give not reproduced. With the peak in the culling window (1b)
they are 7 (8) and 15, which would give partially reproduced. With the corrected potential
above -20 mV at any rise the counts are 5 (6) and 18, and with
spikethresh keeping the slow
events 6 (7) and 20; each would give partially reproduced. In 2-mV bins the 33 values give 7, 8,
10, 4, 2 and 2 neurons, against 7, 8, 11, 3, 2 and 2 read from Figure 3B.
total-variability-population¶
The mean total variability index across neurons is 4.6 mV.
for label, df in S.items():
if label not in SUMMARY_READINGS:
print(f"{label:42s} <q>neurons = {mse(df.q)}")
print("\nBy anesthetic (main):", main.groupby("anesthetic").q.agg(["mean", "count"]).round(2).to_dict("index"))
main <q>neurons = 4.58 +- 0.47 (n = 33)
spikes: fast rise above -20 mV (6) <q>neurons = 5.35 +- 0.54 (n = 33)
spikes: above -20 mV, any rise <q>neurons = 4.43 +- 0.43 (n = 33)
spikes: spikethresh, slow events kept <q>neurons = 5.11 +- 0.51 (n = 33)
peak in culling window (1b) <q>neurons = 4.22 +- 0.44 (n = 33)
peak-to-trough over trace (2b) <q>neurons = 4.58 +- 0.47 (n = 33)
variance over N (3) <q>neurons = 4.29 +- 0.45 (n = 33)
no negative-f rule in PSP culling (4b) <q>neurons = 4.58 +- 0.47 (n = 33)
cullings intersected <q>neurons = 4.58 +- 0.47 (n = 33)
heights re-measured after culling <q>neurons = 4.58 +- 0.47 (n = 33)
By anesthetic (main): {'P': {'mean': 3.69, 'count': 13}, 'P_with_Diaz': {'mean': 6.15, 'count': 3}, 'U': {'mean': 4.99, 'count': 17}}
Result: reproduced (resting on the spike-rule assumption). neurons is 4.58 +- 0.47 mV (SE,
n = 33), against 4.6 +- 0.5. With the claims file's spike rule it is 5.35 mV, outside the
tolerance (4.1 to 5.1); with the peak searched in the culling window 4.22, and with the variance
over N 4.29, both within it. With the corrected potential above -20 mV at any rise it is 4.43
mV, within; with
spikethresh keeping the slow events 5.11 mV, just outside.
lfp-culling-reduction¶
LFP-based culling reduces the mean variability index by 35%, from 4.8 to 3.1 mV, across the 29 neurons with an LFP.
print("Neurons with the LFP flag:", len(lfp_neurons))
for label, df in S.items():
d = df[df.lfp]
print(f"{label:42s} before {mse(d.q)}; after LFP-based culling {mse(d.q_lfp)};"
f" fall {100 * (1 - d.q_lfp.mean() / d.q.mean()):.0f}%")
d = main[main.lfp]
print("\nTones per neuron before and after LFP-based culling (main):")
print(d[["n_tones", "n_tones_lfp", "q", "q_lfp"]].round(2).T.to_string())
Neurons with the LFP flag: 29 main before 4.80 +- 0.52 (n = 29); after LFP-based culling 3.19 +- 0.43 (n = 29); fall 33% spikes: fast rise above -20 mV (6) before 5.41 +- 0.59 (n = 29); after LFP-based culling 3.75 +- 0.56 (n = 29); fall 31% spikes: above -20 mV, any rise before 4.63 +- 0.47 (n = 29); after LFP-based culling 3.13 +- 0.40 (n = 29); fall 32% spikes: spikethresh, slow events kept before 5.22 +- 0.56 (n = 29); after LFP-based culling 3.41 +- 0.48 (n = 29); fall 35% peak in culling window (1b) before 4.43 +- 0.48 (n = 29); after LFP-based culling 2.99 +- 0.41 (n = 29); fall 33% peak-to-trough over trace (2b) before 4.80 +- 0.52 (n = 29); after LFP-based culling 3.12 +- 0.36 (n = 29); fall 35% variance over N (3) before 4.50 +- 0.50 (n = 29); after LFP-based culling 2.94 +- 0.40 (n = 29); fall 35% no negative-f rule in PSP culling (4b) before 4.80 +- 0.52 (n = 29); after LFP-based culling 3.19 +- 0.43 (n = 29); fall 33% cullings intersected before 4.80 +- 0.52 (n = 29); after LFP-based culling 3.19 +- 0.43 (n = 29); fall 33% heights re-measured after culling before 4.80 +- 0.52 (n = 29); after LFP-based culling 3.07 +- 0.37 (n = 29); fall 36% criteria after culling (5b) before 4.80 +- 0.52 (n = 29); after LFP-based culling 3.04 +- 0.35 (n = 29); fall 37% Tones per neuron before and after LFP-based culling (main): neuron 020701md01b 020901md01 021001md01 021401md01b 021401md02 021401md03 021401md05 021901md01 022301md01 022801md01b 022801md02b 030201md01b 030201md03b 030301md01b 030301md02b 030301md03b 102400md01 102500md02 102600md01 102600md03 102800md01 102800md02 102900md02 102900md04 110700md01 110900md05 111200md01 111400md01 111500md01 n_tones 18.00 27.00 6.00 17.00 4.00 14.00 25.00 32.00 29.00 26.00 29.00 6.00 32.00 20.00 16.00 18.00 32.00 30.00 16.00 12.00 22.00 29.00 30.00 25.00 31.00 12.00 32.00 25.00 28.00 n_tones_lfp 17.00 11.00 3.00 12.00 1.00 3.00 7.00 29.00 18.00 19.00 23.00 4.00 25.00 7.00 1.00 12.00 32.00 9.00 11.00 10.00 10.00 23.00 29.00 9.00 20.00 2.00 27.00 7.00 15.00 q 1.70 10.94 1.21 4.87 5.19 1.20 5.91 8.15 6.30 0.88 2.48 1.23 1.57 7.15 2.20 5.40 11.21 8.90 3.79 6.19 2.51 4.19 5.45 4.79 5.14 5.10 3.76 7.36 4.33 q_lfp 1.04 5.90 0.74 3.66 3.13 0.69 2.23 10.41 2.24 0.69 2.42 1.17 1.24 2.35 0.56 1.31 7.11 4.45 2.61 4.27 1.47 5.33 5.04 5.51 2.60 6.19 2.24 3.89 2.11
Result: reproduced (resting on the spike-rule assumption and on detail 1). Over the 29
neurons with the LFP flag, is 4.80 +- 0.52 mV before LFP-based culling and 3.19 +- 0.43 mV
after it (paper 4.8 +- 0.5 and 3.1 +- 0.1; both within tolerance, the second by 0.01 mV), a fall
of 33%. The peak-to-trough height over the trace (2b, 3.12 mV), the criteria after culling (5b,
3.04 mV) and heights re-measured after culling (3.07 mV) keep both values within tolerance. The
peak searched in the culling window (1b) gives 2.99 mV after culling, and the variance over N
2.94 mV, both outside, which would not give reproduced. With the claims file's spike rule the values are 5.41 and 3.75 mV,
both outside; with the corrected potential above -20 mV at any rise 4.63 and 3.13 mV, both
within; with
spikethresh keeping the slow events 5.22 and 3.41 mV, the second outside.
fig5b-residual-lfp-psp-correlation¶
After LFP-based culling, PSP height still correlates with LFP height across trials for many tones, as in the 60 unculled responses to one tone of the Figure 5 neuron.
def residual_correlations(label):
out = []
for name in lfp_neurons:
for r in PT[label][name]:
if not ac2.included(r, "lfp"):
continue
k = r["keep_lfp"]
h, s, lh = r["h"][k], r["s"][k], r["lfp_h"][k]
rs, ps = stats.pearsonr(s, h)
rl, pl = stats.pearsonr(lh, h)
out.append({"neuron": name, "tone": r["tone"], "n": int(k.sum()), "r_scale": rs, "p_scale": ps,
"r_lfp_height": rl, "p_lfp_height": pl})
return pd.DataFrame(out)
for label in ("main", SPK, SPK_ANY, SPK_SLOW, "peak in culling window (1b)", "heights re-measured after culling"):
rc = residual_correlations(label)
print(f"[{label}] {len(rc)} tones of {rc.neuron.nunique()} neurons")
for col in ("scale", "lfp_height"):
r_, p_ = rc[f"r_{col}"], rc[f"p_{col}"]
print(f" LFP {col:10s}: positive {100 * (r_ > 0).mean():.0f}%, p < 0.05 {100 * (p_ < 0.05).mean():.0f}%"
f" (positive and p < 0.05 {100 * ((r_ > 0) & (p_ < 0.05)).mean():.0f}%), median r {r_.median():.3f};"
f" chance for p < 0.05 is 5%")
rc = residual_correlations("main")
FIG5 = cells.reps.idxmax()
print(f"\nFigure 5 neuron (63 repetitions): {FIG5}")
print(rc[rc.neuron == FIG5].round(3).to_string(index=False))
print("Tones of", FIG5, "with 60 trials left:", rc[(rc.neuron == FIG5) & (rc.n == 60)].tone.tolist())
per_neuron = rc.groupby("neuron").r_scale.median()
print("Neurons whose median residual correlation (scale factor) is positive:", int((per_neuron > 0).sum()), "of", len(per_neuron))
[main] 396 tones of 29 neurons LFP scale : positive 81%, p < 0.05 32% (positive and p < 0.05 30%), median r 0.334; chance for p < 0.05 is 5% LFP lfp_height: positive 75%, p < 0.05 26% (positive and p < 0.05 23%), median r 0.265; chance for p < 0.05 is 5%
[spikes: fast rise above -20 mV (6)] 405 tones of 29 neurons LFP scale : positive 79%, p < 0.05 31% (positive and p < 0.05 29%), median r 0.309; chance for p < 0.05 is 5% LFP lfp_height: positive 75%, p < 0.05 24% (positive and p < 0.05 22%), median r 0.248; chance for p < 0.05 is 5%
[spikes: above -20 mV, any rise] 401 tones of 29 neurons LFP scale : positive 80%, p < 0.05 31% (positive and p < 0.05 29%), median r 0.327; chance for p < 0.05 is 5% LFP lfp_height: positive 75%, p < 0.05 24% (positive and p < 0.05 22%), median r 0.257; chance for p < 0.05 is 5%
[spikes: spikethresh, slow events kept] 401 tones of 29 neurons LFP scale : positive 79%, p < 0.05 32% (positive and p < 0.05 30%), median r 0.324; chance for p < 0.05 is 5% LFP lfp_height: positive 74%, p < 0.05 26% (positive and p < 0.05 23%), median r 0.256; chance for p < 0.05 is 5%
[peak in culling window (1b)] 394 tones of 29 neurons LFP scale : positive 80%, p < 0.05 32% (positive and p < 0.05 30%), median r 0.326; chance for p < 0.05 is 5% LFP lfp_height: positive 75%, p < 0.05 26% (positive and p < 0.05 24%), median r 0.265; chance for p < 0.05 is 5%
[heights re-measured after culling] 396 tones of 29 neurons LFP scale : positive 81%, p < 0.05 32% (positive and p < 0.05 30%), median r 0.334; chance for p < 0.05 is 5% LFP lfp_height: positive 75%, p < 0.05 26% (positive and p < 0.05 23%), median r 0.265; chance for p < 0.05 is 5%
Figure 5 neuron (63 repetitions): 020701md01b
neuron tone n r_scale p_scale r_lfp_height p_lfp_height
020701md01b 38134 62 0.555 0.000 0.448 0.000
020701md01b 10172 62 0.750 0.000 0.740 0.000
020701md01b 7498 61 0.470 0.000 0.464 0.000
020701md01b 16909 53 0.306 0.026 0.259 0.061
020701md01b 31118 60 0.699 0.000 0.686 0.000
020701md01b 15275 56 0.762 0.000 0.705 0.000
020701md01b 20722 54 0.283 0.038 0.158 0.254
020701md01b 28110 61 0.743 0.000 0.705 0.000
020701md01b 42214 62 0.388 0.002 0.416 0.001
020701md01b 6773 62 0.605 0.000 0.590 0.000
020701md01b 11260 62 0.566 0.000 0.577 0.000
020701md01b 9188 61 0.639 0.000 0.588 0.000
020701md01b 8300 55 0.609 0.000 0.418 0.001
020701md01b 25393 58 0.528 0.000 0.518 0.000
020701md01b 34448 61 0.553 0.000 0.483 0.000
020701md01b 22939 61 0.421 0.001 0.425 0.001
020701md01b 5527 60 0.651 0.000 0.587 0.000
Tones of 020701md01b with 60 trials left: [31118, 5527]
Neurons whose median residual correlation (scale factor) is positive: 26 of 29
Result: reproduced. Over 396 included tones of the 29 neurons, after LFP-based culling, PSP height correlates positively with the LFP scale factor in 81% of tones and with p < 0.05 in 32% (30% positive and significant; 5% expected by chance), above the criterion's half and quarter. With the LFP peak height in place of the scale factor, 75% and 26%. The Figure 5 neuron (020701md01b, 63 repetitions) has two tones with 60 trials left after culling (31,118 and 5,527 Hz), with r = 0.70 and 0.65 to the scale factor (p < 0.001). With the claims file's spike rule, 405 tones, 79% positive and 31% with p < 0.05; with the two further spike rules 401 tones, 80% and 31%, and 79% and 32%; with heights re-measured after culling the same as the main reading.
fig7a-lfp-culling-estimation¶
LFP-based culling and estimation together reduce the variability index by 40% on average, to 2.9 mV across 29 neurons, reduce it in 26 of 29 neurons, by as much as 71%, and have little effect when total variability is already small.
rows = []
for label, df in S.items():
d = df[df.lfp].dropna(subset=["q_lfp_est"])
red = 1 - d.q_lfp_est / d.q
rows.append({"reading": label, "n": len(d), "<q LFP-culled&est>": mse(d.q_lfp_est),
"reduced": int((red > 0).sum()), "largest %": 100 * red.max(),
"(a) 1 - est/total %": 100 * (1 - d.q_lfp_est.mean() / d.q.mean()),
"(b) mean of reductions %": 100 * red.mean()})
print(pd.DataFrame(rows).set_index("reading").round(1).to_string())
d = main[main.lfp]
show = pd.DataFrame({"total": d.q, "LFP-culled&est": d.q_lfp_est, "reduction %": 100 * (1 - d.q_lfp_est / d.q),
"uncorrected": d.q_lfp_est_uncorrected, "tones": d.n_tones_lfp})
print(show.sort_values("total").round(2).to_string())
print("\nNeurons with total <q>tones below 1.5 mV (shown, not scored):")
print(show[show.total < 1.5].round(2).to_string())
print("Largest reductions (main):", show["reduction %"].nlargest(3).round(1).to_dict())
# Per tone, the neuron with the largest reduction: total q over its included tones, and the
# values after LFP-based culling and estimation for the tones that keep 10 trials.
top = show["reduction %"].idxmax()
tt = tone_table("main", top)
inc = [r["tone"] for r in PT["main"][top] if ac2.included(r)]
kept = [r["tone"] for r in PT["main"][top] if ac2.included(r, "lfp")]
cols = ["n", "mean", "q", "n_lfp", "q_lfp", "q_lfp_est", "q_lfp_shuf"]
print(f"\n{top}: {len(inc)} included tones, {len(kept)} keep 10 trials after LFP-based culling")
print(tt.loc[inc, cols].assign(kept=[t in kept for t in inc]).round(3).to_string())
print(f"<q>tones over all included tones {tt.loc[inc, 'q'].mean():.3f}; over the kept tones {tt.loc[kept, 'q'].mean():.3f};"
f" after culling and estimation {show.loc[top, 'LFP-culled&est']:.3f} mV")
n <q LFP-culled&est> reduced largest % (a) 1 - est/total % (b) mean of reductions %
reading
main 29 2.69 +- 0.33 (n = 29) 26 80.0 43.9 40.1
spikes: fast rise above -20 mV (6) 29 3.22 +- 0.48 (n = 29) 28 80.0 40.5 40.5
spikes: above -20 mV, any rise 29 2.68 +- 0.31 (n = 29) 27 80.0 42.2 38.8
spikes: spikethresh, slow events kept 29 2.85 +- 0.37 (n = 29) 29 80.0 45.3 43.2
peak in culling window (1b) 29 2.53 +- 0.31 (n = 29) 27 80.0 42.8 38.5
peak-to-trough over trace (2b) 29 2.66 +- 0.30 (n = 29) 26 80.0 44.6 38.7
variance over N (3) 29 2.48 +- 0.30 (n = 29) 26 80.2 44.9 40.9
no negative-f rule in PSP culling (4b) 29 2.67 +- 0.32 (n = 29) 26 80.1 44.3 40.4
cullings intersected 29 2.68 +- 0.32 (n = 29) 26 80.1 44.1 40.2
heights re-measured after culling 29 2.62 +- 0.30 (n = 29) 26 80.5 45.3 40.6
criteria after culling (5b) 29 2.59 +- 0.29 (n = 29) 26 73.9 46.0 41.1
total LFP-culled&est reduction % uncorrected tones
neuron
022801md01b 0.88 0.70 20.58 0.64 19
021401md03 1.20 0.61 49.09 0.55 3
021001md01 1.21 0.74 39.17 0.67 3
030201md01b 1.23 1.13 8.61 1.02 4
030201md03b 1.57 1.21 22.94 1.12 25
020701md01b 1.70 0.69 59.61 0.67 17
030301md02b 2.20 0.59 72.95 0.53 1
022801md02b 2.48 2.41 3.01 2.14 23
102800md01 2.51 1.50 40.19 1.40 10
111200md01 3.76 2.13 43.39 1.96 27
102600md01 3.79 2.35 38.09 2.08 11
102800md02 4.19 4.86 -16.01 4.31 23
111500md01 4.33 1.36 68.69 1.12 15
102900md04 4.79 5.19 -8.25 4.54 9
021401md01b 4.87 3.17 34.95 2.86 12
110900md05 5.10 5.69 -11.48 5.23 2
110700md01 5.14 1.34 73.86 1.12 20
021401md02 5.19 3.09 40.47 2.69 1
030301md03b 5.40 1.08 80.02 0.98 12
102900md02 5.45 4.40 19.37 4.03 29
021401md05 5.91 1.89 68.01 1.67 7
102600md03 6.19 3.21 48.12 2.79 10
022301md01 6.30 1.79 71.54 1.69 18
030301md01b 7.15 2.24 68.66 2.06 7
111400md01 7.36 4.20 42.97 3.85 7
021901md01 8.15 6.57 19.36 5.54 29
102500md02 8.90 4.28 51.90 3.81 9
020901md01 10.94 5.14 52.99 4.72 11
102400md01 11.21 4.44 60.45 3.86 32
Neurons with total <q>tones below 1.5 mV (shown, not scored):
total LFP-culled&est reduction % uncorrected tones
neuron
021001md01 1.21 0.74 39.17 0.67 3
021401md03 1.20 0.61 49.09 0.55 3
022801md01b 0.88 0.70 20.58 0.64 19
030201md01b 1.23 1.13 8.61 1.02 4
Largest reductions (main): {'030301md03b': 80.0, '110700md01': 73.9, '030301md02b': 73.0}
030301md03b: 18 included tones, 12 keep 10 trials after LFP-based culling
n mean q n_lfp q_lfp q_lfp_est q_lfp_shuf kept
tone
4074 16 3.068 0.929 12 1.047 0.569 0.982 True
38134 16 7.411 1.500 13 1.214 0.827 1.071 True
10172 16 4.995 1.253 14 1.410 1.172 1.285 True
7498 16 7.520 0.814 16 0.814 0.699 0.782 True
3681 16 3.075 24.786 0 NaN NaN NaN False
31118 16 6.471 4.458 12 1.393 0.536 1.342 True
20722 16 4.199 0.879 9 1.197 1.079 1.091 False
28110 16 6.334 0.893 14 0.892 0.497 0.810 True
42214 16 4.610 2.920 15 2.959 2.784 2.716 True
6773 16 7.075 4.024 14 1.795 0.876 1.670 True
46731 16 10.018 1.414 15 1.426 1.333 1.327 True
11260 16 3.547 1.279 15 1.300 1.265 1.247 True
2000 16 3.752 1.985 4 0.265 0.013 0.194 False
9188 16 5.834 0.470 13 0.570 0.512 0.514 True
2214 16 4.132 35.795 1 NaN NaN NaN False
25393 16 5.024 7.335 5 0.812 0.746 0.648 False
6119 16 3.180 5.745 0 NaN NaN NaN False
22939 16 5.123 0.776 10 0.919 0.686 0.794 True
<q>tones over all included tones 5.403; over the kept tones 1.727; after culling and estimation 1.080 mV
Result: partially reproduced (resting on the spike-rule assumption and on detail 5).
<q LFP-culled&estimation> is 2.69 +- 0.33 mV over 29 neurons (2.9 +- 0.3), and the procedures
reduce tones in 26 of 29 neurons, as in the paper. The average reduction is 43.9% by reading
(a) and 40.1% by reading (b), so the 40% holds. The largest reduction is 80.0% (030301md03b),
outside 71 +- 7.1%. In that neuron LFP-based culling leaves 12 of its 18 included tones; four of the six
it drops have q of 5.7 to 35.8 mV, so
tones is 5.40 mV over all included tones
and 1.73 mV over the 12 kept, before estimation brings it to 1.08 mV. In the four neurons with
total
tones below 1.5 mV the procedures lower it by 8.6 to 49.1%; shown, not scored. With the
tone criteria applied after culling (5b) the value is 2.59 mV, 26 neurons are reduced and the
largest reduction is 73.9%, which would give reproduced; with the claims file's spike rule the
value is 3.22 mV, outside, which would give not reproduced. With the corrected potential above
-20 mV at any rise (2.68 mV, 27 reduced, largest 80.0%), with
spikethresh keeping the slow
events (2.85 mV, 29 reduced, 80.0%) and with heights re-measured after culling (2.62 mV, 26,
80.5%) the outcome is partially reproduced.
psp-culling-reduction¶
PSP-based culling alone reduces the mean variability index to 2.2 mV across 33 neurons, and leaves it at most about 3 mV in 32 of 33 neurons.
for label, df in S.items():
d = df.dropna(subset=["q_psp"])
print(f"{label:42s} <q PSP-culled> {mse(d.q_psp)}; at most 3.3 mV in {int((d.q_psp <= 3.3).sum())} of {len(d)};"
f" at most 3.0 mV in {int((d.q_psp <= 3.0).sum())}")
print("\nPer neuron (main):", main.q_psp.sort_values().round(2).to_dict())
main <q PSP-culled> 2.19 +- 0.31 (n = 33); at most 3.3 mV in 30 of 33; at most 3.0 mV in 30
spikes: fast rise above -20 mV (6) <q PSP-culled> 2.72 +- 0.35 (n = 33); at most 3.3 mV in 25 of 33; at most 3.0 mV in 22
spikes: above -20 mV, any rise <q PSP-culled> 2.42 +- 0.32 (n = 33); at most 3.3 mV in 29 of 33; at most 3.0 mV in 28
spikes: spikethresh, slow events kept <q PSP-culled> 2.54 +- 0.34 (n = 33); at most 3.3 mV in 28 of 33; at most 3.0 mV in 23
peak in culling window (1b) <q PSP-culled> 2.09 +- 0.30 (n = 33); at most 3.3 mV in 31 of 33; at most 3.0 mV in 31
peak-to-trough over trace (2b) <q PSP-culled> 2.22 +- 0.31 (n = 33); at most 3.3 mV in 30 of 33; at most 3.0 mV in 30
variance over N (3) <q PSP-culled> 2.02 +- 0.28 (n = 33); at most 3.3 mV in 32 of 33; at most 3.0 mV in 30
no negative-f rule in PSP culling (4b) <q PSP-culled> 2.69 +- 0.26 (n = 33); at most 3.3 mV in 24 of 33; at most 3.0 mV in 21
cullings intersected <q PSP-culled> 2.19 +- 0.31 (n = 33); at most 3.3 mV in 30 of 33; at most 3.0 mV in 30
heights re-measured after culling <q PSP-culled> 2.18 +- 0.29 (n = 33); at most 3.3 mV in 30 of 33; at most 3.0 mV in 30
criteria after culling (5b) <q PSP-culled> 1.92 +- 0.13 (n = 32); at most 3.3 mV in 32 of 32; at most 3.0 mV in 30
Per neuron (main): {'022801md01b': 0.72, '021001md01': 0.86, '021401md03': 0.87, '030301md03b': 0.95, '102800md01': 0.95, '020701md01b': 1.1, '030201md01b': 1.13, '030201md03b': 1.2, '100500md01': 1.31, '102600md03': 1.4, '100300md02': 1.42, '030301md02b': 1.47, '021401md01b': 1.59, '100300md01': 1.63, '022801md02b': 1.67, '102900md04': 1.8, '022301md01': 1.87, '111200md01': 1.92, '021901md01': 2.08, '030301md01b': 2.11, '020901md01': 2.15, '111500md01': 2.22, '110700md01': 2.39, '021401md05': 2.49, '111400md01': 2.57, '102900md02': 2.8, '102400md01': 2.83, '102600md01': 2.88, '100500md02b': 2.96, '102500md02': 2.97, '102800md02': 3.38, '110900md05': 3.6, '021401md02': 10.98}
Result: partially reproduced (resting on the spike-rule assumption and on details 1 and 5).
is 2.19 +- 0.31 mV over 33 neurons (2.2 +- 0.3), but 30 of 33 neurons are at or
below 3.3 mV, outside 32 +- 1; the three above are 3.38, 3.60 and 10.98 mV. The readings move
this across the tolerance: with the peak searched in the culling window (1b) 2.09 mV and 31
neurons, with the tone criteria after culling (5b) 1.92 mV and 32 of 32, both of which would
give reproduced; without the negative scale factor rule (4b) 2.69 mV and 24, and with the
claims file's spike rule 2.72 mV and 25, both of which would give not reproduced. With the
variance over N, the claims file's variant of detail 3, 2.02 mV and 32. With the corrected
potential above -20 mV at any rise 2.42 mV and 29, with
spikethresh keeping the slow events
2.54 mV and 28, and with heights re-measured after culling 2.18 mV and 30, partially reproduced
each time.
fig4c-psp-culling-can-increase¶
PSP-based culling reduces the tone-averaged variability for all but one neuron, but increases the variability for at least one tone in half (18/33) of the neurons; in the Figure 4 neuron, it raises q for a 9.2-kHz tone from 0.47 to 0.57 mV.
rows = []
for label, df in S.items():
pt = PT[label] if label in PT else PT["main"]
crit = SUMMARY_READINGS.get(label, {}).get("criteria", "before")
up = 0
for name, tones in pt.items():
if any(ac2.included(r, None, crit) and ac2.included(r, "psp", crit) and r["q_psp"] > r["q"] for r in tones):
up += 1
d = df.dropna(subset=["q_psp"])
rows.append({"reading": label, "<q> falls": int((d.q_psp < d.q).sum()), "of": len(d),
"neurons with a tone whose q rises": up})
print(pd.DataFrame(rows).set_index("reading").to_string())
print("Neurons whose <q>tones rises (main):", main[main.q_psp >= main.q][["q", "q_psp"]].round(2).to_dict("index"))
# Figure 4 neuron: one tone where PSP-based culling removes one trace and LFP-based culling removes
# that trace and three others (Figure 4B legend).
cand = []
for name in lfp_neurons:
for r in PT["main"][name]:
if "keep_lfp" not in r:
continue
rp, rl = ~r["keep_psp"], ~r["keep_lfp"]
if rp.sum() == 1 and rl.sum() == 4 and (rp & rl).sum() == 1:
cand.append({"neuron": name, "tone": r["tone"], "n": r["n"]})
cand = pd.DataFrame(cand)
print(f"\nTones fitting the Figure 4B legend: {len(cand)} in {cand.neuron.nunique() if len(cand) else 0} neurons",
cand.groupby("neuron").tone.apply(list).to_dict() if len(cand) else "")
pairs = []
for name, tones in PT["main"].items():
r = next(r for r in tones if r["tone"] == 9188)
if "q_psp" in r and "q" in r:
pairs.append({"neuron": name, "q_9188": r["q"], "q_9188_psp": r["q_psp"], "n": r["n"], "n_psp": r["n_psp"],
"mean": r["mean"], "fits_4B": name in set(cand.neuron) if len(cand) else False})
pairs = pd.DataFrame(pairs).set_index("neuron")
pairs["distance"] = np.hypot(pairs.q_9188 - 0.47, pairs.q_9188_psp - 0.57)
print("9,188-Hz tone, closest five to (0.47, 0.57):")
print(pairs.nsmallest(5, "distance").round(3).to_string())
print("9,188-Hz tone in the neurons fitting Figure 4B:")
print(pairs[pairs.fits_4B].round(3).to_string())
# The Figure 4 panels are labelled 030301md03b (31 kHz in 4A and 4B, 9.2 kHz in 4C). For that
# neuron: the rms errors of each 31,118-Hz trial against the culling thresholds (as Figure 4B),
# the traces (as Figure 4A and 4C) and the 9,188-Hz q under each reading.
FIG4 = "030301md03b"
fig, ax = plt.subplots(1, 4, figsize=(18, 3.6))
for j, tone in enumerate((31118, 9188)):
row = next(r for r in PT["main"][FIG4] if r["tone"] == tone)
wc = N[FIG4]["r"][row["trials"]].astype(float)
lf = N[FIG4]["lfp"][row["trials"]].astype(float)
_, rms_p, thr_p = ac2.scale_and_error(wc)
_, rms_l, thr_l = ac2.scale_and_error(lf)
out_p, out_l = ~row["keep_psp"], ~row["keep_lfp"]
print(f"{FIG4}, {tone} Hz: {row['n']} trials; PSP-based culling removes {int(out_p.sum())}, LFP-based culling"
f" {int(out_l.sum())}, both {int((out_p & out_l).sum())}; q {row['q']:.3f}, after PSP-based culling {row['q_psp']:.3f} mV")
print(" PSP rms / threshold:", np.round(rms_p, 2).tolist(), f"/ {thr_p:.2f}; LFP rms / threshold:",
np.round(rms_l, 3).tolist(), f"/ {thr_l:.3f}")
a = ax[2 * j]
a.plot(tms, wc[~out_p].T, color="k", lw=0.6)
a.plot(tms, wc[out_p].T, color="tab:red", lw=1)
a.plot(tms, wc.mean(0), color="tab:green", lw=2)
a.set(xlabel="ms from onset", ylabel="mV", title=f"{FIG4} {tone} Hz (red: PSP-culled)")
b = ax[2 * j + 1]
b.scatter(rms_l, rms_p, c=np.where(out_p, "tab:red", "k"), s=12)
b.axvline(thr_l, ls="--", color="k")
b.axhline(thr_p, ls="--", color="k")
b.set(xlabel="rms error of LFP (mV)", ylabel="rms error of PSP (mV)", title=f"{tone} Hz")
plt.tight_layout()
plt.show()
for label in PT:
row = next(r for r in PT[label][FIG4] if r["tone"] == 9188)
r31 = next(r for r in PT[label][FIG4] if r["tone"] == 31118)
p31, l31 = ~r31["keep_psp"], ~r31["keep_lfp"]
print(f"[{label}] {FIG4} 9,188 Hz: q {row['q']:.3f} -> {row.get('q_psp', np.nan):.3f} mV after PSP-based culling ({row['n']} -> {row.get('n_psp', 0)} trials);"
f" 31,118 Hz: PSP-based culling removes {int(p31.sum())}, LFP-based {int(l31.sum())}, both {int((p31 & l31).sum())}")
<q> falls of neurons with a tone whose q rises
reading
main 32 33 16
spikes: fast rise above -20 mV (6) 32 33 15
spikes: above -20 mV, any rise 32 33 17
spikes: spikethresh, slow events kept 32 33 16
peak in culling window (1b) 32 33 16
peak-to-trough over trace (2b) 32 33 16
variance over N (3) 32 33 16
no negative-f rule in PSP culling (4b) 32 33 20
cullings intersected 32 33 16
heights re-measured after culling 32 33 18
criteria after culling (5b) 32 32 15
Neurons whose <q>tones rises (main): {'021401md02': {'q': 5.19, 'q_psp': 10.98}}
Tones fitting the Figure 4B legend: 11 in 8 neurons {'021001md01': [6773], '021401md03': [5527], '021401md05': [6773], '022301md01': [15275], '030301md01b': [38134], '110700md01': [16909, 25393], '110900md05': [28110, 6773, 11260], '111400md01': [4993]}
9,188-Hz tone, closest five to (0.47, 0.57):
q_9188 q_9188_psp n n_psp mean fits_4B distance
neuron
030301md03b 0.470 0.502 16 15 5.834 False 0.068
100500md02b 0.549 0.549 14 14 3.627 False 0.082
021001md01 0.464 0.464 13 13 2.569 True 0.106
021401md03 0.407 0.407 10 10 5.322 True 0.174
102800md01 0.347 0.371 14 13 5.920 False 0.234
9,188-Hz tone in the neurons fitting Figure 4B:
q_9188 q_9188_psp n n_psp mean fits_4B distance
neuron
021001md01 0.464 0.464 13 13 2.569 True 0.106
021401md03 0.407 0.407 10 10 5.322 True 0.174
021401md05 2.337 2.337 14 14 18.838 True 2.570
022301md01 5.522 1.822 28 22 9.890 True 5.205
030301md01b 5.597 1.316 14 8 3.422 True 5.181
110700md01 2.452 2.074 11 10 15.870 True 2.488
110900md05 1.485 1.134 13 11 2.905 True 1.162
111400md01 7.984 2.468 12 10 11.722 True 7.750
030301md03b, 31118 Hz: 16 trials; PSP-based culling removes 2, LFP-based culling 4, both 2; q 4.458, after PSP-based culling 1.294 mV
PSP rms / threshold: [0.26, 0.57, 0.89, 0.61, 1.29, 1.98, 0.99, 1.37, 0.76, 0.91, 8.48, 1.17, 1.53, 1.15, 1.09, 0.76] / 2.29; LFP rms / threshold: [0.002, 0.002, 0.003, 0.003, 0.003, 0.005, 0.003, 0.002, 0.003, 0.004, 0.006, 0.004, 0.003, 0.004, 0.002, 0.006] / 0.004
030301md03b, 9188 Hz: 16 trials; PSP-based culling removes 1, LFP-based culling 3, both 1; q 0.470, after PSP-based culling 0.502 mV
PSP rms / threshold: [0.8, 0.67, 0.99, 1.19, 0.94, 0.86, 1.04, 0.58, 1.25, 0.98, 1.23, 0.72, 0.77, 11.27, 1.14, 0.91] / 1.92; LFP rms / threshold: [0.011, 0.004, 0.004, 0.005, 0.003, 0.003, 0.003, 0.004, 0.006, 0.003, 0.003, 0.004, 0.004, 0.009, 0.003, 0.01] / 0.007
[main] 030301md03b 9,188 Hz: q 0.470 -> 0.502 mV after PSP-based culling (16 -> 15 trials); 31,118 Hz: PSP-based culling removes 2, LFP-based 4, both 2 [spikes: fast rise above -20 mV (6)] 030301md03b 9,188 Hz: q 0.470 -> 0.502 mV after PSP-based culling (16 -> 15 trials); 31,118 Hz: PSP-based culling removes 2, LFP-based 4, both 2 [spikes: above -20 mV, any rise] 030301md03b 9,188 Hz: q 0.470 -> 0.502 mV after PSP-based culling (16 -> 15 trials); 31,118 Hz: PSP-based culling removes 2, LFP-based 4, both 2 [spikes: spikethresh, slow events kept] 030301md03b 9,188 Hz: q 0.470 -> 0.502 mV after PSP-based culling (16 -> 15 trials); 31,118 Hz: PSP-based culling removes 2, LFP-based 4, both 2 [peak in culling window (1b)] 030301md03b 9,188 Hz: q 0.470 -> 0.502 mV after PSP-based culling (16 -> 15 trials); 31,118 Hz: PSP-based culling removes 2, LFP-based 4, both 2 [peak-to-trough over trace (2b)] 030301md03b 9,188 Hz: q 0.470 -> 0.502 mV after PSP-based culling (16 -> 15 trials); 31,118 Hz: PSP-based culling removes 2, LFP-based 4, both 2 [variance over N (3)] 030301md03b 9,188 Hz: q 0.440 -> 0.469 mV after PSP-based culling (16 -> 15 trials); 31,118 Hz: PSP-based culling removes 2, LFP-based 4, both 2 [no negative-f rule in PSP culling (4b)] 030301md03b 9,188 Hz: q 0.470 -> 0.502 mV after PSP-based culling (16 -> 15 trials); 31,118 Hz: PSP-based culling removes 1, LFP-based 4, both 1 [cullings intersected] 030301md03b 9,188 Hz: q 0.470 -> 0.502 mV after PSP-based culling (16 -> 15 trials); 31,118 Hz: PSP-based culling removes 2, LFP-based 4, both 2 [heights re-measured after culling] 030301md03b 9,188 Hz: q 0.470 -> 0.502 mV after PSP-based culling (16 -> 15 trials); 31,118 Hz: PSP-based culling removes 2, LFP-based 4, both 2
Result: partially reproduced (resting on the spike-rule assumption and on measuring
heights once before culling). PSP-based culling lowers tones in 32 of 33 neurons (the
exception is 021401md02, whose one remaining tone rises from 5.19 to 10.98 mV), as the paper
says. q rises for at least one tone in 16 of 33 neurons, outside 18 +- 1.8. Readings 1b, 2b, 3
and the intersected cullings give 16, as does
spikethresh keeping the slow events; the claims file's spike rule and the criteria after
culling (5b) 15, and without the negative scale factor rule in PSP culling (4b) 20, all outside.
With heights re-measured after culling the count is 18, and with the corrected potential above
-20 mV at any rise 17; both are within tolerance, with 32 of 33 neurons falling, and would give
reproduced.
The Figure 4 neuron. The search the claims file sets finds 11 tones in 8 neurons that fit the Figure 4B legend, and in none of them does the 9,188-Hz tone give 0.47 and 0.57 mV; the closest 9,188-Hz pair in any neuron is 0.470 to 0.502 mV, in 030301md03b. The panels of Figure 4 are labelled 030301md03b, so the paper names this neuron (asked of a reader, as the claims file reads it as not named). Its 9,188-Hz tone gives 0.470 mV before PSP-based culling, matching 0.47, and 0.502 mV after it, below 0.57 +- 0.057; the culled trial is the one aberrant trace that peaks about 10 ms after onset, as in Figure 4C. At 31,118 Hz (Figure 4A and 4B) PSP-based culling removes two of 16 trials: the trace that falls far below baseline, the paper's aberrant trial, and one with a negative scale factor. LFP-based culling removes four, including both. Without the negative scale factor rule in PSP culling (4b) it removes only the aberrant trace, and LFP-based culling removes that trace and three others, which is what the Figure 4B legend says. That fits reading 4b, while the Methods cite trials 3 and 8 of Figure 5C, a PSP-based culling panel, as trials with negative scale factors, which fits reading 4a. Reading 4a is kept.
fig7b-private-best-estimate¶
PSP-based culling followed by LFP-based culling and estimation reduces variability by as much as 85%, to a best estimate of private variability of 1.6 mV across 28 neurons, comparable and low for all neurons.
rows = []
for label, df in S.items():
d = df.dropna(subset=["q_both_est"])
red = 1 - d.q_both_est / d.q
rows.append({"reading": label, "<q best>": mse(d.q_both_est), "neurons": len(d), "largest reduction %": 100 * red.max(),
"min": d.q_both_est.min(), "max": d.q_both_est.max(), "SD": d.q_both_est.std(ddof=1)})
print(pd.DataFrame(rows).set_index("reading").round(2).to_string())
d = main.dropna(subset=["q_both_est"])
print("\nNeuron with an LFP but no estimate (main):", list(set(lfp_neurons) - set(d.index)))
print("Best estimate per neuron (main):", d.q_both_est.sort_values().round(2).to_dict())
print("Tones kept per neuron:", d.n_tones_both.to_dict())
print("Best estimate within 0.9-3.3 mV:", int(d.q_both_est.between(0.9, 3.3).sum()), "of", len(d),
"; total <q>tones SD across these neurons %.2f, best estimate SD %.2f" % (d.q.std(ddof=1), d.q_both_est.std(ddof=1)))
<q best> neurons largest reduction % min max SD
reading
main 1.48 +- 0.17 (n = 28) 28 84.59 0.31 4.52 0.88
spikes: fast rise above -20 mV (6) 2.01 +- 0.31 (n = 28) 28 84.59 0.31 6.96 1.64
spikes: above -20 mV, any rise 1.56 +- 0.17 (n = 28) 28 84.59 0.31 4.52 0.90
spikes: spikethresh, slow events kept 1.72 +- 0.24 (n = 28) 28 84.59 0.31 4.97 1.25
peak in culling window (1b) 1.42 +- 0.16 (n = 28) 28 84.59 0.31 4.52 0.83
peak-to-trough over trace (2b) 1.57 +- 0.17 (n = 28) 28 84.51 0.32 4.52 0.88
variance over N (3) 1.36 +- 0.15 (n = 28) 28 84.81 0.28 4.14 0.81
no negative-f rule in PSP culling (4b) 1.78 +- 0.22 (n = 28) 28 83.72 0.30 5.92 1.19
cullings intersected 1.51 +- 0.21 (n = 28) 28 87.78 0.30 5.94 1.09
heights re-measured after culling 1.51 +- 0.17 (n = 28) 28 84.88 0.31 4.64 0.92
criteria after culling (5b) 1.47 +- 0.14 (n = 28) 28 83.45 0.31 3.34 0.75
Neuron with an LFP but no estimate (main): ['021401md02']
Best estimate per neuron (main): {'021401md03': 0.31, '020701md01b': 0.54, '030301md02b': 0.6, '022801md01b': 0.63, '102800md01': 0.68, '111500md01': 0.77, '030301md03b': 0.83, '030201md01b': 0.93, '021001md01': 0.95, '030201md03b': 1.09, '022301md01': 1.14, '021401md01b': 1.14, '102900md04': 1.18, '110700md01': 1.29, '021901md01': 1.35, '111200md01': 1.38, '030301md01b': 1.41, '022801md02b': 1.43, '102600md03': 1.45, '102500md02': 1.7, '021401md05': 1.86, '102400md01': 1.95, '020901md01': 2.02, '102600md01': 2.2, '111400md01': 2.54, '102900md02': 2.61, '102800md02': 2.8, '110900md05': 4.52}
Tones kept per neuron: {'020701md01b': 17, '020901md01': 6, '021001md01': 4, '021401md01b': 7, '021401md03': 1, '021401md05': 7, '021901md01': 14, '022301md01': 18, '022801md01b': 19, '022801md02b': 21, '030201md01b': 5, '030201md03b': 26, '030301md01b': 3, '030301md02b': 1, '030301md03b': 12, '102400md01': 29, '102500md02': 6, '102600md01': 10, '102600md03': 4, '102800md01': 4, '102800md02': 20, '102900md02': 25, '102900md04': 3, '110700md01': 21, '110900md05': 3, '111200md01': 25, '111400md01': 6, '111500md01': 11}
Best estimate within 0.9-3.3 mV: 20 of 28 ; total <q>tones SD across these neurons 2.85, best estimate SD 0.88
Result: reproduced (resting on the spike-rule assumption). The best estimate of private
variability, PSP-based culling followed by LFP-based culling on the trials it keeps and
estimation, is 1.48 +- 0.17 mV over 28 neurons (1.6 +- 0.2; 28 exact), and the largest
reduction from total variability is 84.6% (85 +- 8.5). Deciding both cullings on the full
ensemble gives 1.51 mV and 87.8%. The readings of details 1, 2, 4 and 5 give 1.42 to 1.78 mV and
heights re-measured after culling 1.51 mV, all within tolerance; the variance over N, the claims
file's variant of detail 3, gives 1.36 mV, just outside. With the claims file's spike rule the
value is 2.01 mV, outside; with the corrected potential above -20 mV at any rise 1.56 mV and with
spikethresh keeping the slow events 1.72 mV, both within. The values are low for most neurons
(0.31 to 2.80 mV) except 110900md05 (4.52 mV, three tones); their SD across neurons is 0.88 mV
against 2.85 mV for total variability.
culling-controls¶
Random culling of as many traces as LFP-based culling removes leaves the variability index unchanged (4.5 against 4.6 mV), while matching the PSP-based or both cullings lowers it (3.67 and 2.98 mV) because tones with fewer than 10 non-aberrant trials drop out; keeping only those tones without culling gives 2.9 mV. PSP-based culling removes 29% of trials and LFP-based culling 38%.
def removed_mask(r, kind):
if kind == "psp":
return ~r["keep_psp"]
if "keep_lfp" not in r:
return None
return ~r["keep_lfp"] if kind == "lfp" else ~r["keep_both"]
def random_control(pt, kind, names, seed, min_mean=ac2.MIN_MEAN, n_sims=20, ddof=1, min_left=ac2.MIN_TRIALS):
"""Remove at random, per neuron and tone, as many trials as the named culling removed, from the
tones included before culling; keep the tones with at least `min_left` trials left (10: the tone
criteria applied again; 2: not applied again) and take <q>neurons."""
rng = np.random.default_rng(seed)
sims, dropped = [], []
for _ in range(n_sims):
per_neuron, n_drop = [], 0
for name in names:
qs = []
for r in pt[name]:
if not ac2.included(r, None, min_mean=min_mean):
continue
rem = removed_mask(r, kind)
k = int(rem.sum()) if rem is not None else 0
keep = rng.permutation(r["n"])[k:]
if len(keep) >= max(min_left, 2):
qs.append(ac2.q_index(r["h"][keep], ddof))
else:
n_drop += 1
if qs:
per_neuron.append(np.mean(qs))
sims.append((np.mean(per_neuron), len(per_neuron)))
dropped.append(n_drop)
v = np.array([s[0] for s in sims])
return {"mean": v.mean(), "SD over sims": v.std(ddof=1), "SE over sims": v.std(ddof=1) / np.sqrt(n_sims),
"neurons": int(np.median([s[1] for s in sims])), "tones dropped": float(np.mean(dropped))}
all_names, lfp_names = list(N), list(lfp_neurons)
def controls(label, min_mean=ac2.MIN_MEAN, min_left=ac2.MIN_TRIALS):
pt = PT[label]
out = {
"LFP-matched, 29 neurons": random_control(pt, "lfp", lfp_names, 11, min_mean, min_left=min_left),
"LFP-matched, 33 (4 unculled)": random_control(pt, "lfp", all_names, 12, min_mean, min_left=min_left),
"PSP-matched, 33 neurons": random_control(pt, "psp", all_names, 13, min_mean, min_left=min_left),
"PSP-matched, 29 neurons": random_control(pt, "psp", lfp_names, 14, min_mean, min_left=min_left),
"both-matched, 29 neurons": random_control(pt, "both", lfp_names, 15, min_mean, min_left=min_left),
}
return pd.DataFrame(out).T
for label in ("main", SPK, SPK_ANY, SPK_SLOW, "cullings intersected"):
print(f"[{label}] tone criteria applied again after random culling (claims file)")
print(controls(label).round(3).to_string())
df = S[label]
print(f"Second control: unculled <q>tones over tones that keep 10 trials after both cullings: {mse(df.q_on_both_tones)}")
print("[main] tone criteria not applied again (tones included before culling all kept, if 2 trials remain)")
print(controls("main", min_left=2).round(3).to_string())
print(f"Included tones before culling: {sum(ac2.included(r) for name in all_names for r in PT['main'][name])} (33 neurons),"
f" {sum(ac2.included(r) for name in lfp_names for r in PT['main'][name])} (29)")
# Tones per neuron that the LFP-matched random control drops (fewer than 10 left), which are the
# tones that LFP-based culling itself drops, as it removes the same number; <q>tones over the rest.
drop = []
for name in lfp_names:
inc = [r for r in PT["main"][name] if ac2.included(r)]
left = [r for r in inc if r["n"] - int((~r["keep_lfp"]).sum()) >= ac2.MIN_TRIALS]
drop.append({"neuron": name, "included": len(inc), "dropped": len(inc) - len(left),
"<q> all included": np.mean([r["q"] for r in inc]),
"<q> tones left": np.mean([r["q"] for r in left]) if left else np.nan})
drop = pd.DataFrame(drop).set_index("neuron")
print("\nLFP-matched control, tones dropped per neuron (main):")
print(drop.round(2).T.to_string())
print(f"Neurons losing every tone: {int((drop.dropped == drop.included).sum())}; mean of <q> over the tones left {drop['<q> tones left'].mean():.2f}"
f" (n = {drop['<q> tones left'].notna().sum()}) against {drop['<q> all included'].mean():.2f} over all included tones")
def fraction_removed(pt, kind, names):
pooled_rem = pooled_all = 0
per = []
for name in names:
rem_n = all_n = 0
for r in pt[name]:
if not ac2.included(r):
continue
rem = removed_mask(r, kind)
rem_n += int(rem.sum())
all_n += r["n"]
pooled_rem += rem_n
pooled_all += all_n
if all_n:
per.append(rem_n / all_n)
return 100 * pooled_rem / pooled_all, 100 * np.mean(per)
for label in ("main", SPK, SPK_ANY, SPK_SLOW, "peak-to-trough over trace (2b)", "no negative-f rule in PSP culling (4b)"):
pt = PT[label]
f_psp33 = fraction_removed(pt, "psp", all_names)
f_psp29 = fraction_removed(pt, "psp", lfp_names)
f_lfp = fraction_removed(pt, "lfp", lfp_names)
f_both = fraction_removed(pt, "both", lfp_names)
print(f"[{label}] trials removed, pooled / mean of neurons: PSP-based {f_psp33[0]:.1f}% / {f_psp33[1]:.1f}% (33 neurons),"
f" {f_psp29[0]:.1f}% / {f_psp29[1]:.1f}% (29); LFP-based {f_lfp[0]:.1f}% / {f_lfp[1]:.1f}% (29); both {f_both[0]:.1f}% / {f_both[1]:.1f}%")
[main] tone criteria applied again after random culling (claims file)
mean SD over sims SE over sims neurons tones dropped LFP-matched, 29 neurons 3.819 0.106 0.024 29.0 247.0 LFP-matched, 33 (4 unculled) 3.801 0.138 0.031 33.0 247.0 PSP-matched, 33 neurons 3.737 0.234 0.052 33.0 157.0 PSP-matched, 29 neurons 3.861 0.128 0.029 29.0 139.0 both-matched, 29 neurons 2.997 0.097 0.022 28.0 315.0 Second control: unculled <q>tones over tones that keep 10 trials after both cullings: 2.96 +- 0.32 (n = 28) [spikes: fast rise above -20 mV (6)] tone criteria applied again after random culling (claims file)
mean SD over sims SE over sims neurons tones dropped LFP-matched, 29 neurons 4.361 0.192 0.043 29.0 248.0 LFP-matched, 33 (4 unculled) 4.485 0.172 0.039 33.0 248.0 PSP-matched, 33 neurons 3.743 3.423 0.766 33.0 156.0 PSP-matched, 29 neurons 4.831 1.106 0.247 29.0 135.0 both-matched, 29 neurons 3.984 0.136 0.030 28.0 310.0 Second control: unculled <q>tones over tones that keep 10 trials after both cullings: 3.90 +- 0.57 (n = 28) [spikes: above -20 mV, any rise] tone criteria applied again after random culling (claims file)
mean SD over sims SE over sims neurons tones dropped LFP-matched, 29 neurons 3.684 0.076 0.017 29.0 246.0 LFP-matched, 33 (4 unculled) 3.601 0.073 0.016 33.0 246.0 PSP-matched, 33 neurons 3.919 0.367 0.082 33.0 148.0 PSP-matched, 29 neurons 3.988 0.112 0.025 29.0 134.0 both-matched, 29 neurons 3.097 0.073 0.016 28.0 309.0 Second control: unculled <q>tones over tones that keep 10 trials after both cullings: 3.07 +- 0.35 (n = 28) [spikes: spikethresh, slow events kept] tone criteria applied again after random culling (claims file)
mean SD over sims SE over sims neurons tones dropped LFP-matched, 29 neurons 4.110 0.120 0.027 29.0 248.0 LFP-matched, 33 (4 unculled) 4.137 0.108 0.024 33.0 248.0 PSP-matched, 33 neurons 5.088 3.470 0.776 33.0 159.0 PSP-matched, 29 neurons 4.336 0.345 0.077 29.0 138.0 both-matched, 29 neurons 3.468 0.091 0.020 28.0 313.0 Second control: unculled <q>tones over tones that keep 10 trials after both cullings: 3.41 +- 0.40 (n = 28) [cullings intersected] tone criteria applied again after random culling (claims file)
mean SD over sims SE over sims neurons tones dropped LFP-matched, 29 neurons 3.819 0.106 0.024 29.0 247.0 LFP-matched, 33 (4 unculled) 3.801 0.138 0.031 33.0 247.0 PSP-matched, 33 neurons 3.737 0.234 0.052 33.0 157.0 PSP-matched, 29 neurons 3.861 0.128 0.029 29.0 139.0 both-matched, 29 neurons 2.988 0.079 0.018 28.0 327.0 Second control: unculled <q>tones over tones that keep 10 trials after both cullings: 2.95 +- 0.32 (n = 28) [main] tone criteria not applied again (tones included before culling all kept, if 2 trials remain)
mean SD over sims SE over sims neurons tones dropped LFP-matched, 29 neurons 4.832 0.965 0.216 29.0 59.0 LFP-matched, 33 (4 unculled) 4.722 1.277 0.286 33.0 59.0 PSP-matched, 33 neurons 2.298 8.522 1.906 33.0 15.0 PSP-matched, 29 neurons 4.276 1.068 0.239 29.0 15.0 both-matched, 29 neurons 3.586 4.679 1.046 29.0 42.0 Included tones before culling: 737 (33 neurons), 643 (29) LFP-matched control, tones dropped per neuron (main): neuron 020701md01b 020901md01 021001md01 021401md01b 021401md02 021401md03 021401md05 021901md01 022301md01 022801md01b 022801md02b 030201md01b 030201md03b 030301md01b 030301md02b 030301md03b 102400md01 102500md02 102600md01 102600md03 102800md01 102800md02 102900md02 102900md04 110700md01 110900md05 111200md01 111400md01 111500md01 included 18.0 27.00 6.00 17.00 4.00 14.00 25.00 32.00 29.00 26.00 29.00 6.00 32.00 20.00 16.00 18.00 32.00 30.0 16.00 12.00 22.00 29.00 30.00 25.00 31.00 12.00 32.00 25.00 28.00 dropped 1.0 16.00 3.00 5.00 3.00 11.00 18.00 3.00 11.00 7.00 6.00 2.00 7.00 13.00 15.00 6.00 0.00 21.0 5.00 2.00 12.00 6.00 1.00 16.00 11.00 10.00 5.00 18.00 13.00 <q> all included 1.7 10.94 1.21 4.87 5.19 1.20 5.91 8.15 6.30 0.88 2.48 1.23 1.57 7.15 2.20 5.40 11.21 8.9 3.79 6.19 2.51 4.19 5.45 4.79 5.14 5.10 3.76 7.36 4.33 <q> tones left 1.3 6.80 1.21 4.00 3.13 0.69 2.52 7.46 3.29 0.73 2.40 1.31 1.35 3.12 2.79 1.73 11.21 7.1 3.04 5.41 2.12 4.28 5.52 5.66 3.86 6.19 3.82 3.83 2.11 Neurons losing every tone: 0; mean of <q> over the tones left 3.72 (n = 29) against 4.80 over all included tones [main] trials removed, pooled / mean of neurons: PSP-based 24.1% / 19.2% (33 neurons), 25.7% / 20.9% (29); LFP-based 37.8% / 33.3% (29); both 44.6% / 40.0% [spikes: fast rise above -20 mV (6)] trials removed, pooled / mean of neurons: PSP-based 23.8% / 19.2% (33 neurons), 25.6% / 20.8% (29); LFP-based 37.9% / 33.4% (29); both 44.5% / 39.9% [spikes: above -20 mV, any rise] trials removed, pooled / mean of neurons: PSP-based 23.6% / 18.7% (33 neurons), 25.1% / 20.5% (29); LFP-based 37.6% / 33.2% (29); both 44.2% / 39.7% [spikes: spikethresh, slow events kept] trials removed, pooled / mean of neurons: PSP-based 23.8% / 19.3% (33 neurons), 25.6% / 20.9% (29); LFP-based 37.8% / 33.4% (29); both 44.5% / 40.1% [peak-to-trough over trace (2b)] trials removed, pooled / mean of neurons: PSP-based 22.5% / 17.7% (33 neurons), 24.1% / 19.3% (29); LFP-based 36.5% / 31.5% (29); both 42.6% / 37.9% [no negative-f rule in PSP culling (4b)] trials removed, pooled / mean of neurons: PSP-based 20.1% / 15.3% (33 neurons), 21.6% / 16.7% (29); LFP-based 37.8% / 33.3% (29); both 43.7% / 38.8%
Result: not reproduced. Two of the four q values and one of the two fractions match. Random
culling matched to both cullings gives 2.997 mV (2.98 +- 0.02) and the second control 2.96 +-
0.32 mV over 28 neurons (2.9 +- 0.3). Random culling matched to LFP-based culling gives 3.82 mV
over the 29 LFP neurons (3.80 over 33), not the 4.5 +- 0.1 mV of the paper. It removes as many
trials per tone as LFP-based culling, so it drops the same tones: 247 of the 643 included tones
of the 29 neurons fall below 10 trials, and tones over the tones left averages 3.72 mV
against 4.80 over all included tones. Without the tone criteria applied again after random
culling, a reading that may be the paper's (asked of a reader), the LFP-matched control gives
4.83 mV (SD over simulations 0.97), within reach of 4.5, while the PSP-matched one becomes
unstable (2.30 mV, SD 8.5, from tones left with few trials). Matched to PSP-based
culling, with the criteria applied again, it gives 3.737 mV over 33 neurons (3.67 +- 0.02; the
miss, 0.067, is more than three SEs) and 3.86 over 29. PSP-based culling removes 24.1% of the
trials of included tones (pooled; 19.2% as a mean of neurons), outside 29 +- 2.9%, and LFP-based
culling 37.8% (38 +- 3.8). With the claims file's spike rule the controls give 4.36, 3.74 and
3.98 mV, the corrected potential above -20 mV at any rise 3.68, 3.92 and 3.10 mV, and
spikethresh keeping the slow events 4.11, 5.09 and 3.47 mV; the outcome is the same. The
controls use the heights measured on the full ensemble under every reading.
fig7d-large-responses¶
For tones with mean PSPs of at least 15 mV, total variability is lower (2.9 mV, 19 neurons) than for all tones, but the private variability is the same as for all tones (best estimate 1.5 mV, LFP-based bound 1.9 mV, 16 neurons).
LARGE = ac2.LARGE_MEAN - 1e-9 # "at least 15 mV": the criterion is "> min_mean"
SL = {}
for label in READINGS:
SL[label] = ac2.summarise(PT[label], N, min_mean=LARGE)
SL["criteria after culling (5b)"] = ac2.summarise(PT["main"], N, criteria="after", min_mean=LARGE)
rows = []
for label, df in SL.items():
tot = df.q.dropna()
rows.append({"reading": label, "<q> total": mse(tot), "<q LFP-culled&est>": mse(df.q_lfp_est[df.lfp]),
"<q best>": mse(df.q_both_est)})
print(pd.DataFrame(rows).set_index("reading").to_string())
dl = SL["main"]
print("\nNeurons with a tone of at least 15 mV (main):", int(dl.q.notna().sum()), list(dl.index[dl.q.notna()]))
print("Neurons with a best estimate on such tones:", list(dl.index[dl.q_both_est.notna()]))
print("Tones of at least 15 mV per neuron:", dl.n_tones[dl.n_tones > 0].to_dict())
<q> total <q LFP-culled&est> <q best>
reading
main 2.92 +- 0.36 (n = 19) 1.62 +- 0.23 (n = 14) 1.43 +- 0.20 (n = 14)
spikes: fast rise above -20 mV (6) 6.13 +- 1.26 (n = 21) 4.45 +- 1.23 (n = 15) 3.75 +- 0.97 (n = 15)
spikes: above -20 mV, any rise 3.00 +- 0.35 (n = 19) 1.99 +- 0.27 (n = 14) 1.71 +- 0.23 (n = 14)
spikes: spikethresh, slow events kept 5.50 +- 1.29 (n = 21) 3.22 +- 0.98 (n = 15) 2.64 +- 0.83 (n = 15)
peak in culling window (1b) 2.76 +- 0.33 (n = 19) 1.57 +- 0.22 (n = 14) 1.39 +- 0.20 (n = 14)
peak-to-trough over trace (2b) 2.92 +- 0.36 (n = 19) 1.63 +- 0.23 (n = 14) 1.47 +- 0.21 (n = 14)
variance over N (3) 2.72 +- 0.35 (n = 19) 1.50 +- 0.22 (n = 14) 1.32 +- 0.18 (n = 14)
no negative-f rule in PSP culling (4b) 2.92 +- 0.36 (n = 19) 1.61 +- 0.23 (n = 14) 1.62 +- 0.23 (n = 14)
cullings intersected 2.92 +- 0.36 (n = 19) 1.62 +- 0.23 (n = 14) 1.38 +- 0.21 (n = 14)
heights re-measured after culling 2.92 +- 0.36 (n = 19) 1.65 +- 0.25 (n = 14) 1.45 +- 0.20 (n = 14)
criteria after culling (5b) 2.92 +- 0.36 (n = 19) 1.81 +- 0.29 (n = 15) 1.52 +- 0.19 (n = 16)
Neurons with a tone of at least 15 mV (main): 19 ['020901md01', '021401md05', '022301md01', '022801md01b', '022801md02b', '030201md03b', '030301md01b', '030301md02b', '100300md02', '100500md01', '100500md02b', '102400md01', '102500md02', '102800md02', '102900md02', '110700md01', '111200md01', '111400md01', '111500md01']
Neurons with a best estimate on such tones: ['021401md05', '022301md01', '022801md01b', '022801md02b', '030201md03b', '030301md01b', '102400md01', '102500md02', '102800md02', '102900md02', '110700md01', '111200md01', '111400md01', '111500md01']
Tones of at least 15 mV per neuron: {'020901md01': 2, '021401md05': 10, '022301md01': 4, '022801md01b': 3, '022801md02b': 16, '030201md03b': 3, '030301md01b': 3, '030301md02b': 4, '100300md02': 10, '100500md01': 4, '100500md02b': 1, '102400md01': 19, '102500md02': 3, '102800md02': 3, '102900md02': 7, '110700md01': 10, '111200md01': 13, '111400md01': 9, '111500md01': 23}
Result: reproduced (resting on the spike-rule assumption and on details 1, 3 and 4). For
tones with a mean PSP of at least 15 mV, total is 2.92 +- 0.36 mV over 19 neurons (2.9 +-
0.4; 19 exact), the LFP-based bound 1.62 +- 0.23 mV (1.9 +- 0.3) and the best estimate 1.43 +-
0.20 mV (1.5 +- 0.2). The LFP-based bound is within tolerance by 0.02 mV. The peak searched in
the culling window (1b) gives 1.57 mV and the variance over N 1.50 mV, both below 1.6, which
would not give reproduced; without the negative scale factor rule in PSP culling (4b) it is 1.61
mV, within by 0.01, and with heights re-measured after culling 1.65 mV. Both private values come
from 14 neurons, not the paper's 16; with the 3 mV and 15 mV criteria applied after culling
(detail 5b) there are 15 and 16, with values 1.81 and 1.52 mV, also within. With the claims
file's spike rule total
is 6.13 mV over 21 neurons, which would give not reproduced, as
would
spikethresh keeping the slow events (5.50 mV over 21). With the corrected potential above
-20 mV at any rise the values are 3.00, 1.99 and 1.71 mV, the last 0.01 mV outside.
large-responses-controls¶
For tones of at least 15 mV, random culling and keeping only tones that survive both cullings give 3.2 mV, so the private estimate of 1.5 mV is not an undersampling effect.
for label in ("main", SPK, SPK_ANY, SPK_SLOW):
c = random_control(PT[label], "both", lfp_names, 16, LARGE)
print(f"[{label}] random culling matched to both, tones of at least 15 mV: {c['mean']:.3f} mV"
f" (SE over 20 simulations {c['SE over sims']:.3f}, SD {c['SD over sims']:.3f}, {c['neurons']} neurons)")
print(f"[{label}] second control, unculled <q> over tones of at least 15 mV that keep 10 trials after both cullings:"
f" {mse(SL[label].q_on_both_tones)}")
[main] random culling matched to both, tones of at least 15 mV: 2.422 mV (SE over 20 simulations 0.022, SD 0.098, 14 neurons) [main] second control, unculled <q> over tones of at least 15 mV that keep 10 trials after both cullings: 2.40 +- 0.39 (n = 14) [spikes: fast rise above -20 mV (6)] random culling matched to both, tones of at least 15 mV: 4.910 mV (SE over 20 simulations 0.066, SD 0.297, 15 neurons) [spikes: fast rise above -20 mV (6)] second control, unculled <q> over tones of at least 15 mV that keep 10 trials after both cullings: 4.96 +- 1.08 (n = 15)
[spikes: above -20 mV, any rise] random culling matched to both, tones of at least 15 mV: 2.485 mV (SE over 20 simulations 0.019, SD 0.087, 14 neurons) [spikes: above -20 mV, any rise] second control, unculled <q> over tones of at least 15 mV that keep 10 trials after both cullings: 2.51 +- 0.38 (n = 14)
[spikes: spikethresh, slow events kept] random culling matched to both, tones of at least 15 mV: 3.962 mV (SE over 20 simulations 0.030, SD 0.133, 15 neurons) [spikes: spikethresh, slow events kept] second control, unculled <q> over tones of at least 15 mV that keep 10 trials after both cullings: 3.88 +- 1.03 (n = 15)
Result: not reproduced. For tones of at least 15 mV, random culling matched to both cullings
gives 2.42 mV (SE over 20 simulations 0.02) and the second control 2.40 +- 0.39 mV over 14
neurons, against 3.20 +- 0.05 and 3.2 +- 0.5. The two controls agree with each other and are
above the 1.43 mV best estimate, as the paper argues, but both values are outside tolerance
(with the claims file's spike rule, 4.91 and 4.96 mV, with the corrected potential above -20 mV
at any rise 2.49 and 2.51 mV, and with spikethresh keeping the slow events 3.96 and 3.88 mV,
all outside).
cv-large-responses¶
In coefficient of variation, both the total and the private variability of the large (at least 15 mV) responses are lower than for the full data set (total 0.37 against 0.73; private 0.27 against 0.41, significantly lower), unlike the variability index.
rows = []
for label in S:
a, l = S[label], SL[label]
pa, pl = a.cv_both_est.dropna(), l.cv_both_est.dropna()
common = pl.index.intersection(pa.index)
rows.append({"reading": label, "total >= 15 mV": mse(l.cv, 3), "total all": mse(a.cv, 3),
"private >= 15 mV": mse(pl, 3), "private all": mse(pa, 3),
"(a) Student's t p": stats.ttest_ind(pl, pa).pvalue,
"(b) paired t p": stats.ttest_rel(pl[common], pa[common]).pvalue, "paired n": len(common),
"mean paired difference": np.mean(pl[common] - pa[common])})
print(pd.DataFrame(rows).set_index("reading").round(4).to_string())
total >= 15 mV total all private >= 15 mV private all (a) Student's t p (b) paired t p paired n mean paired difference reading main 0.371 +- 0.027 (n = 19) 0.734 +- 0.042 (n = 33) 0.265 +- 0.021 (n = 14) 0.383 +- 0.026 (n = 28) 0.0055 0.0000 14 -0.0808 spikes: fast rise above -20 mV (6) 0.491 +- 0.052 (n = 21) 0.757 +- 0.043 (n = 33) 0.374 +- 0.047 (n = 15) 0.412 +- 0.029 (n = 28) 0.4751 0.3795 15 -0.0195 spikes: above -20 mV, any rise 0.376 +- 0.027 (n = 19) 0.721 +- 0.042 (n = 33) 0.286 +- 0.021 (n = 14) 0.389 +- 0.026 (n = 28) 0.0144 0.0012 14 -0.0687 spikes: spikethresh, slow events kept 0.468 +- 0.054 (n = 21) 0.749 +- 0.043 (n = 33) 0.321 +- 0.044 (n = 15) 0.396 +- 0.028 (n = 28) 0.1362 0.0812 15 -0.0432 peak in culling window (1b) 0.363 +- 0.026 (n = 19) 0.704 +- 0.040 (n = 33) 0.261 +- 0.021 (n = 14) 0.379 +- 0.026 (n = 28) 0.0056 0.0000 14 -0.0736 peak-to-trough over trace (2b) 0.371 +- 0.027 (n = 19) 0.734 +- 0.042 (n = 33) 0.268 +- 0.021 (n = 14) 0.402 +- 0.026 (n = 28) 0.0015 0.0001 14 -0.0927 variance over N (3) 0.358 +- 0.027 (n = 19) 0.709 +- 0.041 (n = 33) 0.254 +- 0.020 (n = 14) 0.367 +- 0.025 (n = 28) 0.0057 0.0000 14 -0.0769 no negative-f rule in PSP culling (4b) 0.371 +- 0.027 (n = 19) 0.734 +- 0.042 (n = 33) 0.282 +- 0.023 (n = 14) 0.424 +- 0.032 (n = 28) 0.0052 0.0000 14 -0.1009 cullings intersected 0.371 +- 0.027 (n = 19) 0.734 +- 0.042 (n = 33) 0.259 +- 0.021 (n = 14) 0.385 +- 0.031 (n = 28) 0.0085 0.0000 14 -0.0798 heights re-measured after culling 0.371 +- 0.027 (n = 19) 0.734 +- 0.042 (n = 33) 0.267 +- 0.021 (n = 14) 0.385 +- 0.027 (n = 28) 0.0058 0.0001 14 -0.0838 criteria after culling (5b) 0.371 +- 0.027 (n = 19) 0.734 +- 0.042 (n = 33) 0.278 +- 0.020 (n = 16) 0.380 +- 0.024 (n = 28) 0.0058 0.0000 16 -0.0718
Result: partially reproduced (resting on the spike-rule assumption and on details 2, 4 and
5). Total CV is 0.371 +- 0.027 for tones of at least 15 mV (19 neurons) and 0.734 +- 0.042 for
all tones (33), matching 0.37 and 0.73. Private CV is 0.265 +- 0.021 (14 neurons) against 0.383
+- 0.026 (28): the first is within 0.27 +- 0.01 and the second 0.007 below 0.41 +- 0.02. The
private CV of large responses is significantly lower under both readings of the test (Student's
t, p = 0.0055; paired over 14 neurons, p < 0.0001). The criterion's reproduced outcome needs all
four values, so the claim holds in part. With the peak-to-trough height over the whole trace (2b)
the private values are 0.268 and 0.402 and all four hold; with the negative scale factor rule
left out of PSP culling (4b) 0.282 and 0.424, and with the criteria after culling (5b) 0.278 and
0.380, one out each time. With the claims file's spike rule the large-response total CV is 0.491
and the private difference is not significant (p = 0.475 and 0.380); spikethresh keeping the
slow events gives 0.468 and p = 0.136 and 0.081. With the corrected potential above -20 mV at any
rise the private values are 0.286 and 0.389 and the difference is significant (p = 0.014), and
with heights re-measured after culling 0.267 and 0.385; one or two values out, partially
reproduced.
lfp-cv-not-different¶
The LFP recorded with the more variable neurons is not more variable than that recorded with the less variable neurons (CV 0.66 against 0.74, P = 0.5), so differences in total variability are not explained by gross state differences seen in the LFP.
def lfp_cv(label, measure):
out = {}
for name in lfp_neurons:
vals = []
for r in PT[label][name]:
if not ac2.included(r):
continue
x = r["lfp_h"] if measure == "height" else r["s"]
vals.append(np.std(x, ddof=1) / np.mean(x))
out[name] = np.mean(vals)
return pd.Series(out)
for label in ("main", SPK, SPK_ANY, SPK_SLOW):
q29 = S[label].q[lfp_neurons].sort_values(ascending=False)
more, less = q29.index[:14], q29.index[14:]
for measure in ("height", "scale factor"):
cv = lfp_cv(label, measure)
t, p = stats.ttest_ind(cv[more], cv[less])
print(f"[{label}; LFP {measure}] more variable (14): {mse(cv[more])}; less variable (15): {mse(cv[less])};"
f" Student's t = {t:.2f}, p = {p:.3f}; difference {cv[more].mean() - cv[less].mean():.3f}")
cvh = lfp_cv("main", "height")
print("\nLFP CV per neuron (height, main):", cvh.round(2).to_dict())
print("Mean LFP height sign check: tones whose mean LFP height is negative:",
sum(np.mean(r["lfp_h"]) < 0 for name in lfp_neurons for r in PT["main"][name] if ac2.included(r)))
[main; LFP height] more variable (14): 0.61 +- 0.05 (n = 14); less variable (15): 0.57 +- 0.06 (n = 15); Student's t = 0.43, p = 0.674; difference 0.035 [main; LFP scale factor] more variable (14): 0.89 +- 0.08 (n = 14); less variable (15): 0.77 +- 0.09 (n = 15); Student's t = 0.95, p = 0.350; difference 0.116 [spikes: fast rise above -20 mV (6); LFP height] more variable (14): 0.55 +- 0.04 (n = 14); less variable (15): 0.62 +- 0.07 (n = 15); Student's t = -0.81, p = 0.428; difference -0.066 [spikes: fast rise above -20 mV (6); LFP scale factor] more variable (14): 0.84 +- 0.08 (n = 14); less variable (15): 0.83 +- 0.10 (n = 15); Student's t = 0.12, p = 0.905; difference 0.015 [spikes: above -20 mV, any rise; LFP height] more variable (14): 0.65 +- 0.06 (n = 14); less variable (15): 0.53 +- 0.05 (n = 15); Student's t = 1.61, p = 0.118; difference 0.128
[spikes: above -20 mV, any rise; LFP scale factor] more variable (14): 0.93 +- 0.08 (n = 14); less variable (15): 0.73 +- 0.08 (n = 15); Student's t = 1.66, p = 0.108; difference 0.196
[spikes: spikethresh, slow events kept; LFP height] more variable (14): 0.56 +- 0.04 (n = 14); less variable (15): 0.61 +- 0.07 (n = 15); Student's t = -0.67, p = 0.511; difference -0.054
[spikes: spikethresh, slow events kept; LFP scale factor] more variable (14): 0.84 +- 0.08 (n = 14); less variable (15): 0.82 +- 0.10 (n = 15); Student's t = 0.15, p = 0.883; difference 0.018
LFP CV per neuron (height, main): {'020701md01b': 0.51, '020901md01': 0.8, '021001md01': 0.56, '021401md01b': 1.07, '021401md02': 1.03, '021401md03': 0.68, '021401md05': 0.62, '021901md01': 0.6, '022301md01': 0.75, '022801md01b': 0.48, '022801md02b': 0.34, '030201md01b': 0.46, '030201md03b': 0.66, '030301md01b': 0.76, '030301md02b': 1.18, '030301md03b': 0.58, '102400md01': 0.49, '102500md02': 0.57, '102600md01': 0.37, '102600md03': 0.34, '102800md01': 0.44, '102800md02': 0.37, '102900md02': 0.38, '102900md04': 0.53, '110700md01': 0.48, '110900md05': 0.68, '111200md01': 0.59, '111400md01': 0.41, '111500md01': 0.34}
Mean LFP height sign check: tones whose mean LFP height is negative: 0
Result: reproduced. Ranking the 29 LFP neurons by total tones, the LFP height CV is 0.61
+- 0.05 for the 14 more variable and 0.57 +- 0.06 for the 15 less variable neurons (paper 0.66
+- 0.06 and 0.74 +- 0.09); Student's t-test finds no difference (t = 0.43, p = 0.674;
difference 0.035). The less variable group's CV is below the paper's tolerance. The LFP scale
factor CV gives 0.89 against 0.77 (p = 0.350), and the claims file's spike rule
p = 0.428 (0.55 against 0.62). Unlike the
paper's slight opposite trend, the more variable neurons have the slightly higher LFP CV here.
With the corrected potential above -20 mV at any rise the LFP height CVs are 0.65 and 0.53 (p =
0.118), and with
spikethresh keeping the slow events 0.56 and 0.61 (p = 0.511); neither
finds a difference.
anesthesia-comparison¶
Total variability does not differ significantly between the three anesthesia protocols (pentobarbital 3.7, pentobarbital with diazepam 6.1, urethane 5.0 mV), and private variability differs slightly but significantly between pentobarbital alone (1.1 mV) and urethane (2.2 mV).
rows = []
for label, df in S.items():
row = {"reading": label}
for an in ("P", "P_with_Diaz", "U"):
g = df[df.anesthetic == an]
row[f"total {an}"] = mse(g.q)
row[f"private {an}"] = mse(g.q_both_est)
groups = [df.q[df.anesthetic == an] for an in ("P", "P_with_Diaz", "U")]
row["total ANOVA p"] = stats.f_oneway(*groups).pvalue
row["total Kruskal-Wallis p"] = stats.kruskal(*groups).pvalue
pp, pu = df.q_both_est[df.anesthetic == "P"].dropna(), df.q_both_est[df.anesthetic == "U"].dropna()
row["private P vs U t p"] = stats.ttest_ind(pp, pu).pvalue
row["private P vs U Mann-Whitney p"] = stats.mannwhitneyu(pp, pu).pvalue
rows.append(row)
an_table = pd.DataFrame(rows).set_index("reading")
print(an_table.iloc[:, :6].to_string())
print(an_table.iloc[:, 6:].round(4).to_string())
total P private P total P_with_Diaz private P_with_Diaz total U private U
reading
main 3.69 +- 0.70 (n = 13) 1.01 +- 0.13 (n = 12) 6.15 +- 2.81 (n = 3) 1.37 +- 0.33 (n = 3) 4.99 +- 0.58 (n = 17) 1.93 +- 0.29 (n = 13)
spikes: fast rise above -20 mV (6) 3.70 +- 0.70 (n = 13) 1.05 +- 0.14 (n = 12) 6.47 +- 2.84 (n = 3) 1.54 +- 0.32 (n = 3) 6.42 +- 0.66 (n = 17) 3.01 +- 0.53 (n = 13)
spikes: above -20 mV, any rise 3.68 +- 0.69 (n = 13) 1.05 +- 0.14 (n = 12) 6.38 +- 2.82 (n = 3) 1.54 +- 0.32 (n = 3) 4.66 +- 0.43 (n = 17) 2.03 +- 0.28 (n = 13)
spikes: spikethresh, slow events kept 3.70 +- 0.70 (n = 13) 1.05 +- 0.14 (n = 12) 6.15 +- 2.81 (n = 3) 1.37 +- 0.33 (n = 3) 6.00 +- 0.65 (n = 17) 2.43 +- 0.41 (n = 13)
peak in culling window (1b) 3.63 +- 0.68 (n = 13) 1.01 +- 0.13 (n = 12) 5.49 +- 2.27 (n = 3) 1.37 +- 0.33 (n = 3) 4.45 +- 0.57 (n = 17) 1.82 +- 0.27 (n = 13)
peak-to-trough over trace (2b) 3.69 +- 0.70 (n = 13) 1.14 +- 0.15 (n = 12) 6.15 +- 2.81 (n = 3) 1.48 +- 0.44 (n = 3) 4.99 +- 0.58 (n = 17) 1.98 +- 0.28 (n = 13)
variance over N (3) 3.46 +- 0.66 (n = 13) 0.93 +- 0.12 (n = 12) 5.93 +- 2.74 (n = 3) 1.27 +- 0.31 (n = 3) 4.63 +- 0.55 (n = 17) 1.77 +- 0.26 (n = 13)
no negative-f rule in PSP culling (4b) 3.69 +- 0.70 (n = 13) 1.09 +- 0.15 (n = 12) 6.15 +- 2.81 (n = 3) 1.37 +- 0.33 (n = 3) 4.99 +- 0.58 (n = 17) 2.51 +- 0.37 (n = 13)
cullings intersected 3.69 +- 0.70 (n = 13) 0.98 +- 0.13 (n = 12) 6.15 +- 2.81 (n = 3) 1.24 +- 0.34 (n = 3) 4.99 +- 0.58 (n = 17) 2.06 +- 0.37 (n = 13)
heights re-measured after culling 3.69 +- 0.70 (n = 13) 1.01 +- 0.13 (n = 12) 6.15 +- 2.81 (n = 3) 1.34 +- 0.32 (n = 3) 4.99 +- 0.58 (n = 17) 2.01 +- 0.30 (n = 13)
criteria after culling (5b) 3.69 +- 0.70 (n = 13) 1.08 +- 0.15 (n = 12) 6.15 +- 2.81 (n = 3) 1.40 +- 0.32 (n = 3) 4.99 +- 0.58 (n = 17) 1.84 +- 0.23 (n = 13)
total ANOVA p total Kruskal-Wallis p private P vs U t p private P vs U Mann-Whitney p
reading
main 0.2512 0.3538 0.0094 0.0083
spikes: fast rise above -20 mV (6) 0.0393 0.0318 0.0023 0.0025
spikes: above -20 mV, any rise 0.1993 0.3248 0.0064 0.0060
spikes: spikethresh, slow events kept 0.0832 0.0719 0.0058 0.0071
peak in culling window (1b) 0.4580 0.5510 0.0168 0.0114
peak-to-trough over trace (2b) 0.2512 0.3538 0.0165 0.0208
variance over N (3) 0.2467 0.3787 0.0098 0.0071
no negative-f rule in PSP culling (4b) 0.2512 0.3538 0.0024 0.0036
cullings intersected 0.2512 0.3538 0.0142 0.0036
heights re-measured after culling 0.2512 0.3538 0.0063 0.0060
criteria after culling (5b) 0.2512 0.3538 0.0112 0.0180
Result: reproduced (resting on the spike-rule assumption and on details 1, 3, 4 and 5). Total
is 3.69 +- 0.70 (pentobarbital, 13), 6.15 +- 2.81 (with diazepam, 3) and 4.99 +- 0.58 mV
(urethane, 17), all within tolerance; neither the one-way ANOVA (p = 0.251) nor Kruskal-Wallis
(p = 0.354) finds a difference. The private best estimate is 1.01 +- 0.13 (12), 1.37 +- 0.33 (3)
and 1.93 +- 0.29 mV (13), with the paper's group sizes (12, 3, 13), all within tolerance, two of
them by 0.01 and 0.03 mV; pentobarbital alone and urethane differ (Student's t p = 0.0094;
Mann-Whitney p = 0.0083). The margins are small and the readings move them across: deciding both
cullings on the full ensemble gives 0.98 mV for pentobarbital, the peak in the culling window
(1b) 1.82 and the criteria after culling (5b) 1.84 mV for urethane, and leaving the negative scale
factor rule out of PSP culling (4b) 2.51 mV for urethane, and the variance over N 0.93 mV for
pentobarbital, each of which would give partially reproduced. With the peak-to-trough height
over the trace (2b) and with heights re-measured after culling (1.01, 1.34 and 2.01 mV; t p =
0.0063) every value holds. With the claims file's spike rule the urethane total is 6.42 mV and
the ANOVA finds a difference (p = 0.039), which would give not reproduced. With the corrected
potential above -20 mV at any rise every value holds and the tests agree (ANOVA p = 0.199,
private t p = 0.0064); with
spikethresh keeping the slow events the urethane total is 6.00 mV,
outside, which would not give reproduced.
fig8-spontaneous-veto¶
The size of spontaneous activity before the stimulus is anticorrelated with the evoked PSP, in the LFP (mean correlation -0.14, 29 neurons) and in the whole-cell record (-0.22, 33 neurons), the whole-cell correlation being negative in 30 of 33 neurons.
The 100 ms before each tone come from the sweeps (pre_wc, pre_lfp). For the first tone of a
sweep this window overlaps the prepulse in neurons whose prepulse ends after 100 ms; those trials
are left out (counted below).
def spont_size(x, channel, reading):
if reading == "a":
return x.max(1) - x.min(1)
dev = x - x.mean(1, keepdims=True)
return dev.max(1) if channel == "wc" else np.abs(dev).max(1)
def veto(label, channel, size_reading, corr_reading, window="100 ms", drop_spiking=False):
"""Correlation of the spontaneous event before each tone with the PSP height, per neuron.
window "100 ms": the 100 ms before onset; "-100 to -15 ms": without the 15-ms baseline of the
PSP height. drop_spiking: leave out trials with a spike in the window (the neuron's spikethresh
above the trace's 15-ms baseline)."""
out, dropped = {}, 0
names = all_names if channel == "wc" else lfp_names
stop = ac2.PRESTIM if window == "100 ms" else ac2.PRESTIM - ac2.PRE
for name in names:
n = N[name]
pre = n["pre_wc"] if channel == "wc" else n["pre_lfp"]
per_tone, pooled_x, pooled_h = [], [], []
for r in PT[label][name]:
if not ac2.included(r):
continue
ok = ~n["prestim_overlaps_prepulse"][r["trials"]]
if drop_spiking:
wc = n["pre_wc"][r["trials"], :stop] - n["rbaseline"][r["trials"], None]
ok &= wc.max(1) <= n["params"]["spikethresh"]
dropped += int((~ok).sum())
x = spont_size(pre[r["trials"], :stop].astype(float), channel, size_reading)[ok]
h = r["h"][ok]
if len(h) >= 3 and np.ptp(x) > 0:
per_tone.append(stats.pearsonr(x, h)[0])
pooled_x.append(x)
pooled_h.append(h / h.mean())
if per_tone:
out[name] = np.mean(per_tone) if corr_reading == "a" else stats.pearsonr(
np.concatenate(pooled_x), np.concatenate(pooled_h))[0]
return pd.Series(out), dropped
def veto_rows(labels, windows, spiking):
rows = []
for label in labels:
for window in windows:
for ds in spiking:
for sr in ("a", "b"):
for cr in ("a", "b"):
wc, dropped = veto(label, "wc", sr, cr, window, ds)
lf, _ = veto(label, "lfp", sr, cr, window, ds)
rows.append({"reading": label, "window": window, "spiking trials out": ds, "size": sr, "correlation": cr,
"LFP mean r": lf.mean(), "LFP SE": lf.std(ddof=1) / np.sqrt(len(lf)), "LFP n": len(lf),
"LFP negative": int((lf < 0).sum()),
"WC mean r": wc.mean(), "WC SE": wc.std(ddof=1) / np.sqrt(len(wc)), "WC n": len(wc),
"WC negative": int((wc < 0).sum()), "trials left out": dropped})
return pd.DataFrame(rows)
print(veto_rows(("main", SPK, SPK_ANY, SPK_SLOW), ("100 ms",), (False,)).round(3).to_string(index=False))
print("\nControls for the overlap with the PSP baseline and for spikes in the window (main):")
print(veto_rows(("main",), ("100 ms", "-100 to -15 ms"), (False, True)).round(3).to_string(index=False))
wc, _ = veto("main", "wc", "a", "a")
print("\nWhole-cell correlation per neuron (main, readings a, a):", wc.round(3).to_dict())
print(f"Null: with PSP heights shuffled within tones the correlation is 0 in expectation; per-tone Pearson r"
f" averaged over tones of {np.mean([len(r['h']) for r in PT['main'][FIG5] if ac2.included(r)]):.0f} trials (Figure 5 neuron).")
print("Figure 5 neuron:", FIG5, "whole-cell r (a, a) = %.3f" % wc[FIG5])
reading window spiking trials out size correlation LFP mean r LFP SE LFP n LFP negative WC mean r WC SE WC n WC negative trials left out
main 100 ms False a a -0.168 0.030 29 25 -0.238 0.040 33 29 257
main 100 ms False a b -0.155 0.028 29 24 -0.193 0.031 33 29 257
main 100 ms False b a -0.167 0.028 29 25 -0.223 0.040 33 28 257
main 100 ms False b b -0.149 0.026 29 25 -0.210 0.038 33 28 257
spikes: fast rise above -20 mV (6) 100 ms False a a -0.156 0.029 29 24 -0.224 0.039 33 28 261
spikes: fast rise above -20 mV (6) 100 ms False a b -0.148 0.027 29 23 -0.179 0.031 33 27 261
spikes: fast rise above -20 mV (6) 100 ms False b a -0.155 0.028 29 23 -0.212 0.039 33 28 261
spikes: fast rise above -20 mV (6) 100 ms False b b -0.143 0.026 29 24 -0.197 0.038 33 28 261
spikes: above -20 mV, any rise 100 ms False a a -0.169 0.028 29 26 -0.236 0.040 33 29 261
spikes: above -20 mV, any rise 100 ms False a b -0.161 0.027 29 26 -0.187 0.031 33 29 261
spikes: above -20 mV, any rise 100 ms False b a -0.167 0.027 29 26 -0.223 0.040 33 28 261
spikes: above -20 mV, any rise 100 ms False b b -0.155 0.026 29 27 -0.205 0.038 33 29 261
spikes: spikethresh, slow events kept 100 ms False a a -0.162 0.030 29 24 -0.230 0.040 33 28 258
spikes: spikethresh, slow events kept 100 ms False a b -0.151 0.028 29 23 -0.181 0.031 33 27 258
spikes: spikethresh, slow events kept 100 ms False b a -0.161 0.029 29 24 -0.218 0.040 33 28 258
spikes: spikethresh, slow events kept 100 ms False b b -0.146 0.026 29 24 -0.201 0.038 33 28 258
Controls for the overlap with the PSP baseline and for spikes in the window (main):
reading window spiking trials out size correlation LFP mean r LFP SE LFP n LFP negative WC mean r WC SE WC n WC negative trials left out main 100 ms False a a -0.168 0.030 29 25 -0.238 0.040 33 29 257 main 100 ms False a b -0.155 0.028 29 24 -0.193 0.031 33 29 257 main 100 ms False b a -0.167 0.028 29 25 -0.223 0.040 33 28 257 main 100 ms False b b -0.149 0.026 29 25 -0.210 0.038 33 28 257 main 100 ms True a a -0.167 0.030 29 25 -0.239 0.040 33 29 298 main 100 ms True a b -0.154 0.028 29 24 -0.197 0.032 33 29 298 main 100 ms True b a -0.166 0.028 29 25 -0.224 0.040 33 28 298 main 100 ms True b b -0.148 0.026 29 25 -0.222 0.040 33 28 298 main -100 to -15 ms False a a -0.158 0.028 29 25 -0.220 0.040 33 29 257 main -100 to -15 ms False a b -0.141 0.028 29 24 -0.172 0.029 33 29 257 main -100 to -15 ms False b a -0.157 0.028 29 25 -0.218 0.041 33 28 257 main -100 to -15 ms False b b -0.135 0.027 29 24 -0.202 0.038 33 28 257 main -100 to -15 ms True a a -0.157 0.028 29 25 -0.219 0.040 33 29 298 main -100 to -15 ms True a b -0.140 0.027 29 24 -0.176 0.030 33 29 298 main -100 to -15 ms True b a -0.155 0.028 29 25 -0.218 0.041 33 28 298 main -100 to -15 ms True b b -0.134 0.027 29 24 -0.216 0.040 33 28 298
Whole-cell correlation per neuron (main, readings a, a): {'020701md01b': -0.451, '020901md01': -0.331, '021001md01': 0.007, '021401md01b': -0.367, '021401md02': 0.559, '021401md03': 0.279, '021401md05': -0.279, '021901md01': -0.523, '022301md01': -0.313, '022801md01b': -0.262, '022801md02b': -0.282, '030201md01b': 0.099, '030201md03b': -0.256, '030301md01b': -0.15, '030301md02b': -0.29, '030301md03b': -0.455, '100300md01': -0.2, '100300md02': -0.439, '100500md01': -0.199, '100500md02b': -0.013, '102400md01': -0.466, '102500md02': -0.376, '102600md01': -0.09, '102600md03': -0.137, '102800md01': -0.064, '102800md02': -0.126, '102900md02': -0.346, '102900md04': -0.442, '110700md01': -0.444, '110900md05': -0.324, '111200md01': -0.317, '111400md01': -0.405, '111500md01': -0.448}
Null: with PSP heights shuffled within tones the correlation is 0 in expectation; per-tone Pearson r averaged over tones of 63 trials (Figure 5 neuron).
Figure 5 neuron: 020701md01b whole-cell r (a, a) = -0.451
Result: reproduced. Under readings (a) size as maximum minus minimum and (b) one correlation
over all trials with heights normalized by the tone's mean, the mean correlation is -0.155 +-
0.028 for the LFP (29 neurons; -0.14 +- 0.02) and -0.193 +- 0.031 for the whole cell (33;
-0.22 +- 0.03), negative in 29 of 33 neurons (30 +- 3). Under readings (b, b) -0.149, -0.210 and
28 also hold. Under (a, a) the LFP mean, -0.168, is outside its tolerance. All combinations, with
either spike rule, give negative means for both channels (-0.143 to -0.168 LFP, -0.179 to -0.238
whole cell). The 15-ms baseline of the PSP height lies inside the 100-ms window, so a
depolarization just before onset lowers the height by construction. With the window ending 15 ms
before onset, so that it does not overlap the baseline, the correlations remain (-0.134 to -0.158
LFP, -0.172 to -0.220 whole cell), and under (a, a) and (b, b) all three values hold (-0.158,
-0.220, 29; -0.135, -0.202, 28); under (a, b) the whole-cell mean, -0.172, is outside. With the
corrected potential above -20 mV at any rise, (b, b) gives -0.155, -0.205 and 29, and with
spikethresh keeping the slow events -0.146, -0.201 and 28, all within; under (a, b) the
whole-cell means (-0.187, -0.181) are outside. Leaving out trials with a spike in the window changes the means by
at most 0.014. 257 trials whose window overlaps the prepulse are left out.
overfitting-correction-example¶
For the Figure 5 neuron, the shuffle control attributes 0.02 mV to overfitting of the linear estimator, and the corrected value after both cullings and estimation is 0.57 mV.
for label in ("main", SPK, SPK_ANY, SPK_SLOW, "peak in culling window (1b)", "peak-to-trough over trace (2b)", "variance over N (3)",
"no negative-f rule in PSP culling (4b)", "cullings intersected", "heights re-measured after culling"):
r = S[label].loc[FIG5]
print(f"[{label}] {FIG5}: <q LFP-culled,PSP-culled> {r.q_both:.3f}; shuffled estimation {r.q_both_shuf:.3f};"
f" overfitting {r.q_both - r.q_both_shuf:.3f}; estimation uncorrected {r.q_both_est_uncorrected:.3f};"
f" corrected {r.q_both_est:.3f}; tones {int(r.n_tones_both)}")
print("Correction per tone (detail 7) gives the same tone average as per neuron: mean(est + q - shuf) = mean(est) + mean(q) - mean(shuf).")
# The shuffle reduction for a regression with one predictor is about q / (N - 1) per tone.
tones = [t for t in PT["main"][FIG5] if ac2.included(t, "both")]
print("Trials left per tone after both cullings:", [int(t["n_both"]) for t in tones])
[main] 020701md01b: <q LFP-culled,PSP-culled> 0.757; shuffled estimation 0.743; overfitting 0.014; estimation uncorrected 0.524; corrected 0.537; tones 17 [spikes: fast rise above -20 mV (6)] 020701md01b: <q LFP-culled,PSP-culled> 0.757; shuffled estimation 0.743; overfitting 0.014; estimation uncorrected 0.524; corrected 0.537; tones 17 [spikes: above -20 mV, any rise] 020701md01b: <q LFP-culled,PSP-culled> 0.757; shuffled estimation 0.743; overfitting 0.014; estimation uncorrected 0.524; corrected 0.537; tones 17 [spikes: spikethresh, slow events kept] 020701md01b: <q LFP-culled,PSP-culled> 0.757; shuffled estimation 0.743; overfitting 0.014; estimation uncorrected 0.524; corrected 0.537; tones 17 [peak in culling window (1b)] 020701md01b: <q LFP-culled,PSP-culled> 0.757; shuffled estimation 0.743; overfitting 0.014; estimation uncorrected 0.524; corrected 0.537; tones 17 [peak-to-trough over trace (2b)] 020701md01b: <q LFP-culled,PSP-culled> 1.118; shuffled estimation 1.085; overfitting 0.033; estimation uncorrected 0.859; corrected 0.892; tones 18 [variance over N (3)] 020701md01b: <q LFP-culled,PSP-culled> 0.743; shuffled estimation 0.730; overfitting 0.013; estimation uncorrected 0.514; corrected 0.528; tones 17 [no negative-f rule in PSP culling (4b)] 020701md01b: <q LFP-culled,PSP-culled> 0.986; shuffled estimation 0.971; overfitting 0.015; estimation uncorrected 0.627; corrected 0.643; tones 17 [cullings intersected] 020701md01b: <q LFP-culled,PSP-culled> 0.759; shuffled estimation 0.746; overfitting 0.013; estimation uncorrected 0.525; corrected 0.538; tones 17 [heights re-measured after culling] 020701md01b: <q LFP-culled,PSP-culled> 0.756; shuffled estimation 0.742; overfitting 0.014; estimation uncorrected 0.522; corrected 0.536; tones 17 Correction per tone (detail 7) gives the same tone average as per neuron: mean(est + q - shuf) = mean(est) + mean(q) - mean(shuf). Trials left per tone after both cullings: [52, 57, 58, 46, 59, 51, 53, 59, 51, 57, 57, 58, 54, 55, 59, 59, 54]
Result: reproduced (resting on details 2 and 4). In the Figure 5 neuron (020701md01b), after
both cullings tones is 0.757 mV (0.76), the shuffled estimation gives 0.743 mV (0.74), the
estimation 0.524 mV uncorrected (0.55) and 0.537 mV corrected (0.57); all four within 10%. The
overfitting share is 0.014 mV against the paper's 0.02. The values are the same with the claims
file's spike rule (no spike trials in this neuron) and nearly so with the cullings decided on the
full ensemble (0.759, 0.746, 0.525, 0.538). With the peak-to-trough height over the whole trace
(2b) they are 1.118, 1.085, 0.859 and 0.892 mV, and without the negative scale factor rule in PSP
culling (4b) 0.986, 0.971, 0.627 and 0.643 mV; in both the corrected value is outside, which
would give not reproduced. Per-tone and per-neuron corrections give the same tone average (the
identity printed above), so detail 7's per-tone reading needs no run of its own. The two further
spike rules and the peak in the culling window give the main values, the variance over N 0.743,
0.730, 0.514 and 0.528 mV and heights re-measured after culling 0.756, 0.742, 0.522 and 0.536 mV,
all within.
private-shared-ranges¶
The private component of the variability is small, about 1 to 3 mV, and relatively homogeneous across neurons, while the shared component is often much larger and more heterogeneous, ranging from about 0 to 10 mV.
rows = []
for label in S:
d = S[label].dropna(subset=["q_both_est"])
priv = d.q_both_est
for total_reading, tot in (("all included tones", d.q), ("tones kept for the private estimate", d.q_on_both_tones)):
shared = tot - priv
rows.append({"reading": label, "total over": total_reading, "private in 0.9-3.3": f"{int(priv.between(0.9, 3.3).sum())} of {len(priv)}",
"%": 100 * priv.between(0.9, 3.3).mean(), "private min": priv.min(), "private max": priv.max(),
"private SD": priv.std(ddof=1), "shared min": shared.min(), "shared max": shared.max(),
"shared within -1 to 11": bool(shared.between(-1, 11).all()), "shared SD": shared.std(ddof=1)})
print(pd.DataFrame(rows).set_index(["reading", "total over"]).round(2).to_string())
d = S["main"].dropna(subset=["q_both_est"])
print("\nPrivate values outside 0.9-3.3 mV (main):", d.q_both_est[~d.q_both_est.between(0.9, 3.3)].round(2).to_dict())
private in 0.9-3.3 % private min private max private SD shared min shared max shared within -1 to 11 shared SD
reading total over
main all included tones 20 of 28 71.43 0.31 4.52 0.88 0.25 9.26 True 2.60
tones kept for the private estimate 20 of 28 71.43 0.31 4.52 0.88 0.10 5.31 True 1.38
spikes: fast rise above -20 mV (6) all included tones 16 of 28 57.14 0.31 6.96 1.64 0.25 8.96 True 2.46
tones kept for the private estimate 16 of 28 57.14 0.31 6.96 1.64 0.10 7.00 True 1.75
spikes: above -20 mV, any rise all included tones 20 of 28 71.43 0.31 4.52 0.90 0.25 8.91 True 2.29
tones kept for the private estimate 20 of 28 71.43 0.31 4.52 0.90 0.10 5.23 True 1.44
spikes: spikethresh, slow events kept all included tones 18 of 28 64.29 0.31 4.97 1.25 0.25 9.19 True 2.58
tones kept for the private estimate 18 of 28 64.29 0.31 4.97 1.25 0.10 5.25 True 1.47
peak in culling window (1b) all included tones 20 of 28 71.43 0.31 4.52 0.83 0.25 8.47 True 2.40
tones kept for the private estimate 20 of 28 71.43 0.31 4.52 0.83 0.10 4.41 True 1.30
peak-to-trough over trace (2b) all included tones 21 of 28 75.00 0.32 4.52 0.88 0.25 9.27 True 2.57
tones kept for the private estimate 21 of 28 75.00 0.32 4.52 0.88 0.09 5.32 True 1.38
variance over N (3) all included tones 18 of 28 64.29 0.28 4.14 0.81 0.25 8.74 True 2.50
tones kept for the private estimate 18 of 28 64.29 0.28 4.14 0.81 0.09 5.00 True 1.33
no negative-f rule in PSP culling (4b) all included tones 20 of 28 71.43 0.30 5.92 1.19 -0.82 8.91 True 2.55
tones kept for the private estimate 20 of 28 71.43 0.30 5.92 1.19 0.08 4.46 True 1.29
cullings intersected all included tones 18 of 28 64.29 0.30 5.94 1.09 -0.84 9.25 True 2.71
tones kept for the private estimate 18 of 28 64.29 0.30 5.94 1.09 0.10 5.30 True 1.32
heights re-measured after culling all included tones 20 of 28 71.43 0.31 4.64 0.92 0.24 9.17 True 2.58
tones kept for the private estimate 20 of 28 71.43 0.31 4.64 0.92 0.10 5.22 True 1.35
criteria after culling (5b) all included tones 21 of 28 75.00 0.31 3.34 0.75 0.18 9.26 True 2.55
tones kept for the private estimate 21 of 28 75.00 0.31 3.34 0.75 0.10 5.31 True 1.38
Private values outside 0.9-3.3 mV (main): {'020701md01b': 0.54, '021401md03': 0.31, '022801md01b': 0.63, '030301md02b': 0.6, '030301md03b': 0.83, '102800md01': 0.68, '110900md05': 4.52, '111500md01': 0.77}
Result: partially reproduced. Two of the three parts hold. Shared variability (total minus private) ranges from 0.25 to 9.26 mV with total over all included tones, and 0.10 to 5.31 mV with total over the tones kept for the private estimate, both within -1 to 11 mV, and its SD across neurons (2.60 or 1.38 mV) exceeds that of private variability (0.88 mV). But only 20 of 28 neurons (71%) have a private value between 0.9 and 3.3 mV, short of the 80% the criterion sets: seven are below 0.9 mV (0.31 to 0.83) and one is above (4.52 mV). Every reading gives 57% to 75%, so the outcome is the same under all of them.
quantal-comparison¶
In 14 of 33 neurons at least one tone elicits responses whose variability index is within a factor of 2 of the 0.38 mV expected from stochastic release alone.
rows = []
for label in ("main", SPK, SPK_ANY, SPK_SLOW, "peak in culling window (1b)", "variance over N (3)"):
for ref in (0.38, 0.39):
qs = {name: [r["q"] for r in tones if ac2.included(r)] for name, tones in PT[label].items()}
rows.append({"reading": label, "reference": ref,
"q <= 2 x ref": sum(any(q <= 2 * ref for q in v) for v in qs.values()),
"ref/2 <= q <= 2 x ref": sum(any(ref / 2 <= q <= 2 * ref for q in v) for v in qs.values())})
print(pd.DataFrame(rows).to_string(index=False))
mins = pd.Series({name: min(r["q"] for r in tones if ac2.included(r)) for name, tones in PT["main"].items()})
print("Smallest included-tone q per neuron (main), sorted:", mins.sort_values().round(2).tolist())
reading reference q <= 2 x ref ref/2 <= q <= 2 x ref
main 0.38 14 14
main 0.39 15 15
spikes: fast rise above -20 mV (6) 0.38 14 14
spikes: fast rise above -20 mV (6) 0.39 15 15
spikes: above -20 mV, any rise 0.38 14 14
spikes: above -20 mV, any rise 0.39 15 15
spikes: spikethresh, slow events kept 0.38 14 14
spikes: spikethresh, slow events kept 0.39 15 15
peak in culling window (1b) 0.38 14 14
peak in culling window (1b) 0.39 15 15
variance over N (3) 0.38 17 17
variance over N (3) 0.39 20 20
Smallest included-tone q per neuron (main), sorted: [0.14, 0.23, 0.3, 0.33, 0.35, 0.35, 0.4, 0.4, 0.41, 0.41, 0.46, 0.47, 0.48, 0.5, 0.77, 0.79, 0.81, 0.82, 0.83, 0.85, 0.89, 0.92, 0.95, 1.02, 1.57, 1.59, 1.63, 1.64, 1.78, 1.79, 2.46, 2.84, 3.02]
Result: reproduced. 14 of 33 neurons have at least one included tone with q within a factor of 2 of 0.38 mV, under both readings of "within a factor of 2" (the smallest q is 0.14 mV, but that neuron has another tone within range). With 0.39 mV the count is 15, also within 14 +- 1.4. The claims file's spike rule, the two further spike rules and the peak in the culling window give the same counts; with the variance over N, the claims file's variant of detail 3, 17 and 20.
mepsp-size¶
Miniature EPSPs recorded in vivo after TTX have a mean size of 0.39 mV with a CV of 0.36 across 6 neurons, consistent with in vitro estimates.
The paper gives no detection method. Two standard detectors are run (claims file readings a and
b): (a) template matching (Clements and Bekkers 1997; template rise 0.5 ms, decay 5 ms), with a
threshold on the ratio of the fitted scale to its standard error and the fitted scale as the
size; (b) a threshold on the derivative of the trace smoothed over 1 ms, in noise SDs of that
derivative, with the size taken as the peak in the 10 ms after onset minus the mean of the 2 ms
before it. The traces are digitized in steps of 0.049 mV, so most differences between samples
are 0; the noise SD is the SD after clipping at 3 SD, which does not collapse (the median-based
SD is printed beside it). Each detector is also run on the sign-inverted trace: mEPSPs are
depolarizing, so the events found there estimate the false positives at that threshold. Sweeps
are those of each recording's stable epoch (first_trial to last_trial), without their first
0.3 s (prepulse); which file is trial k is not documented, so three mappings are run. One
recording, 030401md02h, is divided by its documented gain of 4.
# The TTX recordings: sampling, sweeps used under each mapping, digitization and noise.
rows = []
for name, rec in D["ttx"].items():
fs = 1.0 / rec["dx"]
used = {order: sorted(ac2.mini_sweeps(rec, order)) for order in ("numeric", "dir")}
x = np.concatenate(list(ac2.mini_sweeps(rec, "numeric").values()))
d = np.diff(x)
dd = np.diff(ac2._smooth(x, fs))
rows.append({"recording": name, "fs_hz": fs, "files": len(rec["ep"]), "stimulus files": rec["n_stim_files"],
"first_trial": int(rec["params"]["first_trial"]), "last_trial": int(rec["params"]["last_trial"]),
"files used (numeric)": f"ep_{used['numeric'][0]}..ep_{used['numeric'][-1]}",
"files used (dir order)": ",".join(map(str, used["dir"])), "seconds": len(x) / fs, "gain": rec["gain"],
"median mV": np.median(x), "zero differences": (d == 0).mean(), "step mV": np.abs(d[d != 0]).min(),
"noise SD (clipped)": ac2.noise_sd(dd), "noise SD (MAD)": ac2.noise_sd(dd, "mad")})
ttx = pd.DataFrame(rows).set_index("recording")
print(ttx.round(4).to_string())
fs_hz files stimulus files first_trial last_trial files used (numeric) files used (dir order) seconds gain median mV zero differences step mV noise SD (clipped) noise SD (MAD) recording 030401md02h 4000.0 18 18 1 18 ep_0..ep_17 0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17 138.6 4.0 2.6611 0.1394 0.0122 0.0229 0.0226 050301md01 4000.0 24 24 1 8 ep_0..ep_7 0,1,10,11,12,13,14,15 61.6 1.0 -60.2539 0.4400 0.0488 0.0260 0.0181 050401md02 4000.0 5 5 1 1 ep_0..ep_0 0 7.7 1.0 -70.8008 0.4226 0.0488 0.0329 0.0362 050501md01 4000.0 37 37 12 21 ep_11..ep_20 2,19,20,21,22,23,24,25,26,27 77.0 1.0 -61.9629 0.5151 0.0488 0.0134 0.0181 050501md04 4000.0 81 81 1 36 ep_0..ep_35 0,1,2,3,4,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40 277.2 1.0 -58.0566 0.5151 0.0488 0.0215 0.0181 050801md01c 4000.0 44 44 5 26 ep_4..ep_25 2,3,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31 169.4 1.0 -53.9062 0.6849 0.0488 0.0100 0.0181
# Each detector and threshold, on the trace and on the inverted trace, per recording and mapping.
KS = (3, 4, 5, 6, 7, 8)
def detect(x, fs, det, k):
return ac2.minis_template(x, fs, k) if det == "template" else ac2.minis_derivative(x, fs, k)
rows, EVENTS = [], {}
for order in ("numeric", "dir", "all"):
for name, rec in D["ttx"].items():
fs = 1.0 / rec["dx"]
sweeps = ac2.mini_sweeps(rec, order)
secs = sum(len(s) for s in sweeps.values()) / fs
for det in ("template", "derivative"):
for k in KS:
pos = {i: detect(s, fs, det, k) for i, s in sweeps.items()}
n_neg = sum(len(detect(-s, fs, det, k)["onsets"]) for s in sweeps.values())
size = np.concatenate([e["scale" if det == "template" else "size"] for e in pos.values()])
t_on = np.concatenate([e["onsets"] / fs + ac2.MINI_SKIP_S for e in pos.values()])
EVENTS[order, name, det, k] = pos
rows.append({"mapping": order, "recording": name, "detector": det, "k": k, "seconds": secs,
"events": len(size), "rate_hz": len(size) / secs, "inverted_rate_hz": n_neg / secs,
"inverted/trace": n_neg / len(size) if len(size) else np.nan,
"mean_mV": size.mean() if len(size) else np.nan,
"cv": size.std(ddof=1) / size.mean() if len(size) > 1 else np.nan,
"in 50 ms after a tone": np.mean((t_on - 0.2) % 0.5 < 0.05) if len(size) else np.nan})
minis = pd.DataFrame(rows)
print("Numeric mapping, per recording:")
print(minis[minis.mapping == "numeric"].drop(columns="mapping").round(3).to_string(index=False))
print("\nMean over recordings (paper: 0.39 +- 0.08 mV, CV 0.36 +- 0.03, n = 6); 'in 50 ms after a tone' is 0.10 by chance:")
summ = minis.groupby(["mapping", "detector", "k"]).agg(
mean_size=("mean_mV", "mean"), se_size=("mean_mV", lambda x: x.std(ddof=1) / np.sqrt(x.notna().sum())),
mean_cv=("cv", "mean"), recordings=("mean_mV", "count"), rate_hz=("rate_hz", "mean"),
inverted_over_trace=("inverted/trace", "mean"), after_tone=("in 50 ms after a tone", "mean"))
print(summ.round(3).to_string())
Numeric mapping, per recording:
recording detector k seconds events rate_hz inverted_rate_hz inverted/trace mean_mV cv in 50 ms after a tone
030401md02h template 3 138.6 3382 24.401 22.172 0.909 0.398 0.597 0.091
030401md02h template 4 138.6 3175 22.908 21.003 0.917 0.460 0.578 0.090
030401md02h template 5 138.6 2629 18.968 17.100 0.901 0.538 0.545 0.088
030401md02h template 6 138.6 1729 12.475 11.190 0.897 0.637 0.493 0.102
030401md02h template 7 138.6 792 5.714 5.390 0.943 0.701 0.443 0.100
030401md02h template 8 138.6 185 1.335 1.328 0.995 0.702 0.423 0.124
030401md02h derivative 3 138.6 1097 7.915 6.537 0.826 0.486 0.509 0.111
030401md02h derivative 4 138.6 235 1.696 0.801 0.472 0.648 0.390 0.128
030401md02h derivative 5 138.6 65 0.469 0.072 0.154 0.778 0.337 0.138
030401md02h derivative 6 138.6 15 0.108 0.007 0.067 1.043 0.229 0.067
030401md02h derivative 7 138.6 5 0.036 0.000 0.000 1.264 0.105 0.000
030401md02h derivative 8 138.6 2 0.014 0.000 0.000 1.223 0.025 0.000
050301md01 template 3 61.6 1556 25.260 20.747 0.821 0.396 0.461 0.093
050301md01 template 4 61.6 1112 18.052 13.052 0.723 0.451 0.414 0.091
050301md01 template 5 61.6 686 11.136 6.672 0.599 0.519 0.386 0.099
050301md01 template 6 61.6 387 6.282 2.646 0.421 0.587 0.378 0.085
050301md01 template 7 61.6 190 3.084 0.795 0.258 0.675 0.371 0.074
050301md01 template 8 61.6 90 1.461 0.227 0.156 0.789 0.340 0.067
050301md01 derivative 3 61.6 551 8.945 4.919 0.550 0.483 0.504 0.085
050301md01 derivative 4 61.6 168 2.727 0.601 0.220 0.659 0.468 0.077
050301md01 derivative 5 61.6 69 1.120 0.081 0.072 0.884 0.367 0.101
050301md01 derivative 6 61.6 38 0.617 0.032 0.053 1.058 0.298 0.105
050301md01 derivative 7 61.6 21 0.341 0.000 0.000 1.163 0.291 0.048
050301md01 derivative 8 61.6 13 0.211 0.000 0.000 1.336 0.223 0.000
050401md02 template 3 7.7 88 11.429 6.494 0.568 0.281 0.386 0.080
050401md02 template 4 7.7 36 4.675 1.818 0.389 0.319 0.407 0.056
050401md02 template 5 7.7 15 1.948 0.649 0.333 0.298 0.268 0.000
050401md02 template 6 7.7 8 1.039 0.130 0.125 0.310 0.340 0.000
050401md02 template 7 7.7 4 0.519 0.000 0.000 0.374 0.316 0.000
050401md02 template 8 7.7 2 0.260 0.000 0.000 0.304 0.078 0.000
050401md02 derivative 3 7.7 43 5.584 6.104 1.093 0.532 0.196 0.140
050401md02 derivative 4 7.7 9 1.169 1.558 1.333 0.590 0.150 0.111
050401md02 derivative 5 7.7 1 0.130 0.130 1.000 0.609 NaN 0.000
050401md02 derivative 6 7.7 0 0.000 0.000 NaN NaN NaN NaN
050401md02 derivative 7 7.7 0 0.000 0.000 NaN NaN NaN NaN
050401md02 derivative 8 7.7 0 0.000 0.000 NaN NaN NaN NaN
050501md01 template 3 77.0 2108 27.377 22.805 0.833 0.173 0.439 0.095
050501md01 template 4 77.0 1789 23.234 15.844 0.682 0.192 0.389 0.099
050501md01 template 5 77.0 1313 17.052 8.857 0.519 0.214 0.341 0.100
050501md01 template 6 77.0 750 9.740 3.974 0.408 0.242 0.310 0.092
050501md01 template 7 77.0 347 4.506 1.247 0.277 0.271 0.292 0.072
050501md01 template 8 77.0 122 1.584 0.260 0.164 0.301 0.298 0.090
050501md01 derivative 3 77.0 619 8.039 1.935 0.241 0.213 0.424 0.115
050501md01 derivative 4 77.0 150 1.948 0.065 0.033 0.283 0.335 0.120
050501md01 derivative 5 77.0 46 0.597 0.000 0.000 0.361 0.238 0.043
050501md01 derivative 6 77.0 21 0.273 0.000 0.000 0.385 0.219 0.048
050501md01 derivative 7 77.0 6 0.078 0.000 0.000 0.408 0.165 0.000
050501md01 derivative 8 77.0 2 0.026 0.000 0.000 0.456 0.024 0.000
050501md04 template 3 277.2 4296 15.498 13.431 0.867 0.241 0.377 0.098
050501md04 template 4 277.2 2123 7.659 5.877 0.767 0.288 0.329 0.089
050501md04 template 5 277.2 963 3.474 2.038 0.587 0.331 0.306 0.090
050501md04 template 6 277.2 338 1.219 0.541 0.444 0.379 0.310 0.118
050501md04 template 7 277.2 122 0.440 0.144 0.328 0.424 0.297 0.172
050501md04 template 8 277.2 43 0.155 0.032 0.209 0.436 0.237 0.163
050501md04 derivative 3 277.2 1381 4.982 3.110 0.624 0.301 0.489 0.127
050501md04 derivative 4 277.2 203 0.732 0.152 0.207 0.486 0.432 0.113
050501md04 derivative 5 277.2 71 0.256 0.007 0.028 0.599 0.437 0.085
050501md04 derivative 6 277.2 22 0.079 0.000 0.000 0.771 0.449 0.136
050501md04 derivative 7 277.2 11 0.040 0.000 0.000 0.892 0.496 0.182
050501md04 derivative 8 277.2 5 0.018 0.000 0.000 1.164 0.467 0.000
050801md01c template 3 169.4 4928 29.091 22.916 0.788 0.186 0.447 0.098
050801md01c template 4 169.4 4389 25.909 15.118 0.584 0.201 0.402 0.100
050801md01c template 5 169.4 3690 21.783 7.851 0.360 0.215 0.363 0.101
050801md01c template 6 169.4 2776 16.387 2.887 0.176 0.233 0.333 0.097
050801md01c template 7 169.4 1749 10.325 0.626 0.061 0.254 0.315 0.090
050801md01c template 8 169.4 857 5.059 0.041 0.008 0.273 0.300 0.088
050801md01c derivative 3 169.4 3280 19.362 0.632 0.033 0.202 0.402 0.117
050801md01c derivative 4 169.4 782 4.616 0.000 0.000 0.253 0.321 0.133
050801md01c derivative 5 169.4 163 0.962 0.000 0.000 0.300 0.255 0.166
050801md01c derivative 6 169.4 157 0.927 0.000 0.000 0.300 0.257 0.166
050801md01c derivative 7 169.4 27 0.159 0.000 0.000 0.343 0.226 0.148
050801md01c derivative 8 169.4 5 0.030 0.000 0.000 0.363 0.107 0.200
Mean over recordings (paper: 0.39 +- 0.08 mV, CV 0.36 +- 0.03, n = 6); 'in 50 ms after a tone' is 0.10 by chance:
mean_size se_size mean_cv recordings rate_hz inverted_over_trace after_tone
mapping detector k
all derivative 3 0.944 0.545 0.713 6 8.863 0.652 0.119
4 1.252 0.680 0.708 6 1.972 0.490 0.125
5 1.641 0.881 0.640 6 0.540 0.357 0.156
6 2.255 1.288 0.537 6 0.275 0.422 0.163
7 2.662 1.312 0.488 6 0.118 0.236 0.065
8 1.622 0.455 0.417 5 0.059 0.112 0.040
template 3 0.704 0.363 0.668 6 21.983 0.892 0.097
4 0.803 0.408 0.630 6 16.756 0.821 0.093
5 0.882 0.427 0.593 6 12.098 0.715 0.092
6 0.948 0.433 0.555 6 7.774 0.616 0.096
7 1.062 0.488 0.520 6 4.106 0.514 0.096
8 1.094 0.484 0.510 6 1.605 0.483 0.099
dir derivative 3 0.422 0.077 0.531 6 9.186 0.641 0.116
4 0.540 0.089 0.440 6 2.125 0.470 0.121
5 0.631 0.098 0.440 6 0.559 0.279 0.107
6 0.747 0.158 0.296 5 0.266 0.057 0.117
7 0.862 0.190 0.263 5 0.108 0.027 0.173
8 1.015 0.214 0.221 4 0.056 0.057 0.036
template 3 0.332 0.055 0.603 6 21.977 0.801 0.096
4 0.381 0.066 0.571 6 17.058 0.690 0.089
5 0.425 0.080 0.515 6 12.572 0.584 0.082
6 0.471 0.087 0.497 6 8.249 0.465 0.089
7 0.521 0.090 0.447 6 4.381 0.381 0.098
8 0.525 0.097 0.351 6 1.698 0.341 0.098
numeric derivative 3 0.370 0.061 0.421 6 9.138 0.561 0.116
4 0.487 0.074 0.350 6 2.148 0.378 0.114
5 0.588 0.093 0.327 6 0.589 0.209 0.089
6 0.712 0.160 0.290 5 0.334 0.024 0.104
7 0.814 0.189 0.257 5 0.109 0.000 0.076
8 0.909 0.206 0.169 5 0.050 0.000 0.040
template 3 0.279 0.040 0.451 6 22.176 0.798 0.093
4 0.318 0.048 0.420 6 17.073 0.677 0.088
5 0.353 0.059 0.368 6 12.394 0.550 0.080
6 0.398 0.071 0.361 6 7.857 0.412 0.082
7 0.450 0.080 0.339 6 4.098 0.311 0.085
8 0.467 0.092 0.279 6 1.642 0.255 0.089
# Main reading: the derivative detector at, for each recording, the lowest k at which the inverted
# trace gives at most 10% as many events as the trace (a false-positive bound chosen here); 5% and
# 20% are shown. A recording with no such k has no events counted.
def by_fp(order, bound, det="derivative"):
m = minis[(minis.mapping == order) & (minis.detector == det) & (minis.events > 1)]
ok = m[m["inverted/trace"] <= bound].sort_values("k").groupby("recording").head(1).set_index("recording")
return ok[["k", "events", "rate_hz", "inverted/trace", "mean_mV", "cv"]]
for order in ("numeric", "dir", "all"):
for bound in (0.05, 0.10, 0.20):
for det in ("derivative", "template"):
t = by_fp(order, bound, det)
print(f"[{order}; {det}; inverted/trace <= {bound:.2f}] {len(t)} recordings: mean size {t.mean_mV.mean():.3f}"
f" +- {t.mean_mV.std(ddof=1) / np.sqrt(len(t)):.3f} mV, mean CV {t.cv.mean():.3f}" if len(t) > 1 else
f"[{order}; {det}; inverted/trace <= {bound:.2f}] {len(t)} recordings")
MAIN_MINI = by_fp("numeric", 0.10)
print("\nMain reading per recording:")
print(MAIN_MINI.round(3).to_string())
rest = MAIN_MINI.drop(index="030401md02h")
print(f"Main reading without 030401md02h (offset and gain questioned): {len(rest)} recordings, mean size"
f" {rest.mean_mV.mean():.3f} +- {rest.mean_mV.std(ddof=1) / np.sqrt(len(rest)):.3f} mV, mean CV {rest.cv.mean():.3f}")
# Event-triggered averages and size histograms under the main reading.
fig, ax = plt.subplots(2, 6, figsize=(18, 5.5))
ranges = []
for j, name in enumerate(D["ttx"]):
rec = D["ttx"][name]
fs = 1.0 / rec["dx"]
ax[0, j].set_title(name, fontsize=9)
if name not in MAIN_MINI.index:
ax[0, j].text(0.5, 0.5, "no k passes", ha="center", transform=ax[0, j].transAxes)
continue
k = int(MAIN_MINI.loc[name, "k"])
sweeps = ac2.mini_sweeps(rec, "numeric")
seg = []
for i, e in EVENTS["numeric", name, "derivative", k].items():
x = ac2._smooth(sweeps[i], fs)
for o in e["onsets"]:
if o >= 40 and o + 120 < len(x):
seg.append(x[o - 40:o + 120] - x[o - 8:o].mean())
seg = np.array(seg)
t = (np.arange(-40, 120)) / fs * 1000
ax[0, j].plot(t, seg.mean(0), "k")
ax[0, j].fill_between(t, np.percentile(seg, 25, 0), np.percentile(seg, 75, 0), color="0.8")
ax[0, j].set(xlabel="ms from onset", ylabel="mV")
sizes = np.concatenate([e["size"] for e in EVENTS["numeric", name, "derivative", k].values()])
ax[1, j].hist(sizes, bins=20, color="0.4")
ax[1, j].set(xlabel=f"size (mV), k = {k}", ylabel="events")
ranges.append((name, len(seg), seg.mean(0).min(), seg.mean(0).max(), sizes.min(), sizes.max()))
plt.tight_layout()
plt.show()
for name, n_ev, lo, hi, smin, smax in ranges:
print(f"{name}: {n_ev} events averaged, average {lo:.2f} to {hi:.2f} mV; sizes {smin:.2f} to {smax:.2f} mV (all within the axes)")
[numeric; derivative; inverted/trace <= 0.05] 5 recordings: mean size 0.702 +- 0.220 mV, mean CV 0.314
[numeric; template; inverted/trace <= 0.05] 2 recordings: mean size 0.323 +- 0.051 mV, mean CV 0.308
[numeric; derivative; inverted/trace <= 0.10] 5 recordings: mean size 0.602 +- 0.164 mV, mean CV 0.354
[numeric; template; inverted/trace <= 0.10] 2 recordings: mean size 0.314 +- 0.060 mV, mean CV 0.315
[numeric; derivative; inverted/trace <= 0.20] 5 recordings: mean size 0.549 +- 0.134 mV, mean CV 0.376
[numeric; template; inverted/trace <= 0.20] 4 recordings: mean size 0.408 +- 0.128 mV, mean CV 0.328
[dir; derivative; inverted/trace <= 0.05] 5 recordings: mean size 0.771 +- 0.217 mV, mean CV 0.271
[dir; template; inverted/trace <= 0.05] 2 recordings: mean size 0.325 +- 0.050 mV, mean CV 0.307
[dir; derivative; inverted/trace <= 0.10] 5 recordings: mean size 0.726 +- 0.195 mV, mean CV 0.296
[dir; template; inverted/trace <= 0.10] 3 recordings: mean size 0.508 +- 0.196 mV, mean CV 0.426
[dir; derivative; inverted/trace <= 0.20] 5 recordings: mean size 0.613 +- 0.155 mV, mean CV 0.404
[dir; template; inverted/trace <= 0.20] 3 recordings: mean size 0.466 +- 0.195 mV, mean CV 0.513
[all; derivative; inverted/trace <= 0.05] 3 recordings: mean size 1.520 +- 0.815 mV, mean CV 0.379
[all; template; inverted/trace <= 0.05] 1 recordings
[all; derivative; inverted/trace <= 0.10] 4 recordings: mean size 0.752 +- 0.195 mV, mean CV 0.606
[all; template; inverted/trace <= 0.10] 1 recordings
[all; derivative; inverted/trace <= 0.20] 4 recordings: mean size 0.606 +- 0.137 mV, mean CV 0.631
[all; template; inverted/trace <= 0.20] 2 recordings: mean size 0.674 +- 0.441 mV, mean CV 0.551
Main reading per recording:
k events rate_hz inverted/trace mean_mV cv
recording
050801md01c 3 3280 19.362 0.033 0.202 0.402
050501md01 4 150 1.948 0.033 0.283 0.335
050301md01 5 69 1.120 0.072 0.884 0.367
050501md04 5 71 0.256 0.028 0.599 0.437
030401md02h 6 15 0.108 0.067 1.043 0.229
Main reading without 030401md02h (offset and gain questioned): 4 recordings, mean size 0.492 +- 0.156 mV, mean CV 0.385
030401md02h: 15 events averaged, average -0.06 to 0.99 mV; sizes 0.62 to 1.52 mV (all within the axes) 050301md01: 69 events averaged, average -0.02 to 0.82 mV; sizes 0.41 to 2.08 mV (all within the axes) 050501md01: 149 events averaged, average -0.01 to 0.23 mV; sizes 0.12 to 0.65 mV (all within the axes) 050501md04: 70 events averaged, average -0.03 to 0.55 mV; sizes 0.27 to 1.97 mV (all within the axes) 050801md01c: 3272 events averaged, average -0.00 to 0.15 mV; sizes 0.01 to 0.67 mV (all within the axes)
Result: partially reproduced (resting on the assumed detector and threshold). The traces are digitized in steps of 0.049 mV (0.012 mV for 030401md02h after its gain of 4), and 14% to 68% of the differences between samples are 0. Template matching at thresholds of 3 to 5 finds up to 29 events per s, and on the inverted traces 55% to 80% as many (mean over the six recordings), so most of those events are noise. Its mean size, 0.28 to 0.35 mV at those thresholds, does not measure mEPSPs. The derivative detector finds fewer events on the inverted traces. Its main reading takes, per recording, the lowest threshold at which the inverted trace gives at most 10% as many events (k = 3 to 6). That gives a mean size of 0.602 +- 0.164 mV over five recordings, outside 0.39 +- 0.08, and a mean CV of 0.354, within 0.36 +- 0.03. 050401md02 has one stable sweep (7.7 s), and no threshold keeps its false positives below the bound. The event-triggered averages of the main reading rise within about 2 ms and decay over 10 to 20 ms, peaking at 0.15 to 0.99 mV. The detected events are not locked to the tones played during the recordings: 6.7% to 12.0% of them fall in the 50 ms after a tone, against 10% by chance. The result depends on the detector, threshold and mapping of sweeps. With a fixed threshold, the derivative detector gives 0.370 mV with CV 0.421 at k = 3 (mean holds) and 0.487 mV with CV 0.350 at k = 4 (CV holds). False-positive bounds of 5% and 20% give 0.702 and 0.549 mV (CV 0.314 and 0.376). The documentation's file order (dir) gives 0.726 mV and CV 0.296 under the main rule. Under no reading with the false positives bounded do both values hold. The mean over recordings rests on the threshold each recording gets: event rates under the main reading range from 0.11 Hz (030401md02h, 15 events, 1.04 mV) to 19 Hz (050801md01c, 0.20 mV, sizes down to 0.01 mV). Without 030401md02h, whose offset and gain are questioned (data check), the four others give 0.492 mV and CV 0.385, both outside. The derivative detector with a false-positive bound is adopted as an assumption, because it is the only one of the claims file's two readings whose events are not mostly noise here.
recording-properties¶
Each tone was presented 18 times on average (10 to 63) per neuron; rest potential ranged from -58 to -92 mV (median -70) after junction potential correction, and input resistance from 19 to 141 MOhm (median 64).
props = []
for name, n in N.items():
st = ac2.stable(n)
sweeps = np.unique(n["sweep"][st])
ok = ~ac2.spikes(n, "cell") & st
trials_per_tone = pd.Series(n["tone"][ok]).value_counts()
dv = np.median(n["prepulse_dv"][sweeps])
props.append({"neuron": name, "reps": int(n["params"]["last_cycle"] - n["params"]["first_cycle"] + 1),
"mean trials per tone (no spikes)": trials_per_tone.reindex(ac2.TONE_ORDER).fillna(0).mean(),
"rest, sweeps (raw)": np.median(n["sweep_rest"][sweeps]),
"rest, 15-ms baselines (raw)": np.median(n["rbaseline"][st]),
"prepulse pA": n["info"]["prepulse_height"], "prepulse dV mV": dv,
"Rin MOhm": -dv / n["info"]["prepulse_height"] * 1000})
props = pd.DataFrame(props).set_index("neuron")
print(props.round(2).to_string())
r = props.reps
print(f"\nRepetitions: mean {r.mean():.2f} +- {r.std(ddof=1) / np.sqrt(len(r)):.2f} (SE), range {r.min()} to {r.max()}")
v = props["mean trials per tone (no spikes)"]
print(f"Variant, mean trials per tone without spike trials: mean {v.mean():.2f}, range {v.min():.1f} to {v.max():.1f}")
for jp in (12, 10):
for col in ("rest, sweeps (raw)", "rest, 15-ms baselines (raw)"):
x = props[col] - jp
print(f"Rest potential, {col.split(',')[1].split('(')[0].strip()}, corrected by -{jp} mV: median {x.median():.1f},"
f" range {x.max():.1f} to {x.min():.1f}; at or below -55 mV: {int((x <= -55).sum())} of {len(x)}")
x = props["rest, sweeps (raw)"] + 10
print(f"Documentation's direction (+10 mV, 'more depolarized'): median {x.median():.1f}, range {x.max():.1f} to {x.min():.1f}")
rin = props["Rin MOhm"]
print(f"Input resistance (last 20% of the prepulse against the {int(N[all_names[0]]['info']['prepulse_baselength'])} ms before it):"
f" median {rin.median():.1f}, range {rin.min():.1f} to {rin.max():.1f} MOhm")
# Variant: the step response includes the voltage drop across the access resistance (no bridge
# balance, "no on-line series resistance compensation"); subtract the recorded access resistance.
rin2 = rin - cells.R_access.astype(float)
print(f"Variant, minus the recorded access resistance: median {rin2.median():.1f}, range {rin2.min():.1f} to {rin2.max():.1f} MOhm;"
f" within 17.1-155.1: {int(rin2.between(17.1, 155.1).sum())} of {rin2.notna().sum()}")
print("Prepulse sizes by anesthetic:", props.groupby(cells.anesthetic)["prepulse pA"].unique().to_dict())
reps mean trials per tone (no spikes) rest, sweeps (raw) rest, 15-ms baselines (raw) prepulse pA prepulse dV mV Rin MOhm
neuron
020701md01b 63 63.00 -51.27 -51.25 100.0 -5.82 58.21
020901md01 34 34.00 -54.81 -54.66 100.0 -10.83 108.29
021001md01 13 12.97 -70.50 -70.45 100.0 -17.20 172.00
021401md01b 26 26.00 -53.59 -53.65 100.0 -13.79 137.94
021401md02 11 11.00 -46.26 -46.39 100.0 -6.98 69.79
021401md03 10 9.97 -64.82 -65.04 100.0 -23.75 237.46
021401md05 14 13.91 -51.33 -51.38 100.0 -15.63 156.27
021901md01 27 26.72 -61.71 -62.16 100.0 -18.39 183.86
022301md01 28 27.00 -54.93 -55.28 100.0 -43.83 438.26
022801md01b 20 20.00 -59.75 -59.49 100.0 -7.41 74.05
022801md02b 11 11.00 -72.02 -72.28 100.0 -7.51 75.10
030201md01b 20 18.94 -60.33 -60.19 100.0 -20.49 204.88
030201md03b 17 17.00 -50.54 -50.62 100.0 -5.91 59.14
030301md01b 14 14.00 -68.24 -68.12 100.0 -10.13 101.27
030301md02b 12 12.00 -54.20 -53.97 100.0 -12.40 124.04
030301md03b 16 15.88 -51.76 -51.65 100.0 -9.66 96.62
100300md01 16 16.00 -47.97 -49.29 200.0 -26.00 129.98
100300md02 10 9.41 -56.64 -59.69 200.0 -20.00 100.01
100500md01 11 10.53 -64.70 -67.18 200.0 -60.11 300.55
100500md02b 14 13.84 -62.50 -63.10 200.0 -40.44 202.21
102400md01 17 16.84 -45.65 -50.92 200.0 -30.20 151.02
102500md02 44 42.94 -52.37 -55.90 200.0 -31.10 155.49
102600md01 12 12.00 -56.82 -57.98 200.0 -24.68 123.40
102600md03 12 11.94 -70.01 -69.42 200.0 -25.64 128.20
102800md01 14 14.00 -58.84 -58.98 200.0 -55.45 277.24
102800md02 16 15.03 -54.81 -56.37 200.0 -17.47 87.33
102900md02 17 16.28 -57.19 -59.80 200.0 -45.25 226.23
102900md04 11 10.09 -43.27 -46.89 200.0 -40.17 200.85
110700md01 12 11.59 -54.69 -57.70 200.0 -52.41 262.05
110900md05 13 12.84 -58.96 -58.96 100.0 -9.29 92.93
111200md01 21 20.41 -61.77 -66.22 100.0 -20.03 200.32
111400md01 12 10.84 -54.93 -55.89 100.0 -26.14 261.43
111500md01 10 9.34 -75.87 -79.97 100.0 -13.66 136.61
Repetitions: mean 18.12 +- 1.93 (SE), range 10 to 63
Variant, mean trials per tone without spike trials: mean 17.80, range 9.3 to 63.0
Rest potential, sweeps, corrected by -12 mV: median -68.6, range -55.3 to -87.9; at or below -55 mV: 33 of 33
Rest potential, 15-ms baselines, corrected by -12 mV: median -70.0, range -58.4 to -92.0; at or below -55 mV: 33 of 33
Rest potential, sweeps, corrected by -10 mV: median -66.6, range -53.3 to -85.9; at or below -55 mV: 32 of 33
Rest potential, 15-ms baselines, corrected by -10 mV: median -68.0, range -56.4 to -90.0; at or below -55 mV: 33 of 33
Documentation's direction (+10 mV, 'more depolarized'): median -46.6, range -33.3 to -65.9
Input resistance (last 20% of the prepulse against the 5 ms before it): median 137.9, range 58.2 to 438.3 MOhm
Variant, minus the recorded access resistance: median 65.1, range 19.1 to 140.6 MOhm; within 17.1-155.1: 33 of 33
Prepulse sizes by anesthetic: {'P': array([100.]), 'P_with_Diaz': array([100.]), 'U': array([200., 100.])}
Result: reproduced (resting on the assumption that input resistance excludes the access resistance). Repetitions per tone: 18.12 +- 1.93 (SE), 10 to 63, as in the paper. Rest potential from the 15-ms baselines, corrected by -12 mV: median -70.0 mV, range -58.4 to -92.0, as in the paper; from the sweeps outside the prepulse, median -68.6 and range -55.3 to -87.9 (within the widened range). Input resistance from the prepulse is 58.2 to 438.3 MOhm (median 137.9) as measured; minus the recorded access resistance, which the uncompensated step includes, it is 19.1 to 140.6 MOhm (median 65.1), as in the paper (19 to 141, median 64). The step response alone does not match.
access-resistance-depth-no-effect¶
Response variability is not correlated with access resistance (4.2 against 4.9 mV on either side of the 60 MOhm median) and does not vary systematically with recording depth.
ra = cells.R_access.astype(float)
depth = cells.depth_um.astype(float)
med = ra.median()
print(f"Access resistance: median {med:.1f} MOhm, range {ra.min():.0f} to {ra.max():.0f}; neurons at the median: {list(ra.index[ra == med])}")
print(f"Depth (Sutter travel / sqrt 2): {depth.notna().sum()} of 33 recorded, {depth.min():.0f} to {depth.max():.0f} um (paper: 180 to 834)")
for label in ("main", SPK, SPK_ANY, SPK_SLOW):
q = S[label].q
for placement, above in (("median cell with the upper group", ra >= med), ("median cell with the lower group", ra > med)):
print(f"[{label}; {placement}] above median {mse(q[above])}; below {mse(q[~above])}")
for var, x in (("access resistance", ra), ("depth", depth)):
ok = x.notna()
pr, pp = stats.pearsonr(x[ok], q[ok])
sr, sp = stats.spearmanr(x[ok], q[ok])
print(f"[{label}] {var}: Pearson r = {pr:.3f}, p = {pp:.3f}; Spearman rho = {sr:.3f}, p = {sp:.3f} (n = {ok.sum()})")
Access resistance: median 60.0 MOhm, range 10 to 320; neurons at the median: ['021401md05'] Depth (Sutter travel / sqrt 2): 25 of 33 recorded, 180 to 834 um (paper: 180 to 834) [main; median cell with the upper group] above median 4.88 +- 0.71 (n = 17); below 4.27 +- 0.63 (n = 16) [main; median cell with the lower group] above median 4.81 +- 0.75 (n = 16); below 4.36 +- 0.60 (n = 17) [main] access resistance: Pearson r = 0.087, p = 0.631; Spearman rho = 0.155, p = 0.390 (n = 33) [main] depth: Pearson r = -0.273, p = 0.187; Spearman rho = -0.338, p = 0.098 (n = 25) [spikes: fast rise above -20 mV (6); median cell with the upper group] above median 5.70 +- 0.78 (n = 17); below 4.99 +- 0.75 (n = 16) [spikes: fast rise above -20 mV (6); median cell with the lower group] above median 5.68 +- 0.83 (n = 16); below 5.04 +- 0.71 (n = 17) [spikes: fast rise above -20 mV (6)] access resistance: Pearson r = 0.079, p = 0.661; Spearman rho = 0.133, p = 0.462 (n = 33) [spikes: fast rise above -20 mV (6)] depth: Pearson r = -0.264, p = 0.203; Spearman rho = -0.315, p = 0.125 (n = 25) [spikes: above -20 mV, any rise; median cell with the upper group] above median 4.50 +- 0.57 (n = 17); below 4.36 +- 0.66 (n = 16) [spikes: above -20 mV, any rise; median cell with the lower group] above median 4.41 +- 0.60 (n = 16); below 4.45 +- 0.63 (n = 17) [spikes: above -20 mV, any rise] access resistance: Pearson r = 0.041, p = 0.821; Spearman rho = 0.104, p = 0.564 (n = 33) [spikes: above -20 mV, any rise] depth: Pearson r = -0.255, p = 0.218; Spearman rho = -0.295, p = 0.152 (n = 25) [spikes: spikethresh, slow events kept; median cell with the upper group] above median 5.22 +- 0.73 (n = 17); below 4.99 +- 0.75 (n = 16) [spikes: spikethresh, slow events kept; median cell with the lower group] above median 5.18 +- 0.77 (n = 16); below 5.04 +- 0.71 (n = 17) [spikes: spikethresh, slow events kept] access resistance: Pearson r = 0.001, p = 0.995; Spearman rho = 0.031, p = 0.864 (n = 33) [spikes: spikethresh, slow events kept] depth: Pearson r = -0.263, p = 0.204; Spearman rho = -0.338, p = 0.098 (n = 25)
Result: reproduced (resting on the spike-rule assumption). The median access resistance is
60 MOhm (10 to 320). is 4.88 +- 0.71 mV for the 17 neurons at or above the median and 4.27
+- 0.63 mV for the 16 below it, both within the tolerances of 4.2 +- 0.7 and 4.9 +- 0.7. The
values fit the paper's numbers better the other way round: here the higher access resistance
group is the more variable, where the paper has it the less variable. Neither Pearson (r =
0.087, p = 0.631) nor Spearman (rho = 0.155, p = 0.390) finds a correlation with access
resistance, nor with depth over the 25 neurons with a recorded depth (r = -0.273, p = 0.187; rho
= -0.338, p = 0.098). With the claims file's spike rule the groups are 5.70 and 4.99 mV, the
first outside its tolerance, and no correlation is significant (p = 0.125 to 0.661), which would
not give reproduced. With the corrected potential above -20 mV at any rise the groups are 4.50
and 4.36 mV, both within, and no correlation is significant; with
spikethresh keeping the slow
events the upper group is 5.22 mV, outside.
Summary¶
| Claim | Tested | Outcome |
|---|---|---|
| fig2-stereotyped-shape | yes | partially reproduced |
| fig3a-reliable-neuron | yes | reproduced |
| fig3b-total-variability-distribution | yes | reproduced (assumed) |
| total-variability-population | yes | reproduced (assumed) |
| lfp-culling-reduction | yes | reproduced (assumed) |
| fig5b-residual-lfp-psp-correlation | yes | reproduced |
| fig7a-lfp-culling-estimation | yes | partially reproduced (assumed) |
| psp-culling-reduction | yes | partially reproduced (assumed) |
| fig4c-psp-culling-can-increase | yes | partially reproduced (assumed) |
| fig7b-private-best-estimate | yes | reproduced (assumed) |
| culling-controls | yes | not reproduced |
| fig7d-large-responses | yes | reproduced (assumed) |
| large-responses-controls | yes | not reproduced |
| cv-large-responses | yes | partially reproduced (assumed) |
| lfp-cv-not-different | yes | reproduced |
| anesthesia-comparison | yes | reproduced (assumed) |
| fig8-spontaneous-veto | yes | reproduced |
| overfitting-correction-example | yes | reproduced (assumed) |
| private-shared-ranges | yes | partially reproduced |
| quantal-comparison | yes | reproduced |
| mepsp-size | yes | partially reproduced (assumed) |
| recording-properties | yes | reproduced (assumed) |
| access-resistance-depth-no-effect | yes | reproduced (assumed) |
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 links no code. Trials are each neuron's
stable epoch as given in the data set's cell_specific_parameters.m (its length matches the
paper's repetitions: mean 18.1, 10 to 63). Spike trials are found with the data set's
per-neuron threshold, an assumption. The claims file's rule (a fast rise, taken here as above 10
mV/ms, to above -20 mV after the junction potential correction) misses slow depolarizations to
about -10 mV, which may be regenerative events or large PSPs, broad spikes that cross -20 mV with
a rise below 10 mV/ms, and narrow spikes that peak below -20 mV; it raises from 4.58 to 5.35
mV. Without the rise criterion the -20 mV rule gives 4.43 mV, and
spikethresh keeping the slow
events 5.11 mV. All four rules are shown for every claim. The peak of the mean PSP is
searched over the whole trace after onset; the peak-to-trough height for culling is taken over the
culling window; negative scale factors are removed in both cullings; the tone criteria use the
mean before culling; heights are measured once before culling, and re-measuring them after
culling is shown as a reading (of the claims it is run for, it changes only fig4c, to
reproduced; the random culling controls keep the heights of the full ensemble). LFP-based culling follows
PSP-based culling and is run on the trials it keeps, as the Figure 6 legend says; deciding both
cullings on the full ensemble is shown as a reading. The overfitting correction averages 20
shuffles. Each open detail is run under the readings the claims file lists, and the tables above
show where they change an outcome. Input resistance is taken as the prepulse response minus the
recorded access resistance. mEPSPs are detected with a threshold on the derivative, set per
recording by a bound on the events found in the inverted trace; template matching and fixed
thresholds are shown. The panels drawn as Figures 2, 3A, 3B and 4 match the paper's by eye; the
paper labels its Figure 2 and 3A neuron 022801md01b, the one identified here, and its Figure 4
neuron 030301md03b.
What was found about the data set. It holds the 33 neurons, the 29 LFP flags, the anesthetic
groups and the six TTX recordings the paper describes, and with the per-neuron stable epochs and
spike thresholds of cell_specific_parameters.m the paper's population values follow closely
( 4.58, 4.80 with LFP, 2.19 after PSP-based culling, 3.19 after LFP-based culling, 1.48 for
the best private estimate). Those files are not described for these neurons. The documentation
gives a junction potential of about 10 mV and says the true potentials are about 10 mV more
depolarized than recorded; the paper subtracts 12 mV, and only that reproduces its rest
potentials. It calls the prepulse height a voltage in mV, while the values are currents in pA.
Depth is missing for 8 of 33 neurons. The TTX recording 030401md02h, after its documented
division by 4, sits near +2.6 mV, which is not a membrane potential. Which sweep file of a TTX
recording is trial k of its
first_trial and last_trial is not documented, and the
documentation's own reading code lists the files in name order, which picks different sweeps
from the numeric order in four of the six recordings.