Buccino et al. (2024) on dandiset 000294¶
Paper. Buccino AP, Damart T, Bartram J, Mandge D, Xue X, Zbili M, Gänswein T, Jaquier A,
Emmenegger V, Markram H, Hierlemann A, Van Geit W (2024). A multimodal fitting approach to construct
single-neuron models with patch-clamp and high-density microelectrode arrays. Neural Computation
36(7):1286-1331. doi:10.1162/neco_a_01672. The claims file
(claims/000294/10.1162-neco_a_01672.yaml) was written from the version of record (CC BY 4.0), main
text, Methods, Tables 1 to 5 and the appendix; the paper has no separate supplement.
Dandiset. DANDI:000294, version 0.260528.0854 (doi:10.48324/dandi.000294/0.260528.0854). The paper cites the dandiset without a version. This is the only published version, and the most recent; the draft has the same two files (same sizes and blobs), so pinning it loses nothing.
Outside inputs. The paper's code repository,
alejoe91/multimodalfitting at commit
24e8e6168b2d5be45c551bcc08b62fde187d01c7, which the paper names for its scripts, model files,
morphologies and extracted features. It supplies the ground-truth model and its template (claims on
data elsewhere) and the reconstructed morphologies. It is fetched by gtmodel.py; nothing from it
enters the analyses of the dandiset's recordings except method parameters, cited where used.
Claims.
| Claim | Data | Tested here |
|---|---|---|
| recordings-two-cells-runs | dandiset | tested |
| ecode-protocol-structure | dandiset | tested |
| holding-potential | dandiset | tested |
| extracellular-spikes-above-50uv | dandiset | tested |
| mea-electrodes-and-rates | dandiset | tested |
| culture-age | dandiset | blocked: the files and the archive give no culture age, plating date or DIV |
| fig6-ais-from-dendrite | elsewhere (repository morphologies) | tested |
| gt-rheobase | elsewhere (repository model) | tested |
| fig3b-ap-initiates-distal-ais | elsewhere (repository model) | tested |
| fig3a-largest-channel-at-ais | elsewhere (repository template) | tested |
| table1-feature-counts | elsewhere (repository feature files and evaluators) | tested |
| fig4a-strategy-balance, gt-intracellular-scores, fig4c-hyperdepol-only-single, fig4b-intra-no-difference, fig4b-extra-better-than-soma, fig4d-templates-better-than-soma, fig5-ap-distribution | elsewhere | blocked: too large (10 seeds of 600 generations per strategy; estimate below) |
| fig78a-soma-intra-all-extra, experimental-intracellular-scores | with inputs | blocked: too large (estimate below) |
| fig78b-validation-spiking, fig78c-cell2-worse-extracellular, fig10-cell2-soma-underestimated | with inputs | blocked: depend on fig78a-soma-intra-all-extra |
import hashlib
import time
import warnings
from collections import Counter
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import multimodal as mm
pd.set_option("display.width", 200)
pd.set_option("display.max_columns", 40)
pd.set_option("display.max_rows", 300)
warnings.filterwarnings("ignore", category=RuntimeWarning)
started = time.time()
urls = mm.asset_urls()
meta = {s: mm.metadata(u) for s, u in urls.items()}
runs = {s: mm.runs_in_file(u) for s, u in urls.items()}
patch = mm.load_patch(urls, runs)
print(f"dandiset {mm.DANDISET_ID} version {mm.DANDISET_VERSION}: {len(urls)} files; "
f"patch data of {sum(len(v) for v in patch.values())} sweeps loaded in {time.time() - started:.0f} s")
dandiset 000294 version 0.260528.0854: 2 files; patch data of 540 sweeps loaded in 61 s
Data check¶
Each check of the description against the values is stated under its output.
for s, m in meta.items():
print(f"sub-{s} ({mm.CELLS[s]}): identifier {m['identifier']}, subject_id {m['subject_id']}, "
f"species {m['species']}, '{m['subject_description']}'")
print(f" session start {m['session_start_time']}; '{m['session_description']}'")
print(f" subject fields in the file: {m['subject_fields']}")
print(f" runs in the series names: {runs[s]}; MEA series descriptions: "
f"{[v['description'] for v in m['mea'].values()]}")
print(f" confocal maximum projection {m['confocal_shape']}")
print("\nThe paper's code names cell 1 cell1_211006_3148 and cell 2 cell1_211011_3436 (claims file).")
sub-3148 (cell 1): identifier cell1_211006_3148, subject_id 3148, species Rattus norvegicus, 'Rat cortical embryonic cell culture' session start 2021-10-06T11:24:06+02:00; 'Simultaneous patch-clamp/HD-MEA recording using E-CODE protocols for cell cell1_211006_3148 for 5 experimental runs.' subject fields in the file: ['description', 'species', 'subject_id'] runs in the series names: ['run1', 'run2', 'run3', 'run4', 'run5']; MEA series descriptions: ['HD-MEA extracellular recording for run1', 'HD-MEA extracellular recording for run2', 'HD-MEA extracellular recording for run3', 'HD-MEA extracellular recording for run4', 'HD-MEA extracellular recording for run5'] confocal maximum projection (7745, 11396) sub-3436 (cell 2): identifier cell1_211011_3436, subject_id 3436, species Rattus norvegicus, 'Rat cortical embryonic cell culture' session start 2021-10-11T11:58:48+02:00; 'Simultaneous patch-clamp/HD-MEA recording using E-CODE protocols for cell cell1_211011_3436 for 4 experimental runs.' subject fields in the file: ['description', 'species', 'subject_id'] runs in the series names: ['run1', 'run2', 'run3', 'run4']; MEA series descriptions: ['HD-MEA extracellular recording for run3', 'HD-MEA extracellular recording for run4', 'HD-MEA extracellular recording for run5', 'HD-MEA extracellular recording for run6'] confocal maximum projection (7705, 13234) The paper's code names cell 1 cell1_211006_3148 and cell 2 cell1_211011_3436 (claims file).
import json
import urllib.request
from dandi.dandiapi import DandiAPIClient
with DandiAPIClient() as client:
raw = client.get_dandiset(mm.DANDISET_ID, mm.DANDISET_VERSION).get_raw_metadata()
print("archive summary:", {k: raw["assetsSummary"].get(k) for k in ("numberOfFiles", "numberOfSubjects", "variableMeasured")})
for r in raw.get("relatedResource", []):
print("related resource:", r.get("relation"), r.get("url"), "|", r.get("name"))
if "doi.org/" in r.get("url", ""):
doi = r["url"].split("doi.org/")[1]
title = json.load(urllib.request.urlopen(f"https://api.crossref.org/works/{doi}", timeout=60))["message"]["title"][0]
print(f" Crossref title of {doi}: {title}")
print("paper: doi:10.1162/neco_a_01672")
archive summary: {'numberOfFiles': 2, 'numberOfSubjects': 2, 'variableMeasured': ['ElectricalSeries', 'ProcessingModule', 'CurrentClampSeries', 'CurrentClampStimulusSeries', 'ElectrodeGroup']}
related resource: dcite:IsDescribedBy https://doi.org/10.1162/neco_a_01676 | A Multi-Modal Fitting Approach to Construct Single-Neuron Models with Patch-Clamp and High-Density Microelectrode Arrays
Crossref title of 10.1162/neco_a_01676: Energy Complexity of Convolutional Neural Networks related resource: dcite:IsDocumentedBy https://github.com/alejoe91/multimodalfitting/blob/main/notebooks/experimental/0-export-to-nwb.ipynb | NWB conversion notebook paper: doi:10.1162/neco_a_01672
The archive's metadata name the paper by title but link it to doi:10.1162/neco_a_01676, which Crossref gives as another article; the paper is doi:10.1162/neco_a_01672.
Each file holds one neuron; the identifiers are the code's cell names (cell 1 is sub-3148, cell 2 is
sub-3436). Cell 2's runs are named run1 to run4, while the descriptions of its MEA series say run3 to
run6: the export renamed the original runs 3 to 6, as the code's run_map does. The subject has no
age, date of birth or culture date, so the days in vitro cannot be read from the files.
rows = []
for s, sweeps in patch.items():
for d in sweeps:
dt = np.diff(d["t"])
rows.append(dict(cell=mm.CELLS[s], run=d["run"], protocol=d["protocol"], sweep=d["sweep"], n=len(d["t"]),
dt_us=round(float(np.median(dt)) * 1e6, 3), dt_spread_us=float(np.ptp(dt)) * 1e6,
same_t=d["same_t"], v_min=float(d["v"].min()), v_max=float(d["v"].max()),
i_min=float(d["i"].min()), i_max=float(d["i"].max()),
t0=float(d["t"][0]), t1=float(d["t"][-1])))
sw = pd.DataFrame(rows)
print("Sweeps per cell, run and protocol:")
print(sw.pivot_table(index=["cell", "run"], columns="protocol", values="sweep", aggfunc="count").to_string())
print("\nSampling interval of the patch timestamps (us) by protocol (median; spread within a sweep):")
print(sw.groupby(["cell", "protocol"]).agg(dt_us=("dt_us", "median"), dt_min=("dt_us", "min"), dt_max=("dt_us", "max"),
spread_max=("dt_spread_us", "max")).to_string())
print("\nStimulus and response share their timestamps in every sweep:", bool(sw.same_t.all()))
print("\nRanges by cell:")
print(sw.groupby("cell").agg(v_min=("v_min", "min"), v_max=("v_max", "max"), i_min=("i_min", "min"),
i_max=("i_max", "max")).round(1).to_string())
for s, m in meta.items():
print(f"{mm.CELLS[s]}: response unit '{m['patch_attrs']['data_unit']}', conversion {m['patch_attrs']['data_conversion']}, "
f"gain {m['patch_attrs']['gain']}; stimulus unit '{m['stim_attrs']['data_unit']}', conversion {m['stim_attrs']['data_conversion']}")
Sweeps per cell, run and protocol:
protocol apwaveform firepattern hyperdepol idrest idthres iv poscheops sahp
cell run
cell 1 run1 6 2 4 11 21 11 1 4
run2 6 2 4 11 21 11 1 4
run3 6 2 4 11 21 11 1 4
run4 6 2 4 11 21 11 1 4
run5 6 2 4 11 21 11 1 4
cell 2 run1 6 2 4 11 21 11 1 4
run2 6 2 4 11 21 11 1 4
run3 6 2 4 11 21 11 1 4
run4 6 2 4 11 21 11 1 4
Sampling interval of the patch timestamps (us) by protocol (median; spread within a sweep):
dt_us dt_min dt_max spread_max
cell protocol
cell 1 apwaveform 8.0 8.0 8.0 0.000000e+00
firepattern 49.0 49.0 49.0 0.000000e+00
hyperdepol 49.0 49.0 49.0 0.000000e+00
idrest 49.0 49.0 49.0 5.684342e-08
idthres 50.0 50.0 50.0 1.421085e-08
iv 49.0 49.0 49.0 0.000000e+00
poscheops 49.0 49.0 49.0 0.000000e+00
sahp 49.0 49.0 49.0 0.000000e+00
cell 2 apwaveform 8.0 8.0 8.0 0.000000e+00
firepattern 49.0 49.0 49.0 0.000000e+00
hyperdepol 49.0 49.0 49.0 0.000000e+00
idrest 49.0 49.0 49.0 0.000000e+00
idthres 50.0 50.0 50.0 8.881784e-10
iv 49.0 49.0 49.0 7.105427e-09
poscheops 49.0 49.0 49.0 0.000000e+00
sahp 49.0 49.0 49.0 0.000000e+00
Stimulus and response share their timestamps in every sweep: True
Ranges by cell:
v_min v_max i_min i_max
cell
cell 1 -90.3 46.3 -183.0 435.8
cell 2 -100.3 42.0 -520.9 1019.7
cell 1: response unit 'volts', conversion 1.0, gain 0.001; stimulus unit 'amperes', conversion 1.0
cell 2: response unit 'volts', conversion 1.0, gain 0.001; stimulus unit 'amperes', conversion 1.0
The responses are stated in volts with conversion 1, but range from about -100 to +46: they are in millivolts. The stimuli are stated in amperes, but the holding levels and steps (tens to hundreds) are picoamperes: the IDthres steps of cell 1 start at 60 above holding, which is 50% of the rheobase scale of a cortical neuron in pA, not in nA or A. The analyses use mV and pA. The timestamps are regular within every sweep. Their interval is 50 us for IDthres (20 kHz), 49 us for the other protocols (20.41 kHz) and 8 us for APWaveform (125 kHz). The stimulus step durations measured on these timestamps are the nominal ones (3600 ms, 1350 ms, 270 ms, below), so the timestamps are consistent with the protocols, and the intervals are taken as the real sampling intervals.
for s, m in meta.items():
e = m["electrodes"]
print(f"{mm.CELLS[s]}: electrodes table {len(e['id'])} rows; columns {sorted(e)}")
print(f" x, y, z all NaN: {np.isnan(e['x']).all()}, {np.isnan(e['y']).all()}, {np.isnan(e['z']).all()}; "
f"rel_x {e['rel_x'].min()} to {e['rel_x'].max()} um, rel_y {e['rel_y'].min()} to {e['rel_y'].max()} um")
pitch = np.unique(np.round(np.diff(np.unique(e["rel_x"])), 2))
print(f" smallest spacing of rel_x values {pitch.min()} um; gain_to_uV {np.unique(e['gain_to_uV'])}; "
f"electrode numbers {int(e['electrode'].min())} to {int(e['electrode'].max())} of 26,400")
for name, v in m["mea"].items():
print(f" {name}: {v['shape']} {v['dtype']} chunks {v['chunks']}, conversion {v['conversion']:.4g} {v['unit']}, "
f"{v['n_electrodes']} electrodes, t {v['t0']:.2f} to {v['t1']:.2f} s, rate {v['rate']:.2f} Hz")
cell 1: electrodes table 678 rows; columns ['channel_name', 'electrode', 'filtering', 'gain_to_uV', 'group_name', 'id', 'imp', 'location', 'offset_to_uV', 'rel_x', 'rel_y', 'x', 'y', 'z'] x, y, z all NaN: True, True, True; rel_x 1907.5 to 2432.5 um, rel_y 1242.5 to 1767.5 um smallest spacing of rel_x values 17.5 um; gain_to_uV [3.1471252]; electrode numbers 15733 to 22354 of 26,400 ElectricalSeries_run1: (3915033, 678) uint16 chunks (51966, 9), conversion 3.147e-06 volts, 678 electrodes, t 0.00 to 195.75 s, rate 20000.00 Hz ElectricalSeries_run2: (3915352, 678) uint16 chunks (51966, 9), conversion 3.147e-06 volts, 678 electrodes, t 257.00 to 452.77 s, rate 20000.00 Hz ElectricalSeries_run3: (3915455, 678) uint16 chunks (51975, 9), conversion 3.147e-06 volts, 678 electrodes, t 526.00 to 721.77 s, rate 20000.00 Hz ElectricalSeries_run4: (3915409, 678) uint16 chunks (51966, 9), conversion 3.147e-06 volts, 678 electrodes, t 828.00 to 1023.77 s, rate 20000.00 Hz ElectricalSeries_run5: (3915173, 678) uint16 chunks (51966, 9), conversion 3.147e-06 volts, 678 electrodes, t 1061.00 to 1256.76 s, rate 20000.00 Hz cell 2: electrodes table 675 rows; columns ['channel_name', 'electrode', 'filtering', 'gain_to_uV', 'group_name', 'id', 'imp', 'location', 'offset_to_uV', 'rel_x', 'rel_y', 'x', 'y', 'z'] x, y, z all NaN: True, True, True; rel_x 297.5 to 822.5 um, rel_y 1032.5 to 1557.5 um smallest spacing of rel_x values 17.5 um; gain_to_uV [3.1471252]; electrode numbers 13001 to 19622 of 26,400 ElectricalSeries_run1: (3916458, 675) uint16 chunks (52218, 9), conversion 3.147e-06 volts, 675 electrodes, t 0.00 to 195.82 s, rate 20000.00 Hz ElectricalSeries_run2: (3914532, 675) uint16 chunks (52191, 9), conversion 3.147e-06 volts, 675 electrodes, t 240.00 to 435.73 s, rate 20000.00 Hz ElectricalSeries_run3: (3915387, 675) uint16 chunks (52200, 9), conversion 3.147e-06 volts, 675 electrodes, t 489.00 to 684.77 s, rate 20000.00 Hz ElectricalSeries_run4: (3915454, 675) uint16 chunks (52200, 9), conversion 3.147e-06 volts, 675 electrodes, t 779.00 to 974.77 s, rate 20000.00 Hz
The MEA series are 10-bit counts (values up to 1023) with a conversion of 3.147 uV per count, the
gain_to_uV of the electrodes table, at 20 kHz. Electrode positions are in rel_x and rel_y on a
17.5 um grid (Methods 5.6.2); x, y and z are empty.
# Patch and MEA of a run share one clock: patch sweeps fall inside the run's MEA series.
for s, m in meta.items():
for r in runs[s]:
sub = sw[(sw.cell == mm.CELLS[s]) & (sw.run == r)]
mea = m["mea"][f"ElectricalSeries_{r}"]
print(f"{mm.CELLS[s]} {r}: patch {sub.t0.min():.2f} to {sub.t1.max():.2f} s; MEA {mea['t0']:.2f} to {mea['t1']:.2f} s")
# The icephys grouping tables
for s, m in meta.items():
idx = np.r_[0, m["seq_index"]]
members = set(m["seq_members"].tolist())
missing = sorted(set(range(len(m["ic_protocol_name"]))) - members)
miss_prot = Counter(m["ic_protocol_name"][missing])
print(f"\n{mm.CELLS[s]}: {len(m['ic_protocol_name'])} intracellular recordings; {len(m['seq_type'])} sequential recordings "
f"(stimulus types {dict(Counter(m['seq_type']))}); {m['n_repetitions']} repetitions; "
f"{len(missing)} recordings in no sequential recording, by protocol {dict(miss_prot)}")
print(" table columns protocol_name and run agree with the series names:",
bool((sw[sw.cell == mm.CELLS[s]].sort_values("sweep").protocol.values == m["ic_protocol_name"]).all()))
# Recordings deposited twice: checksums of every patch sweep across both files
hashes = Counter()
for s, sweeps in patch.items():
for d in sweeps:
hashes[hashlib.sha1(d["v"].tobytes()).hexdigest()] += 1
print(f"\nPatch sweeps: {sum(hashes.values())}, distinct voltage traces: {len(hashes)}")
cell 1 run1: patch 0.00 to 195.75 s; MEA 0.00 to 195.75 s
cell 1 run2: patch 257.00 to 452.77 s; MEA 257.00 to 452.77 s
cell 1 run3: patch 526.00 to 721.77 s; MEA 526.00 to 721.77 s
cell 1 run4: patch 828.00 to 1023.77 s; MEA 828.00 to 1023.77 s
cell 1 run5: patch 1061.00 to 1256.76 s; MEA 1061.00 to 1256.76 s
cell 2 run1: patch 0.00 to 195.82 s; MEA 0.00 to 195.82 s
cell 2 run2: patch 240.00 to 435.73 s; MEA 240.00 to 435.73 s
cell 2 run3: patch 489.00 to 684.77 s; MEA 489.00 to 684.77 s
cell 2 run4: patch 779.00 to 974.77 s; MEA 779.00 to 974.77 s
cell 1: 300 intracellular recordings; 35 sequential recordings (stimulus types {np.str_('idthres'): 5, np.str_('firepattern'): 5, np.str_('iv'): 5, np.str_('idrest'): 5, np.str_('apwaveform'): 5, np.str_('hyperdepol'): 5, np.str_('sahp'): 5}); 5 repetitions; 35 recordings in no sequential recording, by protocol {np.str_('firepattern'): 5, np.str_('iv'): 5, np.str_('idrest'): 5, np.str_('apwaveform'): 5, np.str_('hyperdepol'): 5, np.str_('sahp'): 5, np.str_('poscheops'): 5}
table columns protocol_name and run agree with the series names: True
cell 2: 240 intracellular recordings; 28 sequential recordings (stimulus types {np.str_('idthres'): 4, np.str_('firepattern'): 4, np.str_('iv'): 4, np.str_('idrest'): 4, np.str_('apwaveform'): 4, np.str_('hyperdepol'): 4, np.str_('sahp'): 4}); 4 repetitions; 28 recordings in no sequential recording, by protocol {np.str_('firepattern'): 4, np.str_('iv'): 4, np.str_('idrest'): 4, np.str_('apwaveform'): 4, np.str_('hyperdepol'): 4, np.str_('sahp'): 4, np.str_('poscheops'): 4}
table columns protocol_name and run agree with the series names: True
Patch sweeps: 540, distinct voltage traces: 540
Within each run the patch sweeps fall inside the MEA series, on one clock that starts at the first
run. The grouping tables are incomplete: in each run the first sweep of every protocol after the
first, and the PosCheops sweep, belong to no sequential recording (the export loop of
0-export-to-nwb.ipynb starts a new group without its first sweep). The protocol_name and run
columns are complete and agree with the series names, so the analyses use these. No voltage trace
is deposited twice; the MEA runs are compared below by checksums of their raw segments.
# Denoising of every run (code's remove_blank_channels), with checksums of the raw segments.
t_mea = time.time()
blank, mea_hash = {}, {}
for s, m in meta.items():
for r in runs[s]:
name = f"ElectricalSeries_{r}"
blank[(s, r)], mea_hash[(s, r)] = mm.blank_channels(urls[s], name, m["mea"][name]["shape"][0])
dn = pd.DataFrame([dict(cell=mm.CELLS[s], run=r, channels=len(b), blank=int(b.sum()), kept=int((~b).sum()))
for (s, r), b in blank.items()])
print(dn.to_string(index=False))
allh = [h for v in mea_hash.values() for h in v]
print(f"\nRaw MEA segments compared: {len(allh)}, distinct: {len(set(allh))} ({time.time() - t_mea:.0f} s)")
cell run channels blank kept cell 1 run1 678 118 560 cell 1 run2 678 126 552 cell 1 run3 678 118 560 cell 1 run4 678 117 561 cell 1 run5 678 117 561 cell 2 run1 675 130 545 cell 2 run2 675 130 545 cell 2 run3 675 132 543 cell 2 run4 675 130 545 Raw MEA segments compared: 90, distinct: 90 (58 s)
Every run of cell 1 has 678 channels and cell 2 675, but many are constant (blank) in some second of the run; these are the channels the code removes before computing templates. No raw MEA segment is repeated between runs or files.
# Clock check: spike-triggered MEA waveforms at the patch spike times, and at times shifted by 100 ms.
def run_spikes(s, r, include_poscheops=False):
out = []
for d in patch[s]:
if d["run"] == r and (include_poscheops or d["protocol"] != "poscheops"):
out.append(mm.peak_times(d["t"], d["v"]))
return np.concatenate(out)
mea_t = {}
def mea_timestamps(s, r):
if (s, r) not in mea_t:
mea_t[(s, r)] = mm.open_file(urls[s])["acquisition"][f"ElectricalSeries_{r}"]["timestamps"][:]
return mea_t[(s, r)]
clock_rows = []
for s in urls:
r = runs[s][0]
name = f"ElectricalSeries_{r}"
info = meta[s]["mea"][name]
sp = run_spikes(s, r)
for shift in (0.0, 0.1, -0.1):
wf, used, wide = mm.mea_waveforms(urls[s], name, mea_timestamps(s, r), sp + shift, info["conversion"] * 1e6,
info["chunks"][0])
T = mm.template_from_waveforms(wf, wide, keep=~blank[(s, r)])
a = mm.amplitudes(T["median"], T["max_channel"])
clock_rows.append(dict(cell=mm.CELLS[s], run=r, shift_ms=shift * 1e3, spikes=len(used),
max_channel_neg_peak_uV=round(a["neg_peak"], 1), max_channel_ptp_uV=round(a["ptp"], 1)))
print(pd.DataFrame(clock_rows).to_string(index=False))
cell run shift_ms spikes max_channel_neg_peak_uV max_channel_ptp_uV cell 1 run1 0.0 190 189.2 244.1 cell 1 run1 100.0 190 23.3 27.4 cell 1 run1 -100.0 190 23.4 27.3 cell 2 run1 0.0 217 159.2 192.2 cell 2 run1 100.0 217 19.5 22.7 cell 2 run1 -100.0 217 23.5 26.7
With the patch spike times, the template has a spike of 189.2 uV (cell 1) and 159.2 uV (cell 2) on its largest channel; with the same times shifted by 100 ms in either direction, the largest value over more than 500 channels is 19.5 to 23.5 uV, the noise left after averaging. So the patch and MEA series of a run share their clock, at least to well within the 2 ms search window of the centring. This is shown for the first run of each cell; the other runs are checked the same way by their templates below, and the alignment of single spikes within long sweeps is shown there.
recordings-two-cells-runs¶
The experimental data are simultaneous patch-clamp and HD-MEA recordings of two cultured neurons, each stimulated with the eCode protocols repeated in four to six runs of about 200 s.
Criterion: two neurons, each with 4 to 6 runs, and the median run duration (first to last timestamp of a run's patch series) within 200 +/- 20 s.
dur = sw.groupby(["cell", "run"]).agg(first=("t0", "min"), last=("t1", "max"))
dur["duration_s"] = dur["last"] - dur["first"]
dur["mea_s"] = [meta[s]["mea"][f"ElectricalSeries_{r}"]["t1"] - meta[s]["mea"][f"ElectricalSeries_{r}"]["t0"]
for s in urls for r in runs[s]]
print(dur.round(2).to_string())
n_runs = dur.groupby("cell").size()
med = float(dur.duration_s.median())
print(f"\nneurons (NWB files): {len(urls)}; runs per neuron: {n_runs.to_dict()}")
print(f"median run duration (patch series): {med:.1f} s (range {dur.duration_s.min():.1f} to {dur.duration_s.max():.1f}); "
f"MEA series {dur.mea_s.min():.1f} to {dur.mea_s.max():.1f} s")
print("runs per cell within 4 to 6:", bool(((n_runs >= 4) & (n_runs <= 6)).all()), "; duration within 180 to 220 s:", 180 <= med <= 220)
first last duration_s mea_s
cell run
cell 1 run1 0.0 195.75 195.75 195.75
run2 257.0 452.77 195.77 195.77
run3 526.0 721.77 195.77 195.77
run4 828.0 1023.77 195.77 195.77
run5 1061.0 1256.76 195.76 195.76
cell 2 run1 0.0 195.82 195.82 195.82
run2 240.0 435.73 195.73 195.73
run3 489.0 684.77 195.77 195.77
run4 779.0 974.77 195.77 195.77
neurons (NWB files): 2; runs per neuron: {'cell 1': 5, 'cell 2': 4}
median run duration (patch series): 195.8 s (range 195.7 to 195.8); MEA series 195.7 to 195.8 s
runs per cell within 4 to 6: True ; duration within 180 to 220 s: True
Result: reproduced. The dandiset holds two neurons, cell 1 with 5 runs and cell 2 with 4 (the runs the paper's code used: runs 1 to 5 of cell 1 and 3 to 6 of cell 2, renamed). The median run lasts 195.8 s from the first to the last patch sample (195.7 to 195.8 s, the same as the MEA series), within 200 +/- 20 s.
ecode-protocol-structure¶
Each run applies the eCode protocols with the sweep numbers, step durations and amplitudes (in percent of the estimated rheobase) that the Methods list.
Criterion: in every run, all eight protocols with the Methods' sweep counts and durations (within 10%) and every step amplitude within 10% of its stated percentage of the IDthres scale (sweep k of IDthres at 50 + 4k %); IV under reading (a) 11 steps from -140 to 60% or (b) 9 steps from -140 to 20%. Amplitudes are measured above the holding level. For a stated 0% (IV, reading a) a relative tolerance is not defined; 10% of the IV increment (2% of the scale) is used (deviation).
EXPECTED = { # Methods 5.2: (sweeps, [(duration ms, % per sweep), ...] per phase)
"idthres": (21, [(270, [50 + 4 * k for k in range(21)])]),
"firepattern": (2, [(3600, [120, 200])]),
"iv_a": (11, [(3000, list(range(-140, 61, 20)))]),
"iv_b": (9, [(3000, list(range(-140, 21, 20)))]),
"idrest": (11, [(1350, list(range(50, 301, 25)))]),
"apwaveform": (6, [(50, list(range(200, 351, 30)))]),
"hyperdepol": (4, [(450, [-40, -80, -120, -160]), (270, [100] * 4)]),
"sahp": (4, [(250, [40] * 4), (225, [150, 200, 250, 300]), (450, [40] * 4)]),
}
def scale_of(s, r):
ids = [d for d in patch[s] if d["run"] == r and d["protocol"] == "idthres"]
amps = [mm.segments(d)[1][2] for d in ids]
return float(np.median([a / ((50 + 4 * k) / 100) for k, a in enumerate(amps)])), amps
def ok_amp(meas, exp):
return abs(meas - exp) <= (0.1 * abs(exp) if exp != 0 else 2.0)
rows, detail = [], []
for s in urls:
for r in runs[s]:
scale, _ = scale_of(s, r)
by = {}
for d in patch[s]:
if d["run"] == r:
by.setdefault(d["protocol"], []).append(d)
for prot, key_list in [("idthres", ["idthres"]), ("firepattern", ["firepattern"]), ("iv", ["iv_a", "iv_b"]),
("idrest", ["idrest"]), ("apwaveform", ["apwaveform"]), ("hyperdepol", ["hyperdepol"]),
("sahp", ["sahp"])]:
sweeps = sorted(by.get(prot, []), key=lambda d: d["sweep"])
phases = [[seg for seg in mm.segments(d) if abs(seg[2]) > 0.02 * scale] for d in sweeps]
for key in key_list:
n_exp, exp_phases = EXPECTED[key]
count_ok = len(sweeps) == n_exp
dur_ok = amp_ok = count_ok
if count_ok:
for k, ph in enumerate(phases):
if len(ph) != len(exp_phases):
# a 0% step has no phase above holding
if not (len(ph) == 0 and all(e[1][k] == 0 for e in exp_phases)):
dur_ok = amp_ok = False
continue
for (start, dms, lev), (edur, eperc) in zip(ph, exp_phases):
dur_ok &= abs(dms - edur) <= 0.1 * edur
amp_ok &= ok_amp(100 * lev / scale, eperc[k])
detail.append(dict(cell=mm.CELLS[s], run=r, protocol=key, sweep=k, duration_ms=round(dms, 1),
expected_ms=edur, percent=round(100 * lev / scale, 1), expected_percent=eperc[k]))
rows.append(dict(cell=mm.CELLS[s], run=r, protocol=key, scale_pA=round(scale, 1), sweeps=len(sweeps),
expected_sweeps=n_exp, count_ok=count_ok, durations_ok=bool(dur_ok), amplitudes_ok=bool(amp_ok)))
# PosCheops: three ramps to 300% with rise and fall of 4, 2, 1.33 s
pc = by.get("poscheops", [])
rp = mm.ramps(pc[0]) if len(pc) == 1 else []
exp_r = [4000, 2000, 1333]
r_ok = len(rp) == 3 and all(abs(a[1] - e) <= 0.1 * e and abs(a[2] - e) <= 0.1 * e for a, e in zip(rp, exp_r))
a_ok = len(rp) == 3 and all(ok_amp(100 * a[3] / scale, 300) for a in rp)
rows.append(dict(cell=mm.CELLS[s], run=r, protocol="poscheops", scale_pA=round(scale, 1), sweeps=len(pc),
expected_sweeps=1, count_ok=len(pc) == 1, durations_ok=bool(r_ok), amplitudes_ok=bool(a_ok)))
for k, a in enumerate(rp):
detail.append(dict(cell=mm.CELLS[s], run=r, protocol="poscheops", sweep=k, duration_ms=round(a[1], 1),
expected_ms=exp_r[k], percent=round(100 * a[3] / scale, 1), expected_percent=300,
fall_ms=round(a[2], 1)))
ec = pd.DataFrame(rows)
det = pd.DataFrame(detail)
print("IDthres scale (100%) per run, pA:")
print(ec.drop_duplicates(["cell", "run"])[["cell", "run", "scale_pA"]].to_string(index=False))
print("\nPer protocol, every run:")
print(ec.groupby(["cell", "protocol"])[["count_ok", "durations_ok", "amplitudes_ok"]].all().to_string())
print("\nMeasured durations and percentages, range over runs and sweeps:")
print(det.groupby(["protocol", "expected_ms"]).agg(dur_min=("duration_ms", "min"), dur_max=("duration_ms", "max")).to_string())
det["err_percent"] = det.percent - det.expected_percent
print("\nLargest deviation from the stated percentage, by protocol (percentage points):")
print(det.groupby(["cell", "protocol"]).err_percent.agg(lambda x: x.abs().max()).round(2).to_string())
IDthres scale (100%) per run, pA:
cell run scale_pA
cell 1 run1 120.4
cell 1 run2 120.5
cell 1 run3 120.5
cell 1 run4 120.7
cell 1 run5 120.8
cell 2 run1 301.7
cell 2 run2 301.7
cell 2 run3 301.2
cell 2 run4 301.3
Per protocol, every run:
count_ok durations_ok amplitudes_ok
cell protocol
cell 1 apwaveform True True True
firepattern True True True
hyperdepol True True True
idrest True True True
idthres True True True
iv_a True True True
iv_b False False False
poscheops True True True
sahp True True True
cell 2 apwaveform True True True
firepattern True True True
hyperdepol True True True
idrest True True True
idthres True True True
iv_a True True True
iv_b False False False
poscheops True True True
sahp True True True
Measured durations and percentages, range over runs and sweeps:
dur_min dur_max
protocol expected_ms
apwaveform 50 49.9 49.9
firepattern 3600 3599.9 3599.9
hyperdepol 270 269.9 269.9
450 450.0 450.0
idrest 1350 1349.9 1350.0
idthres 270 269.9 270.0
iv_a 3000 2999.9 3000.0
poscheops 1333 1307.3 1322.0
2000 1964.5 1987.9
4000 3924.5 3975.1
sahp 225 224.9 225.0
250 249.9 250.0
450 450.0 450.0
Largest deviation from the stated percentage, by protocol (percentage points):
cell protocol
cell 1 apwaveform 1.0
firepattern 0.3
hyperdepol 0.6
idrest 1.2
idthres 0.4
iv_a 0.7
poscheops 0.9
sahp 1.2
cell 2 apwaveform 0.3
firepattern 0.1
hyperdepol 0.3
idrest 0.3
idthres 0.3
iv_a 0.3
poscheops 0.5
sahp 0.3
ok_a = ec[ec.protocol != "iv_b"][["count_ok", "durations_ok", "amplitudes_ok"]].all()
ok_b = ec[ec.protocol != "iv_a"][["count_ok", "durations_ok", "amplitudes_ok"]].all()
print("reading (a), IV -140 to 60%:", ok_a.to_dict())
print("reading (b), IV -140 to 20%:", ok_b.to_dict())
iv_pc = det[det.protocol == "iv_a"].groupby("sweep").percent.median()
print("IV, median percent per sweep over runs and cells:", iv_pc.round(1).tolist())
fig, ax = plt.subplots(1, 2, figsize=(11, 3.5))
for s in urls:
d = [x for x in patch[s] if x["run"] == runs[s][0] and x["protocol"] == "sahp"][-1]
ax[0].plot((d["t"] - d["t"][0]) * 1e3, d["i"] - mm.holding_current(d), label=mm.CELLS[s])
d = [x for x in patch[s] if x["run"] == runs[s][0] and x["protocol"] == "poscheops"][0]
ax[1].plot((d["t"] - d["t"][0]), d["i"] - mm.holding_current(d), label=mm.CELLS[s])
ax[0].set(title="sAHP 300% sweep, first run", xlabel="ms", ylabel="pA above holding")
ax[1].set(title="PosCheops, first run", xlabel="s")
ax[0].legend(); plt.tight_layout(); plt.show()
print("All plotted values are inside the axes (automatic limits).")
reading (a), IV -140 to 60%: {'count_ok': True, 'durations_ok': True, 'amplitudes_ok': True}
reading (b), IV -140 to 20%: {'count_ok': False, 'durations_ok': False, 'amplitudes_ok': False}
IV, median percent per sweep over runs and cells: [-140.1, -120.3, -100.0, -80.2, -60.0, -40.0, -19.9, 19.9, 40.0, 60.0]
All plotted values are inside the axes (automatic limits).
Result: reproduced. In every run of both cells all eight protocols are present with the Methods' sweep counts (21, 2, 11, 11, 6, 4, 4 and 1), step and ramp durations within 10% and amplitudes within 10% of their stated percentages of the IDthres scale (120.4 to 120.8 pA for cell 1, 301.2 to 301.7 pA for cell 2); the largest deviation is 1.2 percentage points. IV runs from -140 to 60% in 11 steps, reading (a); reading (b), 9 steps to 20%, fails on the sweep count. The PosCheops ramps reach 300% and rise in 3925 to 3975, 1965 to 1988 and 1307 to 1322 ms (measured from 5 pA above holding), within 10% of 4, 2 and 1.33 s.
holding-potential¶
The holding current was set so that the read-out membrane potential was about -70 mV, which is about -84 mV after liquid junction potential correction.
Criterion: per sweep, the median membrane potential before the stimulus onset; median over sweeps
per cell; reading (a) uncorrected, -70 +/- 7 mV; (b) corrected, -84 +/- 8.4 mV; reproduced when both
cells match under either reading. The onset is taken from the stimulus itself (the first change of
the injected current). The code's ton (400 ms for IDrest) is 50 ms after the onset in these files,
so it is shown as a variant, though its window includes 50 ms of the step.
hp = []
for s in urls:
for d in patch[s]:
on = mm.onset_ms(d)
tt = (d["t"] - d["t"][0]) * 1e3
pre = d["v"][tt < on]
pre_code = d["v"][tt < mm.CODE_TON_MS[d["protocol"]]]
hp.append(dict(cell=mm.CELLS[s], run=d["run"], protocol=d["protocol"], onset_ms=round(on, 1),
v_pre=float(np.median(pre)), v_pre_code_ton=float(np.median(pre_code)),
i_hold=mm.holding_current(d)))
hp = pd.DataFrame(hp)
print("Stimulus onset by protocol (ms from sweep start), code's ton:")
print(hp.groupby("protocol").onset_ms.agg(["min", "max"]).assign(code_ton=lambda x: x.index.map(mm.CODE_TON_MS)).to_string())
summ = hp.groupby("cell").agg(v_pre_median=("v_pre", "median"), v_pre_min=("v_pre", "min"), v_pre_max=("v_pre", "max"),
v_code_ton=("v_pre_code_ton", "median"), i_hold_pA=("i_hold", "median"), sweeps=("v_pre", "size"))
print("\nPer cell (mV, pA):")
print(summ.round(2).to_string())
print("\nBy run:")
print(hp.groupby(["cell", "run"]).v_pre.median().round(2).to_string())
for cell, v in summ.v_pre_median.items():
print(f"{cell}: {v:.2f} mV; reading (a) |v + 70| <= 7: {abs(v + 70) <= 7}; reading (b) |v + 84| <= 8.4: {abs(v + 84) <= 8.4}")
Stimulus onset by protocol (ms from sweep start), code's ton:
min max code_ton
protocol
apwaveform 298.2 300.1 350
firepattern 450.0 450.1 500
hyperdepol 350.0 350.0 400
idrest 350.0 350.1 400
idthres 349.9 350.0 400
iv 350.0 350.0 400
poscheops 974.0 1021.8 1000
sahp 350.0 350.0 400
Per cell (mV, pA):
v_pre_median v_pre_min v_pre_max v_code_ton i_hold_pA sweeps
cell
cell 1 -69.54 -76.38 -63.70 -69.48 12.24 300
cell 2 -69.12 -73.45 -64.66 -69.04 -36.72 240
By run:
cell run
cell 1 run1 -70.95
run2 -69.92
run3 -68.44
run4 -67.09
run5 -66.59
cell 2 run1 -69.42
run2 -68.37
run3 -68.19
run4 -70.00
cell 1: -69.54 mV; reading (a) |v + 70| <= 7: True; reading (b) |v + 84| <= 8.4: False
cell 2: -69.12 mV; reading (a) |v + 70| <= 7: True; reading (b) |v + 84| <= 8.4: False
Result: reproduced. The median pre-stimulus potential is -69.54 mV for cell 1 and -69.12 mV for
cell 2 (median holding currents 12.24 and -36.72 pA), within -70 +/- 7 mV, reading (a); run medians
range from -71.0 to -66.6 mV. Neither cell is within -84 +/- 8.4 mV, so reading (b) fails; the files do
not say whether the voltages are corrected for the junction potential. The code's ton window gives the same values to within
0.1 mV.
extracellular-spikes-above-50uv¶
Only cells whose extracellular spikes were at least 50 uV were kept, so both recorded cells show extracellular spikes of 50 uV or more.
Criterion: the template of each cell by the template method; reading (a) the absolute negative peak and (b) the peak-to-peak on the largest channel; reproduced when both cells reach 50 uV under either reading. Templates are computed for every run (Methods 5.7.1), the run with the most channels after denoising is used, and the others are shown.
t_tpl = time.time()
templates, trow = {}, []
for s in urls:
for r in runs[s]:
name = f"ElectricalSeries_{r}"
info = meta[s]["mea"][name]
sp = run_spikes(s, r)
wf, used, wide = mm.mea_waveforms(urls[s], name, mea_timestamps(s, r), sp, info["conversion"] * 1e6, info["chunks"][0])
T = mm.template_from_waveforms(wf, wide, keep=~blank[(s, r)])
T["used_t"] = mea_timestamps(s, r)[used]
templates[(s, r)] = T
a, am = mm.amplitudes(T["median"], T["max_channel"]), mm.amplitudes(T["mean"], T["max_channel"])
ptp = np.ptp(T["median"], axis=0)
trow.append(dict(cell=mm.CELLS[s], run=r, patch_spikes=len(sp), waveforms=T["n_waveforms"], kept=T["n_kept"],
channels_after_denoising=len(T["channels"]), channels_ptp_ge_5uV=int((ptp >= 5).sum()),
max_channel_electrode=meta[s]["electrodes"]["channel_name"][T["channels"][T["max_channel"]]],
neg_peak_uV=round(a["neg_peak"], 1), ptp_uV=round(a["ptp"], 1),
neg_peak_mean_uV=round(am["neg_peak"], 1), ptp_mean_uV=round(am["ptp"], 1),
shift_median_ms=float(np.median(T["shifts_ms"]))))
del wf
tp = pd.DataFrame(trow)
print(tp.to_string(index=False))
print(f"({time.time() - t_tpl:.0f} s)")
chosen = tp.loc[tp.groupby("cell").channels_after_denoising.idxmax()]
print("\nRun used per cell (most channels after denoising; first run on ties):")
print(chosen[["cell", "run", "channels_after_denoising", "neg_peak_uV", "ptp_uV", "neg_peak_mean_uV", "ptp_mean_uV"]].to_string(index=False))
for _, c in chosen.iterrows():
print(f"{c.cell}: (a) negative peak {c.neg_peak_uV} uV >= 50: {c.neg_peak_uV >= 50}; (b) peak-to-peak {c.ptp_uV} uV >= 50: {c.ptp_uV >= 50}")
code_chan = {"3148": "ch885", "3436": "ch384"} # the channel 0-export-to-nwb.ipynb plots for each cell
for s in urls:
names = meta[s]["electrodes"]["channel_name"]
k = [i for i, n in enumerate(names) if n.split()[0] == code_chan[s]][0]
r = chosen[chosen.cell == mm.CELLS[s]].run.iloc[0]
T = templates[(s, r)]
if k in T["channels"]:
j = int(np.where(T["channels"] == k)[0][0])
a = mm.amplitudes(T["median"], j)
print(f"{mm.CELLS[s]}: the export notebook's channel {code_chan[s]} ({names[k]}): negative peak {a['neg_peak']:.1f} uV, ptp {a['ptp']:.1f} uV")
else:
print(f"{mm.CELLS[s]}: the export notebook's channel {code_chan[s]} is removed as blank")
print(f"Mean template minus median template, largest over runs: negative peak {(tp.neg_peak_mean_uV - tp.neg_peak_uV).max():.1f} uV, "
f"peak-to-peak {(tp.ptp_mean_uV - tp.ptp_uV).max():.1f} uV")
print(f"\nOver all runs: negative peak {tp.neg_peak_uV.min()} to {tp.neg_peak_uV.max()} uV; peak-to-peak {tp.ptp_uV.min()} to {tp.ptp_uV.max()} uV")
cell run patch_spikes waveforms kept channels_after_denoising channels_ptp_ge_5uV max_channel_electrode neg_peak_uV ptp_uV neg_peak_mean_uV ptp_mean_uV shift_median_ms cell 1 run1 190 190 174 560 287 ch885 elec19044 189.2 244.1 202.5 256.8 -0.25 cell 1 run2 184 184 177 552 283 ch885 elec19044 202.3 261.6 214.8 273.6 -0.25 cell 1 run3 186 186 177 560 274 ch885 elec19044 206.3 264.3 219.9 278.0 -0.30 cell 1 run4 208 208 204 561 261 ch885 elec19044 217.0 274.0 225.5 281.2 -0.45 cell 1 run5 216 216 209 561 256 ch885 elec19044 218.0 276.4 225.5 283.2 -0.55 cell 2 run1 217 217 204 545 267 ch641 elec16311 159.2 192.2 161.3 194.3 -0.30 cell 2 run2 157 157 150 545 290 ch641 elec16311 164.5 200.3 165.5 200.9 -0.25 cell 2 run3 133 133 126 543 353 ch641 elec16311 166.4 201.7 167.2 201.5 -0.25 cell 2 run4 125 125 115 545 368 ch641 elec16311 164.6 197.8 164.7 197.9 -0.30 (157 s) Run used per cell (most channels after denoising; first run on ties): cell run channels_after_denoising neg_peak_uV ptp_uV neg_peak_mean_uV ptp_mean_uV cell 1 run4 561 217.0 274.0 225.5 281.2 cell 2 run1 545 159.2 192.2 161.3 194.3 cell 1: (a) negative peak 217.0 uV >= 50: True; (b) peak-to-peak 274.0 uV >= 50: True cell 2: (a) negative peak 159.2 uV >= 50: True; (b) peak-to-peak 192.2 uV >= 50: True cell 1: the export notebook's channel ch885 (ch885 elec19044): negative peak 217.0 uV, ptp 274.0 uV cell 2: the export notebook's channel ch384 (ch384 elec16312): negative peak 156.8 uV, ptp 190.6 uV Mean template minus median template, largest over runs: negative peak 13.6 uV, peak-to-peak 13.7 uV Over all runs: negative peak 159.2 to 218.0 uV; peak-to-peak 192.2 to 276.4 uV
fig, axs = plt.subplots(1, 2, figsize=(11, 4))
for ax, (_, c) in zip(axs, chosen.iterrows()):
s = [k for k, v in mm.CELLS.items() if v == c.cell][0]
T = templates[(s, c.run)]
tt = np.arange(T["median"].shape[0]) / 20.0 - mm.MS_BEFORE
ch = T["max_channel"]
ax.plot(tt, T["median"][:, ch], "k", label="median (Methods)")
ax.plot(tt, T["mean"][:, ch], "C1--", label="mean (code)")
ax.axhline(-50, color="gray", lw=0.8)
ax.set(title=f"{c.cell} {c.run}, largest channel", xlabel="ms from extracellular trough", ylabel="uV")
ax.legend()
print(f"{c.cell}: plotted range {T['median'][:, ch].min():.1f} to {T['median'][:, ch].max():.1f} uV, all inside the axes")
plt.tight_layout(); plt.show()
cell 1: plotted range -217.0 to 56.9 uV, all inside the axes cell 2: plotted range -159.2 to 33.0 uV, all inside the axes
Patch and MEA clocks within sweeps. The clock check above uses whole runs. Here each spike's centring shift (the offset of the MEA trough on the largest channel from the patch peak, searched within +/- 2 ms) is set against the spike's time in its sweep, for every run. A patch clock that ran at a different rate from the MEA clock would move late spikes (up to 3.6 s into a firepattern sweep) away from the early ones, or onto the +/- 2 ms edge of the search.
al = []
for (s, r), T in templates.items():
sws = [d for d in patch[s] if d["run"] == r and d["protocol"] != "poscheops"]
t0 = np.array([d["t"][0] for d in sws])
t1 = np.array([d["t"][-1] for d in sws])
for tm, sh in zip(T["used_t"], T["shifts_ms"]):
k = np.where((tm >= t0) & (tm <= t1))[0]
if len(k):
al.append(dict(cell=mm.CELLS[s], run=r, protocol=sws[k[0]]["protocol"], t_in_sweep_s=tm - t0[k[0]], shift_ms=sh))
al = pd.DataFrame(al)
al["time_in_sweep"] = pd.cut(al.t_in_sweep_s, [0, 0.5, 1.0, 2.0, 4.5], labels=["0-0.5 s", "0.5-1 s", "1-2 s", "2-4.5 s"])
print(f"spikes placed in a sweep: {len(al)} of {sum(len(T['used_t']) for T in templates.values())}")
print("\nCentring shift (ms) by time in the sweep, all runs:")
print(al.groupby(["cell", "time_in_sweep"], observed=True).shift_ms.describe(percentiles=[0.1, 0.5, 0.9])
[["count", "10%", "50%", "90%", "min", "max"]].round(3).to_string())
edge = (al.shift_ms.abs() >= mm.CENTER_MS - 0.05)
print(f"\nshifts at the +/- {mm.CENTER_MS} ms edge of the search: {int(edge.sum())} of {len(al)}")
slopes = []
for (cell, r), g in al.groupby(["cell", "run"]):
b = np.polyfit(g.t_in_sweep_s, g.shift_ms, 1)[0]
late = g[g.t_in_sweep_s > 1.0].shift_ms
slopes.append(dict(cell=cell, run=r, spikes=len(g), slope_ms_per_s=round(b, 4), median_shift_ms=round(g.shift_ms.median(), 3),
median_shift_after_1s_ms=round(late.median(), 3) if len(late) else np.nan, spikes_after_1s=len(late),
latest_spike_s=round(g.t_in_sweep_s.max(), 2)))
print("\nPer run: least-squares slope of the shift against time in the sweep")
print(pd.DataFrame(slopes).to_string(index=False))
print("A clock running at 20.41 kHz where 20 kHz was meant would drift by 20.4 ms per s of sweep.")
spikes placed in a sweep: 1616 of 1616 Centring shift (ms) by time in the sweep, all runs:
count 10% 50% 90% min max
cell time_in_sweep
cell 1 0-0.5 s 144.0 -0.535 -0.35 -0.250 -0.70 -0.20
0.5-1 s 349.0 -0.550 -0.35 -0.250 -0.60 -0.15
1-2 s 364.0 -0.550 -0.35 -0.250 -0.65 -0.15
2-4.5 s 127.0 -0.550 -0.30 -0.200 -0.60 -0.15
cell 2 0-0.5 s 225.0 -0.400 -0.30 -0.200 -0.45 -0.15
0.5-1 s 234.0 -0.300 -0.30 -0.215 -0.35 -0.15
1-2 s 144.0 -0.350 -0.30 -0.215 -0.35 -0.20
2-4.5 s 29.0 -0.310 -0.25 -0.200 -0.35 -0.20
shifts at the +/- 2.0 ms edge of the search: 0 of 1616
Per run: least-squares slope of the shift against time in the sweep
cell run spikes slope_ms_per_s median_shift_ms median_shift_after_1s_ms spikes_after_1s latest_spike_s
cell 1 run1 190 0.0131 -0.25 -0.25 96 3.99
cell 1 run2 184 0.0201 -0.25 -0.25 93 4.04
cell 1 run3 186 0.0055 -0.30 -0.30 94 4.03
cell 1 run4 208 0.0025 -0.45 -0.45 102 4.03
cell 1 run5 216 -0.0030 -0.55 -0.55 106 3.95
cell 2 run1 217 0.0019 -0.30 -0.30 68 3.83
cell 2 run2 157 0.0085 -0.25 -0.25 41 3.91
cell 2 run3 133 0.0086 -0.25 -0.30 35 4.01
cell 2 run4 125 0.0233 -0.30 -0.30 29 3.74
A clock running at 20.41 kHz where 20 kHz was meant would drift by 20.4 ms per s of sweep.
Result: reproduced. On the run with the most channels after denoising, the largest channel's template reaches a negative peak of 217.0 uV (peak-to-peak 274.0 uV) for cell 1, run 4, and 159.2 uV (192.2 uV) for cell 2, run 1, at least 50 uV under both readings. Every run of both cells is above 50 uV under both readings (negative peaks 159.2 to 218.0 uV), and the code's mean template is higher by up to 13.6 uV (negative peak) and 13.7 uV (peak-to-peak). The channel the paper's export notebook plots for cell 1 (ch885) is the largest channel found here; for cell 2 it plots ch384, an adjacent electrode (156.8 uV). The patch and MEA clocks stay aligned within sweeps: all 1616 spikes have their MEA trough 0.15 to 0.70 ms before the patch peak, none at the edge of the search, with medians of -0.35 to -0.25 ms in every window of the sweep, also 2 to 4.5 s in, and slopes of -0.003 to 0.023 ms per s of sweep, against the 20.4 ms per s that a 20 kHz clock read as 20.41 kHz would give.
mea-electrodes-and-rates¶
The HD-MEA recorded more than 500 electrodes near the pipette simultaneously, out of at most 1024, at 20 kHz, and the patch signal was digitized at 20 kHz.
Criterion: per cell, more than 500 and at most 1024 recorded electrodes, an MEA rate of 20 kHz, and a patch rate of 20 kHz for every protocol (APWaveform at 20 or 200 kHz); expected rates with tolerance 0.
er = []
for s, m in meta.items():
for name, v in m["mea"].items():
er.append(dict(cell=mm.CELLS[s], series=name, electrodes=v["n_electrodes"], mea_rate_Hz=round(v["rate"], 3)))
er = pd.DataFrame(er)
print(er.to_string(index=False))
pr = sw.groupby(["cell", "protocol"]).dt_us.median().rename("dt_us").reset_index()
pr["patch_rate_Hz"] = (1e6 / pr.dt_us).round(1)
print("\nPatch sampling rate from the timestamps:")
print(pr.pivot(index="protocol", columns="cell", values="patch_rate_Hz").to_string())
el_ok = bool(((er.electrodes > 500) & (er.electrodes <= 1024)).all())
mea_ok = bool((er.mea_rate_Hz.round() == 20000).all())
non_ap = pr[pr.protocol != "apwaveform"]
patch_ok = bool((non_ap.patch_rate_Hz.round() == 20000).all())
ap_ok = bool(pr[pr.protocol == "apwaveform"].patch_rate_Hz.round().isin([20000, 200000]).all())
print(f"\nelectrodes in (500, 1024]: {el_ok}; MEA 20 kHz: {mea_ok}; patch 20 kHz for every protocol but APWaveform: {patch_ok}; "
f"APWaveform at 20 or 200 kHz: {ap_ok}")
print("Protocols at 20 kHz:", sorted(non_ap[non_ap.patch_rate_Hz.round() == 20000].protocol.unique()))
cell series electrodes mea_rate_Hz cell 1 ElectricalSeries_run1 678 20000.0 cell 1 ElectricalSeries_run2 678 20000.0 cell 1 ElectricalSeries_run3 678 20000.0 cell 1 ElectricalSeries_run4 678 20000.0 cell 1 ElectricalSeries_run5 678 20000.0 cell 2 ElectricalSeries_run1 675 20000.0 cell 2 ElectricalSeries_run2 675 20000.0 cell 2 ElectricalSeries_run3 675 20000.0 cell 2 ElectricalSeries_run4 675 20000.0 Patch sampling rate from the timestamps: cell cell 1 cell 2 protocol apwaveform 125000.0 125000.0 firepattern 20408.2 20408.2 hyperdepol 20408.2 20408.2 idrest 20408.2 20408.2 idthres 20000.0 20000.0 iv 20408.2 20408.2 poscheops 20408.2 20408.2 sahp 20408.2 20408.2 electrodes in (500, 1024]: True; MEA 20 kHz: True; patch 20 kHz for every protocol but APWaveform: False; APWaveform at 20 or 200 kHz: False Protocols at 20 kHz: ['idthres']
Result: partially reproduced. Each cell's MEA recorded 678 (cell 1) and 675 (cell 2) electrodes at 20 kHz, as stated. The patch timestamps give 20 kHz only for IDthres; the other protocols are at 20.41 kHz (49 us intervals) and APWaveform at 125 kHz (8 us), neither 20 nor 200 kHz. The intervals agree with the nominal stimulus durations, so they are the sampling of the files as deposited.
culture-age¶
The cultured cells were recorded between days in vitro 18 and 24.
Not tested: the subject has only an id, a species and a description ("Rat cortical embryonic cell culture"); neither the files nor the archive's metadata give an age, a plating date or days in vitro (fields printed in the data check). The session dates (2021-10-06 and 2021-10-11) are 5 days apart, which fits either range of the criterion but does not give the age.
fig6-ais-from-dendrite¶
In both cells the axon initial segment does not leave the soma but a dendrite (an axon-bearing dendrite).
Data elsewhere: the reconstructed morphologies cell_models/<cell>/morphology_corrected.swc of the
repository at the pinned commit. Criterion: the AIS is the path labelled "apic" (type 4, Methods
5.7.2); reproduced when, in both cells, the first AIS point's parent is a dendrite point (type 3),
not the soma (type 1).
import gtmodel as gm
root = gm.fetch_repo()
for s, code_name in mm.CODE_NAMES.items():
a = gm.read_swc(root / "cell_models" / code_name / "morphology_corrected.swc")
ids, typ, par = a[:, 0].astype(int), a[:, 1].astype(int), a[:, 6].astype(int)
pos = {i: k for k, i in enumerate(ids)}
counts = dict(Counter(typ.tolist()))
roots = [k for k in np.where(typ == 4)[0] if par[k] == -1 or typ[pos[par[k]]] != 4]
print(f"{mm.CELLS[s]} ({code_name}): points by type {counts}")
for k in roots:
p = pos[par[k]]
soma = np.where(typ == 1)[0][0]
print(f" AIS starts at point {ids[k]}; its parent {ids[p]} has type {typ[p]}; "
f"distance of the AIS start from the soma point {np.linalg.norm(a[k, 2:5] - a[soma, 2:5]):.1f} um")
ax_roots = [k for k in np.where(typ == 2)[0] if typ[pos[par[k]]] != 2]
print(f" axon (type 2) starts from a point of type {[int(typ[pos[par[k]]]) for k in ax_roots]}")
cell 1 (cell1_211006_3148): points by type {1: 1, 3: 61472, 4: 184, 2: 18455}
AIS starts at point 53831; its parent 42786 has type 3; distance of the AIS start from the soma point 42.9 um
axon (type 2) starts from a point of type [4]
cell 2 (cell1_211011_3436): points by type {1: 1, 3: 22496, 4: 192, 2: 56230}
AIS starts at point 1072; its parent 205 has type 3; distance of the AIS start from the soma point 39.1 um
axon (type 2) starts from a point of type [4]
Result: reproduced. In both reconstructions the single AIS path (type 4) starts from a dendrite point (type 3), not from the soma, and the axon (type 2) starts from the end of the AIS.
gt-rheobase¶
The ground-truth model's rheobase, estimated with a 270 ms step, is 150 pA.
Data elsewhere: the released hay_ais model of the repository. Criterion: the lowest amplitude of 270 ms steps that elicits a spike, with the step series of the ground-truth notebook (from 0.15 nA by 0.005 nA). Because that series starts at the stated value, a series from 0.10 nA is also run.
t_sim = time.time()
rheo_nb = gm.rheobase()
rheo_low = gm.rheobase(step_min=0.10)
print("notebook series (nA, spikes):", rheo_nb)
print("series from 0.10 nA (nA, spikes):", rheo_low)
print(f"rheobase: notebook series {rheo_nb[-1][0] * 1000:.0f} pA; lower series {rheo_low[-1][0] * 1000:.0f} pA; "
f"within 150 +/- 15 pA: {abs(rheo_nb[-1][0] * 1000 - 150) <= 15}, {abs(rheo_low[-1][0] * 1000 - 150) <= 15} "
f"({time.time() - t_sim:.0f} s)")
p300 = gm.protocol_json()["IDrest_300"]["stimuli"][1]
print(f"IDrest 300% in the repository's protocol file: {p300['amp']:.3f} nA, {p300['thresh_perc']:.0f}% -> 100% = {p300['amp'] / 3 * 1000:.0f} pA")
import contextlib, io
cell_lfpy, sim_lfpy = gm.lfpy_cell(gm.probe())
with contextlib.redirect_stdout(io.StringIO()):
code_rheo = {smin: gm.package().ecode.compute_rheobase_for_model(cell_lfpy, sim=sim_lfpy, step_min=smin, step_increment=0.005,
isolate=False)[0]
for smin in (0.15, 0.10)}
print("the code's compute_rheobase_for_model on the LFPy cell model (first step with exactly one spike), nA:",
{k: round(float(v), 4) for k, v in code_rheo.items()})
Warning: no DISPLAY environment variable. --No graphics will be displayed.
notebook series (nA, spikes): [(0.15, 1)] series from 0.10 nA (nA, spikes): [(0.1, 0), (0.105, 0), (0.11, 0), (0.115, 0), (0.12, 0), (0.125, 0), (0.13, 0), (0.135, 0), (0.14, 1)] rheobase: notebook series 150 pA; lower series 140 pA; within 150 +/- 15 pA: True, True (13 s) IDrest 300% in the repository's protocol file: 0.450 nA, 300% -> 100% = 150 pA
the code's compute_rheobase_for_model on the LFPy cell model (first step with exactly one spike), nA: {0.15: 0.15, 0.1: 0.14}
Result: reproduced. With the ground-truth notebook's series the first step, 150 pA, fires one
spike. Starting lower, the model first fires at 140 pA, also within 150 +/- 15 pA. The code's own
compute_rheobase_for_model on the LFPy cell model gives the same two values. The repository's
IDrest 300% step of 0.45 nA implies the same 150 pA scale.
fig3b-ap-initiates-distal-ais¶
In the ground-truth model the action potential starts in the distal AIS and travels back toward the soma.
Criterion: IDrest 300% on the released model, membrane potential at the soma, the middle and the end of the AIS, cut-outs averaged around the somatic peaks; readings (a) the time the averaged AP crosses -20 mV and (b) the time of its peak; reproduced when under both the distal AIS comes first, then the middle, then the soma.
The run is the code's: 7-AIS-test.ipynb's protocol (define_protocols with the AIS recordings) on
the LFPy cell model and LFPy simulator, and its get_peak_cutout (3 ms before, 10 ms after, first
and last spikes left out). That function interpolates each trace on a 20 kHz grid, so the averaged
cut-outs have 0.05 ms samples. The paper does not say at what resolution the order is read, so the
times are given at the sample (first sample at or above -20 mV, largest sample), and between samples
(linear interpolation for the crossing, a parabola through the three samples around the maximum for
the peak). As a check below the code's grid, each spike is also timed on a fixed-step NEURON run
(dt 0.005 ms) of the whole protocol.
t_sim = time.time()
resp300, prots300 = gm.code_idrest300()
cut = gm.code_cutouts(resp300)
st_ = resp300["IDrest_300.soma.v"]
n_pk = len(gm.peak_times(np.asarray(st_["time"]), np.asarray(st_["voltage"])))
t_soma, v_soma = cut["soma.v"]
t_ref = t_soma[np.argmax(v_soma)]
order = []
for name in ("ais_end_v.v", "ais_mid_v.v", "soma.v"):
t, v = cut[name]
tp = t - t_ref
k = int(np.where(v >= -20)[0][0])
lin = tp[k - 1] + (-20 - v[k - 1]) / (v[k] - v[k - 1]) * (tp[k] - tp[k - 1])
j = int(np.argmax(v))
a, b, c = v[j - 1], v[j], v[j + 1]
par = tp[j] + 0.5 * (a - c) / (a - 2 * b + c) * (tp[1] - tp[0])
order.append(dict(location=name.split(".")[0], cross_sample_ms=round(float(tp[k]), 3), cross_linear_ms=round(float(lin), 4),
peak_sample_ms=round(float(tp[j]), 3), peak_parabola_ms=round(float(par), 4), peak_mV=round(float(b), 1)))
order = pd.DataFrame(order)
print(f"code's cut-outs: {n_pk} somatic spikes, {n_pk - 2} averaged; grid {np.diff(t_soma)[0]:.3f} ms; "
"times relative to the somatic peak sample")
print(order.to_string(index=False))
def strict(col):
return bool(order[col].is_monotonic_increasing and order[col].is_unique)
for col in ("cross_sample_ms", "cross_linear_ms", "peak_sample_ms", "peak_parabola_ms"):
print(f"{col}: distal AIS, then middle, then soma: {strict(col)}; distal minus middle "
f"{order[col][0] - order[col][1]:+.4f} ms; middle minus soma {order[col][1] - order[col][2]:+.4f} ms")
sim5 = gm.fixed_step_simulator(0.005)
r5 = gm.run_step(p300["amp"], p300["delay"], p300["duration"], p300["totduration"], ais=True, sim=sim5)
pt = pd.DataFrame(gm.per_spike_timing(r5), columns=["cross_end", "cross_mid", "cross_soma", "peak_end", "peak_mid", "peak_soma"])
print(f"\nfixed step 0.005 ms: {len(pt) + 2} somatic spikes, {len(pt)} timed (first and last left out); ms from each somatic peak")
print(pt.describe().loc[["mean", "min", "max"]].round(4).to_string())
for kind in ("cross", "peak"):
d1 = pt[f"{kind}_end"] - pt[f"{kind}_mid"]
d2 = pt[f"{kind}_mid"] - pt[f"{kind}_soma"]
print(f"{kind}: distal minus middle {d1.min():+.4f} to {d1.max():+.4f} ms; middle minus soma {d2.min():+.4f} to {d2.max():+.4f} ms; "
f"order distal, middle, soma in {int(((d1 < 0) & (d2 < 0)).sum())} of {len(pt)} spikes")
print(f"({time.time() - t_sim:.0f} s)")
fig, ax = plt.subplots(figsize=(6, 3.5))
for name in ("soma.v", "ais_mid_v.v", "ais_end_v.v"):
t, v = cut[name]
ax.plot(t - t_ref, v, ".-", ms=3, label=name.split(".")[0])
ax.axhline(-20, color="gray", lw=0.8)
ax.set(xlim=(-1, 1), ylim=(-80, 40), xlabel="ms from somatic peak sample", ylabel="mV", title="code's averaged cut-outs (20 kHz)")
ax.legend(); plt.show()
w = (t_soma - t_ref >= -1) & (t_soma - t_ref <= 1)
print(f"plotted window -1 to 1 ms of the cut-out; values there {min(cut[n][1][w].min() for n in cut):.1f} to "
f"{max(cut[n][1][w].max() for n in cut):.1f} mV, inside the axes")
code's cut-outs: 14 somatic spikes, 12 averaged; grid 0.050 ms; times relative to the somatic peak sample
location cross_sample_ms cross_linear_ms peak_sample_ms peak_parabola_ms peak_mV
ais_end_v -0.35 -0.3996 -0.25 -0.2558 29.8
ais_mid_v -0.35 -0.3954 -0.25 -0.2356 33.2
soma -0.20 -0.2246 0.00 -0.0145 27.4
cross_sample_ms: distal AIS, then middle, then soma: False; distal minus middle +0.0000 ms; middle minus soma -0.1500 ms
cross_linear_ms: distal AIS, then middle, then soma: True; distal minus middle -0.0042 ms; middle minus soma -0.1708 ms
peak_sample_ms: distal AIS, then middle, then soma: False; distal minus middle +0.0000 ms; middle minus soma -0.2500 ms
peak_parabola_ms: distal AIS, then middle, then soma: True; distal minus middle -0.0202 ms; middle minus soma -0.2211 ms
fixed step 0.005 ms: 14 somatic spikes, 12 timed (first and last left out); ms from each somatic peak
cross_end cross_mid cross_soma peak_end peak_mid peak_soma
mean -0.3765 -0.3711 -0.2047 -0.2696 -0.2492 -0.0133
min -0.4219 -0.4165 -0.2497 -0.3150 -0.2950 -0.0550
max -0.3268 -0.3217 -0.1538 -0.2200 -0.2000 0.0400
cross: distal minus middle -0.0070 to -0.0051 ms; middle minus soma -0.1679 to -0.1621 ms; order distal, middle, soma in 12 of 12 spikes
peak: distal minus middle -0.0250 to -0.0200 ms; middle minus soma -0.2400 to -0.2300 ms; order distal, middle, soma in 12 of 12 spikes
(88 s)
plotted window -1 to 1 ms of the cut-out; values there -66.8 to 33.2 mV, inside the axes
Result: reproduced, resting on an assumption about the time resolution. On the code's own cut-outs (12 spikes averaged, 0.05 ms samples) the distal and middle AIS reach -20 mV in the same sample (-0.35 ms) and peak in the same sample (-0.25 ms), so at the sample the two tie under both readings, and the soma comes 0.15 and 0.25 ms later. Between samples the distal AIS leads the middle by 0.0042 ms at the crossing and 0.0202 ms at the peak. A fixed-step run (dt 0.005 ms) times each of the 12 spikes on its own: the distal AIS leads the middle by 0.0051 to 0.0070 ms at the crossing and 0.020 to 0.025 ms at the peak, and the middle leads the soma by 0.162 to 0.168 and 0.23 to 0.24 ms, in 12 of 12 spikes under both readings. The order of the claim holds whenever the resolution can separate the two AIS points; the outcome takes the times between samples (an assumption: the paper does not say how the order is read).
fig3a-largest-channel-at-ais¶
On the simulated 80-electrode MEA, the largest extracellular signal is on the electrode closest to the AIS.
Data elsewhere: the ground-truth template cell_models/hay_ais/fitting/extracellular/template.npy
with its locations.npy (20 by 4 electrodes, 50 um pitch). Criterion: the channel with the largest
peak-to-peak amplitude is the electrode nearest the AIS midpoint in the plane; the largest absolute
negative peak is shown too.
The AIS is placed in the frame of the simulation: the segment midpoints of the instantiated
LFPyCellModel, where LFPy puts the soma midpoint at the origin. To show that the repository's
template was simulated in this frame, the template is recomputed here from the IDrest 300% run above
with the code's calculate_eap on the code's probe, and again with the probe moved by a few um.
t_sim = time.time()
ext = root / "cell_models" / "hay_ais" / "fitting" / "extracellular"
eap, loc, fs_gt = np.load(ext / "template.npy"), np.load(ext / "locations.npy"), float(np.load(ext / "fs.npy"))
print(f"repository template {eap.shape} (channels, samples) at {fs_gt:.0f} Hz; electrodes x {loc[:, 0].min()} to {loc[:, 0].max()}, "
f"y {loc[:, 1].min()} to {loc[:, 1].max()} um, z {np.unique(loc[:, 2])}")
pr0 = gm.probe()
print("the code's planar probe has the template's electrode positions:", bool(np.allclose(pr0.positions, loc)))
def compare(e):
p1, p0 = np.ptp(e, axis=1), np.ptp(eap, axis=1)
return dict(largest_ptp_channel=int(np.argmax(p1)), corr_waveforms=round(float(np.corrcoef(e.ravel(), eap.ravel())[0, 1]), 5),
corr_ptp=round(float(np.corrcoef(p1, p0)[0, 1]), 5),
rel_rms_difference=round(float(np.sqrt(np.mean((e - eap) ** 2)) / np.sqrt(np.mean(eap ** 2))), 4))
eap_here = gm.template(resp300, prots300)
rows = [dict(probe_shift_um="none", **compare(eap_here))]
for dx, dy in ((2, 0), (-2, 0), (5, 0), (-5, 0), (0, 5), (0, -5)):
pr = gm.probe()
pr.move([dx, dy, 0])
r_, p_ = gm.code_idrest300(electrode=pr, ais=False)
rows.append(dict(probe_shift_um=f"({dx}, {dy})", **compare(gm.template(r_, p_))))
print("\nTemplate recomputed with the code (IDrest 300%, calculate_eap) against the repository's file:")
print(pd.DataFrame(rows).to_string(index=False))
mids = gm.lfpy_segment_midpoints()
ais = mids["ais"]
n_ais = len(ais)
ais_mid = ais[n_ais // 2 - 1:n_ais // 2 + 1].mean(axis=0) if n_ais % 2 == 0 else ais[n_ais // 2]
print(f"\nLFPy frame: soma midpoint {mids['_somapos'].round(2)}; AIS {n_ais} segments, midpoints from "
f"{ais[0, :2].round(1)} to {ais[-1, :2].round(1)} um (x, y); AIS midpoint {ais_mid[:2].round(1)}")
d_ais = np.linalg.norm(loc[:, :2] - ais_mid[:2], axis=1)
ptp = np.ptp(eap, axis=1)
neg = -eap.min(axis=1)
near = np.argsort(d_ais)[:3]
print("electrodes nearest the AIS midpoint (channel, position, distance um):",
[(int(k), loc[k, :2].tolist(), round(float(d_ais[k]), 1)) for k in near])
top = np.argsort(-ptp)[:3]
print("largest peak-to-peak (channel, position, fraction of the largest, distance to the AIS midpoint um):",
[(int(k), loc[k, :2].tolist(), round(float(ptp[k] / ptp.max()), 3), round(float(d_ais[k]), 1)) for k in top])
print(f"largest negative peak: channel {int(np.argmax(neg))} at {loc[np.argmax(neg), :2].tolist()}")
print("largest peak-to-peak on the nearest electrode:", int(np.argmax(ptp)) == int(near[0]),
"; largest negative peak on the nearest electrode:", int(np.argmax(neg)) == int(near[0]))
d_seg = np.linalg.norm(loc[:, None, :2] - ais[None, :, :2], axis=2)
print("electrode nearest each AIS segment midpoint:", sorted(set(np.argmin(d_seg, axis=0).tolist())))
pts = gm.section_points()
soma3d = [v for k, v in pts.items() if k.endswith("soma[0]")][0]
ais3d = [v for k, v in pts.items() if k.endswith("ais[0]")][0]
mid_v = ais3d[len(ais3d) // 2] - soma3d.mean(axis=0)
d_v = np.linalg.norm(loc[:, :2] - mid_v[:2], axis=1)
print(f"\nvariant (round 1): AIS midpoint from the 3D points relative to the mean soma point, {mid_v[:2].round(1)}; "
f"nearest electrode {int(np.argmin(d_v))} ({np.sort(d_v)[0]:.1f} um; next {np.sort(d_v)[1]:.1f} um)")
print(f"({time.time() - t_sim:.0f} s)")
repository template (80, 160) (channels, samples) at 40000 Hz; electrodes x -75.0 to 75.0, y -175.0 to 775.0 um, z [20.] the code's planar probe has the template's electrode positions: True
Template recomputed with the code (IDrest 300%, calculate_eap) against the repository's file:
probe_shift_um largest_ptp_channel corr_waveforms corr_ptp rel_rms_difference
none 43 0.99992 0.99999 0.0127
(2, 0) 24 0.99668 0.99702 0.0828
(-2, 0) 43 0.99799 0.99770 0.0658
(5, 0) 24 0.98975 0.98879 0.1428
(-5, 0) 43 0.98930 0.98770 0.1541
(0, 5) 43 0.98794 0.98923 0.1607
(0, -5) 44 0.98759 0.99116 0.1582
LFPy frame: soma midpoint [0. 0. 0.]; AIS 8 segments, midpoints from [11.8 -1. ] to [ 13.8 -28.3] um (x, y); AIS midpoint [ 12.6 -13.7] electrodes nearest the AIS midpoint (channel, position, distance um): [(43, [25.0, -25.0], 16.8), (23, [-25.0, -25.0], 39.3), (44, [25.0, 25.0], 40.6)] largest peak-to-peak (channel, position, fraction of the largest, distance to the AIS midpoint um): [(43, [25.0, -25.0], 1.0, 16.8), (24, [-25.0, 25.0], 0.939, 54.0), (44, [25.0, 25.0], 0.878, 40.6)] largest negative peak: channel 43 at [25.0, -25.0] largest peak-to-peak on the nearest electrode: True ; largest negative peak on the nearest electrode: True electrode nearest each AIS segment midpoint: [43]
variant (round 1): AIS midpoint from the 3D points relative to the mean soma point, [ 1.7 -23.4]; nearest electrode 43 (23.3 um; next 26.8 um) (42 s)
Result: reproduced. The template recomputed with the code on the LFPy cell model matches the repository's file (correlation 0.99992, RMS difference 1.3% of the template), and moving the probe by 2 to 5 um raises the difference to 6.6 to 16.1%, so the file was simulated in the LFPy frame read here. In that frame the AIS runs from (11.8, -1.0) to (13.8, -28.3) um, its midpoint is (12.6, -13.7) um, and the nearest electrode is channel 43 at (25, -25) um, 16.8 um away, against 39.3 um for the next; it is also the nearest electrode to every AIS segment. Channel 43 has the largest peak-to-peak and the largest negative peak. The second-largest peak-to-peak is 0.939 of the largest (channel 24, 54.0 um from the AIS midpoint), so the largest channel itself is sensitive to the placement: a probe moved by 2 um in x makes channel 24 the largest.
table1-feature-counts¶
For the ground-truth model, the soma strategy uses 76 intracellular features; all adds 11 extracellular features (87), sections adds 33 (109), and single adds 77 (153).
Data elsewhere: the repository's feature and protocol files for hay_ais. The objectives are counted
on the evaluator the code builds for each strategy (create_evaluator, one objective per feature,
define_fitness_calculator), as the criterion says; the feature files features_BPO_<strategy>.json
are counted as a check. The Methods' 79 features are shown as a variant.
counts = {}
for strategy in ("soma", "all", "sections", "single"):
ev = gm.evaluator(strategy)
objs = ev.fitness_calculator.objectives
extra = sum("MEA" in o.name for o in objs)
f = gm.features_json(name=f"features_BPO_{strategy}.json")
f_intra = sum(len(v) for p in f.values() for loc_, v in p.items() if loc_ != "MEA")
f_extra = sum(len(v) for p in f.values() for loc_, v in p.items() if loc_ == "MEA")
counts[strategy] = (len(objs) - extra, extra, len(objs))
print(f"{strategy}: evaluator objectives {len(objs)} (intracellular {len(objs) - extra}, extracellular {extra}); "
f"feature file {f_intra} + {f_extra}")
expected = {"soma": 76, "all": 87, "sections": 109, "single": 153}
print("totals match Table 1:", {k: counts[k][2] == v for k, v in expected.items()},
"; intracellular 76 for every strategy:", all(c[0] == 76 for c in counts.values()))
fj = gm.features_json(name="features.json")
train = ("IDrest_150", "IDrest_250", "IDrest_300", "IV_-100", "IV_-20", "APWaveform_290")
n79 = sum(len(v) for p, d in fj.items() if p in train for v in d.values())
print(f"variant: features.json, the six training protocols: {n79} features")
soma: evaluator objectives 76 (intracellular 76, extracellular 0); feature file 76 + 0
all: evaluator objectives 87 (intracellular 76, extracellular 11); feature file 76 + 11
sections: evaluator objectives 109 (intracellular 76, extracellular 33); feature file 76 + 33
single: evaluator objectives 153 (intracellular 76, extracellular 77); feature file 76 + 77
totals match Table 1: {'soma': True, 'all': True, 'sections': True, 'single': True} ; intracellular 76 for every strategy: True
variant: features.json, the six training protocols: 79 features
Result: reproduced. The code's evaluators have 76 intracellular objectives for every strategy and 0, 11, 33 and 77 extracellular ones, totals 76, 87, 109 and 153, as in Table 1; the feature files give the same counts. The Methods' 79 is the count in features.json for the six training protocols, before the code drops the three sag features of IV -20%.
The model fits (Figures 4, 5, 7, 8 and 10)¶
All remaining claims need the optimizations: 10 seeds per strategy, 600 generations of 20 offspring
each (scripts/run_optimizations.py), for four strategies on the ground-truth model and two
strategies on each experimental cell, with extracellular signals simulated for the extracellular
strategies. No optimization results are in the repository. The cost is estimated from one
evaluation of the six training protocols on the released ground-truth model: intracellular only with
the NEURON simulator, and with the code's evaluator for the all strategy on the LFPy simulator,
which adds the extracellular signals. The experimental models are taken to cost the same, so the
estimate is a rough lower bound.
t0 = time.time()
prot = gm.protocol_json()
per = {}
for name in train:
st = prot[name]["stimuli"][1]
t1 = time.time()
gm.run_step(st["amp"], st["delay"], st["duration"], st["totduration"])
per[name] = time.time() - t1
t_eval = time.time() - t0
print("seconds per training protocol:", {k: round(v, 2) for k, v in per.items()}, f"; one intracellular evaluation {t_eval:.1f} s")
ev_all = gm.evaluator("all")
pv = gm.release_params(ev_all.param_names)
t1 = time.time()
ev_all.run_protocols(ev_all.fitness_protocols.values(), param_values=pv)
t_eval_x = time.time() - t1
print(f"one evaluation of the all strategy with the code's evaluator (LFPy, extracellular included): {t_eval_x:.1f} s")
per_run = 600 * 20
cost_gt = 10 * per_run * (t_eval + 3 * t_eval_x)
cost_exp = 2 * 10 * per_run * (t_eval + t_eval_x)
for label, n, c in (("ground truth, 4 strategies", 10 * 4 * per_run, cost_gt),
("experimental, 2 cells x 2 strategies", 10 * 2 * 2 * per_run, cost_exp)):
print(f"{label}: {n} evaluations, at least {c / 3600:.0f} CPU hours, {c / 3600 / 4:.0f} hours on the 4 cores of this runner")
print(f"each estimate above 300 CPU hours: {min(cost_gt, cost_exp) / 3600 > 300}; above 60 hours on 4 cores: {min(cost_gt, cost_exp) / 3600 / 4 > 60}")
print(f"Notebook limit: 0.5 hours; whole author run: 5 hours. ({time.time() - started:.0f} s since the start of the notebook)")
seconds per training protocol: {'IDrest_150': 3.08, 'IDrest_250': 3.66, 'IDrest_300': 3.74, 'IV_-100': 0.44, 'IV_-20': 0.41, 'APWaveform_290': 0.83} ; one intracellular evaluation 12.2 s
one evaluation of the all strategy with the code's evaluator (LFPy, extracellular included): 20.4 s ground truth, 4 strategies: 480000 evaluations, at least 2446 CPU hours, 611 hours on the 4 cores of this runner experimental, 2 cells x 2 strategies: 480000 evaluations, at least 2171 CPU hours, 543 hours on the 4 cores of this runner each estimate above 300 CPU hours: True; above 60 hours on 4 cores: True Notebook limit: 0.5 hours; whole author run: 5 hours. (544 s since the start of the notebook)
Result: the fits are blocked as too large. The printed estimates, from one evaluation of the six training protocols intracellularly and one with the extracellular signals, are each above 300 CPU hours, for the ground-truth fits and for the experimental fits, and above 60 hours on this runner's 4 cores, far beyond the 30-minute notebook limit and the run's 5 hours. The exact figures vary with the runner's load; the record gives those of the run it cites. Every claim on the fits, and those that use their solutions (Figures 7B, 7C/8C and 10), is blocked.
What can be computed without the fits¶
Some parts of the blocked claims need no fitted model. They are computed here; none decides a claim, because each claim compares fitted models with these values.
- fig4c-hyperdepol-only-single: the ground truth's spike count in the 270 ms depolarizing phase of each HyperDepol sweep (eFEL, -20 mV), from the validation run of the code's evaluator.
- fig4b-intra-no-difference and fig4b-extra-better-than-soma: the number of validation features,
counted as 4-optimization-performance.ipynb counts them (one entry per feature name of the
validation evaluator; the extracellular features of the firepattern 120% template on 80
electrodes,
compute_extra_features). - fig78b-validation-spiking: the range over runs of the recorded spike count in each validation sweep, reading (a)'s reference.
t_val = time.time()
ev_val = gm.evaluator(None, validation=True)
objs = ev_val.fitness_calculator.objectives
pv = gm.release_params(ev_val.param_names)
resp_val = ev_val.run_protocols(ev_val.fitness_protocols.values(), param_values=pv)
print(f"validation protocols: {sorted(ev_val.fitness_protocols)}; simulated in {time.time() - t_val:.0f} s")
feat_names = {o.features[0].name for o in objs}
print(f"validation evaluator: {len(objs)} objectives, {len(feat_names)} distinct feature names (the notebook's count; claim: 86)")
kw = gm.package().utils.get_extra_kwargs()
eap_val = gm.package().utils.calculate_eap(responses=resp_val, protocols=ev_val.fitness_protocols, protocol_name="firepattern_120", **kw)
xf = gm.package().efeatures_extraction.compute_extra_features(eap_val, fs=kw["fs"], upsample=kw["upsample"])
n_x = sum(len(v) for v in xf.values())
n_fin = sum(int(np.isfinite(np.asarray(v, float)).sum()) for v in xf.values())
print(f"extracellular validation features: {len(xf)} kinds x {eap_val.shape[0]} electrodes = {n_x}; finite {n_fin} "
f"(not finite by kind: { {k: int((~np.isfinite(np.asarray(v, float))).sum()) for k, v in xf.items() if (~np.isfinite(np.asarray(v, float))).any()} }); claim: 796")
hd = []
for k, r in resp_val.items():
if k.startswith("HyperDepol") and k.endswith("soma.v"):
tt, vv = np.asarray(r["time"]), np.asarray(r["voltage"])
pk = gm.peak_times(tt, vv)
hd.append(dict(sweep=k.split(".")[0], spikes_hyperpolarizing_phase=int(((pk >= 250.0) & (pk < 700.0)).sum()),
spikes_depolarizing_phase=int(((pk >= 700.0) & (pk < 970.0)).sum()), spikes_whole_sweep=len(pk)))
print("\nground truth, HyperDepol (eFEL peak times at -20 mV; hyperpolarizing phase 250 to 700 ms, depolarizing phase 700 to 970 ms, protocols.json tmid and toff):")
print(pd.DataFrame(hd).sort_values("sweep").to_string(index=False))
val = []
for s in urls:
for r in runs[s]:
for prot_ in ("firepattern", "hyperdepol", "sahp", "poscheops"):
sws = sorted([d for d in patch[s] if d["run"] == r and d["protocol"] == prot_], key=lambda d: d["sweep"])
for k, d in enumerate(sws):
val.append(dict(cell=mm.CELLS[s], run=r, protocol=prot_, sweep=k, spikes=len(mm.peak_times(d["t"], d["v"]))))
val = pd.DataFrame(val)
print("\nrecorded spike count per validation sweep (whole sweep), range over runs:")
print(val.groupby(["cell", "protocol", "sweep"]).spikes.agg(["min", "max", "size"]).rename(columns={"size": "runs"})
.unstack("cell").to_string())
print(f"({time.time() - t_val:.0f} s)")
print(f"whole notebook: {(time.time() - started) / 60:.1f} min")
validation protocols: ['HyperDepol_-120', 'HyperDepol_-160', 'HyperDepol_-40', 'HyperDepol_-80', 'PosCheops_300', 'firepattern_120', 'firepattern_200', 'sAHP_150', 'sAHP_200', 'sAHP_250', 'sAHP_300']; simulated in 93 s
validation evaluator: 110 objectives, 86 distinct feature names (the notebook's count; claim: 86)
extracellular validation features: 11 kinds x 80 electrodes = 880; finite 800 (not finite by kind: {'recovery_slope': 80}); claim: 796
ground truth, HyperDepol (eFEL peak times at -20 mV; hyperpolarizing phase 250 to 700 ms, depolarizing phase 700 to 970 ms, protocols.json tmid and toff):
sweep spikes_hyperpolarizing_phase spikes_depolarizing_phase spikes_whole_sweep
HyperDepol_-120 0 0 0
HyperDepol_-160 0 0 0
HyperDepol_-40 0 1 1
HyperDepol_-80 0 1 1
recorded spike count per validation sweep (whole sweep), range over runs:
min max runs
cell cell 1 cell 2 cell 1 cell 2 cell 1 cell 2
protocol sweep
firepattern 0 12 0 15 1 5 4
1 33 14 37 30 5 4
hyperdepol 0 0 0 0 0 5 4
1 0 0 0 0 5 4
2 0 0 0 0 5 4
3 0 0 0 0 5 4
poscheops 0 76 52 90 80 5 4
sahp 0 1 0 1 2 5 4
1 2 3 3 5 5 4
2 3 7 4 8 5 4
3 4 8 4 10 5 4
(94 s)
whole notebook: 10.6 min
These values are references only. The ground truth fires no spike in the depolarizing phase of HyperDepol -160% and one at -40%, as fig4c-hyperdepol-only-single expects of it. The validation evaluator has 86 distinct intracellular features, the number fig4b-intra-no-difference gives. The extracellular features number 880, 800 of them finite, within the 796 +/- 80 of fig4b-extra-better-than-soma. In the recordings no HyperDepol sweep has a spike in any run of either cell. The claims stay blocked, since each compares fitted models with these values.
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 | Data | Outcome |
|---|---|---|
| recordings-two-cells-runs | dandiset | reproduced |
| ecode-protocol-structure | dandiset | reproduced |
| holding-potential | dandiset | reproduced |
| extracellular-spikes-above-50uv | dandiset | reproduced |
| mea-electrodes-and-rates | dandiset | partially reproduced |
| culture-age | dandiset | blocked: metadata insufficient |
| fig6-ais-from-dendrite | elsewhere | reproduced |
| gt-rheobase | elsewhere | reproduced |
| fig3b-ap-initiates-distal-ais | elsewhere | reproduced (assumed: time resolution) |
| fig3a-largest-channel-at-ais | elsewhere | reproduced |
| table1-feature-counts | elsewhere | reproduced |
| fig4a, gt-intracellular-scores, fig4c, fig4b (two), fig4d, fig5 | elsewhere | blocked: too large |
| fig78a, experimental-intracellular-scores, fig78b, fig78c, fig10 | with inputs | blocked: too large |
How the analysis differs from the paper's. The template follows Methods 5.7.1 with the code's
parameters, on the filtered MEA series of the dandiset, after the code's removal of blank channels;
the filter is applied forward and backward. The holding potential is read before the stimulus onset
found in the stimulus, 50 ms before the code's ton. The ground-truth model is simulated with
BluePyOpt 1.14.26 in place of 1.13.88: the AIS timing, the template check and the evaluators use its
LFPy simulator as the paper's notebooks do, and the rheobase and the fixed-step AIS check its NEURON
simulator. The AIS order is read between the samples of the code's 0.05 ms cut-outs (an
assumption), where at the sample the distal and middle AIS tie. Table 1 is counted on the code's
evaluators.
What was found about the dataset. The patch voltages are in mV and the currents in pA, though
the files say volts and amperes. Most protocols are sampled at 20.41 kHz and APWaveform at 125 kHz,
where the paper says 20 kHz and about 200 kHz. The icephys grouping tables leave out the first sweep
of most protocols and every PosCheops sweep. Many MEA channels are constant for part of a run (117 to
132 per run). In the files, the stimulus starts 50 ms before the code's ton timings. No culture age is given. The dandiset's metadata cite the paper as
doi:10.1162/neco_a_01676, not 10.1162/neco_a_01672.