Eyherabide, Rokem, Herz & Samengo (2008) on CRCNS ia-1¶
Paper: H. G. Eyherabide, A. Rokem, A. V. M. Herz and I. Samengo (2008), "Burst firing is a neural
code in an insect auditory system", Frontiers in Computational Neuroscience,
doi:10.3389/neuro.10.003.2008. The claims file
(claims/crcns-ia-1/10.3389-neuro.10.003.2008.yaml) was written from the version of record in
PubMed Central (PMC2525941), with its Methods and the legends of Figures 1 to 10 and Table 1; the
paper has no supplement and links no code.
Dataset: CRCNS data set ia-1, doi:10.6080/K0BG2KWB, read from
the public bucket crcnsarchive (ia-1/crcns-ia1.zip), streamed and opened in place.
Version: CRCNS data sets are not versioned, so the version is recorded as unversioned and the
file read is pinned by its MD5 sum in the data set's checksums.md5: crcns-ia1.zip
f0740514fabc6995fe54b2a361b6c763 (the helper checks the object's ETag against it before reading;
the object was last modified 2026-02-11). This is the only published state of the data set, so it
is the one pinned. Only Data1 is used: the data descriptor names it as the data of this paper;
Data2 (2004) is unpublished and used other stimuli.
| Claim | Tested here |
|---|---|
| data-sessions-and-trials | tested |
| table1-sessions-per-protocol | tested |
| firing-rate-uniform | tested |
| threshold-population | blocked: the thresholds are not in the data set |
| table1-bursting-sessions | tested |
| table1-isolated-spikes | tested |
| table1-burst-length | tested |
| fig5-burst-probability-vs-stimulus | tested |
| fig3-correlation-shape | tested (the data part; the model panels need the thresholds) |
| method-correlation-error | tested |
| model-matches-classification | blocked: the model needs the thresholds s_TH |
| fig6-bursts-less-correlated | tested |
| fig7ab-information-per-burst-increases | tested |
| fig7cd-occurrence-decreases | tested |
| fig7ef-jitter-independent-of-n | tested |
| fig7gh-information-rate | tested |
| direct-method-comparison | tested |
| fig8-nbta-differ | tested |
| fig8-latency-decreases | tested |
| fig8e-sigma-minimum | tested |
| sigma-n-below-sta-sigma | tested |
| fig9-nbta-not-concatenated-1bta | tested (data only; the model sessions need the thresholds) |
| fig9c-population-nbta | tested |
| fig10-burst-n-vs-height | tested |
| firing-rate-no-effect-on-limit | tested |
| model-reproduces-coding-trends | blocked: depends on model-matches-classification |
Two analyses draw random numbers, the shifted-onset baseline of Figure 7 and the size-matched control of sigma-n-below-sta-sigma; each has its own seeded generator.
Revised after the first verification round. The main population is now selection (c), which also
leaves out sessions with fewer than 98 trials and sessions sorted off-line (see "Session
selection"); the first round's population, selection (a), is shown beside it. The central maximum
of the burst-triggered averages is searched before the burst onset; a first-peak reading with the
refractory plateau left out, a shifted-onset baseline for the information measures, a size-matched
control for sigma_n, sensitivity to the minimum number of bursts, and checks of the trial blocks and
of the firing_rates file were added.
import collections
import itertools
import re
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from scipy import stats
import ia1
pd.set_option("display.width", 200)
pd.set_option("display.max_columns", 40)
pd.set_option("display.max_rows", 200)
print("data set", ia1.CRCNS_ID, ia1.VERSION, "doi", ia1.DATASET_DOI, "zip md5", ia1.ZIP_MD5)
z = ia1.open_zip()
sessions = ia1.load_sessions(z)
cells = ia1.load_singlecells(z)
shapes = {name: ia1.load_stimulus_shape(z, name) for name in sorted({s.stimulus for s in sessions})}
cell_of = {e: i for i, row in enumerate(cells) for e in row}
data set ia-1 unversioned doi 10.6080/K0BG2KWB zip md5 f0740514fabc6995fe54b2a361b6c763
Data check¶
What was loaded, and the checks of the description against the values. The data descriptor
says: spike times in microseconds from the stimulus trigger; trials of the first 1000 ms of the
stimulus, 98 to 533 per stimulus; 700 ms pauses; two interleaved stimuli per experiment;
cell.dat gives the stimulus names and intensities (the maximal intensity, mean plus 3 SD);
stimulus files at 20 kHz, 10 s long; singlecells.txt lists the experiments made on one cell.
rows = []
for s in sessions:
r = ia1.trial_rates(s.trials)
rows.append(dict(session=s.name, exp=s.exp, date=s.exp[:8], cell=cell_of.get(s.exp, -1),
stimulus=s.stimulus, protocol=s.protocol, cutoff=s.cutoff, intensity=s.intensity,
header_intensity=float(s.header["intensity (dB)"]) if "intensity (dB)" in s.header else np.nan,
desired_rate=float(s.header["desired firing rate (Hz)"]) if "desired firing rate (Hz)" in s.header else np.nan,
pause_ms=float(s.header["pause between stimuli (msec)"]) if "pause between stimuli (msec)" in s.header else np.nan,
header_stim=s.header.get("stimulus file", "").split("/")[-1].removesuffix(".dat"),
offline=s.offline, n_trials=s.n_trials, rate=r.mean(), rate_sd=r.std(ddof=1),
rate_full=ia1.trial_rates(s.trials, 0, 1).mean(), rate_sd_full=ia1.trial_rates(s.trials, 0, 1).std(ddof=1)))
D = pd.DataFrame(rows)
print("Data1 experiments:", D.exp.nunique(), " sessions:", len(D), " offline-sorted sessions (no OEL header):", D.offline.sum())
print("experiments not listed in singlecells.txt:", sorted(set(D.exp) - set(cell_of)))
print("cells (rows of singlecells.txt, plus experiments not listed):", len(cells) + len(set(D.exp) - set(cell_of)),
" recording dates:", D.date.nunique())
print("\nSessions by stimulus file (cell.dat):")
print(D.groupby("stimulus").agg(sessions=("session", "size"), protocol=("protocol", "first")).to_string())
hdr = D[D.header_stim != ""]
print("\nHeader stimulus agrees with cell.dat in", (hdr.header_stim.str.replace("gauss_", "") == hdr.stimulus.str.replace("gauss_", "")).sum(),
"of", len(hdr), "sessions with a stimulus line in the header")
d = (D.header_intensity - D.intensity).dropna()
print("header intensity minus cell.dat intensity (dB): n =", len(d), " median |diff|", round(d.abs().median(), 2),
" max |diff|", round(d.abs().max(), 2), " sessions over 1 dB:", list(D.session[(D.header_intensity - D.intensity).abs() > 1]))
print("desired firing rate in headers:", D.desired_rate.value_counts(dropna=False).to_dict())
print("pause between stimuli (ms) in headers:", D.pause_ms.value_counts(dropna=False).to_dict())
allsp = np.concatenate([np.concatenate(s.trials) for s in sessions])
print("\nspike times (s): min", allsp.min(), " max", allsp.max(), " fraction after 1 s (in the pause):", round((allsp > 1).mean(), 4))
print("trials per session: min", D.n_trials.min(), " max", D.n_trials.max())
print("rate over 200-1000 ms (Hz): min", round(D.rate.min(), 1), " max", round(D.rate.max(), 1))
unsorted = sum(int((np.diff(t) <= 0).any()) for s in sessions for t in s.trials)
print("trials with unsorted or repeated spike times:", unsorted)
isi = np.concatenate([np.diff(t) for s in sessions for t in s.trials])
print("shortest ISI (ms):", round(isi.min() * 1e3, 2), " 0.1 ms quantile of ISIs (ms):", round(np.quantile(isi, 0.001) * 1e3, 2))
for name, sh in shapes.items():
print(f"stimulus {name}: {len(sh)} samples in the first 1.05 s at 20 kHz, dB re its 10 s maximum, range {sh.min():.1f} to {sh.max():.1f}, SD over 0-1 s {sh[:20000].std():.2f} dB")
# Duplicate recordings: identical spike trains or identical first trials.
seen = collections.defaultdict(list)
for s in sessions:
seen[np.concatenate(s.trials).round(6).tobytes()].append(s.name)
seen[("first", s.trials[0].round(6).tobytes())].append(s.name)
print("sessions with identical spike times or identical first trials:", [v for v in seen.values() if len(v) > 1])
# Trial blocks: spikes after the 1 s stimulus fall in the pause. A block holding two trials (a
# missing header) would restart its spike times near 0, so its times would not be sorted (counted
# above), and would hold about twice the spike count of the session's other blocks.
rows = []
for s in sessions:
n01 = np.array([np.sum(t < 1) for t in s.trials])
mx = max(t.max() for t in s.trials if len(t))
rows.append(dict(session=s.name, offline=s.offline, max_t=mx, late_blocks=sum(len(t) > 0 and t.max() > 1.7 for t in s.trials),
rate_0_1=n01.mean(), rate_pause=np.mean([np.sum((t >= 1) & (t < 1.7)) for t in s.trials]) / 0.7,
rate_after_1_7=np.mean([np.sum(t >= 1.7) for t in s.trials]) / max(1e-9, mx - 1.7),
doubled=int((n01 > 1.5 * np.median(n01)).sum())))
blk = pd.DataFrame(rows)
lb = blk[blk.max_t > 1.2]
print("\nsessions with spikes later than 1.2 s:", len(lb), "(off-line sorted:", int(lb.offline.sum()), "of", int(blk.offline.sum()), "); with blocks reaching past 1.7 s:",
int((lb.late_blocks > 0).sum()), "; latest spike", round(blk.max_t.max(), 3), "s; latest spike in on-line sessions", round(blk.max_t[~blk.offline].max(), 3), "s")
print(" in those sessions, rate 0-1 s median", round(lb.rate_0_1.median(), 1), "Hz; 1-1.7 s median", round(lb.rate_pause.median(), 1),
"Hz (max", round(lb.rate_pause.max(), 1), "); after 1.7 s median", round(lb.rate_after_1_7.median(), 1), "Hz (max", round(lb.rate_after_1_7.max(), 1), ")")
print(" blocks with more than 1.5 times the session's median spike count in 0-1 s:", int(blk.doubled.sum()), "in", list(blk.session[blk.doubled > 0]),
"; none over 2 times:", all(max(np.sum(t < 1) for t in s.trials) < 2 * np.median([np.sum(t < 1) for t in s.trials]) for s in sessions))
# Data1/firing_rates: one record per spike file, its number of trials and a rate with a +/- value.
frt = z.read(ia1.ROOT + "firing_rates").decode("latin1")
FR = pd.DataFrame([dict(session=f"{e}_sig{g}", fr_stim=st_, fr_n=int(n), fr_rate=float(r), fr_pm=float(pm)) for e, g, st_, n, r, pm in
re.findall(r"\./\S+/(\S+)_sig(\d)_spikes\.dat\s*\n(\S+)\s*\n\(n=(\d+)\) firing rate = ([\d.]+) \+/- ([\d.]+)", frt)])
FRM = FR.merge(D[["session", "n_trials", "rate_full", "rate_sd", "offline", "protocol"]], on="session")
print("\nfiring_rates: records", len(FR), "of", len(D), "spike files; trials equal to the blocks parsed here in", int((FRM.fr_n == FRM.n_trials).sum()),
"; differing:", FRM.loc[FRM.fr_n != FRM.n_trials, ["session", "fr_n", "n_trials"]].values.tolist())
print(" rate / rate over 0-1 s: median", round((FRM.fr_rate / FRM.rate_full).median(), 3), "(on-line sessions", round((FRM.fr_rate / FRM.rate_full)[~FRM.offline].median(), 3),
"; off-line", round((FRM.fr_rate / FRM.rate_full)[FRM.offline].median(), 3), ")")
print(" +/- value: median", round(FR.fr_pm.median(), 1), "Hz, range", FR.fr_pm.min(), "to", FR.fr_pm.max(), "; trial-to-trial SD of the rate here: median",
round(D.rate_sd.median(), 1), "Hz; seven-protocol sessions without a record:", sorted(set(D.session[D.protocol.notna()]) - set(FR.session)))
Data1 experiments: 75 sessions: 150 offline-sorted sessions (no OEL header): 28
experiments not listed in singlecells.txt: ['03-01-15-aa', '03-06-27-ae']
cells (rows of singlecells.txt, plus experiments not listed): 45 recording dates: 25
Sessions by stimulus file (cell.dat):
sessions protocol
stimulus
gaussModality_co200 14 NaN
gauss_st12_co200 8 3.0
gauss_st3_co200 9 2.0
gauss_st6_co100 11 5.0
gauss_st6_co200 54 1.0
gauss_st6_co25 9 4.0
gauss_st6_co400 8 6.0
gauss_st6_co800 9 7.0
lmd_co200 14 NaN
st6_co1600 3 NaN
st6_co200 7 1.0
st6_co25 1 4.0
st6_co50 3 NaN
Header stimulus agrees with cell.dat in 116 of 116 sessions with a stimulus line in the header
header intensity minus cell.dat intensity (dB): n = 122 median |diff| 0.04 max |diff| 6.58 sessions over 1 dB: ['03-04-11-aa_sig2', '03-04-16-ae_sig2', '03-04-23-ae_sig1']
desired firing rate in headers: {100.0: 118, nan: 28, 90.0: 4}
pause between stimuli (ms) in headers: {700.0: 116, nan: 28, 450.0: 6}
spike times (s): min 0.0 max 2.763 fraction after 1 s (in the pause): 0.0172
trials per session: min 69 max 503
rate over 200-1000 ms (Hz): min 65.7 max 169.1
trials with unsorted or repeated spike times: 0
shortest ISI (ms): 1.0 0.1 ms quantile of ISIs (ms): 1.6
stimulus gaussModality_co200: 21000 samples in the first 1.05 s at 20 kHz, dB re its 10 s maximum, range -39.3 to 0.0, SD over 0-1 s 6.41 dB
stimulus gauss_st12_co200: 21000 samples in the first 1.05 s at 20 kHz, dB re its 10 s maximum, range -72.0 to -0.4, SD over 0-1 s 11.94 dB
stimulus gauss_st3_co200: 21000 samples in the first 1.05 s at 20 kHz, dB re its 10 s maximum, range -18.0 to 0.0, SD over 0-1 s 2.76 dB
stimulus gauss_st6_co100: 21000 samples in the first 1.05 s at 20 kHz, dB re its 10 s maximum, range -35.8 to 0.0, SD over 0-1 s 5.62 dB
stimulus gauss_st6_co200: 21000 samples in the first 1.05 s at 20 kHz, dB re its 10 s maximum, range -35.2 to 0.0, SD over 0-1 s 6.15 dB
stimulus gauss_st6_co25: 21000 samples in the first 1.05 s at 20 kHz, dB re its 10 s maximum, range -33.0 to 0.0, SD over 0-1 s 5.83 dB
stimulus gauss_st6_co400: 21000 samples in the first 1.05 s at 20 kHz, dB re its 10 s maximum, range -36.0 to 0.0, SD over 0-1 s 5.91 dB
stimulus gauss_st6_co800: 21000 samples in the first 1.05 s at 20 kHz, dB re its 10 s maximum, range -36.0 to 0.0, SD over 0-1 s 5.91 dB
stimulus lmd_co200: 21000 samples in the first 1.05 s at 20 kHz, dB re its 10 s maximum, range -38.5 to -0.2, SD over 0-1 s 6.85 dB
stimulus st6_co1600: 21000 samples in the first 1.05 s at 20 kHz, dB re its 10 s maximum, range -36.0 to 0.0, SD over 0-1 s 5.96 dB
stimulus st6_co200: 21000 samples in the first 1.05 s at 20 kHz, dB re its 10 s maximum, range -35.2 to 0.0, SD over 0-1 s 6.15 dB
stimulus st6_co25: 21000 samples in the first 1.05 s at 20 kHz, dB re its 10 s maximum, range -33.0 to 0.0, SD over 0-1 s 5.83 dB
stimulus st6_co50: 21000 samples in the first 1.05 s at 20 kHz, dB re its 10 s maximum, range -36.0 to -1.2, SD over 0-1 s 5.75 dB
sessions with identical spike times or identical first trials: []
sessions with spikes later than 1.2 s: 25 (off-line sorted: 25 of 28 ); with blocks reaching past 1.7 s: 24 ; latest spike 2.763 s; latest spike in on-line sessions 1.02 s in those sessions, rate 0-1 s median 116.2 Hz; 1-1.7 s median 2.7 Hz (max 25.0 ); after 1.7 s median 0.9 Hz (max 19.7 ) blocks with more than 1.5 times the session's median spike count in 0-1 s: 6 in ['03-04-16-ad_sig1', '03-04-16-ad_sig2'] ; none over 2 times: True firing_rates: records 138 of 150 spike files; trials equal to the blocks parsed here in 136 ; differing: [['03-04-04-aa_sig1', 80, 81], ['03-04-04-aa_sig2', 79, 80]] rate / rate over 0-1 s: median 0.997 (on-line sessions 0.998 ; off-line 0.428 ) +/- value: median 127.9 Hz, range 56.73 to 188.94 ; trial-to-trial SD of the rate here: median 5.8 Hz; seven-protocol sessions without a record: ['03-06-27-ab_sig1', '03-06-27-ac_sig1', '03-06-27-ae_sig1', '03-06-27-af_sig1', '03-06-27-ag_sig1', '03-07-03-aa_sig1']
Trial blocks. Spikes reach 2.76 s, later than one trial and its pause, so a block without a
header could hold two trials. The check above shows that late spikes occur only in off-line
sorted sessions, at the low rate of the pause, and that no block restarts its times or holds
twice the spike count of the others: each block is one trial and the trial counts stand (the data
set's own firing_rates file gives the same counts in 136 of 138 files). The firing_rates file
is not a selection: its +/- value is far above any trial-to-trial SD of the rate.
The stimulus files are named for their SD and cutoff (gauss_st6_co200 is SD 6 dB, cutoff 200
Hz). Their SDs over the first second can be checked against the names (6 dB files near 6 dB).
Clock check. Spikes and stimuli share the trial's clock if the spike-triggered average of the stimulus (STA) shows the upstroke a few ms before the spike, which the paper describes (Figure 1: "responses typically appear after stimulus upstrokes with an delay of 3.4 msec"). The contrast is the same STA computed from the stimulus shifted by 100 ms, which keeps every property of the stimulus and the spike train except their alignment.
P = [s for s in sessions if s.protocol] # sessions in the paper's seven protocols
def stim_db(s):
return shapes[s.stimulus] + s.intensity
def spikes_in(s, t0=ia1.T0, t1=ia1.T1):
return np.concatenate([t[(t >= t0) & (t < t1)] for t in s.trials])
rows = []
for s in P:
st = stim_db(s)
sp = spikes_in(s)
seg, tt = ia1.segments(st, sp)
sta = seg.mean(0) - st[4000:20000].mean()
seg2, _ = ia1.segments(st, sp - 0.1)
sta2 = seg2.mean(0) - st[4000:20000].mean()
rows.append(dict(session=s.name, protocol=s.protocol, sta_peak=sta.max(), lat_ms=-tt[np.argmax(sta[tt <= 0])] * 1e3,
shifted_peak=sta2.max()))
C = pd.DataFrame(rows)
print("STA peak above the mean stimulus (dB): median", round(C.sta_peak.median(), 2), " range", round(C.sta_peak.min(), 2), "to", round(C.sta_peak.max(), 2))
print("with the stimulus shifted by 100 ms: median", round(C.shifted_peak.median(), 2), " range", round(C.shifted_peak.min(), 2), "to", round(C.shifted_peak.max(), 2))
print("sessions whose STA peak exceeds the shifted one:", (C.sta_peak > C.shifted_peak).sum(), "of", len(C))
print("latency of the STA maximum (ms): median", C.lat_ms.median(), " range", C.lat_ms.min(), "to", C.lat_ms.max())
print("sessions where the true STA peak is under 1 dB:", list(C.session[C.sta_peak < 1]))
print("sessions whose STA maximum is at the spike (latency 0 ms, the edge of the search):", list(C.session[C.lat_ms <= 0]))
STA peak above the mean stimulus (dB): median 6.02 range 1.53 to 15.71 with the stimulus shifted by 100 ms: median 0.47 range 0.16 to 1.85 sessions whose STA peak exceeds the shifted one: 116 of 116 latency of the STA maximum (ms): median 5.6499999999999995 range -0.0 to 9.950000000000001 sessions where the true STA peak is under 1 dB: [] sessions whose STA maximum is at the spike (latency 0 ms, the edge of the search): ['03-03-06-aa_sig2']
Session selection and burst identification¶
Kept sessions (claim data-sessions-and-trials). The paper keeps sessions whose trial-to-trial SD
of the firing rate is below 35 Hz and whose mean rate is about 100 Hz. Readings, from the
criterion: (a) the SD condition alone; (b) the SD condition and a target rate of 100 Hz where
the header records one. Per-trial rates over 200 to 1000 ms, and over the whole 1 s as a
variant. Only sessions of the paper's seven protocols (Table 1) can be kept; the data set's other
stimuli (gaussModality_co200, lmd_co200, cutoffs 50 and 1600 Hz) are not in the paper's design
(Methods: "The other signal could be one of six different stimulation protocols").
Readings (a) and (b) keep 115 and 113 of the 116 seven-protocol sessions, against the paper's
86, so the paper removed sessions by a rule it does not state. Reading (c), added after the first
verification, is the SD condition and two exclusions the data suggest: sessions with fewer than
98 trials (the Methods give "N ranging between 98 and 503" for the paper's sessions, so a session
with fewer was not among them) and sessions whose spikes were detected off-line (OFFOEL
header; the data descriptor says off-line detection was done where the on-line detection "was
deemed to be inexact"). The off-line exclusion has no passage in the paper behind it; it is taken
because, with it, the kept sessions per protocol and the bursting sessions per protocol come
close to Table 1 (below). Because it was chosen in part from the counts of Table 1, it does not
decide data-sessions-and-trials, table1-sessions-per-protocol or table1-bursting-sessions, where it
is shown as a variant. For the other claims it is the main population, an assumption; the results
over the 115 sessions of reading (a) are shown beside it.
Bursts, by the Methods (Burst identification), under every combination of the readings of the
open details (claims file, opening comment): peak (a) first local maximum, (b) largest value in
0 < tau < 5 ms, and (c) the first local maximum after C_s first rises significantly above its
refractory plateau (added after the first verification, see ia1.py); minimum (a) first
significant local minimum after the peak, (b) lowest value up to 1.25 / cutoff; error bar N (a)
the number of product terms, (b) the number of trials, and (c) sigma the SD over trials of each
trial's C_s with N the number of trials (added here, see ia1.py).
def kept(reading="a", window="200-1000"):
col = "rate_sd" if window == "200-1000" else "rate_sd_full"
out = []
for s in P:
row = D.loc[D.session == s.name].iloc[0]
ok = row[col] < 35
if reading == "b" and not np.isnan(row.desired_rate):
ok = ok and row.desired_rate == 100
if reading == "c":
ok = ok and not s.offline and s.n_trials >= 98
if ok:
out.append(s)
return out
T1_SESS = [43, 6, 7, 7, 8, 7, 8]
sv = D[D.protocol.notna() & (D.rate_sd < 35)]
rows = []
for label, g in [("SD < 35 Hz (reading a)", sv), (" without off-line sorted sessions", sv[~sv.offline]),
(" without sessions under 98 trials", sv[sv.n_trials >= 98]),
(" without both (reading c)", sv[~sv.offline & (sv.n_trials >= 98)]),
(" without the 450 ms pause sessions", sv[sv.pause_ms != 450]),
(" without off-line and 450 ms pause", sv[~sv.offline & (sv.pause_ms != 450)]),
(" without off-line, 450 ms pause and under 98 trials", sv[~sv.offline & (sv.pause_ms != 450) & (sv.n_trials >= 98)])]:
v = [int((g.protocol == p).sum()) for p in range(1, 8)]
rows.append(dict(selection=label, sessions=len(g), **{f"p{p}": x for p, x in zip(range(1, 8), v)},
within_10pct=sum(abs(x - e) <= 0.1 * e for x, e in zip(v, T1_SESS)),
min_trials=g.n_trials.min(), mean_trials=round(g.n_trials.mean(), 1)))
rows.append(dict(selection="paper (Table 1)", sessions=86, **{f"p{p}": x for p, x in zip(range(1, 8), T1_SESS)}))
print(pd.DataFrame(rows).to_string(index=False))
print("off-line sorted seven-protocol sessions by protocol:", sv[sv.offline].protocol.value_counts().sort_index().to_dict(),
"; sessions under 98 trials:", sv.loc[sv.n_trials < 98, ["session", "n_trials", "offline"]].values.tolist())
KEPT = kept("c")
KEPT_A = kept("a")
POP = {"c": {s.name for s in KEPT}, "a": {s.name for s in KEPT_A}}
COMBOS = list(itertools.product("abc", "ab", "abc"))
MAIN = ("b", "a", "a")
limits = {c: {s.name: ia1.burst_limit(s.trials, s.cutoff, *c)[0] for s in P} for c in COMBOS}
for c in COMBOS:
lim = limits[c]
nb = sum(lim[s.name] > 0 for s in KEPT)
na = sum(lim[s.name] > 0 for s in KEPT_A)
print("peak %s, minimum %s, N %s: bursting %3d of %d kept (c), %3d of %d (a); median limiting ISI %s ms" % (
*c, nb, len(KEPT), na, len(KEPT_A), np.round(np.median([lim[s.name] for s in KEPT if lim[s.name] > 0]) * 1e3, 2) if nb else "-"))
LIM = limits[MAIN]
BURSTING = [s for s in KEPT if LIM[s.name] > 0]
BURSTING_A = [s for s in KEPT_A if LIM[s.name] > 0]
print("\nmain reading", MAIN, ": kept", len(KEPT), " bursting", len(BURSTING), "; reading (a) of the selection: kept", len(KEPT_A), " bursting", len(BURSTING_A))
selection sessions p1 p2 p3 p4 p5 p6 p7 within_10pct min_trials mean_trials
SD < 35 Hz (reading a) 115 60 9 8 10 11 8 9 0.0 69.0 160.3
without off-line sorted sessions 93 49 7 7 7 8 8 7 3.0 80.0 165.2
without sessions under 98 trials 106 56 8 8 8 10 7 9 1.0 98.0 167.0
without both (reading c) 89 47 6 7 7 8 7 7 6.0 98.0 168.7
without the 450 ms pause sessions 109 57 8 8 10 11 7 8 2.0 69.0 161.2
without off-line and 450 ms pause 87 46 6 7 7 8 7 6 6.0 80.0 166.6
without off-line, 450 ms pause and under 98 trials 83 44 5 7 7 8 6 6 4.0 98.0 170.4
paper (Table 1) 86 43 6 7 7 8 7 8 NaN NaN NaN
off-line sorted seven-protocol sessions by protocol: {1.0: 11, 2.0: 2, 3.0: 1, 4.0: 3, 5.0: 3, 7.0: 2} ; sessions under 98 trials: [['03-03-06-aa_sig2', 97, True], ['03-03-06-ab_sig1', 70, True], ['03-03-06-ab_sig2', 69, True], ['03-03-14-aa_sig1', 76, True], ['03-03-14-aa_sig2', 75, True], ['03-04-04-aa_sig1', 81, False], ['03-04-04-aa_sig2', 80, False], ['03-04-25-ah_sig1', 93, False], ['03-04-25-ah_sig2', 92, False]]
peak a, minimum a, N a: bursting 11 of 89 kept (c), 15 of 115 (a); median limiting ISI 2.9 ms
peak a, minimum a, N b: bursting 0 of 89 kept (c), 0 of 115 (a); median limiting ISI - ms
peak a, minimum a, N c: bursting 11 of 89 kept (c), 15 of 115 (a); median limiting ISI 2.9 ms
peak a, minimum b, N a: bursting 11 of 89 kept (c), 15 of 115 (a); median limiting ISI 3.8 ms
peak a, minimum b, N b: bursting 0 of 89 kept (c), 0 of 115 (a); median limiting ISI - ms
peak a, minimum b, N c: bursting 11 of 89 kept (c), 15 of 115 (a); median limiting ISI 3.8 ms
peak b, minimum a, N a: bursting 64 of 89 kept (c), 84 of 115 (a); median limiting ISI 3.55 ms
peak b, minimum a, N b: bursting 0 of 89 kept (c), 0 of 115 (a); median limiting ISI - ms
peak b, minimum a, N c: bursting 64 of 89 kept (c), 84 of 115 (a); median limiting ISI 3.55 ms
peak b, minimum b, N a: bursting 64 of 89 kept (c), 84 of 115 (a); median limiting ISI 4.6 ms
peak b, minimum b, N b: bursting 0 of 89 kept (c), 0 of 115 (a); median limiting ISI - ms
peak b, minimum b, N c: bursting 64 of 89 kept (c), 84 of 115 (a); median limiting ISI 4.6 ms
peak c, minimum a, N a: bursting 61 of 89 kept (c), 82 of 115 (a); median limiting ISI 3.5 ms
peak c, minimum a, N b: bursting 0 of 89 kept (c), 0 of 115 (a); median limiting ISI - ms
peak c, minimum a, N c: bursting 61 of 89 kept (c), 82 of 115 (a); median limiting ISI 3.5 ms
peak c, minimum b, N a: bursting 61 of 89 kept (c), 83 of 115 (a); median limiting ISI 4.4 ms
peak c, minimum b, N b: bursting 0 of 89 kept (c), 0 of 115 (a); median limiting ISI - ms
peak c, minimum b, N c: bursting 61 of 89 kept (c), 83 of 115 (a); median limiting ISI 4.4 ms
main reading ('b', 'a', 'a') : kept 89 bursting 64 ; reading (a) of the selection: kept 115 bursting 84
data-sessions-and-trials¶
Thirty-seven receptor cells were recorded from 23 animals, giving 132 sessions, of which 86 met the firing-rate and stationarity conditions; each stimulus was repeated 98 to 503 times, 172 on average.
Data1 does not record the animal. The recording date is shown as a stand-in (one animal per
day is an assumption the descriptor does not make); cells are the rows of singlecells.txt.
def describe(sel, label):
names = {s.name for s in sel}
d = D[D.session.isin(names)]
# cells: rows of singlecells.txt, plus each experiment it does not list
return dict(population=label, cells=d.cell[d.cell >= 0].nunique() + d.exp[d.cell < 0].nunique(),
dates=d.date.nunique(), sessions=len(d), min_trials=d.n_trials.min(), max_trials=d.n_trials.max(),
mean_trials=round(d.n_trials.mean(), 1), mean_rate_sd=round(d.rate_sd.mean(), 2),
mean_rate_sd_full=round(d.rate_sd_full.mean(), 2))
tab = pd.DataFrame([describe(sessions, "all Data1 sessions"),
describe([s for s in sessions if s.stimulus.replace("gauss_", "").startswith("st")], "Gaussian AM stimuli (any cutoff)"),
describe(P, "seven protocols"),
describe(kept("a"), "kept (a), 200-1000 ms"), describe(kept("b"), "kept (b), 200-1000 ms"),
describe(kept("a", "full"), "kept (a), 0-1000 ms"), describe(kept("b", "full"), "kept (b), 0-1000 ms"),
describe(kept("c"), "kept (c), 200-1000 ms"), describe(kept("c", "full"), "kept (c), 0-1000 ms")])
print(tab.to_string(index=False))
print("\nsessions in the seven protocols with rate SD >= 35 Hz (200-1000 ms):",
list(D.session[D.protocol.notna() & (D.rate_sd >= 35)]))
print("SD of the trial rates over the seven-protocol sessions (Hz): quantiles 0, 0.25, 0.5, 0.75, 0.9, 1:",
np.quantile(D.rate_sd[D.protocol.notna()], [0, .25, .5, .75, .9, 1]).round(1))
expected = dict(cells=(37, 3.7), sessions_recorded=(132, 13.2), kept=(86, 8.6), min_trials=(98, 9.8), max_trials=(503, 50.3),
mean_trials=(172, 17.2), mean_rate_sd=(6, 0.6))
print("\nexpected (value, tolerance):", expected)
population cells dates sessions min_trials max_trials mean_trials mean_rate_sd mean_rate_sd_full
all Data1 sessions 45 25 150 69 503 155.9 7.08 6.66
Gaussian AM stimuli (any cutoff) 37 23 122 69 503 157.2 6.75 6.29
seven protocols 37 23 116 69 503 159.8 6.78 6.30
kept (a), 200-1000 ms 37 23 115 69 503 160.3 6.48 6.05
kept (b), 200-1000 ms 36 23 113 69 503 161.3 6.46 6.04
kept (a), 0-1000 ms 37 23 116 69 503 159.8 6.78 6.30
kept (b), 0-1000 ms 36 23 114 69 503 160.7 6.76 6.29
kept (c), 200-1000 ms 29 18 89 98 503 168.7 6.00 5.72
kept (c), 0-1000 ms 29 18 90 98 503 167.9 6.39 6.04
sessions in the seven protocols with rate SD >= 35 Hz (200-1000 ms): ['03-04-16-ae_sig1']
SD of the trial rates over the seven-protocol sessions (Hz): quantiles 0, 0.25, 0.5, 0.75, 0.9, 1: [ 3. 4.5 5.6 7.4 9.5 40.9]
expected (value, tolerance): {'cells': (37, 3.7), 'sessions_recorded': (132, 13.2), 'kept': (86, 8.6), 'min_trials': (98, 9.8), 'max_trials': (503, 50.3), 'mean_trials': (172, 17.2), 'mean_rate_sd': (6, 0.6)}
Result: partially reproduced. The 116 sessions of the seven protocols come from 37 cells recorded on 23 dates, the paper's 37 cells and 23 animals (if one animal was recorded per date, an assumption, as Data1 does not record the animal). The 122 Gaussian AM sessions (any cutoff) are within 10% of the paper's 132 sessions; all 150 Data1 sessions are not. Under the criterion's readings the selection does not reproduce: the SD < 35 Hz condition keeps 115 of the 116 sessions (113 with the 100 Hz target of reading (b)), against 86, and over them the trials range from 69 (outside 98 +/- 9.8) to 503 and average 160.3 (within 172 +/- 17.2), with a mean trial-to-trial SD of the rate of 6.48 Hz (within 6 +/- 0.6). Reading (c), which also leaves out the sessions with fewer than 98 trials and those sorted off-line, keeps 89 sessions with 98 to 503 trials (mean 168.7) and a mean SD of 6.00 Hz, within tolerance, but they come from 29 cells recorded on 18 dates, outside 37 +/- 3.7 and 23 +/- 2.3; the 37 cells and 23 dates hold only when counted over all seven-protocol sessions. It adds two exclusions the paper does not give as conditions, chosen in part from these counts (its fewest trials are 98 by construction), so it does not decide this claim; it is used as the population of the claims that follow.
table1-sessions-per-protocol¶
The 86 kept sessions divide among the seven stimulus protocols as 43 (6 dB, 200 Hz), 6 (3 dB, 200 Hz), 7 (12 dB, 200 Hz), 7 (6 dB, 25 Hz), 8 (6 dB, 100 Hz), 7 (6 dB, 400 Hz) and 8 (6 dB, 800 Hz).
rows = []
for label, sel in [("kept (a) 200-1000", kept("a")), ("kept (b) 200-1000", kept("b")),
("kept (a) 0-1000", kept("a", "full")), ("kept (b) 0-1000", kept("b", "full")),
("kept (c) 200-1000", kept("c")), ("kept (c) 0-1000", kept("c", "full")), ("seven protocols", P)]:
cnt = collections.Counter(s.protocol for s in sel)
vals = [cnt[p] for p in range(1, 8)]
rows.append(dict(population=label, **{f"p{p}": v for p, v in zip(range(1, 8), vals)},
within_10pct=sum(abs(v - e) <= 0.1 * e for v, e in zip(vals, T1_SESS))))
rows.append(dict(population="paper", **{f"p{p}": v for p, v in zip(range(1, 8), T1_SESS)}, within_10pct=7))
print(pd.DataFrame(rows).to_string(index=False))
exps = collections.defaultdict(list)
for s in P:
exps[s.exp].append(s.protocol)
print("\nexperiments with both stimuli in the protocols, by pair:",
dict(collections.Counter(tuple(sorted(v)) for v in exps.values() if len(v) == 2)))
population p1 p2 p3 p4 p5 p6 p7 within_10pct
kept (a) 200-1000 60 9 8 10 11 8 9 0
kept (b) 200-1000 58 9 8 10 11 8 9 0
kept (a) 0-1000 61 9 8 10 11 8 9 0
kept (b) 0-1000 59 9 8 10 11 8 9 0
kept (c) 200-1000 47 6 7 7 8 7 7 6
kept (c) 0-1000 48 6 7 7 8 7 7 5
seven protocols 61 9 8 10 11 8 9 0
paper 43 6 7 7 8 7 8 7
experiments with both stimuli in the protocols, by pair: {(1, 3): 8, (1, 7): 9, (1, 4): 10, (1, 2): 9, (1, 6): 8, (1, 5): 11}
Result: not reproduced. None of the seven counts is within 10% of Table 1 under the criterion's readings: the kept sessions are 60, 9, 8, 10, 11, 8 and 9 per protocol (reading (a), 200 to 1000 ms) against 43, 6, 7, 7, 8, 7 and 8. Under reading (c) of the selection they are 47, 6, 7, 7, 8, 7 and 7, six of seven within tolerance (800 Hz has 7 against 8), but that reading was chosen in part from these counts. The data set holds every reference session paired with each other protocol (55 experiments) and more reference sessions than the paper counts.
firing-rate-uniform¶
The kept sessions had a mean firing rate of 113 Hz (SD 16 Hz) that did not differ significantly across stimulus conditions (one-way ANOVA, p = 0.58).
rows = []
for label, sel, col in [("kept (c), 200-1000 ms", kept("c"), "rate"), ("kept (c), 0-1000 ms", kept("c", "full"), "rate_full"),
("kept (a), 200-1000 ms", kept("a"), "rate"), ("kept (a), 0-1000 ms", kept("a", "full"), "rate_full"),
("all seven-protocol sessions, 200-1000 ms", P, "rate"), ("all seven-protocol sessions, 0-1000 ms", P, "rate_full"),
("kept (b), 200-1000 ms", kept("b"), "rate")]:
d = D[D.session.isin({s.name for s in sel})]
groups = [g[col].values for _, g in d.groupby("protocol")]
F, p = stats.f_oneway(*groups)
rows.append(dict(population=label, n=len(d), mean_rate=round(d[col].mean(), 1), sd_rate=round(d[col].std(ddof=1), 1),
F=round(F, 2), p=round(p, 3)))
print(pd.DataFrame(rows).to_string(index=False))
d = D[D.session.isin({s.name for s in KEPT})]
print("\nmean rate by protocol (kept (c), 200-1000 ms):", d.groupby("protocol").rate.mean().round(1).to_dict())
population n mean_rate sd_rate F p
kept (c), 200-1000 ms 89 107.1 14.5 0.49 0.814
kept (c), 0-1000 ms 90 111.6 14.2 0.82 0.560
kept (a), 200-1000 ms 115 110.2 16.7 0.49 0.818
kept (a), 0-1000 ms 116 114.7 16.9 0.94 0.471
all seven-protocol sessions, 200-1000 ms 116 110.3 16.6 0.50 0.809
all seven-protocol sessions, 0-1000 ms 116 114.7 16.9 0.94 0.471
kept (b), 200-1000 ms 113 110.4 16.8 0.50 0.805
mean rate by protocol (kept (c), 200-1000 ms): {1.0: 107.5, 2.0: 115.3, 3.0: 104.8, 4.0: 102.4, 5.0: 104.9, 6.0: 107.3, 7.0: 107.2}
Result: reproduced. Over the 89 sessions of selection (c) the mean rate is 107.1 Hz (within 113 +/- 16) with SD 14.5 Hz (within 16 +/- 1.6), and a one-way ANOVA across the seven protocols finds no difference (F = 0.49, p = 0.814). Over the 115 sessions of selection (a) the mean is 110.2 Hz (SD 16.7), F = 0.49, p = 0.818; over all seven-protocol sessions and with whole-trial rates the ANOVA gives p between 0.471 and 0.809. No population shows a significant difference.
threshold-population (blocked)¶
The mean threshold of the cells at their best frequency was 58 dB (SD 14 dB).
Data1 holds spike times, cell.dat (stimulus names and the intensities they were played at) and
the stimuli. It holds no threshold, best-frequency or tuning (f-I) measurement: the listing below
shows every kind of file in Data1. Blocked, data_missing.
kinds = collections.Counter()
for n in z.namelist():
if n.startswith(ia1.ROOT) and not n.endswith("/"):
base = n.split("/")[-1]
kinds[base.split("_")[-1] if "_" in base and n.count("/") > 3 and "stimuli" not in n else ("stimuli/" + base if "stimuli" in n else base)] += 1
print(dict(kinds))
print("cell.dat keys:", sorted({line.split(":")[0].strip() for n in z.namelist() if n.endswith("cell.dat")
for line in z.read(n).decode("latin1").splitlines() if ":" in line}))
{'.DS_Store': 70, '03-01-15-aa_sig1_long.dat': 1, '03-01-15-aa_sig1_spikes.dat': 1, '03-01-15-aa_sig2_long.dat': 1, '03-01-15-aa_sig2_spikes.dat': 1, 'cell.dat': 75, '03-01-15-ab_sig1_long.dat': 1, '03-01-15-ab_sig1_spikes.dat': 1, '03-01-15-ab_sig2_long.dat': 1, '03-01-15-ab_sig2_spikes.dat': 1, '03-01-16-ab_sig1_long.dat': 1, '03-01-16-ab_sig1_spikes.dat': 1, '03-01-16-ab_sig2_long.dat': 1, '03-01-16-ab_sig2_spikes.dat': 1, '03-01-30-ad_sig1_long.dat': 1, '03-01-30-ad_sig1_spikes.dat': 1, '03-01-30-ad_sig2_long.dat': 1, '03-01-30-ad_sig2_spikes.dat': 1, '03-01-30-ae_sig1_long.dat': 1, '03-01-30-ae_sig1_spikes.dat': 1, '03-01-30-ae_sig2_long.dat': 1, '03-01-30-ae_sig2_spikes.dat': 1, '03-02-14-aa_sig1_long.dat': 1, '03-02-14-aa_sig1_spikes.dat': 1, '03-02-14-aa_sig2_long.dat': 1, '03-02-14-aa_sig2_spikes.dat': 1, '03-02-18-ab_sig1_long.dat': 1, '03-02-18-ab_sig1_spikes.dat': 1, '03-02-18-ab_sig2_long.dat': 1, '03-02-18-ab_sig2_spikes.dat': 1, '03-02-25-ac_sig1_long.dat': 1, '03-02-25-ac_sig1_spikes.dat': 1, '03-02-25-ac_sig2_long.dat': 1, '03-02-25-ac_sig2_spikes.dat': 1, '03-02-25-aj_sig1_long.dat': 1, '03-02-25-aj_sig1_spikes.dat': 1, '03-02-25-aj_sig2_long.dat': 1, '03-02-25-aj_sig2_spikes.dat': 1, '03-03-05-ab_sig1_long.dat': 1, '03-03-05-ab_sig1_spikes.dat': 1, '03-03-05-ab_sig2_long.dat': 1, '03-03-05-ab_sig2_spikes.dat': 1, '03-03-05-ac_sig1_long.dat': 1, '03-03-05-ac_sig1_spikes.dat': 1, '03-03-05-ac_sig2_long.dat': 1, '03-03-05-ac_sig2_spikes.dat': 1, '03-03-06-aa_sig1_long.dat': 1, '03-03-06-aa_sig1_spikes.dat': 1, '03-03-06-aa_sig2_long.dat': 1, '03-03-06-aa_sig2_spikes.dat': 1, '03-03-06-ab_sig1_long.dat': 1, '03-03-06-ab_sig1_spikes.dat': 1, '03-03-06-ab_sig2_long.dat': 1, '03-03-06-ab_sig2_spikes.dat': 1, '03-03-14-aa_sig1_long.dat': 1, '03-03-14-aa_sig1_spikes.dat': 1, '03-03-14-aa_sig2_long.dat': 1, '03-03-14-aa_sig2_spikes.dat': 1, '03-03-14-ab_sig1_long.dat': 1, '03-03-14-ab_sig1_spikes.dat': 1, '03-03-14-ab_sig2_long.dat': 1, '03-03-14-ab_sig2_spikes.dat': 1, '03-03-14-ac_sig1_long.dat': 1, '03-03-14-ac_sig1_spikes.dat': 1, '03-03-14-ac_sig2_long.dat': 1, '03-03-14-ac_sig2_spikes.dat': 1, '03-03-14-ad_sig1_long.dat': 1, '03-03-14-ad_sig1_spikes.dat': 1, '03-03-14-ad_sig2_long.dat': 1, '03-03-14-ad_sig2_spikes.dat': 1, '03-03-18-aa_sig1_long.dat': 1, '03-03-18-aa_sig1_spikes.dat': 1, '03-03-18-aa_sig2_long.dat': 1, '03-03-18-aa_sig2_spikes.dat': 1, '03-03-18-ab_sig1_long.dat': 1, '03-03-18-ab_sig1_spikes.dat': 1, '03-03-18-ab_sig2_long.dat': 1, '03-03-18-ab_sig2_spikes.dat': 1, '03-03-18-ac_sig1_long.dat': 1, '03-03-18-ac_sig1_spikes.dat': 1, '03-03-18-ac_sig2_long.dat': 1, '03-03-18-ac_sig2_spikes.dat': 1, '03-03-18-ad_sig1_long.dat': 1, '03-03-18-ad_sig1_spikes.dat': 1, '03-03-18-ad_sig2_long.dat': 1, '03-03-18-ad_sig2_spikes.dat': 1, '03-03-18-ae_sig1_long.dat': 1, '03-03-18-ae_sig1_spikes.dat': 1, '03-03-18-ae_sig2_long.dat': 1, '03-03-18-ae_sig2_spikes.dat': 1, '03-03-21-aa_sig1_long.dat': 1, '03-03-21-aa_sig1_spikes.dat': 1, '03-03-21-aa_sig2_long.dat': 1, '03-03-21-aa_sig2_spikes.dat': 1, '03-03-21-ab_sig1_long.dat': 1, '03-03-21-ab_sig1_spikes.dat': 1, '03-03-21-ab_sig2_long.dat': 1, '03-03-21-ab_sig2_spikes.dat': 1, '03-03-21-ac_sig1_long.dat': 1, '03-03-21-ac_sig1_spikes.dat': 1, '03-03-21-ac_sig2_long.dat': 1, '03-03-21-ac_sig2_spikes.dat': 1, '03-04-02-aa_sig1_long.dat': 1, '03-04-02-aa_sig1_spikes.dat': 1, '03-04-02-aa_sig2_long.dat': 1, '03-04-02-aa_sig2_spikes.dat': 1, '03-04-04-aa_sig1_long.dat': 1, '03-04-04-aa_sig1_spikes.dat': 1, '03-04-04-aa_sig2_long.dat': 1, '03-04-04-aa_sig2_spikes.dat': 1, '03-04-08-aa_sig1_long.dat': 1, '03-04-08-aa_sig1_spikes.dat': 1, '03-04-08-aa_sig2_long.dat': 1, '03-04-08-aa_sig2_spikes.dat': 1, '03-04-08-ab_sig1_long.dat': 1, '03-04-08-ab_sig1_spikes.dat': 1, '03-04-08-ab_sig2_long.dat': 1, '03-04-08-ab_sig2_spikes.dat': 1, '03-04-08-ac_sig1_long.dat': 1, '03-04-08-ac_sig1_spikes.dat': 1, '03-04-08-ac_sig2_long.dat': 1, '03-04-08-ac_sig2_spikes.dat': 1, '03-04-11-aa_sig1_long.dat': 1, '03-04-11-aa_sig1_spikes.dat': 1, '03-04-11-aa_sig2_long.dat': 1, '03-04-11-aa_sig2_spikes.dat': 1, '03-04-11-ab_sig1_long.dat': 1, '03-04-11-ab_sig1_spikes.dat': 1, '03-04-11-ab_sig2_long.dat': 1, '03-04-11-ab_sig2_spikes.dat': 1, '03-04-15-aa_sig1_long.dat': 1, '03-04-15-aa_sig1_spikes.dat': 1, '03-04-15-aa_sig2_long.dat': 1, '03-04-15-aa_sig2_spikes.dat': 1, '03-04-16-aa_sig1_long.dat': 1, '03-04-16-aa_sig1_spikes.dat': 1, '03-04-16-aa_sig2_long.dat': 1, '03-04-16-aa_sig2_spikes.dat': 1, '03-04-16-ab_sig1_long.dat': 1, '03-04-16-ab_sig1_spikes.dat': 1, '03-04-16-ab_sig2_long.dat': 1, '03-04-16-ab_sig2_spikes.dat': 1, '03-04-16-ac_sig1_long.dat': 1, '03-04-16-ac_sig1_spikes.dat': 1, '03-04-16-ac_sig2_long.dat': 1, '03-04-16-ac_sig2_spikes.dat': 1, '03-04-16-ad_sig1_long.dat': 1, '03-04-16-ad_sig1_spikes.dat': 1, '03-04-16-ad_sig2_long.dat': 1, '03-04-16-ad_sig2_spikes.dat': 1, '03-04-16-ae_sig1_long.dat': 1, '03-04-16-ae_sig1_spikes.dat': 1, '03-04-16-ae_sig2_long.dat': 1, '03-04-16-ae_sig2_spikes.dat': 1, '03-04-16-af_sig1_long.dat': 1, '03-04-16-af_sig1_spikes.dat': 1, '03-04-16-af_sig2_long.dat': 1, '03-04-16-af_sig2_spikes.dat': 1, '03-04-16-ag_sig1_long.dat': 1, '03-04-16-ag_sig1_spikes.dat': 1, '03-04-16-ag_sig2_long.dat': 1, '03-04-16-ag_sig2_spikes.dat': 1, '03-04-23-aa_sig1_long.dat': 1, '03-04-23-aa_sig1_spikes.dat': 1, '03-04-23-aa_sig2_long.dat': 1, '03-04-23-aa_sig2_spikes.dat': 1, '03-04-23-ab_sig1_long.dat': 1, '03-04-23-ab_sig1_spikes.dat': 1, '03-04-23-ab_sig2_long.dat': 1, '03-04-23-ab_sig2_spikes.dat': 1, '03-04-23-ac_sig1_long.dat': 1, '03-04-23-ac_sig1_spikes.dat': 1, '03-04-23-ac_sig2_long.dat': 1, '03-04-23-ac_sig2_spikes.dat': 1, '03-04-23-ad_sig1_long.dat': 1, '03-04-23-ad_sig1_spikes.dat': 1, '03-04-23-ad_sig2_long.dat': 1, '03-04-23-ad_sig2_spikes.dat': 1, '03-04-23-ae_sig1_long.dat': 1, '03-04-23-ae_sig1_spikes.dat': 1, '03-04-23-ae_sig2_long.dat': 1, '03-04-23-ae_sig2_spikes.dat': 1, '03-04-23-af_sig1_long.dat': 1, '03-04-23-af_sig1_spikes.dat': 1, '03-04-23-af_sig2_long.dat': 1, '03-04-23-af_sig2_spikes.dat': 1, '03-04-23-ag_sig1_long.dat': 1, '03-04-23-ag_sig1_spikes.dat': 1, '03-04-23-ag_sig2_long.dat': 1, '03-04-23-ag_sig2_spikes.dat': 1, '03-04-23-ah_sig1_long.dat': 1, '03-04-23-ah_sig1_spikes.dat': 1, '03-04-23-ah_sig2_long.dat': 1, '03-04-23-ah_sig2_spikes.dat': 1, '03-04-23-ai_sig1_long.dat': 1, '03-04-23-ai_sig1_spikes.dat': 1, '03-04-23-ai_sig2_long.dat': 1, '03-04-23-ai_sig2_spikes.dat': 1, '03-04-23-aj_sig1_long.dat': 1, '03-04-23-aj_sig1_spikes.dat': 1, '03-04-23-aj_sig2_long.dat': 1, '03-04-23-aj_sig2_spikes.dat': 1, '03-04-25-ac_sig1_long.dat': 1, '03-04-25-ac_sig1_spikes.dat': 1, '03-04-25-ac_sig2_long.dat': 1, '03-04-25-ac_sig2_spikes.dat': 1, '03-04-25-ad_sig1_long.dat': 1, '03-04-25-ad_sig1_spikes.dat': 1, '03-04-25-ad_sig2_long.dat': 1, '03-04-25-ad_sig2_spikes.dat': 1, '03-04-25-ae_sig1_long.dat': 1, '03-04-25-ae_sig1_spikes.dat': 1, '03-04-25-ae_sig2_long.dat': 1, '03-04-25-ae_sig2_spikes.dat': 1, '03-04-25-af_sig1_long.dat': 1, '03-04-25-af_sig1_spikes.dat': 1, '03-04-25-af_sig2_long.dat': 1, '03-04-25-af_sig2_spikes.dat': 1, '03-04-25-ag_sig1_long.dat': 1, '03-04-25-ag_sig1_spikes.dat': 1, '03-04-25-ag_sig2_long.dat': 1, '03-04-25-ag_sig2_spikes.dat': 1, '03-04-25-ah_sig1_long.dat': 1, '03-04-25-ah_sig1_spikes.dat': 1, '03-04-25-ah_sig2_long.dat': 1, '03-04-25-ah_sig2_spikes.dat': 1, '03-04-30-aa_sig1_long.dat': 1, '03-04-30-aa_sig1_spikes.dat': 1, '03-04-30-aa_sig2_long.dat': 1, '03-04-30-aa_sig2_spikes.dat': 1, '03-04-30-ac_sig1_long.dat': 1, '03-04-30-ac_sig1_spikes.dat': 1, '03-04-30-ac_sig2_long.dat': 1, '03-04-30-ac_sig2_spikes.dat': 1, '03-05-02-aa_sig1_long.dat': 1, '03-05-02-aa_sig1_spikes.dat': 1, '03-05-02-aa_sig2_long.dat': 1, '03-05-02-aa_sig2_spikes.dat': 1, '03-05-02-ab_sig1_long.dat': 1, '03-05-02-ab_sig1_spikes.dat': 1, '03-05-02-ab_sig2_long.dat': 1, '03-05-02-ab_sig2_spikes.dat': 1, '03-05-02-ac_sig1_long.dat': 1, '03-05-02-ac_sig1_spikes.dat': 1, '03-05-02-ac_sig2_long.dat': 1, '03-05-02-ac_sig2_spikes.dat': 1, '03-05-02-ad_sig1_long.dat': 1, '03-05-02-ad_sig1_spikes.dat': 1, '03-05-02-ad_sig2_long.dat': 1, '03-05-02-ad_sig2_spikes.dat': 1, '03-05-02-ae_sig1_long.dat': 1, '03-05-02-ae_sig1_spikes.dat': 1, '03-05-02-ae_sig2_long.dat': 1, '03-05-02-ae_sig2_spikes.dat': 1, '03-05-06-aa_sig1_long.dat': 1, '03-05-06-aa_sig1_spikes.dat': 1, '03-05-06-aa_sig2_long.dat': 1, '03-05-06-aa_sig2_spikes.dat': 1, '03-05-06-ab_sig1_long.dat': 1, '03-05-06-ab_sig1_spikes.dat': 1, '03-05-06-ab_sig2_long.dat': 1, '03-05-06-ab_sig2_spikes.dat': 1, '03-05-06-ac_sig1_long.dat': 1, '03-05-06-ac_sig1_spikes.dat': 1, '03-05-06-ac_sig2_long.dat': 1, '03-05-06-ac_sig2_spikes.dat': 1, '03-05-09-ab_sig1_long.dat': 1, '03-05-09-ab_sig1_spikes.dat': 1, '03-05-09-ab_sig2_long.dat': 1, '03-05-09-ab_sig2_spikes.dat': 1, '03-05-09-ac_long.dat': 1, '03-05-09-ac_sig1_long.dat': 1, '03-05-09-ac_sig1_spikes.dat': 1, '03-05-09-ac_sig2_spikes.dat': 1, '03-05-09-ad_sig1_long.dat': 1, '03-05-09-ad_sig1_spikes.dat': 1, '03-05-09-ad_sig2_long.dat': 1, '03-05-09-ad_sig2_spikes.dat': 1, '03-06-27-ab_sig1_long.dat': 1, '03-06-27-ab_sig1_spikes.dat': 1, '03-06-27-ab_sig2_long.dat': 1, '03-06-27-ab_sig2_spikes.dat': 1, '03-06-27-ac_sig1_long.dat': 1, '03-06-27-ac_sig1_spikes.dat': 1, '03-06-27-ac_sig2_long.dat': 1, '03-06-27-ac_sig2_spikes.dat': 1, '03-06-27-ae_sig1_long.dat': 1, '03-06-27-ae_sig1_spikes.dat': 1, '03-06-27-ae_sig2_long.dat': 1, '03-06-27-ae_sig2_spikes.dat': 1, '03-06-27-af_sig1_long.dat': 1, '03-06-27-af_sig1_spikes.dat': 1, '03-06-27-af_sig2_long.dat': 1, '03-06-27-af_sig2_spikes.dat': 1, '03-06-27-ag_sig1_long.dat': 1, '03-06-27-ag_sig1_spikes.dat': 1, '03-06-27-ag_sig2_long.dat': 1, '03-06-27-ag_sig2_spikes.dat': 1, '03-07-03-aa_sig1_long.dat': 1, '03-07-03-aa_sig1_spikes.dat': 1, '03-07-03-aa_sig2_long.dat': 1, '03-07-03-aa_sig2_spikes.dat': 1, 'firing_rates': 1, 'singlecells.txt': 1, 'stimuli/.DS_Store': 1, 'stimuli/gauss_st12_co200.dat': 1, 'stimuli/gauss_st3_co200.dat': 1, 'stimuli/gauss_st6_co100.dat': 1, 'stimuli/gauss_st6_co1600.dat': 1, 'stimuli/gauss_st6_co200.dat': 1, 'stimuli/gauss_st6_co25.dat': 1, 'stimuli/gauss_st6_co400.dat': 1, 'stimuli/gauss_st6_co50.dat': 1, 'stimuli/gauss_st6_co800.dat': 1, 'stimuli/gaussModality_co200.dat': 1, 'stimuli/lmd_co200.dat': 1, 'stimuli/stimtodB.m': 1}
cell.dat keys: ['Cell name', 'Cut off frequency', 'Intensity', 'Name', 'same cell as']
table1-bursting-sessions¶
The number of sessions with bursts of n > 1 depends on the stimulus, from none for 3 dB and for 800 Hz to all 7 for 12 dB and 7 of 8 for 100 Hz; 40 of 43 reference sessions, 6 of 7 at 25 Hz and 3 of 7 at 400 Hz.
Counted over the kept sessions, under every combination of the burst-identification readings. The fraction of sessions is shown beside the counts, since the counts depend on how many sessions are kept.
T1_BURST = [40, 0, 7, 6, 7, 3, 0]
tol = [max(0.1 * e, 0) for e in T1_BURST]
def n_gt1_sessions(sel, lim):
cnt = collections.Counter()
for s in sel:
if lim[s.name] > 0 and (ia1.burst_counts(s.trials, lim[s.name])[1:] > 0).any():
cnt[s.protocol] += 1
return [cnt[p] for p in range(1, 8)]
rows = []
for kr in ["c", "a", "b"]:
sel = kept(kr)
tot = collections.Counter(s.protocol for s in sel)
for c in COMBOS:
v = n_gt1_sessions(sel, limits[c])
rows.append(dict(kept=kr, readings="".join(c), **{f"p{p}": f"{x}/{tot[p]}" for p, x in zip(range(1, 8), v)},
within=sum(abs(x - e) <= t for x, e, t in zip(v, T1_BURST, tol))))
print(pd.DataFrame(rows).to_string(index=False))
print("paper:", [f"{a}/{b}" for a, b in zip(T1_BURST, T1_SESS)])
kept readings p1 p2 p3 p4 p5 p6 p7 within c aaa 5/47 0/6 2/7 1/7 3/8 0/7 0/7 2 c aab 0/47 0/6 0/7 0/7 0/8 0/7 0/7 2 c aac 5/47 0/6 2/7 1/7 3/8 0/7 0/7 2 c aba 5/47 0/6 2/7 1/7 3/8 0/7 0/7 2 c abb 0/47 0/6 0/7 0/7 0/8 0/7 0/7 2 c abc 5/47 0/6 2/7 1/7 3/8 0/7 0/7 2 c baa 40/47 1/6 7/7 6/7 7/8 3/7 0/7 6 c bab 0/47 0/6 0/7 0/7 0/8 0/7 0/7 2 c bac 40/47 1/6 7/7 6/7 7/8 3/7 0/7 6 c bba 40/47 1/6 7/7 6/7 7/8 3/7 0/7 6 c bbb 0/47 0/6 0/7 0/7 0/8 0/7 0/7 2 c bbc 40/47 1/6 7/7 6/7 7/8 3/7 0/7 6 c caa 37/47 1/6 7/7 6/7 7/8 3/7 0/7 6 c cab 0/47 0/6 0/7 0/7 0/8 0/7 0/7 2 c cac 37/47 1/6 7/7 6/7 7/8 3/7 0/7 6 c cba 37/47 1/6 7/7 6/7 7/8 3/7 0/7 6 c cbb 0/47 0/6 0/7 0/7 0/8 0/7 0/7 2 c cbc 37/47 1/6 7/7 6/7 7/8 3/7 0/7 6 a aaa 6/60 0/9 2/8 2/10 5/11 0/8 0/9 2 a aab 0/60 0/9 0/8 0/10 0/11 0/8 0/9 2 a aac 6/60 0/9 2/8 2/10 5/11 0/8 0/9 2 a aba 6/60 0/9 2/8 2/10 5/11 0/8 0/9 2 a abb 0/60 0/9 0/8 0/10 0/11 0/8 0/9 2 a abc 6/60 0/9 2/8 2/10 5/11 0/8 0/9 2 a baa 52/60 2/9 8/8 9/10 10/11 3/8 0/9 2 a bab 0/60 0/9 0/8 0/10 0/11 0/8 0/9 2 a bac 52/60 2/9 8/8 9/10 10/11 3/8 0/9 2 a bba 52/60 2/9 8/8 9/10 10/11 3/8 0/9 2 a bbb 0/60 0/9 0/8 0/10 0/11 0/8 0/9 2 a bbc 52/60 2/9 8/8 9/10 10/11 3/8 0/9 2 a caa 49/60 3/9 8/8 9/10 10/11 3/8 0/9 2 a cab 0/60 0/9 0/8 0/10 0/11 0/8 0/9 2 a cac 49/60 3/9 8/8 9/10 10/11 3/8 0/9 2 a cba 49/60 3/9 8/8 9/10 10/11 3/8 1/9 1 a cbb 0/60 0/9 0/8 0/10 0/11 0/8 0/9 2 a cbc 49/60 3/9 8/8 9/10 10/11 3/8 1/9 1 b aaa 6/58 0/9 2/8 2/10 5/11 0/8 0/9 2 b aab 0/58 0/9 0/8 0/10 0/11 0/8 0/9 2 b aac 6/58 0/9 2/8 2/10 5/11 0/8 0/9 2 b aba 6/58 0/9 2/8 2/10 5/11 0/8 0/9 2 b abb 0/58 0/9 0/8 0/10 0/11 0/8 0/9 2 b abc 6/58 0/9 2/8 2/10 5/11 0/8 0/9 2 b baa 50/58 2/9 8/8 9/10 10/11 3/8 0/9 2 b bab 0/58 0/9 0/8 0/10 0/11 0/8 0/9 2 b bac 50/58 2/9 8/8 9/10 10/11 3/8 0/9 2 b bba 50/58 2/9 8/8 9/10 10/11 3/8 0/9 2 b bbb 0/58 0/9 0/8 0/10 0/11 0/8 0/9 2 b bbc 50/58 2/9 8/8 9/10 10/11 3/8 0/9 2 b caa 47/58 3/9 8/8 9/10 10/11 3/8 0/9 2 b cab 0/58 0/9 0/8 0/10 0/11 0/8 0/9 2 b cac 47/58 3/9 8/8 9/10 10/11 3/8 0/9 2 b cba 47/58 3/9 8/8 9/10 10/11 3/8 1/9 1 b cbb 0/58 0/9 0/8 0/10 0/11 0/8 0/9 2 b cbc 47/58 3/9 8/8 9/10 10/11 3/8 1/9 1 paper: ['40/43', '0/6', '7/7', '6/7', '7/8', '3/7', '0/8']
Result: not reproduced. Under every selection the criterion supports, at most 2 of the 7 counts are within tolerance: over the 115 sessions of selection (a) the main burst reading (peak b, minimum a, N a) gives 52/60, 2/9, 8/8, 9/10, 10/11, 3/8 and 0/9 against 40/43, 0/6, 7/7, 6/7, 7/8, 3/7 and 0/8, and selection (b) gives 2 of 7 as well. Over selection (c) the main reading gives 40/47, 1/6, 7/7, 6/7, 7/8, 3/7 and 0/7, six of seven equal to the paper's, and the first local maximum read with the refractory plateau left out (peak c) gives 37/47 reference sessions and the same other counts. Selection (c) leaves out the off-line sorted sessions, which no passage of the paper supports, and was chosen in part because these counts then come close to Table 1, so it is a variant here and does not decide the outcome. The outcome recorded after the first verification, partially reproduced, rested on selection (c); the second verification found this circular, and the third verification agreed with the change to not reproduced. Under peak reading (a) 5 of 47 reference sessions burst, and with N the number of trials (b) none does.
table1-isolated-spikes¶
The percentage of isolated spikes among all bursts of the bursting sessions depends on the stimulus (74, 100, 62, 60, 55, 96 and 100% for protocols 1 to 7), and in every protocol isolated spikes are the most probable burst.
T1_ISO = [74, 100, 62, 60, 55, 96, 100]
def protocol_counts(sel, lim):
# Pooled n counts per protocol over the bursting sessions (all sessions where a protocol has none).
out = {}
for p in range(1, 8):
ses = [s for s in sel if s.protocol == p]
bur = [s for s in ses if lim[s.name] > 0]
use = bur if bur else ses
tot = np.zeros(400)
per = []
for s in use:
c = ia1.burst_counts(s.trials, lim[s.name])
tot[:len(c)] += c
per.append(c[0] / c.sum())
out[p] = (np.trim_zeros(tot, "b"), per, len(bur))
return out
rows = []
PC = {}
for kr, sel, c in [("c", KEPT, c) for c in COMBOS] + [("a", KEPT_A, c) for c in COMBOS]:
pc = protocol_counts(sel, limits[c])
PC[(kr, c)] = pc
iso = [100 * pc[p][0][0] / pc[p][0].sum() for p in range(1, 8)]
iso_avg = [100 * np.mean(pc[p][1]) for p in range(1, 8)]
most = [int(np.argmax(pc[p][0]) + 1) for p in range(1, 8)]
rows.append(dict(kept=kr, readings="".join(c), **{f"p{p}": round(v, 1) for p, v in zip(range(1, 8), iso)},
within=sum(abs(v - e) <= 0.1 * e for v, e in zip(iso, T1_ISO)), most_probable_n=most,
session_average=[round(v, 1) for v in iso_avg]))
print(pd.DataFrame(rows).to_string(index=False))
print("paper:", T1_ISO)
kept readings p1 p2 p3 p4 p5 p6 p7 within most_probable_n session_average c aaa 92.2 100.0 53.1 89.6 93.3 100.0 100.0 3 [1, 1, 1, 1, 1, 1, 1] [83.9, 100.0, 51.5, 89.6, 83.7, 100.0, 100.0] c aab 100.0 100.0 100.0 100.0 100.0 100.0 100.0 3 [1, 1, 1, 1, 1, 1, 1] [100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0] c aac 92.2 100.0 53.1 89.6 93.3 100.0 100.0 3 [1, 1, 1, 1, 1, 1, 1] [83.9, 100.0, 51.5, 89.6, 83.7, 100.0, 100.0] c aba 89.3 100.0 51.4 7.5 87.8 100.0 100.0 3 [1, 1, 1, 2, 1, 1, 1] [78.4, 100.0, 50.1, 7.5, 59.7, 100.0, 100.0] c abb 100.0 100.0 100.0 100.0 100.0 100.0 100.0 3 [1, 1, 1, 1, 1, 1, 1] [100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0] c abc 89.3 100.0 51.4 7.5 87.8 100.0 100.0 3 [1, 1, 1, 2, 1, 1, 1] [78.4, 100.0, 50.1, 7.5, 59.7, 100.0, 100.0] c baa 80.7 71.7 66.1 78.0 71.3 95.0 100.0 4 [1, 1, 1, 1, 1, 1, 1] [80.8, 71.7, 66.8, 78.6, 71.3, 95.3, 100.0] c bab 100.0 100.0 100.0 100.0 100.0 100.0 100.0 3 [1, 1, 1, 1, 1, 1, 1] [100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0] c bac 80.7 71.7 66.1 79.0 71.3 95.0 100.0 4 [1, 1, 1, 1, 1, 1, 1] [80.8, 71.7, 66.8, 79.3, 71.3, 95.3, 100.0] c bba 71.5 71.7 61.2 15.8 51.9 94.4 100.0 5 [1, 1, 1, 1, 1, 1, 1] [71.7, 71.7, 61.0, 11.9, 49.7, 94.6, 100.0] c bbb 100.0 100.0 100.0 100.0 100.0 100.0 100.0 3 [1, 1, 1, 1, 1, 1, 1] [100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0] c bbc 71.5 71.7 61.2 15.8 51.9 94.4 100.0 5 [1, 1, 1, 1, 1, 1, 1] [71.7, 71.7, 61.0, 11.9, 49.7, 94.6, 100.0] c caa 81.9 88.6 66.1 80.0 83.6 97.0 100.0 3 [1, 1, 1, 1, 1, 1, 1] [80.7, 88.6, 66.8, 80.0, 72.5, 97.2, 100.0] c cab 100.0 100.0 100.0 100.0 100.0 100.0 100.0 3 [1, 1, 1, 1, 1, 1, 1] [100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0] c cac 81.9 88.6 66.1 80.9 83.6 96.9 100.0 3 [1, 1, 1, 1, 1, 1, 1] [80.7, 88.6, 66.8, 80.7, 72.5, 97.1, 100.0] c cba 73.9 80.6 61.2 15.8 76.7 96.1 100.0 4 [1, 1, 1, 1, 1, 1, 1] [72.3, 80.6, 61.0, 11.9, 56.2, 96.3, 100.0] c cbb 100.0 100.0 100.0 100.0 100.0 100.0 100.0 3 [1, 1, 1, 1, 1, 1, 1] [100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0] c cbc 73.9 80.6 61.2 15.8 76.7 96.1 100.0 4 [1, 1, 1, 1, 1, 1, 1] [72.3, 80.6, 61.0, 11.9, 56.2, 96.3, 100.0] a aaa 90.8 100.0 53.1 87.3 92.4 100.0 100.0 3 [1, 1, 1, 1, 1, 1, 1] [82.4, 100.0, 51.5, 85.7, 83.8, 100.0, 100.0] a aab 100.0 100.0 100.0 100.0 100.0 100.0 100.0 3 [1, 1, 1, 1, 1, 1, 1] [100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0] a aac 90.8 100.0 53.1 87.3 92.4 100.0 100.0 3 [1, 1, 1, 1, 1, 1, 1] [82.4, 100.0, 51.5, 85.7, 83.8, 100.0, 100.0] a aba 88.1 100.0 51.4 8.2 88.1 100.0 100.0 3 [1, 1, 1, 2, 1, 1, 1] [77.9, 100.0, 50.1, 8.4, 64.6, 100.0, 100.0] a abb 100.0 100.0 100.0 100.0 100.0 100.0 100.0 3 [1, 1, 1, 1, 1, 1, 1] [100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0] a abc 88.1 100.0 51.4 8.2 88.1 100.0 100.0 3 [1, 1, 1, 2, 1, 1, 1] [77.9, 100.0, 50.1, 8.4, 64.6, 100.0, 100.0] a baa 79.5 77.7 67.7 78.1 73.5 95.0 100.0 4 [1, 1, 1, 1, 1, 1, 1] [79.2, 80.7, 68.6, 78.5, 73.0, 95.3, 100.0] a bab 100.0 100.0 100.0 100.0 100.0 100.0 100.0 3 [1, 1, 1, 1, 1, 1, 1] [100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0] a bac 79.5 77.7 67.7 78.8 73.5 95.0 100.0 4 [1, 1, 1, 1, 1, 1, 1] [79.2, 80.7, 68.6, 79.0, 73.0, 95.3, 100.0] a bba 71.1 68.4 61.8 14.0 49.4 94.4 100.0 4 [1, 1, 1, 1, 1, 1, 1] [71.0, 65.0, 61.7, 10.6, 48.1, 94.6, 100.0] a bbb 100.0 100.0 100.0 100.0 100.0 100.0 100.0 3 [1, 1, 1, 1, 1, 1, 1] [100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0] a bbc 71.1 68.4 61.8 14.0 49.4 94.4 100.0 4 [1, 1, 1, 1, 1, 1, 1] [71.0, 65.0, 61.7, 10.6, 48.1, 94.6, 100.0] a caa 82.4 87.2 67.7 79.6 85.7 97.0 100.0 3 [1, 1, 1, 1, 1, 1, 1] [80.3, 88.4, 68.6, 79.5, 75.7, 97.2, 100.0] a cab 100.0 100.0 100.0 100.0 100.0 100.0 100.0 3 [1, 1, 1, 1, 1, 1, 1] [100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0] a cac 82.4 87.2 67.7 80.3 85.7 96.9 100.0 3 [1, 1, 1, 1, 1, 1, 1] [80.3, 88.4, 68.6, 79.9, 75.7, 97.1, 100.0] a cba 75.6 75.2 61.8 14.0 79.4 96.1 99.9 4 [1, 1, 1, 1, 1, 1, 1] [72.9, 76.4, 61.7, 10.6, 58.5, 96.3, 99.9] a cbb 100.0 100.0 100.0 100.0 100.0 100.0 100.0 3 [1, 1, 1, 1, 1, 1, 1] [100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0] a cbc 75.6 75.2 61.8 14.0 79.4 96.1 99.9 4 [1, 1, 1, 1, 1, 1, 1] [72.9, 76.4, 61.7, 10.6, 58.5, 96.3, 99.9] paper: [74, 100, 62, 60, 55, 96, 100]
Result: not reproduced. Under the main reading the pooled percentage of isolated spikes over selection (c) is 80.7, 71.7, 66.1, 78.0, 71.3, 95.0 and 100% against 74, 100, 62, 60, 55, 96 and 100%: 4 of 7 within tolerance (protocols 1, 3, 6, 7). At 3 dB the one bursting session gives 71.7% where the paper has no bursts, and 25 Hz and 100 Hz have more isolated spikes than the paper. No combination of readings brings more than 5 of 7 within tolerance (4 over selection (a)). n = 1 is the most probable burst in every protocol under the main reading, as the paper says. The criterion requires every value; 4 of 7 is below the five of seven at which the other Table 1 criteria call a result partial, so the outcome is not reproduced (the first round recorded it as partially reproduced).
table1-burst-length¶
The highest n and the average n depend on the stimulus; the longest bursts (up to 15 spikes, average 2.10 with SD 3.67) occur at 25 Hz, and no bursts occur at 3 dB or at 800 Hz.
EXP_HI = [9, 1, 8, 15, 9, 3, 1]
EXP_AV = [1.3, 1, 1.5, 2.1, 1.7, 1.04, 1]
EXP_SD = {1: 1.3, 3: 1.6, 4: 3.67, 5: 1.8, 6: 0.42}
rows = []
for kr, c in PC:
pc = PC[(kr, c)]
hi = [len(pc[p][0]) for p in range(1, 8)]
av = [ia1.average_n(pc[p][0])[0] for p in range(1, 8)]
sd = [ia1.average_n(pc[p][0])[1] for p in range(1, 8)]
ok = (sum(abs(h - e) <= 0.1 * e for h, e in zip(hi, EXP_HI)) + sum(abs(a - e) <= 0.1 * e for a, e in zip(av, EXP_AV))
+ sum(abs(sd[p - 1] - e) <= 0.1 * e for p, e in EXP_SD.items()))
rows.append(dict(kept=kr, readings="".join(c), highest_n=hi, average_n=[round(a, 2) for a in av], sd_n=[round(x, 2) for x in sd], within_of_19=ok))
print(pd.DataFrame(rows).to_string(index=False))
print("paper: highest", EXP_HI, " average", EXP_AV, " SD", EXP_SD)
kept readings highest_n average_n sd_n within_of_19
c aaa [5, 1, 6, 12, 6, 1, 1] [1.09, 1.0, 1.63, 1.2, 1.09, 1.0, 1.0] [0.31, 0.0, 0.81, 0.79, 0.37, 0.0, 0.0] 6
c aab [1, 1, 1, 1, 1, 1, 1] [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0] [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] 5
c aac [5, 1, 6, 12, 6, 1, 1] [1.09, 1.0, 1.63, 1.2, 1.09, 1.0, 1.0] [0.31, 0.0, 0.81, 0.79, 0.37, 0.0, 0.0] 6
c aba [6, 1, 8, 43, 8, 1, 1] [1.13, 1.0, 1.66, 8.14, 1.23, 1.0, 1.0] [0.39, 0.0, 0.84, 6.83, 0.73, 0.0, 0.0] 6
c abb [1, 1, 1, 1, 1, 1, 1] [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0] [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] 5
c abc [6, 1, 8, 43, 8, 1, 1] [1.13, 1.0, 1.66, 8.14, 1.23, 1.0, 1.0] [0.39, 0.0, 0.84, 6.83, 0.73, 0.0, 0.0] 6
c baa [7, 9, 8, 14, 7, 5, 1] [1.23, 1.39, 1.49, 1.9, 1.39, 1.05, 1.0] [0.51, 0.72, 0.82, 2.22, 0.71, 0.23, 0.0] 8
c bab [1, 1, 1, 1, 1, 1, 1] [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0] [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] 5
c bac [7, 9, 8, 14, 7, 5, 1] [1.23, 1.39, 1.49, 1.87, 1.39, 1.05, 1.0] [0.51, 0.72, 0.82, 2.2, 0.71, 0.23, 0.0] 7
c bba [10, 9, 8, 73, 9, 5, 1] [1.38, 1.39, 1.59, 7.26, 1.78, 1.06, 1.0] [0.71, 0.72, 0.93, 7.6, 1.04, 0.24, 0.0] 8
c bbb [1, 1, 1, 1, 1, 1, 1] [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0] [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] 5
c bbc [10, 9, 8, 73, 9, 5, 1] [1.38, 1.39, 1.59, 7.26, 1.78, 1.06, 1.0] [0.71, 0.72, 0.93, 7.6, 1.04, 0.24, 0.0] 8
c caa [7, 5, 8, 14, 9, 3, 1] [1.22, 1.13, 1.49, 1.83, 1.24, 1.03, 1.0] [0.5, 0.38, 0.82, 2.16, 0.62, 0.17, 0.0] 9
c cab [1, 1, 1, 1, 1, 1, 1] [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0] [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] 5
c cac [7, 5, 8, 14, 9, 3, 1] [1.22, 1.13, 1.49, 1.8, 1.24, 1.03, 1.0] [0.5, 0.38, 0.82, 2.13, 0.62, 0.17, 0.0] 9
c cba [10, 5, 8, 73, 9, 3, 1] [1.36, 1.24, 1.59, 7.26, 1.4, 1.04, 1.0] [0.69, 0.55, 0.93, 7.6, 0.9, 0.2, 0.0] 8
c cbb [1, 1, 1, 1, 1, 1, 1] [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0] [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] 5
c cbc [10, 5, 8, 73, 9, 3, 1] [1.36, 1.24, 1.59, 7.26, 1.4, 1.04, 1.0] [0.69, 0.55, 0.93, 7.6, 0.9, 0.2, 0.0] 8
a aaa [7, 1, 6, 13, 6, 1, 1] [1.11, 1.0, 1.63, 1.32, 1.1, 1.0, 1.0] [0.37, 0.0, 0.81, 1.19, 0.4, 0.0, 0.0] 6
a aab [1, 1, 1, 1, 1, 1, 1] [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0] [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] 5
a aac [7, 1, 6, 13, 6, 1, 1] [1.11, 1.0, 1.63, 1.32, 1.1, 1.0, 1.0] [0.37, 0.0, 0.81, 1.19, 0.4, 0.0, 0.0] 6
a aba [7, 1, 8, 43, 11, 1, 1] [1.14, 1.0, 1.66, 7.76, 1.23, 1.0, 1.0] [0.43, 0.0, 0.84, 6.44, 0.76, 0.0, 0.0] 6
a abb [1, 1, 1, 1, 1, 1, 1] [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0] [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] 5
a abc [7, 1, 8, 43, 11, 1, 1] [1.14, 1.0, 1.66, 7.76, 1.23, 1.0, 1.0] [0.43, 0.0, 0.84, 6.44, 0.76, 0.0, 0.0] 6
a baa [7, 9, 8, 14, 7, 5, 1] [1.25, 1.3, 1.46, 1.85, 1.36, 1.05, 1.0] [0.54, 0.64, 0.8, 2.12, 0.69, 0.23, 0.0] 7
a bab [1, 1, 1, 1, 1, 1, 1] [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0] [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] 5
a bac [7, 9, 8, 14, 7, 5, 1] [1.25, 1.3, 1.46, 1.82, 1.36, 1.05, 1.0] [0.54, 0.64, 0.8, 2.1, 0.69, 0.23, 0.0] 7
a bba [10, 12, 8, 73, 11, 5, 1] [1.39, 1.48, 1.59, 7.38, 1.87, 1.06, 1.0] [0.72, 0.85, 0.93, 7.31, 1.14, 0.24, 0.0] 7
a bbb [1, 1, 1, 1, 1, 1, 1] [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0] [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] 5
a bbc [10, 12, 8, 73, 11, 5, 1] [1.39, 1.48, 1.59, 7.38, 1.87, 1.06, 1.0] [0.72, 0.85, 0.93, 7.31, 1.14, 0.24, 0.0] 7
a caa [7, 5, 8, 14, 9, 3, 1] [1.21, 1.15, 1.46, 1.8, 1.2, 1.03, 1.0] [0.5, 0.43, 0.8, 2.07, 0.57, 0.17, 0.0] 9
a cab [1, 1, 1, 1, 1, 1, 1] [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0] [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] 5
a cac [7, 5, 8, 14, 9, 3, 1] [1.21, 1.15, 1.46, 1.78, 1.2, 1.03, 1.0] [0.5, 0.43, 0.8, 2.05, 0.57, 0.17, 0.0] 9
a cba [10, 12, 8, 73, 11, 3, 2] [1.33, 1.39, 1.59, 7.38, 1.37, 1.04, 1.0] [0.67, 0.82, 0.93, 7.31, 0.89, 0.2, 0.03] 6
a cbb [1, 1, 1, 1, 1, 1, 1] [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0] [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] 5
a cbc [10, 12, 8, 73, 11, 3, 2] [1.33, 1.39, 1.59, 7.38, 1.37, 1.04, 1.0] [0.67, 0.82, 0.93, 7.31, 0.89, 0.2, 0.03] 6
paper: highest [9, 1, 8, 15, 9, 3, 1] average [1.3, 1, 1.5, 2.1, 1.7, 1.04, 1] SD {1: 1.3, 3: 1.6, 4: 3.67, 5: 1.8, 6: 0.42}
Result: not reproduced. Under the main reading over selection (c) 8 of the 19 values are within tolerance (highest n 7, 9, 8, 14, 7, 5, 1 against 9, 1, 8, 15, 9, 3, 1; average n 1.23, 1.39, 1.49, 1.90, 1.39, 1.05, 1.00 against 1.3, 1, 1.5, 2.1, 1.7, 1.04, 1), below the 10 the criterion needs for a partial reproduction, and no combination of readings under either selection gives more than 9 (peak reading c). The longest bursts are at 25 Hz (14 spikes) and 800 Hz has none, as in the paper, but 3 dB has bursts of up to 9 spikes in its one bursting session.
fig5-burst-probability-vs-stimulus¶
The probability of generating isolated spikes is lowest, and high-n bursts most frequent, for large amplitude fluctuations (12 dB) and cutoff frequencies around 100 Hz.
def pn_table(sel, lim, bursting_only):
rows = []
for s in sel:
if bursting_only and lim[s.name] <= 0:
continue
c = ia1.burst_counts(s.trials, lim[s.name]).astype(float)
c = np.r_[c, np.zeros(3)][:3] / c.sum()
rows.append(dict(protocol=s.protocol, P1=c[0], P2=c[1], P3=c[2]))
return pd.DataFrame(rows).groupby("protocol").agg(["mean", "std", "size"])
for kr, sel, c in [("c", KEPT, MAIN), ("c", KEPT, ("b", "b", "a")), ("c", KEPT, ("b", "a", "c")), ("c", KEPT, ("a", "a", "a")),
("c", KEPT, ("c", "a", "a")), ("a", KEPT_A, MAIN)]:
for reading, bo in [("(a) all kept", False), ("(b) bursting only", True)]:
t = pn_table(sel, limits[c], bo)
means = t.xs("mean", axis=1, level=1)
cut = {25: 4, 100: 5, 200: 1, 400: 6, 800: 7}
sdp = {3: 2, 6: 1, 12: 3}
p1c = {k: means.P1.get(v, np.nan) for k, v in cut.items()}
p1s = {k: means.P1.get(v, np.nan) for k, v in sdp.items()}
print(f"kept ({kr}), readings {''.join(c)}, sessions {reading}")
print(" P(n=1) by cutoff:", {k: round(v, 3) for k, v in p1c.items()}, " lowest at", min(p1c, key=p1c.get), "Hz")
print(" P(n=1) by SD:", {k: round(v, 3) for k, v in p1s.items()}, " lowest at", min(p1s, key=p1s.get), "dB")
for col in ["P2", "P3"]:
print(f" {col} by cutoff:", {k: round(means[col].get(v, np.nan), 3) for k, v in cut.items()},
" by SD:", {k: round(means[col].get(v, np.nan), 3) for k, v in sdp.items()})
t = pn_table(KEPT, LIM, False)
print("\nmain reading, all kept (c): mean, SD and number of sessions per protocol")
print(t.round(3).to_string())
kept (c), readings baa, sessions (a) all kept
P(n=1) by cutoff: {25: np.float64(0.816), 100: np.float64(0.749), 200: np.float64(0.837), 400: np.float64(0.98), 800: np.float64(1.0)} lowest at 100 Hz
P(n=1) by SD: {3: np.float64(0.953), 6: np.float64(0.837), 12: np.float64(0.668)} lowest at 12 dB
P2 by cutoff: {25: np.float64(0.058), 100: np.float64(0.175), 200: np.float64(0.138), 400: np.float64(0.02), 800: np.float64(0.0)} by SD: {3: np.float64(0.033), 6: np.float64(0.138), 12: np.float64(0.248)}
P3 by cutoff: {25: np.float64(0.027), 100: np.float64(0.056), 200: np.float64(0.022), 400: np.float64(0.0), 800: np.float64(0.0)} by SD: {3: np.float64(0.011), 6: np.float64(0.022), 12: np.float64(0.061)}
kept (c), readings baa, sessions (b) bursting only
P(n=1) by cutoff: {25: np.float64(0.786), 100: np.float64(0.713), 200: np.float64(0.808), 400: np.float64(0.953), 800: nan} lowest at 100 Hz
P(n=1) by SD: {3: np.float64(0.717), 6: np.float64(0.808), 12: np.float64(0.668)} lowest at 12 dB
P2 by cutoff: {25: np.float64(0.068), 100: np.float64(0.201), 200: np.float64(0.162), 400: np.float64(0.046), 800: nan} by SD: {3: np.float64(0.2), 6: np.float64(0.162), 12: np.float64(0.248)}
P3 by cutoff: {25: np.float64(0.032), 100: np.float64(0.064), 200: np.float64(0.026), 400: np.float64(0.001), 800: nan} by SD: {3: np.float64(0.064), 6: np.float64(0.026), 12: np.float64(0.061)}
kept (c), readings bba, sessions (a) all kept
P(n=1) by cutoff: {25: np.float64(0.245), 100: np.float64(0.56), 200: np.float64(0.759), 400: np.float64(0.977), 800: np.float64(1.0)} lowest at 25 Hz
P(n=1) by SD: {3: np.float64(0.953), 6: np.float64(0.759), 12: np.float64(0.61)} lowest at 12 dB
P2 by cutoff: {25: np.float64(0.085), 100: np.float64(0.257), 200: np.float64(0.177), 400: np.float64(0.023), 800: np.float64(0.0)} by SD: {3: np.float64(0.033), 6: np.float64(0.177), 12: np.float64(0.275)}
P3 by cutoff: {25: np.float64(0.092), 100: np.float64(0.1), 200: np.float64(0.047), 400: np.float64(0.0), 800: np.float64(0.0)} by SD: {3: np.float64(0.011), 6: np.float64(0.047), 12: np.float64(0.076)}
kept (c), readings bba, sessions (b) bursting only
P(n=1) by cutoff: {25: np.float64(0.119), 100: np.float64(0.497), 200: np.float64(0.717), 400: np.float64(0.946), 800: nan} lowest at 25 Hz
P(n=1) by SD: {3: np.float64(0.717), 6: np.float64(0.717), 12: np.float64(0.61)} lowest at 12 dB
P2 by cutoff: {25: np.float64(0.099), 100: np.float64(0.294), 200: np.float64(0.207), 400: np.float64(0.053), 800: nan} by SD: {3: np.float64(0.2), 6: np.float64(0.207), 12: np.float64(0.275)}
P3 by cutoff: {25: np.float64(0.107), 100: np.float64(0.114), 200: np.float64(0.055), 400: np.float64(0.001), 800: nan} by SD: {3: np.float64(0.064), 6: np.float64(0.055), 12: np.float64(0.076)}
kept (c), readings bac, sessions (a) all kept
P(n=1) by cutoff: {25: np.float64(0.822), 100: np.float64(0.749), 200: np.float64(0.837), 400: np.float64(0.98), 800: np.float64(1.0)} lowest at 100 Hz
P(n=1) by SD: {3: np.float64(0.953), 6: np.float64(0.837), 12: np.float64(0.668)} lowest at 12 dB
P2 by cutoff: {25: np.float64(0.056), 100: np.float64(0.175), 200: np.float64(0.138), 400: np.float64(0.02), 800: np.float64(0.0)} by SD: {3: np.float64(0.033), 6: np.float64(0.138), 12: np.float64(0.248)}
P3 by cutoff: {25: np.float64(0.026), 100: np.float64(0.056), 200: np.float64(0.022), 400: np.float64(0.0), 800: np.float64(0.0)} by SD: {3: np.float64(0.011), 6: np.float64(0.022), 12: np.float64(0.061)}
kept (c), readings bac, sessions (b) bursting only
P(n=1) by cutoff: {25: np.float64(0.793), 100: np.float64(0.713), 200: np.float64(0.808), 400: np.float64(0.953), 800: nan} lowest at 100 Hz
P(n=1) by SD: {3: np.float64(0.717), 6: np.float64(0.808), 12: np.float64(0.668)} lowest at 12 dB
P2 by cutoff: {25: np.float64(0.065), 100: np.float64(0.201), 200: np.float64(0.162), 400: np.float64(0.046), 800: nan} by SD: {3: np.float64(0.2), 6: np.float64(0.162), 12: np.float64(0.248)}
P3 by cutoff: {25: np.float64(0.03), 100: np.float64(0.064), 200: np.float64(0.026), 400: np.float64(0.001), 800: nan} by SD: {3: np.float64(0.064), 6: np.float64(0.026), 12: np.float64(0.061)}
kept (c), readings aaa, sessions (a) all kept
P(n=1) by cutoff: {25: np.float64(0.985), 100: np.float64(0.939), 200: np.float64(0.983), 400: np.float64(1.0), 800: np.float64(1.0)} lowest at 100 Hz
P(n=1) by SD: {3: np.float64(1.0), 6: np.float64(0.983), 12: np.float64(0.862)} lowest at 12 dB
P2 by cutoff: {25: np.float64(0.009), 100: np.float64(0.042), 200: np.float64(0.015), 400: np.float64(0.0), 800: np.float64(0.0)} by SD: {3: np.float64(0.0), 6: np.float64(0.015), 12: np.float64(0.102)}
P3 by cutoff: {25: np.float64(0.003), 100: np.float64(0.015), 200: np.float64(0.002), 400: np.float64(0.0), 800: np.float64(0.0)} by SD: {3: np.float64(0.0), 6: np.float64(0.002), 12: np.float64(0.027)}
kept (c), readings aaa, sessions (b) bursting only
P(n=1) by cutoff: {25: np.float64(0.896), 100: np.float64(0.837), 200: np.float64(0.839), 400: nan, 800: nan} lowest at 100 Hz
P(n=1) by SD: {3: nan, 6: np.float64(0.839), 12: np.float64(0.515)} lowest at 3 dB
P2 by cutoff: {25: np.float64(0.063), 100: np.float64(0.112), 200: np.float64(0.142), 400: nan, 800: nan} by SD: {3: nan, 6: np.float64(0.142), 12: np.float64(0.358)}
P3 by cutoff: {25: np.float64(0.02), 100: np.float64(0.04), 200: np.float64(0.019), 400: nan, 800: nan} by SD: {3: nan, 6: np.float64(0.019), 12: np.float64(0.096)}
kept (c), readings caa, sessions (a) all kept
P(n=1) by cutoff: {25: np.float64(0.829), 100: np.float64(0.759), 200: np.float64(0.848), 400: np.float64(0.988), 800: np.float64(1.0)} lowest at 100 Hz
P(n=1) by SD: {3: np.float64(0.981), 6: np.float64(0.848), 12: np.float64(0.668)} lowest at 12 dB
P2 by cutoff: {25: np.float64(0.053), 100: np.float64(0.164), 200: np.float64(0.127), 400: np.float64(0.012), 800: np.float64(0.0)} by SD: {3: np.float64(0.017), 6: np.float64(0.127), 12: np.float64(0.248)}
P3 by cutoff: {25: np.float64(0.025), 100: np.float64(0.053), 200: np.float64(0.021), 400: np.float64(0.0), 800: np.float64(0.0)} by SD: {3: np.float64(0.002), 6: np.float64(0.021), 12: np.float64(0.061)}
kept (c), readings caa, sessions (b) bursting only
P(n=1) by cutoff: {25: np.float64(0.8), 100: np.float64(0.725), 200: np.float64(0.807), 400: np.float64(0.972), 800: nan} lowest at 100 Hz
P(n=1) by SD: {3: np.float64(0.886), 6: np.float64(0.807), 12: np.float64(0.668)} lowest at 12 dB
P2 by cutoff: {25: np.float64(0.061), 100: np.float64(0.187), 200: np.float64(0.162), 400: np.float64(0.028), 800: nan} by SD: {3: np.float64(0.101), 6: np.float64(0.162), 12: np.float64(0.248)}
P3 by cutoff: {25: np.float64(0.029), 100: np.float64(0.06), 200: np.float64(0.026), 400: np.float64(0.0), 800: nan} by SD: {3: np.float64(0.012), 6: np.float64(0.026), 12: np.float64(0.061)}
kept (a), readings baa, sessions (a) all kept
P(n=1) by cutoff: {25: np.float64(0.806), 100: np.float64(0.755), 200: np.float64(0.82), 400: np.float64(0.982), 800: np.float64(1.0)} lowest at 100 Hz
P(n=1) by SD: {3: np.float64(0.957), 6: np.float64(0.82), 12: np.float64(0.686)} lowest at 12 dB
P2 by cutoff: {25: np.float64(0.065), 100: np.float64(0.172), 200: np.float64(0.15), 400: np.float64(0.017), 800: np.float64(0.0)} by SD: {3: np.float64(0.032), 6: np.float64(0.15), 12: np.float64(0.234)}
P3 by cutoff: {25: np.float64(0.03), 100: np.float64(0.056), 200: np.float64(0.025), 400: np.float64(0.0), 800: np.float64(0.0)} by SD: {3: np.float64(0.008), 6: np.float64(0.025), 12: np.float64(0.058)}
kept (a), readings baa, sessions (b) bursting only
P(n=1) by cutoff: {25: np.float64(0.785), 100: np.float64(0.73), 200: np.float64(0.792), 400: np.float64(0.953), 800: nan} lowest at 100 Hz
P(n=1) by SD: {3: np.float64(0.807), 6: np.float64(0.792), 12: np.float64(0.686)} lowest at 12 dB
P2 by cutoff: {25: np.float64(0.072), 100: np.float64(0.189), 200: np.float64(0.173), 400: np.float64(0.046), 800: nan} by SD: {3: np.float64(0.146), 6: np.float64(0.173), 12: np.float64(0.234)}
P3 by cutoff: {25: np.float64(0.033), 100: np.float64(0.062), 200: np.float64(0.029), 400: np.float64(0.001), 800: nan} by SD: {3: np.float64(0.037), 6: np.float64(0.029), 12: np.float64(0.058)}
main reading, all kept (c): mean, SD and number of sessions per protocol
P1 P2 P3
mean std size mean std size mean std size
protocol
1 0.837 0.108 47 0.138 0.087 47 0.022 0.020 47
2 0.953 0.116 6 0.033 0.082 6 0.011 0.026 6
3 0.668 0.136 7 0.248 0.088 7 0.061 0.045 7
4 0.816 0.169 7 0.058 0.030 7 0.027 0.021 7
5 0.749 0.150 8 0.175 0.099 8 0.056 0.035 8
6 0.980 0.027 7 0.020 0.026 7 0.000 0.001 7
7 1.000 0.000 7 0.000 0.000 7 0.000 0.000 7
Result: reproduced. Under the main reading over selection (c) the population average of P(n = 1) is lowest at 100 Hz among the cutoffs (0.749, against 0.816 at 25, 0.837 at 200, 0.980 at 400 and 1.000 at 800 Hz) and lowest at 12 dB among the SDs (0.668, against 0.837 at 6 and 0.953 at 3 dB), with all kept sessions and with bursting sessions only; the same holds over selection (a) (0.755 and 0.686) and with the plateau-excluded peak (c). P(n = 2) and P(n = 3) are highest at 100 Hz and at 12 dB over all kept sessions. With the bounding minimum read as the lowest value (b) the lowest P(n = 1) moves to 25 Hz, where that reading gives very long bursts.
fig3-correlation-shape¶
For a cell tested with all stimuli, larger stimulus SD gives correlation functions with sharper peaks, and faster stimuli lead from multi-modal (25 Hz) to single-peaked (100, 200 Hz) to shallower and broader correlation functions (400, 800 Hz).
singlecells.txt lists one cell recorded in seven experiments, 03-04-23-ad to -aj, which together
hold all seven protocols (six reference sessions, one of each other protocol): this is taken as the
cell of Figure 3. C_s is normalised by the squared mean rate so that sessions compare. "Sharper"
is read as (a) a higher first peak over the baseline (the mean of C_s/rbar^2 over 20 to 40 ms) and
(b) a narrower first peak (full width at half height above that baseline). Peaks are counted on
C_s/rbar^2 smoothed with a 1 ms running mean (a standard choice that keeps the ms-scale shape), as
local maxima within 0 < tau <= 1.25 / cutoff that rise above the neighbouring minima by more than a
prominence threshold; the threshold is not in the paper, so 0.05, 0.1 and 0.2 are shown.
from scipy.signal import find_peaks
fig3_cell = next(row for row in cells if "03-04-23-ad" in row)
F3 = [s for s in P if s.exp in fig3_cell]
rows = []
curves = []
fig, axes = plt.subplots(1, 2, figsize=(12, 4))
for s in F3:
tau, c, sd, count, ntr, sdt = ia1.correlation(s.trials, 0.06)
rb = ia1.binned(s.trials).astype(bool).mean() / ia1.DT
cn = c / rb ** 2
base = cn[(tau >= 0.02) & (tau <= 0.04)].mean()
sm = np.convolve(cn, np.ones(10) / 10, mode="same")
sm[:6] = cn[:6]
lim_i = int(1.25 / s.cutoff / ia1.DT)
ip = 1 + np.argmax(cn[1:51])
height = cn[ip] - base
width = ia1.fwhm(cn, ip, base) * ia1.DT * 1e3
npk = {}
for prom in [0.05, 0.1, 0.2]:
pk, _ = find_peaks(sm[1:lim_i + 1], prominence=prom)
npk[prom] = len(pk)
rows.append(dict(session=s.name, sd=s.sd_cut[0], cutoff=s.cutoff, bursting=LIM[s.name] > 0,
limit_ms=round(LIM[s.name] * 1e3, 1), peak_ms=round(tau[ip] * 1e3, 1), peak_height=round(height, 2),
peak_fwhm_ms=round(width, 1), **{f"peaks_prom{k}": v for k, v in npk.items()}))
curves.append(cn[1:])
ax = axes[0] if s.sd_cut[1] == 200 else axes[1]
ax.plot(tau[1:] * 1e3, cn[1:], lw=0.8, label=f"{s.sd_cut[0]} dB, {s.cutoff} Hz ({s.exp[-2:]})")
F3T = pd.DataFrame(rows)
print(F3T.to_string(index=False))
for ax, t in zip(axes, ["cutoff 200 Hz", "SD 6 dB"]):
ax.set(xlabel="tau (ms)", ylabel="C_s / rbar^2", title=t, ylim=(-1.5, 4))
ax.legend(fontsize=7)
allc = np.concatenate(curves)
print("C_s/rbar^2 for tau > 0: range", allc.min().round(2), "to", allc.max().round(2), "; points outside the axes (-1.5 to 4):", int(((allc < -1.5) | (allc > 4)).sum()))
plt.show()
ref = F3T[(F3T.sd == 6) & (F3T.cutoff == 200)]
print("reference sessions: peak height range", ref.peak_height.min(), "to", ref.peak_height.max(), "; FWHM range", ref.peak_fwhm_ms.min(), "to", ref.peak_fwhm_ms.max(), "ms")
session sd cutoff bursting limit_ms peak_ms peak_height peak_fwhm_ms peaks_prom0.05 peaks_prom0.1 peaks_prom0.2 03-04-23-ad_sig1 6 200 True 4.4 3.5 2.32 0.4 1 1 1 03-04-23-ad_sig2 6 800 False 0.0 4.9 0.40 0.4 0 0 0 03-04-23-ae_sig1 6 200 True 4.1 3.5 2.35 0.4 1 1 1 03-04-23-ae_sig2 12 200 True 4.2 3.3 5.06 0.3 1 1 1 03-04-23-ag_sig1 6 200 True 4.6 3.9 1.40 0.5 1 1 1 03-04-23-ag_sig2 6 400 False 0.0 4.6 0.35 0.2 0 0 0 03-04-23-ah_sig1 6 200 True 4.7 4.2 1.29 0.5 1 1 1 03-04-23-ah_sig2 6 100 True 5.3 4.3 2.59 0.8 2 2 1 03-04-23-ai_sig1 6 200 True 4.7 4.5 1.14 0.5 1 1 1 03-04-23-ai_sig2 6 25 True 4.9 4.2 4.96 0.7 14 5 4 03-04-23-aj_sig1 6 200 True 4.7 4.0 1.75 0.5 1 1 1 03-04-23-aj_sig2 3 200 False 0.0 5.0 0.50 0.4 1 1 0 C_s/rbar^2 for tau > 0: range -1.01 to 5.06 ; points outside the axes (-1.5 to 4): 6
reference sessions: peak height range 1.14 to 2.35 ; FWHM range 0.4 to 0.5 ms
Result: reproduced (resting on the peak-counting threshold). In the cell recorded with all seven protocols (03-04-23-ad to -aj) the first peak of C_s/rbar^2 rises above its baseline by 0.50 at 3 dB, 1.14 to 2.35 in the six 6 dB sessions and 5.06 at 12 dB, so it sharpens with SD under reading (a); its width (b) does not order (0.4, 0.4 to 0.5, 0.3 ms). At 25 Hz C_s has 4 to 14 peaks, at 200 Hz one; at 100 Hz one peak at prominence 0.2 and two at 0.05 and 0.1. The 3 dB, 400 Hz and 800 Hz sessions have no separable first peak, as the paper names for Figures 3A, E and F.
method-correlation-error¶
The temporal average of the estimated error of the correlation function was on average 3.4% (SD 1.5%) of the span of C_s.
Delta = sigma / sqrt(N) at each tau, averaged over tau, over max minus min of C_s, averaged over kept sessions. Readings: sigma and N (a), (b) and (c); tau range (a) 0 to 1.25 / cutoff, (b) the whole range up to 800 ms. Whether tau = 0, where C_s holds the self term (rbar / dt, two orders of magnitude above the rest), enters the span is not stated: both are shown.
rows = []
for s in KEPT_A:
tau, c, sd, count, ntr, sdt = ia1.correlation(s.trials, 0.7999)
L = int(round(1.25 / s.cutoff / ia1.DT))
for nm in "abc":
err = ia1.correlation_error(sd, count, ntr, sdt, nm)
for rng_name, sl in [("(a) 0 to 1.25/cutoff", slice(0, L + 1)), ("(b) 0 to 800 ms", slice(0, None))]:
for z0 in [True, False]:
ss = sl if z0 else slice(1, sl.stop)
rows.append(dict(session=s.name, kept_c=s.name in POP["c"], N=nm, tau_range=rng_name, tau0=z0,
pct=100 * err[ss].mean() / (c[ss].max() - c[ss].min())))
E = pd.DataFrame(rows)
print("kept (c)")
print(E[E.kept_c].groupby(["N", "tau_range", "tau0"]).pct.agg(["mean", "std", "size"]).round(2).to_string())
print("\nkept (a)")
print(E.groupby(["N", "tau_range", "tau0"]).pct.agg(["mean", "std", "size"]).round(2).to_string())
kept (c)
mean std size
N tau_range tau0
a (a) 0 to 1.25/cutoff False 49.84 176.45 89
True 0.07 0.02 89
(b) 0 to 800 ms False 2.63 1.39 89
True 0.16 0.03 89
b (a) 0 to 1.25/cutoff False 4454.51 15773.85 89
True 6.33 1.67 89
(b) 0 to 800 ms False 132.04 72.35 89
True 7.73 1.35 89
c (a) 0 to 1.25/cutoff False 33.96 130.35 89
True 0.06 0.02 89
(b) 0 to 800 ms False 2.61 1.38 89
True 0.15 0.03 89
kept (a)
mean std size
N tau_range tau0
a (a) 0 to 1.25/cutoff False 44.80 165.44 115
True 0.07 0.02 115
(b) 0 to 800 ms False 2.69 1.37 115
True 0.16 0.03 115
b (a) 0 to 1.25/cutoff False 4004.05 14789.74 115
True 6.54 1.73 115
(b) 0 to 800 ms False 135.30 71.45 115
True 7.93 1.46 115
c (a) 0 to 1.25/cutoff False 30.33 119.83 115
True 0.06 0.02 115
(b) 0 to 800 ms False 2.67 1.36 115
True 0.16 0.03 115
Result: reproduced (resting on leaving tau = 0 out of the span). With sigma the SD of the product terms and N their number (reading a), over the whole range of lags and without the zero-lag self term, the error averages 2.63% (SD 1.39%) of the span of C_s over the 89 sessions of selection (c), within 3.4 +/- 1.5%; reading (c) of sigma and N gives 2.61%, and over the 115 sessions of selection (a) 2.69% and 2.67%. With tau = 0 in the span every reading is below 0.2% or, with N the number of trials, 6.3 to 7.9%; over 0 to 1.25 / cutoff without tau = 0 the values are unstable (means 30 to 4500%).
model-matches-classification (blocked)¶
A threshold-linear Poisson model with refractoriness and no free parameters classifies sessions as bursting or not in agreement with the real data in 86% of sessions, and where both are bursting their limiting ISIs differ by less than 1 ms in 81% of cases.
The model thresholds the stimulus at the cell's measured threshold s_TH before filtering (Eq. 2),
and the criterion says the claim is not testable as the paper ran it if the data set does not hold
s_TH. Data1 holds no threshold (see threshold-population). Blocked, data_missing; the dependent
model-reproduces-coding-trends is blocked for the same reason.
fig6-bursts-less-correlated¶
Representing the spike train as n-bursts reduces correlations between time bins; the mean squared Pearson coefficient over t from 200 to 990 ms and tau from 0 to 10 ms is 2.89 times larger (SD 1.49) for spikes than for bursts, on average over the bursting sessions.
Eq. 11 across trials at each (t, tau). Bin width (a) 0.1 ms, (b) 0.4 ms; tau = 0 included or not. The spike variable is 1 when the bin holds a spike. Bursting sessions under the main burst reading; detail (4) reading (b) (bursting sessions that have bursts of n > 1) is shown too.
Result: reproduced. Over the 64 bursting sessions of selection (c) the mean squared Pearson coefficient is larger for spikes than for bursts in every session; the mean ratio is 1.99 (SD 0.92) with 0.1 ms bins and 2.17 (SD 1.86) with 0.4 ms bins when tau = 0 is left out, within 2.89 +/- 1.49, and 1.29 and 1.21 when it is included (outside). Over the 84 bursting sessions of selection (a) the ratios are 1.99 and 2.13 without tau = 0. Every bursting session has bursts of n > 1, so reading (b) of the sessions is the same.
Information per burst (Figure 7): per-session measures¶
For every kept session, under the main burst reading, with the bin widths (a) 0.1 ms and (b) 0.4 ms: rbar_n, I_n(1) in bits per event (reading i; reading ii is T = 0.8 s times it), I'_n in bits/s, and the jitter of each n (readings (a) and (b) of the events). Fits are over the n values present. Sessions of reading (a) of the selection are computed and summarised under both selections.
The plug-in estimate of I_n(1) is biased upwards for rare events, and more so the rarer they are,
so part of a rise of I_n(1) with n can come from the fall in occurrence. As a baseline (not in the
paper), the same measures are computed after shifting each trial's burst onsets by a random
offset, wrapped within 200 to 1000 ms (ia1.circular_shift), which keeps the bursts of each
trial and removes their locking to the stimulus; ten shifts per session, averaged.
T = ia1.T1 - ia1.T0
rng7 = np.random.default_rng(2008)
info, info_sh = {}, {}
for s in KEPT_A:
tr, on, n = ia1.bursts(s.trials, LIM[s.name])
shifted = [ia1.circular_shift(tr, on, s.n_trials, rng7) for _ in range(10)]
for bw in [1e-4, 4e-4]:
info[(s.name, bw)] = ia1.information(ia1.rates_of(on, n, s.n_trials, bw))
sh = [ia1.information(ia1.rates_of(o, n, s.n_trials, bw)) for o in shifted]
info_sh[(s.name, bw)] = {m: {k: np.mean([x[m][k] for x in sh]) for k in sh[0][m]} for m in sh[0]}
def fit_rows(sel, bw, src):
rows = []
for s in sel:
inf = src[(s.name, bw)]
ns = np.array(sorted(inf))
i1 = np.array([inf[n]["i1"] for n in ns])
rbar = np.array([inf[n]["rbar"] for n in ns])
rate = np.array([inf[n]["rate"] for n in ns])
occ = rbar / rbar.sum()
row = dict(session=s.name, kept_c=s.name in POP["c"], protocol=s.protocol, n_max=ns.max(), n_values=len(ns))
for name, y in [("i1", i1), ("i1_ii", i1 * T), ("occ_log2", np.log2(occ)), ("occ_ln", np.log(occ)),
("occ_log10", np.log10(occ)), ("rate", rate)]:
sl, se, pg, p2 = ia1.line_fit(ns, y)
row.update({f"{name}_slope": sl, f"{name}_se": se, f"{name}_pg": pg, f"{name}_p2": p2})
row["total_rate"] = rate.sum()
row["frac_gt1"] = rate[ns > 1].sum() / rate.sum()
rows.append(row)
return pd.DataFrame(rows)
F7 = {bw: fit_rows(BURSTING_A, bw, info) for bw in [1e-4, 4e-4]}
F7S = {bw: fit_rows(BURSTING_A, bw, info_sh) for bw in [1e-4, 4e-4]}
for kr in "ca":
f = F7[1e-4] if kr == "a" else F7[1e-4][F7[1e-4].kept_c]
print(f"kept ({kr}): bursting sessions {len(f)}; with at least 3 n values (a slope with an error):",
int((f.n_values >= 3).sum()), "; with 2:", int((f.n_values == 2).sum()), "; with 1:", int((f.n_values == 1).sum()))
kept (c): bursting sessions 64; with at least 3 n values (a slope with an error): 63 ; with 2: 1 ; with 1: 0 kept (a): bursting sessions 84; with at least 3 n values (a slope with an error): 83 ; with 2: 1 ; with 1: 0
fig7ab-information-per-burst-increases¶
The information carried by each n-burst increases with n; every bursting session has a significantly positive slope of I_n(1) against n, with an average slope of 1.5 bits (SD 0.7 bits) over 59 bursting cells.
Significance (a) slope more than one standard error above zero, (b) one-sided t-test p < 0.05. Sessions with fewer than three n values have no error for the slope and count as not significant.
rows = []
for kr, bw in itertools.product("ca", [1e-4, 4e-4]):
f, fs = F7[bw], F7S[bw]
if kr == "c":
f, fs = f[f.kept_c], fs[fs.kept_c]
for norm, col in [("(i) bits per event", "i1"), ("(ii) Eq. 10 as printed", "i1_ii")]:
for subset, g, gs in [("(a) all bursting", f, fs), ("(b) n_max > 1", f[f.n_max > 1], fs[f.n_max.values > 1])]:
sig_a = (g[f"{col}_slope"] - g[f"{col}_se"] > 0)
sig_b = (g[f"{col}_pg"] < 0.05)
rows.append(dict(kept=kr, bin_ms=bw * 1e3, I1=norm, sessions=subset, n=len(g), mean_slope=round(g[f"{col}_slope"].mean(), 2),
sd_slope=round(g[f"{col}_slope"].std(ddof=1), 2), positive=int((g[f"{col}_slope"] > 0).sum()),
sig_a=int(sig_a.sum()), sig_b=int(sig_b.sum()), shifted_mean_slope=round(gs[f"{col}_slope"].mean(), 2),
slope_minus_shifted=round((g[f"{col}_slope"].values - gs[f"{col}_slope"].values).mean(), 2),
above_shifted=int((g[f"{col}_slope"].values > gs[f"{col}_slope"].values).sum())))
print(pd.DataFrame(rows).to_string(index=False))
kept bin_ms I1 sessions n mean_slope sd_slope positive sig_a sig_b shifted_mean_slope slope_minus_shifted above_shifted c 0.1 (i) bits per event (a) all bursting 64 2.22 0.87 64 63 62 2.59 -0.37 3 c 0.1 (i) bits per event (b) n_max > 1 64 2.22 0.87 64 63 62 2.59 -0.37 3 c 0.1 (ii) Eq. 10 as printed (a) all bursting 64 1.78 0.69 64 63 62 2.07 -0.30 3 c 0.1 (ii) Eq. 10 as printed (b) n_max > 1 64 1.78 0.69 64 63 62 2.07 -0.30 3 c 0.4 (i) bits per event (a) all bursting 64 1.91 0.75 64 63 62 2.23 -0.32 9 c 0.4 (i) bits per event (b) n_max > 1 64 1.91 0.75 64 63 62 2.23 -0.32 9 c 0.4 (ii) Eq. 10 as printed (a) all bursting 64 1.53 0.60 64 63 62 1.78 -0.26 9 c 0.4 (ii) Eq. 10 as printed (b) n_max > 1 64 1.53 0.60 64 63 62 1.78 -0.26 9 a 0.1 (i) bits per event (a) all bursting 84 2.14 0.85 84 83 82 2.49 -0.35 6 a 0.1 (i) bits per event (b) n_max > 1 84 2.14 0.85 84 83 82 2.49 -0.35 6 a 0.1 (ii) Eq. 10 as printed (a) all bursting 84 1.71 0.68 84 83 82 1.99 -0.28 6 a 0.1 (ii) Eq. 10 as printed (b) n_max > 1 84 1.71 0.68 84 83 82 1.99 -0.28 6 a 0.4 (i) bits per event (a) all bursting 84 1.85 0.73 84 83 82 2.15 -0.30 13 a 0.4 (i) bits per event (b) n_max > 1 84 1.85 0.73 84 83 82 2.15 -0.30 13 a 0.4 (ii) Eq. 10 as printed (a) all bursting 84 1.48 0.58 84 83 82 1.72 -0.24 13 a 0.4 (ii) Eq. 10 as printed (b) n_max > 1 84 1.48 0.58 84 83 82 1.72 -0.24 13
Result: partially reproduced. Over the 64 bursting sessions of selection (c), within 59 +/- 5.9, the mean slope of I_n(1) on n is 2.22 bits (SD 0.87) with 0.1 ms bins, just outside 1.5 +/- 0.7, and 1.91 (SD 0.75) with 0.4 ms bins in bits per event, 1.78 and 1.53 for Eq. 10 as printed, within. All 64 slopes are positive; 63 are significant by one standard error and 62 by the t-test (one session has only two n values), so not every session. Over the 84 bursting sessions of selection (a), outside 59 +/- 5.9, the slopes are 2.14, 1.85, 1.71 and 1.48.
The baseline with each trial's burst onsets shifted at random gives a larger mean slope, 2.59 bits with 0.1 ms bins and 2.23 with 0.4 ms bins (selection (c)); the slope of the data exceeds that of its shifted copy in only 3 of 64 sessions (9 with 0.4 ms bins). A shifted train has lost its locking to the stimulus, so its rate is spread over more bins and its plug-in bias is at least that of the data; the comparison shows that the bias of the estimate for rare events is of the order of the rise of I_n(1) with n, not that the whole rise is bias. The criterion fixes the formula, so the outcome rests on it.
fig7cd-occurrence-decreases¶
High-n bursts are rare; in every session the frequency of occurrence decreases significantly with n, with an average slope of -2.6 (SD 0.7).
Slope of the logarithm of rbar_n / sum rbar_n on n; logarithm (a) base 2, (b) natural, (c) base 10. Significance as in fig7ab (negative here). The occurrence does not depend on the bin width.
rows = []
for kr in "ca":
f = F7[1e-4] if kr == "a" else F7[1e-4][F7[1e-4].kept_c]
for col in ["occ_log2", "occ_ln", "occ_log10"]:
for subset, g in [("(a) all bursting", f), ("(b) n_max > 1", f[f.n_max > 1])]:
rows.append(dict(kept=kr, log=col, sessions=subset, n=len(g), mean_slope=round(g[f"{col}_slope"].mean(), 2),
sd_slope=round(g[f"{col}_slope"].std(ddof=1), 2), negative=int((g[f"{col}_slope"] < 0).sum()),
sig_a=int((g[f"{col}_slope"] + g[f"{col}_se"] < 0).sum()), sig_b=int(((1 - g[f"{col}_pg"]) < 0.05).sum())))
print(pd.DataFrame(rows).to_string(index=False))
kept log sessions n mean_slope sd_slope negative sig_a sig_b c occ_log2 (a) all bursting 64 -2.82 1.12 64 63 61 c occ_log2 (b) n_max > 1 64 -2.82 1.12 64 63 61 c occ_ln (a) all bursting 64 -1.96 0.77 64 63 61 c occ_ln (b) n_max > 1 64 -1.96 0.77 64 63 61 c occ_log10 (a) all bursting 64 -0.85 0.34 64 63 61 c occ_log10 (b) n_max > 1 64 -0.85 0.34 64 63 61 a occ_log2 (a) all bursting 84 -2.71 1.09 84 83 81 a occ_log2 (b) n_max > 1 84 -2.71 1.09 84 83 81 a occ_ln (a) all bursting 84 -1.88 0.76 84 83 81 a occ_ln (b) n_max > 1 84 -1.88 0.76 84 83 81 a occ_log10 (a) all bursting 84 -0.82 0.33 84 83 81 a occ_log10 (b) n_max > 1 84 -0.82 0.33 84 83 81
Result: partially reproduced. Over the 64 bursting sessions of selection (c) the log2 occurrence falls with n with a mean slope of -2.82 (SD 1.12), within -2.6 +/- 0.7 (natural log -1.96, log10 -0.85, outside); every slope is negative and 63 of 64 are significant by reading (a), 61 by reading (b), so not every session, as the claim states. Over selection (a): -2.71, with 83 and 81 of 84 significant.
fig7ef-jitter-independent-of-n¶
The jitter of the first spike of a burst is roughly independent of n in 80% of the bursting sessions, and the mean slope of jitter against n is -0.03 ms (SD 0.2 ms).
Jitter per n is the mean over events of the SD of the onset times. Events (a) peaks of the onset
histogram in 0.1 ms bins, between the troughs on either side; (b) onsets within the session's
first-peak domain (the limiting ISI) centred on an event centre, the centres being the histogram
bins taken from the highest down (ia1.jitter_events_b). Slopes over the n values with a jitter (events with at
least two onsets); "not significant" is (a) |slope| below one standard error, (b) two-sided p >=
0.05; sessions with fewer than three such n values have no test and are counted apart.
rows = []
for s in BURSTING_A:
tr, on, n = ia1.bursts(s.trials, LIM[s.name])
row = dict(session=s.name, kept_c=s.name in POP["c"])
for ev in "ab":
ns, js = [], []
for m in np.unique(n):
obt = [on[(n == m) & (tr == i)] for i in range(s.n_trials)]
j = ia1.jitter_events_a(obt) if ev == "a" else ia1.jitter_events_b(obt, LIM[s.name])
if np.isfinite(j):
ns.append(m)
js.append(j * 1e3)
sl, se, pg, p2 = ia1.line_fit(ns, js)
row.update({f"{ev}_slope": sl, f"{ev}_se": se, f"{ev}_p2": p2, f"{ev}_nvals": len(ns),
f"{ev}_j1": js[0] if ns and ns[0] == 1 else np.nan})
rows.append(row)
J = pd.DataFrame(rows)
summ = []
for kr, ev in itertools.product("ca", "ab"):
Jk = J[J.kept_c] if kr == "c" else J
g = Jk[Jk[f"{ev}_nvals"] >= 3]
summ.append(dict(kept=kr, events=ev, sessions_with_test=len(g), no_test=int((Jk[f"{ev}_nvals"] < 3).sum()),
mean_slope_ms=round(Jk[f"{ev}_slope"].mean(), 3), sd_slope_ms=round(Jk[f"{ev}_slope"].std(ddof=1), 3),
pct_not_sig_a=round(100 * (g[f"{ev}_slope"].abs() < g[f"{ev}_se"]).mean(), 1),
pct_not_sig_b=round(100 * (g[f"{ev}_p2"] >= 0.05).mean(), 1),
median_jitter_n1_ms=round(Jk[f"{ev}_j1"].median(), 3)))
print(pd.DataFrame(summ).to_string(index=False))
kept events sessions_with_test no_test mean_slope_ms sd_slope_ms pct_not_sig_a pct_not_sig_b median_jitter_n1_ms c a 63 1 -0.001 0.005 44.4 88.9 0.058 c b 61 3 0.005 0.052 32.8 86.9 0.389 a a 83 1 -0.001 0.005 43.4 91.6 0.058 a b 81 3 0.005 0.050 34.6 87.7 0.389
Result: reproduced (resting on how events are formed). Over the bursting sessions of selection (c), with events read as onsets within the first-peak domain around an event centre (b), the mean slope of jitter on n is 0.005 ms (SD 0.052), within -0.03 +/- 0.2, and 86.9% of the 61 sessions with a test have a slope not significantly different from 0 by the t-test, inside 80 +/- 8. With events between histogram troughs (a) the slope is -0.001 ms but 88.9% are not significant (just outside); with significance read as one standard error, 32.8% and 44.4% (outside). Over selection (a), 87.7% (b) and 91.6% (a). The median jitter of 1-bursts is 0.389 ms under (b) and 0.058 ms under (a), the latter below the paper's 0.2 ms error bar, so reading (a) splits events.
fig7gh-information-rate¶
The information rate of n-bursts decreases with n (average slope -11 bits/s, SD 16 bits/s), yet bursts with n > 1 together carry on average 47% of the total information over the bursting cells.
I'_n = rbar_n I_n(1) (bits/s), the same under readings (i) and (ii) of I_n(1).
rows = []
for kr, bw in itertools.product("ca", [1e-4, 4e-4]):
f, fs = F7[bw], F7S[bw]
if kr == "c":
f, fs = f[f.kept_c], fs[fs.kept_c]
for subset, g, gs in [("(a) all bursting", f, fs), ("(b) n_max > 1", f[f.n_max > 1], fs[f.n_max.values > 1])]:
rows.append(dict(kept=kr, bin_ms=bw * 1e3, sessions=subset, n=len(g), mean_slope=round(g.rate_slope.mean(), 1),
sd_slope=round(g.rate_slope.std(ddof=1), 1), negative=int((g.rate_slope < 0).sum()),
mean_frac_gt1_pct=round(100 * g.frac_gt1.mean(), 1), sd_frac_pct=round(100 * g.frac_gt1.std(ddof=1), 1),
mean_total_bits_s=round(g.total_rate.mean(), 1), shifted_slope=round(gs.rate_slope.mean(), 1),
shifted_total_bits_s=round(gs.total_rate.mean(), 1), shifted_frac_gt1_pct=round(100 * gs.frac_gt1.mean(), 1)))
print(pd.DataFrame(rows).to_string(index=False))
kept bin_ms sessions n mean_slope sd_slope negative mean_frac_gt1_pct sd_frac_pct mean_total_bits_s shifted_slope shifted_total_bits_s shifted_frac_gt1_pct c 0.1 (a) all bursting 64 -46.5 39.0 64 35.5 13.3 241.4 -16.1 106.5 46.8 c 0.1 (b) n_max > 1 64 -46.5 39.0 64 35.5 13.3 241.4 -16.1 106.5 46.8 c 0.4 (a) all bursting 64 -38.1 33.1 64 37.7 13.8 199.8 -3.9 37.5 63.2 c 0.4 (b) n_max > 1 64 -38.1 33.1 64 37.7 13.8 199.8 -3.9 37.5 63.2 a 0.1 (a) all bursting 84 -44.7 36.0 84 36.3 12.7 246.6 -16.0 112.5 47.5 a 0.1 (b) n_max > 1 84 -44.7 36.0 84 36.3 12.7 246.6 -16.0 112.5 47.5 a 0.4 (a) all bursting 84 -36.6 30.6 84 38.6 13.2 203.4 -3.9 40.2 63.8 a 0.4 (b) n_max > 1 84 -36.6 30.6 84 38.6 13.2 203.4 -3.9 40.2 63.8
Result: not reproduced. Over the 64 bursting sessions of selection (c) the information rate falls with n in every session, but the mean slope is -46.5 bits/s (SD 39.0) with 0.1 ms bins and -38.1 (SD 33.1) with 0.4 ms bins, outside -11 +/- 16, and bursts with n > 1 carry 35.5% and 37.7% of the information, outside 47 +/- 4.7. The total rate averages 241.4 and 199.8 bits/s. With the burst onsets shifted at random in each trial the same formula still gives 106.5 and 37.5 bits/s, slopes of -16.1 and -3.9 bits/s and fractions of 46.8% and 63.2%: a large part of these rates is estimator bias. Over selection (a): -44.7 and -36.6 bits/s, 36.3% and 38.6%.
direct-method-comparison¶
Information rates from the burst code are highly correlated with those of the direct method (R = 0.95), averaging 191 against 222 bits/s, the direct method being higher in all sessions but one by 31 bits/s on average, so the burst code keeps about 85% of the information.
Direct method on binary words of 0.4 ms bins, words starting at every bin over 200 to 1000 ms,
entropies NSB-corrected (ia1.nsb_entropy). Readings (a) the rate at 8 bins (3.2 ms), (b) a
linear extrapolation in 1 / word length over 1 to 8 bins. The burst-code rate is I' of
fig7gh (non-bursting sessions as all 1-bursts), with bins (a) 0.1 ms and (b) 0.4 ms. Sessions (a)
all kept, (b) bursting.
rows = []
for s in KEPT_A:
rates = [ia1.direct_information(s.trials, L) for L in range(1, 9)]
sl, ic, *_ = stats.linregress(1 / np.arange(1, 9), rates)
row = dict(session=s.name, kept_c=s.name in POP["c"], bursting=LIM[s.name] > 0, direct_8=rates[-1], direct_extrap=ic)
for bw in [1e-4, 4e-4]:
row[f"burst_{bw * 1e3:.1f}"] = sum(v["rate"] for v in info[(s.name, bw)].values())
rows.append(row)
DM = pd.DataFrame(rows)
summ = []
for dcol in ["direct_8", "direct_extrap"]:
for bcol in ["burst_0.1", "burst_0.4"]:
for lab, g in [("(a) kept (c)", DM[DM.kept_c]), ("(b) bursting (c)", DM[DM.kept_c & DM.bursting]),
("(a) kept (a)", DM), ("(b) bursting (a)", DM[DM.bursting])]:
diff = g[dcol] - g[bcol]
summ.append(dict(direct=dcol, burst=bcol, sessions=lab, n=len(g), R=round(np.corrcoef(g[dcol], g[bcol])[0, 1], 3),
mean_direct=round(g[dcol].mean(), 1), sd_direct=round(g[dcol].std(ddof=1), 1),
mean_burst=round(g[bcol].mean(), 1), sd_burst=round(g[bcol].std(ddof=1), 1),
mean_diff=round(diff.mean(), 1), sd_diff=round(diff.std(ddof=1), 1),
direct_higher=f"{int((diff > 0).sum())}/{len(g)}",
lost_pct=round(100 * (1 - g[bcol].mean() / g[dcol].mean()), 1)))
print(pd.DataFrame(summ).to_string(index=False))
direct burst sessions n R mean_direct sd_direct mean_burst sd_burst mean_diff sd_diff direct_higher lost_pct
direct_8 burst_0.1 (a) kept (c) 89 0.946 198.4 76.6 220.1 64.0 -21.8 26.2 15/89 -11.0
direct_8 burst_0.1 (b) bursting (c) 64 0.920 220.7 71.1 241.4 56.4 -20.7 29.4 15/64 -9.4
direct_8 burst_0.1 (a) kept (a) 115 0.943 204.3 77.2 227.3 66.8 -23.1 26.5 19/115 -11.3
direct_8 burst_0.1 (b) bursting (a) 84 0.916 224.3 72.6 246.6 61.1 -22.3 29.6 19/84 -9.9
direct_8 burst_0.4 (a) kept (c) 89 0.967 198.4 76.6 178.1 65.6 20.3 21.2 74/89 10.2
direct_8 burst_0.4 (b) bursting (c) 64 0.951 220.7 71.1 199.8 59.2 21.0 23.5 52/64 9.5
direct_8 burst_0.4 (a) kept (a) 115 0.965 204.3 77.2 183.6 67.6 20.7 21.4 94/115 10.1
direct_8 burst_0.4 (b) bursting (a) 84 0.950 224.3 72.6 203.4 62.8 20.9 23.5 66/84 9.3
direct_extrap burst_0.1 (a) kept (c) 89 0.945 197.7 79.0 220.1 64.0 -22.4 27.9 15/89 -11.3
direct_extrap burst_0.1 (b) bursting (c) 64 0.922 219.8 74.5 241.4 56.4 -21.6 31.3 15/64 -9.8
direct_extrap burst_0.1 (a) kept (a) 115 0.941 203.5 79.6 227.3 66.8 -23.9 28.1 18/115 -11.7
direct_extrap burst_0.1 (b) bursting (a) 84 0.918 223.2 76.0 246.6 61.1 -23.4 31.3 18/84 -10.5
direct_extrap burst_0.4 (a) kept (c) 89 0.967 197.7 79.0 178.1 65.6 19.6 22.8 73/89 9.9
direct_extrap burst_0.4 (b) bursting (c) 64 0.954 219.8 74.5 199.8 59.2 20.1 25.3 51/64 9.1
direct_extrap burst_0.4 (a) kept (a) 115 0.965 203.5 79.6 183.6 67.6 19.9 22.8 94/115 9.8
direct_extrap burst_0.4 (b) bursting (a) 84 0.953 223.2 76.0 203.4 62.8 19.8 25.0 66/84 8.9
Result: partially reproduced. With 0.4 ms bins for the burst code, over the 89 kept sessions of selection (c), R = 0.967, the direct method gives 198.4 bits/s (8 bins) and the burst code 178.1, a mean difference of 20.3 bits/s, all within tolerance, and the burst code keeps 90% (10.2% lost, outside 14 +/- 1.4). The direct method is higher in 74 of 89 sessions (83%), not in all but one. Over the bursting sessions R = 0.951, 220.7 against 199.8 bits/s, difference 21.0, higher in 52 of 64. Over selection (a) the values are alike (R = 0.965, 204.3 against 183.6, 10.1% lost, 94 of 115). The extrapolated direct rate changes these by under 1 bit/s. With 0.1 ms bins the burst code is above the direct method on average (220.1 against 198.4 bits/s).
Burst-triggered averages: per-session measures¶
For each bursting session (main reading), the stimulus in dB (stimulus file shifted to the
intensity in cell.dat, as stimtodB.m does) around every burst onset in 200 to 1000 ms, from
-25 to +15 ms at 20 kHz (Methods, Estimation of burst-triggered averages). tau_n is the distance
from the nBTA's central maximum to the onset (Eq. 13 and the sentence after it: "The time tau_n
between the maximum of nBTA and tau = 0"), which the Results call "the mean delay between stimulus
upstroke and burst generation" (Qualitative description, paragraph 3). The central maximum is
therefore searched (a) before the onset, -25 to 0 ms, the main reading since revision 1; (b) over
the whole window, -25 to +15 ms, the reading of the first round, shown as a variant.
A session enters the analyses of n = 1 to 4 when it has at least two bursts of each n, the least that gives a variance; the paper says only that eight sessions had "too few bursts" (Results, paragraph 6), so this is an assumption, and the results with at least 20 and at least 50 bursts of each n are shown too.
def p_two_sided_ttests(a, b):
return stats.ttest_ind(a, b, axis=0).pvalue
SEARCH = {"a": (-ia1.PRE, 0.0), "b": None} # where the central maximum is searched
MIN_BURSTS = [2, 20, 50]
BTA = {}
for s in BURSTING_A:
st = stim_db(s)
tr, on, n = ia1.bursts(s.trials, LIM[s.name])
seg = {}
for m in range(1, 5):
if (n == m).sum() >= 2:
seg[m], tt = ia1.segments(st, on[n == m])
sp = spikes_in(s)
sta_seg, tt = ia1.segments(st, sp)
BTA[s.name] = dict(seg=seg, tt=tt, sta_seg=sta_seg, on=on, n=n, tr=tr, mean_stim=st[4000:20000].mean(),
counts=[len(seg[m]) if m in seg else int((n == m).sum()) for m in range(1, 5)])
have4 = [s for s in BURSTING_A if len(BTA[s.name]["seg"]) == 4]
for kr in "ca":
bs = [s for s in BURSTING_A if s.name in POP[kr]]
print(f"kept ({kr}): bursting sessions {len(bs)}; with at least k bursts of each n = 1 to 4:",
{k: sum(min(BTA[s.name]["counts"]) >= k for s in bs) for k in MIN_BURSTS})
cnt4 = np.array([BTA[s.name]["counts"][3] for s in have4 if s.name in POP["c"]])
print("kept (c), sessions with n = 1 to 4: number of 4-bursts, quantiles 0, 0.25, 0.5, 0.75, 1:", np.quantile(cnt4, [0, .25, .5, .75, 1]))
def subset(rows_df, kr, k):
# Rows of sessions in selection kr with at least k bursts of each n = 1 to 4.
ok = {s.name for s in have4 if s.name in POP[kr] and min(BTA[s.name]["counts"]) >= k}
return rows_df[rows_df.session.isin(ok)]
kept (c): bursting sessions 64; with at least k bursts of each n = 1 to 4: {2: 53, 20: 39, 50: 28}
kept (a): bursting sessions 84; with at least k bursts of each n = 1 to 4: {2: 71, 20: 52, 50: 38}
kept (c), sessions with n = 1 to 4: number of 4-bursts, quantiles 0, 0.25, 0.5, 0.75, 1: [ 2. 17. 57. 153. 553.]
fig8-nbta-differ¶
In every bursting cell the burst-triggered averages for n = 1 to 4 differ significantly (in the example cell their central peaks grow higher and wider with n); 50 of the 58 bursting sessions show significantly different central peaks in the time domain, the other eight having too few bursts.
Two-way ANOVAs (frequency bin x n) on the real and the imaginary parts of the Fourier transform of each segment, each at 0.01 / 2; frequency bins (a) all, (b) up to the cutoff. Then at each time point, t-tests of each nBTA against the three others, p < 0.01; a point is marked for n when all three are significant. Central peak: from the minimum before to the minimum after the central maximum, searched (a) before the onset or (b) over the whole window. Readings (a) at least one marked point in the central peak of every n, (b) more than half of the points of every central peak marked. The two-way ANOVAs treat every segment and frequency bin as an observation, as the Methods describe; with hundreds to thousands of segments they reject in nearly every session.
rows = []
for s in have4:
B = BTA[s.name]
tt = B["tt"]
segs = [B["seg"][m] for m in range(1, 5)]
ffts = [np.fft.rfft(g, axis=1) for g in segs]
freqs = np.fft.rfftfreq(len(tt), 1 / ia1.STIM_RATE)
row = dict(session=s.name, protocol=s.protocol, counts=[len(g) for g in segs])
for fr, sel in [("a", slice(None)), ("b", freqs <= s.cutoff)]:
ps = []
for part in [np.real, np.imag]:
gs = [part(f)[:, sel] for f in ffts]
if part is np.imag:
gs = [g[:, 1:] for g in gs] # the imaginary part of the 0 Hz bin is 0
ps.append(min(ia1.two_way_anova(gs)))
row[f"anova_{fr}"] = min(ps) < 0.01 / 2
btas = [g.mean(0) for g in segs]
sig = []
for i in range(4):
x = np.ones(len(tt), bool)
for j in range(4):
if j != i:
x &= p_two_sided_ttests(segs[i], segs[j]) < 0.01
sig.append(x)
for sr in "ab":
marked = []
for i in range(4):
ip, a, b = ia1.central_peak(btas[i], tt, SEARCH[sr])
marked.append(sig[i][a:b + 1].mean())
row[f"peak_any_{sr}"] = all(m > 0 for m in marked)
row[f"peak_most_{sr}"] = all(m > 0.5 for m in marked)
row[f"marked_{sr}"] = [round(m, 2) for m in marked]
rows.append(row)
N8 = pd.DataFrame(rows)
print(subset(N8, "c", 2).to_string(index=False))
out = []
for kr, k, sr in itertools.product("ca", MIN_BURSTS, "ab"):
g = subset(N8, kr, k)
out.append(dict(kept=kr, min_bursts=k, search=sr, sessions=len(g), bursting=sum(s.name in POP[kr] for s in BURSTING_A),
anova_a=int(g.anova_a.sum()), anova_b=int(g.anova_b.sum()),
peak_any=int((g.anova_a & g[f"peak_any_{sr}"]).sum()), peak_most=int((g.anova_a & g[f"peak_most_{sr}"]).sum())))
print()
print(pd.DataFrame(out).to_string(index=False))
session protocol counts anova_a anova_b peak_any_a peak_most_a marked_a peak_any_b peak_most_b marked_b 03-02-25-aj_sig1 1 [12482, 3808, 703, 147] True True True False [0.7, 0.54, 0.4, 0.17] True False [0.7, 0.54, 0.4, 0.17] 03-03-05-ab_sig1 1 [39236, 3669, 364, 8] True True False False [0.5, 0.03, 0.29, 0.0] False False [0.5, 0.03, 0.29, 0.0] 03-03-05-ab_sig2 2 [25529, 7111, 2271, 553] True True True False [0.45, 0.57, 0.04, 0.39] True False [0.45, 0.57, 0.6, 0.75] 03-03-05-ac_sig2 4 [7267, 628, 232, 111] True True True False [0.78, 0.53, 0.47, 0.43] True False [0.78, 0.53, 0.47, 0.64] 03-03-14-ad_sig1 1 [30960, 5345, 427, 35] True True True False [0.46, 0.66, 0.75, 0.75] True False [0.46, 0.66, 0.71, 0.64] 03-03-14-ad_sig2 5 [18897, 6486, 2531, 258] True True True True [0.86, 0.89, 0.8, 0.79] True True [0.86, 0.89, 0.8, 0.79] 03-03-18-ad_sig1 1 [18133, 3820, 704, 102] True True True False [0.62, 0.52, 0.41, 0.17] True False [0.62, 0.52, 0.41, 0.17] 03-03-18-ae_sig1 1 [17271, 4403, 1021, 247] True True True False [0.77, 0.64, 0.46, 0.24] True False [0.77, 0.64, 0.46, 0.24] 03-03-18-ae_sig2 3 [12104, 4529, 2639, 473] True True True True [0.83, 0.51, 0.66, 0.58] True True [0.83, 0.51, 0.66, 0.58] 03-03-21-aa_sig1 1 [7048, 708, 46, 3] True True False False [0.23, 0.0, 0.0, 0.0] False False [0.23, 0.0, 0.0, 0.0] 03-03-21-ab_sig1 1 [12775, 4203, 845, 157] True True True False [0.77, 0.61, 0.57, 0.4] True False [0.77, 0.61, 0.57, 0.4] 03-03-21-ab_sig2 6 [18153, 1183, 43, 4] True True False False [0.0, 0.06, 0.0, 0.13] False False [0.0, 0.06, 0.0, 0.41] 03-03-21-ac_sig1 1 [18790, 8229, 1443, 228] True True True True [0.69, 0.85, 0.75, 0.89] True True [0.69, 0.85, 0.75, 0.75] 03-03-21-ac_sig2 4 [5331, 677, 339, 380] True True True True [0.76, 0.68, 0.63, 0.74] True True [0.76, 0.68, 0.63, 0.74] 03-04-11-ab_sig2 3 [8325, 1305, 159, 20] True True True False [0.39, 0.17, 0.81, 0.22] True False [0.39, 0.17, 0.81, 0.22] 03-04-15-aa_sig2 5 [4316, 1275, 547, 165] True True True True [0.75, 0.73, 0.71, 0.67] True True [0.75, 0.73, 0.71, 0.67] 03-04-16-aa_sig1 1 [5695, 2440, 469, 4] True True False False [0.26, 0.09, 0.4, 0.0] True False [0.26, 0.09, 0.4, 0.3] 03-04-16-ab_sig1 1 [5189, 1910, 399, 3] True True True False [0.26, 0.23, 0.42, 0.04] True False [0.26, 0.23, 0.42, 0.04] 03-04-16-ac_sig1 1 [5180, 1989, 411, 2] True True True False [0.34, 0.34, 0.15, 0.13] True False [0.34, 0.34, 0.15, 0.13] 03-04-16-ac_sig2 4 [7965, 563, 181, 60] True True False False [0.81, 0.31, 0.0, 0.0] False False [0.81, 0.57, 0.0, 0.0] 03-04-16-ae_sig2 3 [2297, 1870, 473, 167] True True True False [0.75, 0.79, 0.38, 0.1] True False [0.75, 0.79, 0.38, 0.1] 03-04-16-af_sig1 1 [7904, 4104, 720, 115] True True True False [0.69, 0.58, 0.71, 0.48] True True [0.69, 0.58, 0.71, 0.79] 03-04-16-af_sig2 3 [5580, 3316, 949, 24] True True False False [0.6, 0.63, 0.43, 0.0] True False [0.6, 0.63, 0.43, 0.72] 03-04-23-aa_sig1 1 [7319, 1390, 295, 48] True True True False [0.42, 0.37, 0.12, 0.14] True False [0.42, 0.37, 0.12, 0.6] 03-04-23-aa_sig2 5 [5074, 1292, 527, 159] True True True False [0.75, 0.72, 0.34, 0.42] True False [0.75, 0.72, 0.34, 0.42] 03-04-23-ac_sig1 1 [6336, 1046, 99, 4] True True False False [0.0, 0.15, 0.06, 0.08] False False [0.0, 0.15, 0.06, 0.21] 03-04-23-ac_sig2 3 [5331, 1099, 160, 31] True True False False [0.51, 0.56, 0.0, 0.56] False False [0.51, 0.56, 0.0, 0.56] 03-04-23-ad_sig1 1 [6039, 1490, 357, 109] True True True False [0.66, 0.68, 0.58, 0.44] True True [0.66, 0.68, 0.58, 0.79] 03-04-23-ae_sig1 1 [6963, 1635, 410, 98] True True True False [0.64, 0.66, 0.5, 0.4] True False [0.64, 0.66, 0.5, 0.76] 03-04-23-ae_sig2 3 [4259, 1660, 188, 125] True True True False [0.73, 0.46, 0.06, 0.6] True False [0.73, 0.46, 0.06, 0.84] 03-04-23-ag_sig1 1 [6972, 1279, 322, 120] True True True False [0.6, 0.48, 0.47, 0.57] True False [0.6, 0.48, 0.47, 0.79] 03-04-23-ah_sig1 1 [7159, 1160, 267, 91] True True True False [0.59, 0.34, 0.27, 0.64] True False [0.59, 0.34, 0.27, 0.62] 03-04-23-ah_sig2 5 [3286, 2279, 605, 283] True True True False [0.92, 0.38, 0.61, 0.49] True False [0.92, 0.38, 0.61, 0.82] 03-04-23-ai_sig1 1 [7541, 936, 138, 36] True True True False [0.55, 0.42, 0.27, 0.32] True False [0.55, 0.42, 0.27, 0.48] 03-04-23-ai_sig2 4 [3496, 445, 350, 243] True True False False [0.75, 0.71, 0.0, 0.0] False False [0.75, 0.71, 0.0, 0.0] 03-04-23-aj_sig1 1 [6683, 1488, 371, 163] True True True False [0.69, 0.39, 0.3, 0.68] True False [0.69, 0.39, 0.3, 0.8] 03-04-25-ad_sig2 3 [5979, 1924, 307, 153] True True True False [0.71, 0.61, 0.49, 0.46] True False [0.71, 0.61, 0.49, 0.85] 03-04-25-af_sig1 1 [9399, 1161, 114, 17] True True True False [0.45, 0.39, 0.2, 0.16] True False [0.45, 0.39, 0.2, 0.14] 03-04-25-ag_sig1 1 [7680, 1652, 221, 27] True True False False [0.5, 0.49, 0.0, 0.03] False False [0.5, 0.49, 0.0, 0.03] 03-04-30-aa_sig1 1 [7520, 950, 120, 9] True True False False [0.16, 0.19, 0.0, 0.0] False False [0.16, 0.19, 0.0, 0.41] 03-04-30-aa_sig2 5 [6034, 1641, 383, 85] True True True False [0.55, 0.46, 0.3, 0.4] True False [0.55, 0.46, 0.3, 0.4] 03-05-02-ab_sig1 1 [6593, 582, 98, 25] True True True False [0.44, 0.59, 0.24, 0.21] True False [0.44, 0.59, 0.32, 0.53] 03-05-06-aa_sig2 5 [8381, 1344, 347, 34] True True True False [0.61, 0.64, 0.47, 0.39] True False [0.61, 0.64, 0.47, 0.39] 03-05-06-ab_sig2 4 [7577, 687, 246, 106] True True True False [0.85, 0.73, 0.48, 0.17] True False [0.85, 0.73, 0.48, 0.48] 03-05-09-ab_sig1 1 [8600, 667, 50, 5] True True False False [0.0, 0.0, 0.0, 0.0] False False [0.0, 0.0, 0.0, 0.16] 03-05-09-ab_sig2 5 [8475, 1111, 237, 28] True True True False [0.48, 0.38, 0.23, 0.27] True False [0.48, 0.38, 0.23, 0.27] 03-05-09-ac_sig2 4 [7613, 336, 129, 57] True True True False [0.76, 0.45, 0.27, 0.33] True False [0.76, 0.45, 0.27, 0.48] 03-06-27-ab_sig1 1 [5376, 1271, 199, 50] True True False False [0.41, 0.47, 0.3, 0.0] True False [0.41, 0.47, 0.3, 0.33] 03-06-27-ac_sig1 1 [6982, 942, 185, 28] True True True False [0.41, 0.26, 0.08, 0.45] True False [0.41, 0.26, 0.08, 0.43] 03-06-27-ae_sig1 1 [4629, 1565, 195, 13] True True True False [0.43, 0.45, 0.53, 0.43] True False [0.43, 0.45, 0.53, 0.51] 03-06-27-af_sig1 1 [5315, 1372, 102, 8] True True True False [0.19, 0.34, 0.31, 0.13] True False [0.19, 0.34, 0.31, 0.38] 03-06-27-ag_sig1 1 [5843, 1528, 122, 17] True True False False [0.33, 0.57, 0.0, 0.0] False False [0.33, 0.57, 0.0, 0.27] 03-07-03-aa_sig1 1 [4414, 1742, 352, 4] True True False False [0.38, 0.08, 0.36, 0.0] True False [0.38, 0.08, 0.36, 0.25] kept min_bursts search sessions bursting anova_a anova_b peak_any peak_most c 2 a 53 64 53 53 38 5 c 2 b 53 64 53 53 42 7 c 20 a 39 64 39 39 33 5 c 20 b 39 64 39 39 35 7 c 50 a 28 64 28 28 25 5 c 50 b 28 64 28 28 26 7 a 2 a 71 84 71 71 54 11 a 2 b 71 84 71 71 60 14 a 20 a 52 84 52 52 44 11 a 20 b 52 84 52 52 48 14 a 50 a 38 84 38 38 35 11 a 50 b 38 84 38 38 36 14
Result: partially reproduced. Of the 64 bursting sessions of selection (c) (58 +/- 5.8 in the paper, so just outside), 53 have at least two bursts of each n = 1 to 4, and all 53 pass the frequency-domain ANOVA, with all frequency bins and with bins up to the cutoff. In the time domain, with the central maximum searched before the onset (main reading), 38 of them have marked points in every central peak (reading a), outside 50 +/- 5, and 5 have most of every central peak marked (reading b); with the maximum searched over the whole window, 42 and 7. With at least 20 bursts of each n, 39 sessions, 33 and 35; with at least 50, 28 sessions, 25 and 26. Over selection (a), 71 sessions all pass the ANOVA and 54 (60 over the whole window) have marked points in every central peak. The ANOVA part of the claim holds; the time-domain count does not.
fig8-latency-decreases¶
The latency between the stimulus upstroke and burst onset decreases with n, so that high or wide deflections produce prompt responses with high-n bursts.
tau_n per session with bursts of n = 1 to 4. Readings (a) tau_n decreasing from n = 1 to 4 (tau_4 < tau_1 with no increase between consecutive n) in the majority of those sessions; (b) the mean over sessions of tau_n decreasing from 1 to 4.
rows = []
for s in have4:
B = BTA[s.name]
row = dict(session=s.name, protocol=s.protocol)
for sr in "ab":
taus = [-B["tt"][ia1.central_peak(B["seg"][m].mean(0), B["tt"], SEARCH[sr])[0]] * 1e3 for m in range(1, 5)]
row.update({f"tau{m}_{sr}": t for m, t in zip(range(1, 5), taus)})
row[f"decreasing_{sr}"] = all(np.diff(taus) <= 0) and taus[3] < taus[0]
row[f"tau4_below_tau1_{sr}"] = taus[3] < taus[0]
rows.append(row)
LT = pd.DataFrame(rows)
print(subset(LT, "c", 2).round(2).to_string(index=False))
out = []
for kr, k, sr in itertools.product("ca", MIN_BURSTS, "ab"):
g = subset(LT, kr, k)
tau = g[[f"tau{m}_{sr}" for m in range(1, 5)]]
out.append(dict(kept=kr, min_bursts=k, search=sr, sessions=len(g), decreasing=int(g[f"decreasing_{sr}"].sum()),
tau4_below_tau1=int(g[f"tau4_below_tau1_{sr}"].sum()), mean_tau_ms=tau.mean().round(2).tolist(),
negative_tau=int((tau.values < 0).any(axis=1).sum()),
wilcoxon_p=stats.wilcoxon(g[f"tau1_{sr}"], g[f"tau4_{sr}"]).pvalue if len(g) > 5 else np.nan))
print()
print(pd.DataFrame(out).to_string(index=False))
# The search before the onset ends at 0 ms, so a maximum found there sits on the edge of the search.
print()
for kr in "ca":
g = subset(LT, kr, 2)
tau = g[[f"tau{m}_a" for m in range(1, 5)]].values
edge = np.isclose(tau, 0)
free, free3 = ~edge.any(axis=1), ~edge[:, :3].any(axis=1)
print(f"kept ({kr}), search before onset: sessions with a tau_n at the 0 ms edge", int(edge.any(axis=1).sum()), "of", len(g),
"; values at the edge", int(edge.sum()), "(per n", edge.sum(0).tolist(), ")")
print(" sessions with no tau_n at the edge:", int(free.sum()), "; their mean tau_n", tau[free].mean(0).round(2).tolist(),
"; decreasing at every step", int(g.decreasing_a.values[free].sum()), "; tau_4 < tau_1", int(g.tau4_below_tau1_a.values[free].sum()))
dec3 = (tau[:, 0] > tau[:, 1]) & (tau[:, 1] > tau[:, 2])
print(" sessions with tau_1 to tau_3 off the edge:", int(free3.sum()), "; tau_1 > tau_2 > tau_3 in", int(dec3[free3].sum()),
"; tau_3 < tau_1 in", int((tau[:, 2] < tau[:, 0])[free3].sum()))
print(" mean tau_n with the values at the edge left out:", [round(tau[~edge[:, m], m].mean(), 2) for m in range(4)])
session protocol tau1_a tau2_a tau3_a tau4_a decreasing_a tau4_below_tau1_a tau1_b tau2_b tau3_b tau4_b decreasing_b tau4_below_tau1_b 03-02-25-aj_sig1 1 4.85 3.85 3.25 3.85 False True 4.85 3.85 3.25 3.85 False True 03-03-05-ab_sig1 1 5.20 4.00 2.90 4.00 False True 5.20 4.00 2.90 4.00 False True 03-03-05-ab_sig2 2 5.35 4.15 -0.00 0.75 False True 5.35 4.15 -0.60 -4.90 True True 03-03-05-ac_sig2 4 3.95 1.95 0.40 -0.00 True True 3.95 1.95 0.40 -1.35 True True 03-03-14-ad_sig1 1 5.00 3.50 -0.00 -0.00 True True 5.00 3.50 -0.50 -4.25 True True 03-03-14-ad_sig2 5 5.30 3.20 1.55 0.45 True True 5.30 3.20 1.55 0.45 True True 03-03-18-ad_sig1 1 5.25 4.30 3.25 4.20 False True 5.25 4.30 3.25 4.20 False True 03-03-18-ae_sig1 1 5.60 4.45 4.15 4.00 True True 5.60 4.45 4.15 4.00 True True 03-03-18-ae_sig2 3 5.60 4.35 4.15 4.35 False True 5.60 4.35 4.15 4.35 False True 03-03-21-aa_sig1 1 6.35 5.20 5.05 5.50 False True 6.35 5.20 5.05 5.50 False True 03-03-21-ab_sig1 1 6.30 5.30 4.60 5.25 False True 6.30 5.30 4.60 5.25 False True 03-03-21-ab_sig2 6 6.30 6.15 6.25 6.45 False False 6.30 6.15 6.25 -1.75 False True 03-03-21-ac_sig1 1 5.90 4.95 3.40 -0.00 True True 5.90 4.95 3.40 -0.55 True True 03-03-21-ac_sig2 4 8.05 4.20 2.15 0.25 True True 8.05 4.20 2.15 0.25 True True 03-04-11-ab_sig2 3 6.85 4.90 -0.00 4.65 False True 6.85 4.90 -0.00 4.65 False True 03-04-15-aa_sig2 5 5.75 3.55 2.25 0.95 True True 5.75 3.55 2.25 0.95 True True 03-04-16-aa_sig1 1 3.50 2.70 1.50 2.60 False True 3.50 2.70 1.50 -2.30 True True 03-04-16-ab_sig1 1 3.40 2.60 1.40 0.10 True True 3.40 2.60 1.40 0.10 True True 03-04-16-ac_sig1 1 3.35 2.55 1.45 0.60 True True 3.35 2.55 1.45 0.60 True True 03-04-16-ac_sig2 4 0.60 -0.00 -0.00 -0.00 True True 0.60 -0.95 -1.90 -2.70 True True 03-04-16-ae_sig2 3 3.45 2.55 1.70 1.85 False True 3.45 2.55 1.70 1.85 False True 03-04-16-af_sig1 1 3.45 2.55 0.95 2.35 False True 3.45 2.55 0.95 -3.10 True True 03-04-16-af_sig2 3 3.50 2.50 1.55 1.80 False True 3.50 2.50 1.55 -2.45 True True 03-04-23-aa_sig1 1 7.70 6.15 6.60 6.40 False True 7.70 6.15 6.60 -1.00 False True 03-04-23-aa_sig2 5 8.20 6.10 4.60 4.00 True True 8.20 6.10 4.60 4.00 True True 03-04-23-ac_sig1 1 8.55 6.85 3.20 -0.00 True True 8.55 6.85 3.20 -0.90 True True 03-04-23-ac_sig2 3 8.25 7.10 6.80 1.35 True True 8.25 7.10 6.80 1.35 True True 03-04-23-ad_sig1 1 6.40 4.95 5.40 5.45 False True 6.40 4.95 5.40 -3.80 False True 03-04-23-ae_sig1 1 6.40 4.80 5.20 5.55 False True 6.40 4.80 5.20 -3.15 False True 03-04-23-ae_sig2 3 6.25 5.25 5.20 -0.00 True True 6.25 5.25 5.20 -0.25 True True 03-04-23-ag_sig1 1 6.70 5.10 5.60 0.60 False True 6.70 5.10 5.60 -3.50 False True 03-04-23-ah_sig1 1 6.90 5.10 5.75 1.20 False True 6.90 5.10 5.75 -3.90 False True 03-04-23-ah_sig2 5 6.85 4.70 3.15 4.65 False True 6.85 4.70 3.15 -4.55 True True 03-04-23-ai_sig1 1 7.25 5.25 5.95 1.35 False True 7.25 5.25 5.95 -4.70 False True 03-04-23-ai_sig2 4 11.85 4.60 -0.00 -0.00 True True 11.85 4.60 -1.15 -1.85 True True 03-04-23-aj_sig1 1 7.00 5.05 5.90 1.10 False True 7.00 5.05 5.90 -4.05 False True 03-04-25-ad_sig2 3 5.15 3.80 3.85 2.90 False True 5.15 3.80 3.85 -2.30 False True 03-04-25-af_sig1 1 6.35 5.00 5.15 5.30 False True 6.35 5.00 5.15 -1.45 False True 03-04-25-ag_sig1 1 6.30 4.95 5.20 5.05 False True 6.30 4.95 5.20 5.05 False True 03-04-30-aa_sig1 1 7.10 5.65 6.35 16.10 False False 7.10 5.65 6.35 -2.95 False True 03-04-30-aa_sig2 5 7.45 5.40 3.85 3.60 True True 7.45 5.40 3.85 3.60 True True 03-05-02-ab_sig1 1 6.00 3.75 4.30 1.00 False True 6.00 3.75 -1.30 -4.65 True True 03-05-06-aa_sig2 5 6.60 4.40 3.10 1.85 True True 6.60 4.40 3.10 1.85 True True 03-05-06-ab_sig2 4 6.10 2.50 0.65 -0.00 True True 6.10 2.50 0.65 -1.35 True True 03-05-09-ab_sig1 1 6.35 4.90 5.45 5.75 False True 6.35 4.90 5.45 -2.75 False True 03-05-09-ab_sig2 5 6.00 4.45 2.85 2.10 True True 6.00 4.45 2.85 2.10 True True 03-05-09-ac_sig2 4 3.65 2.10 1.35 -0.00 True True 3.65 2.10 1.35 -0.75 True True 03-06-27-ab_sig1 1 6.35 5.30 5.55 -0.00 False True 6.35 5.30 5.55 -1.20 False True 03-06-27-ac_sig1 1 6.35 5.15 5.50 -0.00 False True 6.35 5.15 5.50 -1.65 False True 03-06-27-ae_sig1 1 3.85 2.90 1.55 2.80 False True 3.85 2.90 1.55 -2.40 True True 03-06-27-af_sig1 1 4.05 3.00 1.65 3.25 False True 4.05 3.00 1.65 -2.30 True True 03-06-27-ag_sig1 1 4.30 3.20 2.15 2.40 False True 4.30 3.20 2.15 -4.45 True True 03-07-03-aa_sig1 1 4.40 3.55 2.30 3.65 False True 4.40 3.55 2.30 -1.40 True True kept min_bursts search sessions decreasing tau4_below_tau1 mean_tau_ms negative_tau wilcoxon_p c 2 a 53 20 51 [5.8, 4.23, 3.32, 2.67] 0 4.483954e-09 c 2 b 53 30 53 [5.8, 4.22, 3.14, -0.5] 33 2.384565e-10 c 20 a 39 17 39 [6.01, 4.26, 3.33, 2.12] 0 5.239677e-08 c 20 b 39 22 39 [6.01, 4.24, 3.08, -0.36] 23 5.247481e-08 c 50 a 28 13 28 [5.88, 4.08, 3.14, 2.05] 0 3.781414e-06 c 50 b 28 16 28 [5.88, 4.04, 3.01, -0.38] 17 3.781414e-06 a 2 a 71 30 69 [5.75, 4.15, 3.14, 2.46] 0 4.457235e-12 a 2 b 71 43 71 [5.75, 4.13, 2.88, -0.55] 45 2.427312e-13 a 20 a 52 24 52 [5.95, 4.18, 3.17, 2.14] 0 3.495086e-10 a 20 b 52 30 52 [5.95, 4.16, 2.93, -0.38] 31 3.498010e-10 a 50 a 38 18 38 [6.0, 4.12, 3.07, 2.14] 0 7.727609e-08 a 50 b 38 22 38 [6.0, 4.1, 2.94, -0.44] 23 7.727609e-08 kept (c), search before onset: sessions with a tau_n at the 0 ms edge 13 of 53 ; values at the edge 17 (per n [0, 1, 5, 11] ) sessions with no tau_n at the edge: 40 ; their mean tau_n [5.77, 4.33, 3.77, 3.4] ; decreasing at every step 11 ; tau_4 < tau_1 38 sessions with tau_1 to tau_3 off the edge: 48 ; tau_1 > tau_2 > tau_3 in 32 ; tau_3 < tau_1 in 48 mean tau_n with the values at the edge left out: [np.float64(5.8), np.float64(4.32), np.float64(3.67), np.float64(3.37)] kept (a), search before onset: sessions with a tau_n at the 0 ms edge 19 of 71 ; values at the edge 26 (per n [0, 2, 7, 17] ) sessions with no tau_n at the edge: 52 ; their mean tau_n [5.9, 4.32, 3.53, 3.26] ; decreasing at every step 16 ; tau_4 < tau_1 50 sessions with tau_1 to tau_3 off the edge: 64 ; tau_1 > tau_2 > tau_3 in 44 ; tau_3 < tau_1 in 64 mean tau_n with the values at the edge left out: [np.float64(5.75), np.float64(4.27), np.float64(3.48), np.float64(3.24)]
Result: reproduced (reading b). With the central maximum searched before the burst onset, the mean latency over the 53 sessions of selection (c) falls with n: 5.80, 4.23, 3.32 and 2.67 ms for n = 1 to 4; no tau_n is negative, which holds by construction as the search ends at 0 ms. In 13 of the 53 sessions 17 values of tau_n sit on that edge (1 for n = 2, 5 for n = 3, 11 for n = 4), and the mean tau_4 includes them. Leaving those values out, the means are 5.80, 4.32, 3.67 and 3.37 ms, still falling; over the 40 sessions with no value on the edge they are 5.77, 4.33, 3.77 and 3.40 ms, and of the 48 sessions with tau_1 to tau_3 off the edge, 32 have tau_1 > tau_2 > tau_3. 51 of 53 sessions have tau_4 < tau_1 (Wilcoxon p = 4.5e-9), but only 20 of 53 decrease at every step (reading a, not a majority). With the maximum searched over the whole window, as in the first round, the means are 5.80, 4.22, 3.14 and -0.50 ms, with a negative tau_n in 33 sessions, where the 4BTA peaks after the onset. With at least 20 or 50 bursts of each n the means also fall (6.01 to 2.12 and 5.88 to 2.05 ms), and over selection (a) they fall from 5.75 to 2.46 ms.
fig8e-sigma-minimum¶
The standard deviation of the stimuli preceding an n-burst has a minimum at a position nearly the same for all n; over the bursting cells the SD of the minimum's position across n = 1 to 4 averages 0.4 ms (SD 0.7 ms), roughly 0.05 times the inverse cutoff frequency.
sigma_n(tau) by Eq. 14. Window for the minimum (a) -25 to 0 ms, (b) -25 to +15 ms. The SD across
the four n is the sample SD (N - 1), as MATLAB's std; the population SD (N) is shown too. With
few segments sigma_n is noisy and its minimum wanders, so the result is shown for sessions with
at least 2, 20 and 50 bursts of each n.
rows = []
for s in have4:
B = BTA[s.name]
tt = B["tt"]
row = dict(session=s.name, cutoff=s.cutoff)
for w, sel in [("a", tt <= 0), ("b", tt <= tt.max())]:
pos = [tt[sel][np.argmin(B["seg"][m].std(0, ddof=1)[sel])] * 1e3 for m in range(1, 5)]
row[f"pos_{w}"] = pos
row[f"sd_{w}"] = np.std(pos, ddof=1)
row[f"sd0_{w}"] = np.std(pos)
rows.append(row)
SG = pd.DataFrame(rows)
out = []
for kr, k, w in itertools.product("ca", MIN_BURSTS, "ab"):
g = subset(SG, kr, k)
for col in [f"sd_{w}", f"sd0_{w}"]:
out.append(dict(kept=kr, min_bursts=k, window=w, sd="N-1" if col.startswith("sd_") else "N", sessions=len(g),
mean_sd_ms=round(g[col].mean(), 3), sd_of_sd_ms=round(g[col].std(ddof=1), 3), median_ms=round(g[col].median(), 3),
times_cutoff=round((g[col] * g.cutoff * 1e-3).mean(), 4)))
print(pd.DataFrame(out).to_string(index=False))
print("\nkept (c), at least 2 bursts: median position of the minimum (ms), window (a), per n:",
np.median(np.array(subset(SG, "c", 2).pos_a.tolist()), 0))
kept min_bursts window sd sessions mean_sd_ms sd_of_sd_ms median_ms times_cutoff c 2 a N-1 53 2.366 2.442 1.958 0.4005 c 2 a N 53 2.049 2.115 1.696 0.3468 c 2 b N-1 53 2.835 2.704 2.295 0.4765 c 2 b N 53 2.455 2.342 1.988 0.4127 c 20 a N-1 39 2.136 2.423 1.898 0.2994 c 20 a N 39 1.850 2.098 1.644 0.2593 c 20 b N-1 39 2.218 2.540 2.028 0.3029 c 20 b N 39 1.920 2.199 1.756 0.2623 c 50 a N-1 28 1.744 2.359 1.775 0.1863 c 50 a N 28 1.511 2.043 1.538 0.1613 c 50 b N-1 28 1.858 2.543 1.791 0.1911 c 50 b N 28 1.609 2.202 1.551 0.1655 a 2 a N-1 71 2.164 2.404 1.834 0.3640 a 2 a N 71 1.874 2.082 1.589 0.3153 a 2 b N-1 71 2.641 2.678 2.148 0.4346 a 2 b N 71 2.287 2.319 1.860 0.3763 a 20 a N-1 52 1.974 2.467 1.775 0.2836 a 20 a N 52 1.709 2.136 1.538 0.2456 a 20 b N-1 52 2.207 2.641 1.993 0.3050 a 20 b N 52 1.912 2.287 1.726 0.2641 a 50 a N-1 38 1.459 2.119 0.919 0.1511 a 50 a N 38 1.264 1.835 0.796 0.1308 a 50 b N-1 38 1.589 2.280 1.273 0.1558 a 50 b N 38 1.376 1.974 1.102 0.1349 kept (c), at least 2 bursts: median position of the minimum (ms), window (a), per n: [-6.85 -6.3 -4.8 -6.6 ]
Result: not reproduced. The minimum of sigma_n lies near -6 ms for every n (medians -6.85, -6.3, -4.8 and -6.6 ms over the 53 sessions of selection (c)), but the SD of its position across n = 1 to 4 averages 2.37 ms (SD 2.44) in window (a) and 2.84 ms in window (b), outside 0.4 +/- 0.7, and 2.05 and 2.46 ms with the population SD; times the cutoff it is 0.35 to 0.48, against 0.05. Requiring more bursts of each n narrows it but not to the paper's value: 2.14 ms with at least 20 (39 sessions) and 1.74 ms with at least 50 (28 sessions; 0.19 times the cutoff), and 1.46 ms over the 38 such sessions of selection (a).
sigma-n-below-sta-sigma¶
In 98% of bursting cells, for all n, the SD of the stimuli preceding an n-burst is smaller than the SD of the stimuli preceding all spikes, by a ratio of minima of 1.62 on average (SD 0.56).
"Smaller" read as (a) the minimum of sigma_n below the minimum of sigma, (b) sigma_n below sigma at every tau in -25 to +15 ms; n = 1 to 4. The ratio is min sigma over the mean over n of min sigma_n, minima over -25 to +15 ms (and over -25 to 0 ms as a variant).
The minimum over 800 time points of an SD estimated from few segments is biased low, so a smaller minimum for the n-bursts could come from their smaller number alone. As a control (not in the paper), sigma is also computed from random subsets of the segments preceding all spikes, of the same size as each n's, 20 draws per n, and the minimum of sigma_n is compared with the mean minimum of these size-matched sigmas.
rng8 = np.random.default_rng(14)
rows = []
for s in have4:
B = BTA[s.name]
tt = B["tt"]
sig = B["sta_seg"].std(0, ddof=1)
sn = [B["seg"][m].std(0, ddof=1) for m in range(1, 5)]
matched = [np.mean([B["sta_seg"][rng8.choice(len(B["sta_seg"]), len(B["seg"][m]), replace=False)].std(0, ddof=1).min()
for _ in range(20)]) for m in range(1, 5)]
rows.append(dict(session=s.name, a=all(x.min() < sig.min() for x in sn), b=all((x < sig).all() for x in sn),
ratio=sig.min() / np.mean([x.min() for x in sn]),
ratio_pre=sig[tt <= 0].min() / np.mean([x[tt <= 0].min() for x in sn]),
**{f"ratio_n{m}": sig.min() / x.min() for m, x in zip(range(1, 5), sn)},
a_matched=all(x.min() < mm for x, mm in zip(sn, matched)),
ratio_matched=np.mean(matched) / np.mean([x.min() for x in sn]),
matched_ratio_of_sta=sig.min() / np.mean(matched)))
SB = pd.DataFrame(rows)
out = []
for kr, k in itertools.product("ca", MIN_BURSTS):
g = subset(SB, kr, k)
out.append(dict(kept=kr, min_bursts=k, sessions=len(g), a_pct=round(100 * g.a.mean(), 1), b_pct=round(100 * g.b.mean(), 1),
ratio=round(g.ratio.mean(), 3), ratio_sd=round(g.ratio.std(ddof=1), 3), ratio_pre=round(g.ratio_pre.mean(), 3),
**{f"ratio_n{m}": round(g[f"ratio_n{m}"].mean(), 3) for m in range(1, 5)},
a_vs_matched_pct=round(100 * g.a_matched.mean(), 1), ratio_vs_matched=round(g.ratio_matched.mean(), 3),
sta_vs_matched=round(g.matched_ratio_of_sta.mean(), 3)))
print(pd.DataFrame(out).to_string(index=False))
kept min_bursts sessions a_pct b_pct ratio ratio_sd ratio_pre ratio_n1 ratio_n2 ratio_n3 ratio_n4 a_vs_matched_pct ratio_vs_matched sta_vs_matched c 2 53 98.1 1.9 1.596 0.280 1.581 1.269 1.449 1.623 23.528 94.3 1.517 1.056 c 20 39 97.4 2.6 1.562 0.288 1.557 1.278 1.463 1.635 2.376 97.4 1.549 1.008 c 50 28 96.4 3.6 1.571 0.317 1.565 1.314 1.507 1.623 2.185 96.4 1.563 1.005 a 2 71 98.6 1.4 1.604 0.278 1.590 1.283 1.467 1.634 18.560 94.4 1.533 1.049 a 20 52 98.1 1.9 1.558 0.278 1.550 1.294 1.478 1.614 2.265 96.2 1.544 1.009 a 50 38 97.4 2.6 1.581 0.297 1.574 1.329 1.523 1.620 2.160 97.4 1.572 1.006
Result: reproduced. In 98.1% of the 53 sessions of selection (c) the minimum of sigma_n is below the minimum of sigma for all n = 1 to 4 (reading a), within 98 +/- 9.8, and the ratio of the minima averages 1.596 (SD 0.280), within 1.62 +/- 0.56. Read as sigma_n below sigma at every tau (b), it holds in 1.9%. The result is not an effect of the smaller number of n-bursts: sigma from random subsets of the spike-triggered segments of the same size has a minimum about 5% below that of all segments (0.5 to 0.8% with at least 20 or 50 bursts of each n), the minimum of sigma_n is below that size-matched minimum for every n in 94.3% of sessions, and the ratio to it is 1.517. With at least 20 or 50 bursts of each n the values are 97.4% and 1.562, and 96.4% and 1.571. The ratios per n are 1.27, 1.45 and 1.62 for n = 1 to 3; for n = 4 the mean is 23.5, inflated by sessions with few 4-bursts, and 2.38 with at least 20.
fig9-nbta-not-concatenated-1bta¶
The nBTA differs from n 1BTAs shifted by the burst's own ISIs and summed; for a large fraction of cells the difference is significant, the fraction diminishing as n grows; results are similar with shifted STAs.
f_n(t) = sum over the burst's spikes of 1BTA(t - t_i) (Eq. 15). Summing n copies of a stimulus in dB only compares with the nBTA once the mean stimulus is subtracted (as the paper does for Figure 9C); this is taken as the reading of Eq. 15, an assumption, and the average of the shifted copies (the sum over n) is shown as a variant. The 1BTA (or STA) is computed from -60 to +15 ms so that its shifts cover the window. Comparison over the interval between the minima flanking the central maximum of the nBTA (searched before the onset): two-way ANOVAs (frequency bin x group) on the real and imaginary parts of the Fourier transform, at 0.01 / 2, with observations (a) each burst's segment against its own f_n (paired: ANOVA on the differences), (b) the segments and the f_n as two groups. Then t-tests at each time point (p < 0.01) and sessions with more than 70% of the interval significant.
def paired_anova(d):
# Differences segment - f_n per burst (rows) and frequency bin (columns): p of a mean
# difference (group main effect) and of its variation across bins (the interaction).
Nb, F = d.shape
m = d.mean()
mf = d.mean(0)
ss_err = ((d - mf) ** 2).sum()
df_err = (Nb - 1) * F
ms = ss_err / df_err
p_main = stats.f.sf(Nb * F * m ** 2 / ms, 1, df_err)
p_int = stats.f.sf(Nb * ((mf - m) ** 2).sum() / (F - 1) / ms, F - 1, df_err)
return min(p_main, p_int)
rows = []
for s in BURSTING_A:
B = BTA[s.name]
st = stim_db(s)
ms = B["mean_stim"]
tr, on, n = B["tr"], B["on"], B["n"]
long_t0 = 0.06
one, tl = ia1.segments(st, on[n == 1], pre=long_t0)
allsp, _ = ia1.segments(st, spikes_in(s), pre=long_t0)
kernels = {"1BTA": one.mean(0) - ms, "STA": allsp.mean(0) - ms}
i_zero = int(round(long_t0 * ia1.STIM_RATE))
for m in [2, 3, 4]:
if m not in B["seg"]:
continue
idx = np.flatnonzero(n == m)
seg = B["seg"][m] - ms
tt = B["tt"]
ip, a, b = ia1.central_peak(seg.mean(0), tt, SEARCH["a"])
if b - a < 3:
continue
# spike offsets of each burst relative to its onset
trial_spikes = s.trials
for kname, ker in kernels.items():
for scale in ["sum", "mean"]:
fn = np.zeros_like(seg)
for r, k in enumerate(idx):
t = trial_spikes[tr[k]]
j = np.searchsorted(t, on[k])
offs = t[j:j + m] - on[k]
for o in offs:
sh = int(round(o * ia1.STIM_RATE))
fn[r] += ker[i_zero - 500 - sh:i_zero + 300 - sh]
if scale == "mean":
fn[r] /= m
x, y = seg[:, a:b + 1], fn[:, a:b + 1]
X, Y = np.fft.rfft(x, axis=1), np.fft.rfft(y, axis=1)
pa = min(paired_anova(np.real(X - Y)), paired_anova(np.imag(X - Y)[:, 1:]) if X.shape[1] > 2 else 1)
pb = min(min(ia1.two_way_anova([np.real(X), np.real(Y)])), min(ia1.two_way_anova([np.imag(X)[:, 1:], np.imag(Y)[:, 1:]])) if X.shape[1] > 2 else 1)
tp = stats.ttest_ind(x, y, axis=0).pvalue
rows.append(dict(session=s.name, kept_c=s.name in POP["c"], n=m, kernel=kname, scale=scale, bursts=len(idx), anova_a=pa < 0.005, anova_b=pb < 0.005,
t_frac=(tp < 0.01).mean()))
N9 = pd.DataFrame(rows)
summ = []
for (kname, scale, m), g0 in N9.groupby(["kernel", "scale", "n"]):
for kr, g in [("c", g0[g0.kept_c]), ("a", g0)]:
summ.append(dict(kept=kr, kernel=kname, scale=scale, n=m, sessions=len(g), anova_a_pct=round(100 * g.anova_a.mean(), 1),
anova_b_pct=round(100 * g.anova_b.mean(), 1),
t70_a_pct=round(100 * (g.anova_a & (g.t_frac > 0.7)).mean(), 1),
t70_b_pct=round(100 * (g.anova_b & (g.t_frac > 0.7)).mean(), 1)))
print(pd.DataFrame(summ).to_string(index=False))
kept kernel scale n sessions anova_a_pct anova_b_pct t70_a_pct t70_b_pct c 1BTA mean 2 64 100.0 100.0 100.0 100.0 a 1BTA mean 2 84 100.0 100.0 100.0 100.0 c 1BTA mean 3 63 100.0 100.0 93.7 93.7 a 1BTA mean 3 83 100.0 100.0 94.0 94.0 c 1BTA mean 4 53 100.0 100.0 66.0 66.0 a 1BTA mean 4 71 100.0 100.0 71.8 71.8 c 1BTA sum 2 64 100.0 100.0 93.8 93.8 a 1BTA sum 2 84 100.0 100.0 94.0 94.0 c 1BTA sum 3 63 100.0 100.0 84.1 84.1 a 1BTA sum 3 83 100.0 100.0 85.5 85.5 c 1BTA sum 4 53 100.0 98.1 52.8 52.8 a 1BTA sum 4 71 100.0 98.6 54.9 54.9 c STA mean 2 64 100.0 100.0 95.3 95.3 a STA mean 2 84 100.0 100.0 95.2 95.2 c STA mean 3 63 100.0 100.0 93.7 93.7 a STA mean 3 83 100.0 100.0 92.8 92.8 c STA mean 4 53 100.0 100.0 62.3 62.3 a STA mean 4 71 100.0 100.0 66.2 66.2 c STA sum 2 64 100.0 100.0 100.0 100.0 a STA sum 2 84 100.0 100.0 100.0 100.0 c STA sum 3 63 100.0 100.0 82.5 82.5 a STA sum 3 83 100.0 100.0 83.1 83.1 c STA sum 4 53 98.1 94.3 37.7 37.7 a STA sum 4 71 97.2 94.4 39.4 39.4
Result: reproduced (resting on subtracting the mean stimulus before summing the 1BTAs). Over the bursting sessions of selection (c) the nBTA differs from the summed 1BTAs by the ANOVA in 100% of the sessions for n = 2, 3 and 4 with paired observations (a), and in 100%, 100% and 98.1% with grouped ones (b), so the fraction does not grow with n; the share with more than 70% of the interval significant falls from 93.8% (n = 2) to 84.1% and 52.8% (n = 4). With the STA the ANOVA fractions are 100%, 100% and 98.1% (a) and the 70% share falls from 100% to 37.7%; with the average of the shifted copies in place of their sum the results are alike, and so are those over selection (a). An ANOVA that rejects in nearly every session separates little; the criterion specifies it.
fig9c-population-nbta¶
In the seven cells driven at 100 Hz cutoff and 6 dB, the population-average nBTAs show that bursts with more spikes follow higher or wider stimulus deflections, with large spread across cells.
Each session's 1BTA, 2BTA and 3BTA minus its mean stimulus, shifted so that its central maximum sits at the mean tau_n over the sessions (shifting by tau_n, averaging and shifting back), averaged over the bursting protocol-5 sessions with at least two bursts of n = 1 to 3. Width (a) full width at half height, (b) between the minima flanking the central maximum. The central maximum is searched before the onset and, as a variant, over the whole window.
rows = []
fig, ax = plt.subplots(figsize=(6, 4))
for kr, sr in itertools.product("ca", "ab"):
p5 = [s for s in BURSTING_A if s.name in POP[kr] and s.protocol == 5 and all(m in BTA[s.name]["seg"] for m in (1, 2, 3))]
if sr == "a":
print(f"kept ({kr}): protocol-5 sessions kept", sum(s.protocol == 5 and s.name in POP[kr] for s in KEPT_A),
"; bursting with n = 1 to 3:", len(p5), [s.name for s in p5])
cl = D.set_index("session").loc[[s.name for s in p5]]
ids = [f"cell {r.cell}" if r.cell >= 0 else r.exp for r in cl.itertuples()]
print(" their cells (rows of singlecells.txt):", ids, "; distinct:", len(set(ids)))
for m in (1, 2, 3):
tt = BTA[p5[0].name]["tt"]
btas, peaks = [], []
for s in p5:
b = BTA[s.name]["seg"][m].mean(0) - BTA[s.name]["mean_stim"]
btas.append(b)
peaks.append(ia1.central_peak(b, tt, SEARCH[sr])[0])
target = int(round(np.mean(peaks)))
aligned = np.array([np.roll(b, target - p) for b, p in zip(btas, peaks)])
avg = aligned.mean(0)
ip, a, b2 = ia1.central_peak(avg, tt, SEARCH[sr])
rows.append(dict(kept=kr, search=sr, sessions=len(p5), n=m, height_dB=avg[ip], fwhm_ms=ia1.fwhm(avg, ip, 0.0) / 20,
minima_width_ms=(b2 - a) / 20, spread_at_peak_dB=aligned[:, ip].std(ddof=1), tau_ms=-tt[ip] * 1e3))
if kr == "c" and sr == "a":
ax.plot(tt * 1e3, avg, label=f"{m}BTA")
ax.fill_between(tt * 1e3, avg - aligned.std(0, ddof=1), avg + aligned.std(0, ddof=1), alpha=0.2)
ax.set(xlabel="time from burst onset (ms)", ylabel="stimulus minus its mean (dB)", title="kept (c), maximum searched before onset")
ax.legend()
plt.show()
F9C = pd.DataFrame(rows)
print(F9C.round(2).to_string(index=False))
kept (c): protocol-5 sessions kept 8 ; bursting with n = 1 to 3: 7 ['03-03-14-ad_sig2', '03-04-15-aa_sig2', '03-04-23-aa_sig2', '03-04-23-ah_sig2', '03-04-30-aa_sig2', '03-05-06-aa_sig2', '03-05-09-ab_sig2'] their cells (rows of singlecells.txt): ['cell 11', 'cell 21', 'cell 25', 'cell 27', 'cell 31', 'cell 35', 'cell 37'] ; distinct: 7 kept (a): protocol-5 sessions kept 11 ; bursting with n = 1 to 3: 10 ['03-03-14-aa_sig2', '03-03-14-ab_sig2', '03-03-14-ac_sig2', '03-03-14-ad_sig2', '03-04-15-aa_sig2', '03-04-23-aa_sig2', '03-04-23-ah_sig2', '03-04-30-aa_sig2', '03-05-06-aa_sig2', '03-05-09-ab_sig2'] their cells (rows of singlecells.txt): ['cell 10', 'cell 11', 'cell 11', 'cell 11', 'cell 21', 'cell 25', 'cell 27', 'cell 31', 'cell 35', 'cell 37'] ; distinct: 8
kept search sessions n height_dB fwhm_ms minima_width_ms spread_at_peak_dB tau_ms c a 7 1 3.13 5.00 12.55 0.54 6.60 c a 7 2 7.03 5.50 13.60 0.83 4.55 c a 7 3 8.46 7.35 17.80 1.51 3.05 c b 7 1 3.13 5.00 12.55 0.54 6.60 c b 7 2 7.03 5.50 13.60 0.83 4.55 c b 7 3 8.46 7.35 17.80 1.51 3.05 a a 10 1 3.22 5.10 12.70 0.57 6.30 a a 10 2 7.11 5.50 13.50 0.74 4.30 a a 10 3 8.60 7.15 17.15 1.28 2.80 a b 10 1 3.22 5.10 12.70 0.57 6.30 a b 10 2 7.11 5.50 13.50 0.74 4.30 a b 10 3 8.60 7.15 17.15 1.28 2.80
Result: reproduced. Selection (c) holds 7 bursting protocol-5 sessions with bursts of n = 1 to
3, from 7 distinct cells (rows of singlecells.txt; 03-04-23-aa and 03-04-23-ah are different
cells recorded on the same date), the paper's seven cells. Their population-average central peak is 3.13, 7.03 and 8.46 dB above
the mean stimulus for n = 1, 2 and 3, and its width grows from 5.00 to 5.50 to 7.35 ms at half
height (12.55, 13.60, 17.80 ms between minima). The SD across sessions at the peak is 0.54, 0.83
and 1.51 dB. The two searches for the maximum give the same values, and the 10 sessions of
selection (a) give 3.22, 7.11 and 8.60 dB.
fig10-burst-n-vs-height¶
The probability of a burst with n spikes after a stimulus maximum depends on its height, low excursions giving no response or an isolated spike and high ones doublets, triplets or 4-bursts; all bursting cells show this.
Methods (Relating burst probabilities to the height of stimulus excursions), per trial: each n-burst shifted back by tau_n (the maximum of the nBTA searched (a) before the onset, (b) over the whole window); each local maximum of the stimulus in 200 to 1000 ms matched to the burst whose shifted onset is closest within [t0 - T, t0 + T]; a burst matched to more than one maximum kept by the nearest, the other maxima counting as no response. T is the width of the central STA peak, (a) between its flanking minima, (b) full width at half height. P(n | h) in windows of half-width 5% of the span of h, at 41 centres. Segregation (a) the h at which P(n | h) peaks increases with n for n = 0 to the largest n with bursts (up to 4); (b) the mean of h weighted by P(n | h) increases with n. Bursting sessions with bursts of n = 1 and 2 at least.
from scipy.signal import argrelmax
def p_n_given_h(s, Twin, sr):
B = BTA[s.name]
st = stim_db(s)
i0, i1 = int(ia1.T0 * ia1.STIM_RATE), int(ia1.T1 * ia1.STIM_RATE)
mx = argrelmax(st[i0:i1])[0] + i0
tmax, h = mx / ia1.STIM_RATE, st[mx]
tt = B["tt"]
tau = {m: -tt[ia1.central_peak(B["seg"][m].mean(0), tt, SEARCH[sr])[0]] for m in B["seg"]}
resp = []
for i in range(s.n_trials):
k = B["tr"] == i
on, n = B["on"][k], B["n"][k]
sh = on - np.array([tau.get(m, tau[max(tau)]) for m in n])
best = np.full(len(tmax), -1)
for j, t0 in enumerate(tmax):
d = np.abs(sh - t0)
if len(d) and d.min() <= Twin:
best[j] = np.argmin(d)
r = np.zeros(len(tmax), int)
for bi in np.unique(best[best >= 0]):
js = np.flatnonzero(best == bi)
jbest = js[np.argmin(np.abs(tmax[js] - sh[bi]))]
r[jbest] = n[bi]
resp.append(r)
resp = np.concatenate(resp)
matched = int((resp > 0).sum())
hh = np.tile(h, s.n_trials)
span = hh.max() - hh.min()
centres = np.linspace(hh.min(), hh.max(), 41)
nmax = min(4, resp.max())
P = np.array([[np.mean(resp[np.abs(hh - c) <= 0.05 * span] == m) if (np.abs(hh - c) <= 0.05 * span).any() else np.nan
for c in centres] for m in range(nmax + 1)])
occupancy = np.array([(np.abs(hh - c) <= 0.05 * span).sum() for c in centres])
return centres, P, nmax, matched, np.bincount(np.minimum(resp, 4), minlength=5), occupancy
rows = []
example = None
for s, sr in itertools.product(BURSTING_A, "ab"):
B = BTA[s.name]
if not (1 in B["seg"] and 2 in B["seg"]):
continue
tt = B["tt"]
sta = B["sta_seg"].mean(0) - B["mean_stim"]
ip, a, b = ia1.central_peak(sta, tt, SEARCH[sr])
for tw, Tw in [("a", (b - a) / ia1.STIM_RATE), ("b", ia1.fwhm(sta, ip, 0.0) / ia1.STIM_RATE)]:
centres, Pm, nmax, matched, nresp, occ = p_n_given_h(s, Tw, sr)
peak_h = [centres[np.nanargmax(p)] for p in Pm]
mean_h = [np.nansum(centres * p) / np.nansum(p) for p in Pm]
rows.append(dict(session=s.name, kept_c=s.name in POP["c"], search=sr, T=tw, T_ms=Tw * 1e3, nmax=nmax, matched=matched, peak_h=np.round(peak_h, 1),
mean_h=np.round(mean_h, 1), n_resp=nresp.tolist(), peak_at_top=[int(np.nanargmax(p) == len(centres) - 1) for p in Pm], top_occupancy=int(occ[-1]), a=all(np.diff(peak_h) > 0), b=all(np.diff(mean_h) > 0),
a3=all(np.diff(peak_h[:4]) > 0), b3=all(np.diff(mean_h[:4]) > 0)))
if example is None and nmax == 4 and tw == "a" and sr == "a" and s.name in POP["c"]:
example = (s.name, centres, Pm)
H = pd.DataFrame(rows)
out = []
for kr, sr, tw in itertools.product("ca", "ab", "ab"):
g = H[(H.kept_c | (kr == "a")) & (H.search == sr) & (H["T"] == tw)]
out.append(dict(kept=kr, search=sr, T=tw, sessions=len(g), median_T_ms=round(g.T_ms.median(), 2), seg_a=int(g.a.sum()),
seg_a_pct=round(100 * g.a.mean(), 1), seg_b=int(g.b.sum()), seg_b_pct=round(100 * g.b.mean(), 1),
n0_3_a=int(g.a3.sum()), n0_3_b=int(g.b3.sum()), median_matched=int(g.matched.median()),
nmax_peak_at_top_pct=round(100 * np.mean([p[-1] for p in g.peak_at_top]), 1),
median_maxima_in_top_window=int(g.top_occupancy.median())))
print(pd.DataFrame(out).to_string(index=False))
side = H[H.kept_c & (H.search == "a")].pivot(index="session", columns="T", values=["matched", "a", "b", "nmax"])
print("\nkept (c), search (a): sessions where the two T readings match a different number of maxima:",
int((side["matched"]["a"] != side["matched"]["b"]).sum()), "of", len(side), "; where flag (a) differs:",
int((side["a"]["a"] != side["a"]["b"]).sum()), "; flag (b) differs:", int((side["b"]["a"] != side["b"]["b"]).sum()))
print(side.head(15).to_string())
ex = H[H.kept_c & (H.search == "a") & (H["T"] == "a")]
print("\nkept (c), search (a), T (a), first sessions: responses per n (0 to 4+), peak h and mean h per n")
print(ex[["session", "nmax", "n_resp", "peak_h", "mean_h", "peak_at_top", "top_occupancy"]].head(10).to_string(index=False))
name, centres, Pm = example
fig, ax = plt.subplots(figsize=(6, 4))
for m, p in enumerate(Pm):
ax.plot(centres, p, label=f"n = {m}")
ax.set(xlabel="height of the stimulus maximum h (dB)", ylabel="P(n | h)", title=name, ylim=(0, 1))
ax.legend()
plt.show()
/tmp/ipykernel_2316/3297635730.py:52: RuntimeWarning: invalid value encountered in scalar divide mean_h = [np.nansum(centres * p) / np.nansum(p) for p in Pm]
kept search T sessions median_T_ms seg_a seg_a_pct seg_b seg_b_pct n0_3_a n0_3_b median_matched nmax_peak_at_top_pct median_maxima_in_top_window
c a a 64 7.40 7 10.9 15 23.4 34 37 8316 29.7 258
c a b 64 3.05 7 10.9 15 23.4 34 36 8173 29.7 258
c b a 64 7.40 7 10.9 16 25.0 34 37 8316 32.8 258
c b b 64 3.05 7 10.9 16 25.0 34 36 8173 32.8 258
a a a 84 7.40 8 9.5 19 22.6 44 48 8316 28.6 258
a a b 84 3.05 8 9.5 19 22.6 44 48 8173 28.6 258
a b a 84 7.40 8 9.5 21 25.0 44 48 8316 32.1 258
a b b 84 3.05 8 9.5 21 25.0 44 48 8173 32.1 258
kept (c), search (a): sessions where the two T readings match a different number of maxima: 58 of 64 ; where flag (a) differs: 0 ; flag (b) differs: 2
matched a b nmax
T a b a b a b a b
session
03-02-25-ac_sig1 8555 8331 False False True True 4 4
03-02-25-ac_sig2 11520 11152 False False True False 3 3
03-02-25-aj_sig1 16532 16092 False False False False 4 4
03-03-05-ab_sig1 39195 37401 False False False False 4 4
03-03-05-ab_sig2 34668 32596 False False False False 4 4
03-03-05-ac_sig1 9555 9022 True True True True 3 3
03-03-05-ac_sig2 1032 1032 False False False False 4 4
03-03-14-ad_sig1 35500 34911 False False False False 4 4
03-03-14-ad_sig2 24744 24040 False False True True 4 4
03-03-18-ad_sig1 21256 20757 False False False False 4 4
03-03-18-ae_sig1 21795 21211 True True True True 4 4
03-03-18-ae_sig2 20180 20068 False False False False 4 4
03-03-21-aa_sig1 7243 7015 False False False False 4 4
03-03-21-ab_sig1 17205 16860 False False False False 4 4
03-03-21-ab_sig2 19299 18228 False False False True 4 4
kept (c), search (a), T (a), first sessions: responses per n (0 to 4+), peak h and mean h per n
session nmax n_resp peak_h mean_h peak_at_top top_occupancy
03-02-25-ac_sig1 4 [8525, 6149, 2044, 361, 1] [51.8, 64.5, 71.5, 74.0, 74.0] [55.8, 65.4, 70.3, 71.5, 73.5] [0, 0, 0, 1, 1] 280
03-02-25-ac_sig2 3 [23060, 11016, 502, 2, 0] [51.7, 69.7, 72.0, 70.8] [57.4, 68.3, 70.5, 71.6] [0, 0, 0, 0] 420
03-02-25-aj_sig1 4 [13358, 12150, 3524, 688, 170] [21.8, 34.5, 43.4, 44.0, 41.5] [26.3, 34.9, 39.4, 39.6, 38.7] [0, 0, 0, 1, 0] 490
03-03-05-ab_sig1 4 [22171, 35953, 2910, 324, 8] [40.8, 53.5, 62.4, 63.0, 60.5] [44.5, 53.6, 60.2, 61.1, 60.6] [0, 0, 0, 1, 0] 1006
03-03-05-ab_sig2 4 [23564, 24742, 6972, 2247, 707] [44.6, 50.9, 54.5, 54.2, 52.8] [48.3, 50.1, 52.1, 52.7, 52.2] [0, 0, 0, 0, 0] 1506
03-03-05-ac_sig1 3 [6305, 8896, 603, 56, 0] [37.8, 53.7, 59.4, 60.0] [41.9, 50.9, 57.4, 58.1] [0, 0, 0, 1] 260
03-03-05-ac_sig2 4 [0, 441, 191, 141, 259] [52.2, 52.2, 53.4, 54.3, 55.7] [nan, 53.4, 54.0, 54.7, 55.7] [0, 0, 0, 0, 0] 129
03-03-14-ad_sig1 4 [25500, 30158, 4880, 427, 35] [59.9, 74.5, 84.0, 77.0, 72.6] [65.4, 76.3, 80.0, 77.2, 74.2] [0, 0, 1, 0, 0] 1000
03-03-14-ad_sig2 4 [6256, 15475, 6461, 2530, 278] [59.1, 69.6, 75.0, 77.3, 77.3] [63.8, 69.9, 74.5, 75.9, 76.6] [0, 0, 0, 1, 1] 4000
03-03-18-ad_sig1 4 [12782, 16862, 3573, 692, 129] [39.5, 52.2, 58.5, 61.7, 59.2] [43.4, 52.4, 57.2, 57.3, 56.7] [0, 0, 0, 1, 0] 558
Result: not reproduced. Over the 64 bursting sessions of selection (c), with each n-burst shifted back by tau_n found before the onset, the h at which P(n | h) peaks increases from n = 0 to the largest n (up to 4) in 7 sessions (10.9%) and the P-weighted mean h in 15 (23.4%), under both widths T, against all bursting cells. With tau_n searched over the whole window the counts are 7 and 16, and over selection (a) 8 and 19 to 21 of 84. Restricted to n = 0 to 3, a variant outside the criterion, the ordering holds in 34 and 36 to 37 of 64 sessions: the 4-bursts, which are few, most often break it. The two T readings match a different number of maxima in 58 of 64 sessions (a median of 8316 against 8173), yet the peak flag changes in none and the mean flag in 2, one each way: the peak of P(n | h) on the 41-point grid of h does not move when about 2% of the matches change, which is why the counts are equal.
firing-rate-no-effect-on-limit¶
The firing rate has practically no effect on the limiting ISI; a 50 Hz increase in rate shifts it by less than 0.4 ms and raises the average n by less than 25%.
Reading (a), the same cell and stimulus at two target rates 50 Hz apart: the headers record target rates of 100 and 90 Hz only (Data check), so no pair exists. Reading (b): linear regressions over the kept sessions of the limiting ISI (bursting sessions) and of the average n on the session's rate, scaled to 50 Hz; the increase of the average n is relative to its mean.
rows = []
for s in KEPT_A:
c = ia1.burst_counts(s.trials, LIM[s.name])
rows.append(dict(session=s.name, kept_c=s.name in POP["c"], protocol=s.protocol, rate=D.loc[D.session == s.name, "rate"].iloc[0],
limit_ms=LIM[s.name] * 1e3, avg_n=ia1.average_n(c)[0]))
RT = pd.DataFrame(rows)
pairs = D[D.protocol.notna()].groupby(["cell", "stimulus"]).desired_rate.agg(lambda x: sorted(set(x.dropna())))
print("cell and stimulus with more than one target rate:", [k for k, v in pairs.items() if len(v) > 1])
for lab, g in [("kept (c)", RT[RT.kept_c]), ("kept (c), protocol 1 only", RT[RT.kept_c & (RT.protocol == 1)]),
("kept (a)", RT), ("kept (a), protocol 1 only", RT[RT.protocol == 1])]:
gb = g[g.limit_ms > 0]
r1 = stats.linregress(gb.rate, gb.limit_ms)
r2 = stats.linregress(g.rate, g.avg_n)
print(f"{lab}: limiting ISI slope {r1.slope:.4f} ms/Hz (p {r1.pvalue:.3g}, n {len(gb)}) -> shift for +50 Hz {50 * r1.slope:.3f} ms;"
f" average n slope {r2.slope:.5f} per Hz (p {r2.pvalue:.3g}, n {len(g)}) -> +50 Hz changes it by {100 * 50 * r2.slope / g.avg_n.mean():.1f}%"
f" of its mean {g.avg_n.mean():.3f}")
cell and stimulus with more than one target rate: [] kept (c): limiting ISI slope -0.0067 ms/Hz (p 0.44, n 64) -> shift for +50 Hz -0.337 ms; average n slope 0.00386 per Hz (p 0.132, n 89) -> +50 Hz changes it by 15.7% of its mean 1.231 kept (c), protocol 1 only: limiting ISI slope -0.0011 ms/Hz (p 0.914, n 40) -> shift for +50 Hz -0.057 ms; average n slope 0.00450 per Hz (p 0.00169, n 47) -> +50 Hz changes it by 18.9% of its mean 1.193 kept (a): limiting ISI slope -0.0038 ms/Hz (p 0.488, n 84) -> shift for +50 Hz -0.189 ms; average n slope 0.00338 per Hz (p 0.0688, n 115) -> +50 Hz changes it by 13.6% of its mean 1.242 kept (a), protocol 1 only: limiting ISI slope -0.0001 ms/Hz (p 0.981, n 52) -> shift for +50 Hz -0.007 ms; average n slope 0.00494 per Hz (p 5.92e-07, n 60) -> +50 Hz changes it by 20.3% of its mean 1.217
Result: reproduced (reading b). No cell was recorded at two target rates (the headers hold 100 and 90 Hz only), so reading (a) cannot be formed. Across the sessions of selection (c) a 50 Hz higher rate goes with a limiting ISI shorter by 0.337 ms (p = 0.44, 64 bursting sessions), below 0.44 ms, and an average n higher by 15.7% (p = 0.132, 89 sessions), below 27.5%; within protocol 1 alone, 0.057 ms and 18.9%. Over selection (a), 0.189 ms and 13.6%. These are differences between cells, not the within-cell change the paper may have measured.
model-reproduces-coding-trends (blocked)¶
The threshold-linear model with refractoriness reproduces the coding trends of the data, at single-cell and population level, though it is less precise in time than the real cells.
Depends on model-matches-classification, which is blocked (data_missing, no thresholds s_TH).
Summary¶
A failure to reproduce here may be a failure of this attempt, not of the paper: see the reasoning and the differences from the paper's method.
| Claim | Outcome |
|---|---|
| data-sessions-and-trials | partially reproduced |
| table1-sessions-per-protocol | not reproduced |
| firing-rate-uniform | reproduced |
| threshold-population | blocked: data missing |
| table1-bursting-sessions | not reproduced |
| table1-isolated-spikes | not reproduced |
| table1-burst-length | not reproduced |
| fig5-burst-probability-vs-stimulus | reproduced (assumed) |
| fig3-correlation-shape | reproduced (assumed) |
| method-correlation-error | reproduced (assumed) |
| model-matches-classification | blocked: data missing |
| fig6-bursts-less-correlated | reproduced (assumed) |
| fig7ab-information-per-burst-increases | partially reproduced (assumed) |
| fig7cd-occurrence-decreases | partially reproduced (assumed) |
| fig7ef-jitter-independent-of-n | reproduced (assumed) |
| fig7gh-information-rate | not reproduced (assumed) |
| direct-method-comparison | partially reproduced (assumed) |
| fig8-nbta-differ | partially reproduced (assumed) |
| fig8-latency-decreases | reproduced (assumed) |
| fig8e-sigma-minimum | not reproduced (assumed) |
| sigma-n-below-sta-sigma | reproduced (assumed) |
| fig9-nbta-not-concatenated-1bta | reproduced (assumed) |
| fig9c-population-nbta | reproduced (assumed) |
| fig10-burst-n-vs-height | not reproduced (assumed) |
| firing-rate-no-effect-on-limit | reproduced (assumed) |
| model-reproduces-coding-trends | blocked: data missing |
How the analysis differs from the paper's:
- Population. The paper keeps 86 of its 132 sessions by the rate SD (< 35 Hz) and a rate near 100 Hz. In the data set that rule keeps 115 of the 116 sessions of the paper's seven protocols. Leaving out also the 9 sessions with fewer than 98 trials (the paper gives 98 to 503) and the sessions sorted off-line (selection c) keeps 89, with 47, 6, 7, 7, 8, 7 and 7 per protocol against the paper's 43, 6, 7, 7, 8, 7 and 8, and bursting sessions per protocol equal to Table 1's in six of seven, but from 29 cells on 18 dates. Selection (c) is the main population from firing-rate-uniform on, an assumption: the off-line exclusion has no passage in the paper behind it, and as it was chosen in part from Table 1's counts it is only a variant for the claims on those counts. The results over the 115 sessions of selection (a), the population of the first round, are shown beside it.
- Burst identification. The Methods leave open which maximum is the first peak, which minimum bounds it and what sigma and N are in the error bar. The first local maximum (peak reading a) falls on edge effects of the refractory plateau, where C_s has no spike pairs, and leaves 11 of 89 sessions bursting; started after the plateau (peak reading c) it gives 61, close to the 64 of the main reading. With N the number of trials the error bars are so wide that no session bursts. The main reading is the largest value below 5 ms, the first significant minimum after it, and N the number of product terms, under which the error bar is 2.6% of the span of C_s (the paper's 3.4%). Outcomes that rest on bursts are marked as resting on this choice.
- Burst-triggered averages. The central maximum of an nBTA, which gives tau_n, is searched before the burst onset (the paper calls tau_n the delay between stimulus upstroke and burst generation); the first round searched the whole window, which put the 4BTA's maximum after the onset in many sessions. Sessions need at least two bursts of each n; results with at least 20 and 50 are shown.
- Information. I_n(1) and I'_n are plug-in estimates, as the criterion gives them. A baseline with each trial's burst onsets shifted at random gives a larger rise of I_n(1) with n than the data, so the bias of the estimate is of the order of that rise.
- Thresholds. The data set holds no threshold s_TH, so the model of Eq. 2 and the claims on it are blocked.
- Smaller choices are stated in each section: the stimulus is the AM signal in dB from
stimtodB.mwith the intensity ofcell.dat; f_n of Eq. 15 sums 1BTAs after subtracting the mean stimulus; peaks of C_s are counted with a prominence threshold.
What was found about the dataset:
- Data1 holds 150 sessions (75 experiments of two stimuli each) of 45 cells; 116 sessions from
37 cells on 23 dates use the paper's seven protocols, the paper's 37 cells; 28 sessions use the
gaussModalityandlmdstimuli and 6 use cutoffs of 50 and 1600 Hz, which the paper does not mention. Nothing marks the 86 sessions the paper kept. - Data1 records no animal, no threshold and no best frequency;
cell.datgives the intensity played. 28 sessions were sorted off-line and have no header; in 25 of them the trial blocks run on into the pause (spikes up to 2.76 s, at a median 2.7 Hz from 1 to 1.7 s), while no on-line block goes past 1.02 s. Each block is one trial: none is sorted out of order or holds twice the usual spike count.firing_ratesgives each spike file's trial count, which matches the blocks parsed here in 136 of 138 files, and a rate whose +/- value (median 128 Hz) is not the trial-to-trial SD (median 5.8 Hz). - The intensity in three headers differs from
cell.datby more than 1 dB (up to 6.6 dB); the descriptor says 98 to 533 trials, the files hold 69 to 503; six sessions had pauses of 450 ms, not 700. Two experiments (03-01-15-aa, 03-06-27-ae) are missing fromsinglecells.txt. - Spikes and stimuli share the trial clock: the STA peaks 5.6 ms (median) before the spike, at 0 ms in one session (03-03-06-aa sig2), exceeds the STA with the stimulus shifted by 100 ms in all 116 sessions, and almost vanishes with the stimulus shifted by 100 ms.