Reproducing Quattrocolo & Maccaferri (2013) from CRCNS hc-9¶

Paper. Quattrocolo G, Maccaferri G (2013). Novel GABAergic circuits mediating excitation/inhibition of Cajal-Retzius cells in the developing hippocampus. J Neurosci 33(13):5486-5498. doi:10.1523/jneurosci.5680-12.2013. The claims file (claims/crcns-hc-9/) was written from the version of record, which has no supplement and links no code.

Data set. CRCNS hc-9, "Firing pattern of O-LM cells in mouse hippocampal CA1" (doi:10.6080/K0HQ3WTH), read from the public bucket at https://crcnsarchive.s3-us-west-2.amazonaws.com/hc-9/live/. CRCNS data sets are not versioned, so the version is recorded as unversioned. The data are pinned by MD5: the one data file, Maccaferri_data.zip, has MD5 4e3175055681b99ad5be06ba542a93b8, as the data set's checksums.md5 gives it, and the data check below confirms the MD5 of the zip as read. This is the only published copy of the data (the bucket's hc-9/old/ folder holds only earlier documentation), so there is no version to choose between; the paper cites no data set. The description read is version 0.7 (Dec 2, 2016) of live/docs/crcns_hc-9_data_description.pdf.

What the data set holds. Current-clamp responses of three O-LM interneurons to 1 s current steps (15 sweeps each), the firing patterns of the paper's Figure 7. None of the recordings from Cajal-Retzius cells behind the paper's main findings (Figures 1 to 6 and 8 to 12) are in it, and the light-evoked responses of Figure 7C and 7D are not either. The zip is opened in place over range requests; nothing is written to disk.

Claims. One of the 26 claims can be tested.

Claim Tested or blocked
fig7a-olm-sag tested
fig1-spontaneous-amplitude blocked, data missing
fig4-upsc-amplitude blocked, data missing
fig4-upsc-kinetics blocked, data missing
paired-presynaptic-neurogliaform blocked, data missing
upsc-gabazine-block blocked, data missing
fig5-no-gabab-component blocked, data missing
fig6-spontaneous-larger-than-unitary blocked, data missing
fig6-spontaneous-faster-than-unitary blocked, data missing
fig7-chr2-olm-light-firing blocked, data missing
fig8-light-evoked-psc blocked, data missing
fig9-opsc-rise-time blocked, data missing
fig9-opsc-amplitude-decay blocked, data missing
fig10-dhpg-increases-frequency blocked, data missing
fig10-dhpg-amplitude-unchanged blocked, data missing
fig10-reduced-slice-no-effect blocked, data missing
fig11-slm-inhibition blocked, data missing
fig11-slm-inhibition-gabazine blocked, data missing
fig11-slm-excitation blocked, data missing
fig11d-baseline-frequency-dichotomy blocked, data missing
gabazine-baseline-firing blocked, data missing
fig12a-slow-conductance-depolarized blocked, data missing
fig12a-slow-conductance-hyperpolarized blocked, data missing
fig12b-fast-conductance-depolarized blocked, data missing
fig12b-fast-conductance-hyperpolarized blocked, data missing
no-age-dependence blocked, data missing
In [1]:
import io
import itertools
import time
import zipfile

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from scipy import stats

import hc9

pd.set_option("display.width", 200)
pd.set_option("display.max_columns", 40)

started = time.time()
cells = hc9.load_cells()
print(f"loaded {len(cells)} cells from CRCNS {hc9.DATASET} ({hc9.VERSION}), {hc9.ZIP_URL}, "
      f"in {time.time() - started:.1f} s")
loaded 3 cells from CRCNS hc-9 (unversioned), https://crcnsarchive.s3-us-west-2.amazonaws.com/hc-9/live/Maccaferri_data.zip, in 3.3 s

Data check¶

The data set is one zip in the CRCNS bucket. Its MD5 is compared with the data set's own checksums.md5, and the zip's members are listed. Each check of the description against the values is stated under its output.

In [2]:
published = hc9.published_checksums()
print("checksums.md5:")
for k, v in published.items():
    print(f"  {v}  {k}")
md5 = hc9.zip_md5()
print(f"\nMD5 of Maccaferri_data.zip as read now: {md5}")
print("matches checksums.md5:", md5 == published["Maccaferri_data.zip"],
      "; matches the MD5 pinned in hc9.py:", md5 == hc9.EXPECTED_ZIP_MD5)
print("\nMembers of the zip:")
print(hc9.zip_listing().to_string(index=False))
print("\nMD5 of each ATF member:")
for name, c in cells.items():
    print(f"  {c['raw_md5']}  {name}_raw.atf")
checksums.md5:
  32366f0262dc21cb97cb1f1efaade1a7  docs/crcns_hc-9_data_description.docx
  7b36e1831089cad1aa854b0e1bf9b011  docs/crcns_hc-9_data_description.pdf
  23542905cd6ea65f31f492f63fae0bfc  filelist.txt
  4e3175055681b99ad5be06ba542a93b8  Maccaferri_data.zip
MD5 of Maccaferri_data.zip as read now: 4e3175055681b99ad5be06ba542a93b8
matches checksums.md5: True ; matches the MD5 pinned in hc9.py: True

Members of the zip:
                                  member   bytes       date
                        Maccaferri_data/       0 2016-11-23
     Maccaferri_data/1May2012_P3_raw.atf 7388537 2014-05-06
   Maccaferri_data/20Sept2011_P2_raw.atf 7388498 2014-05-06
Maccaferri_data/24October2012_C2_raw.atf 7385994 2014-05-06

MD5 of each ATF member:
  3a1a661c64058505f0d221c3ba309b6b  1May2012_P3_raw.atf
  c4097c125efe088d6666fad6371d7947  20Sept2011_P2_raw.atf
  a4601337699cf96ea8aa3b0413a15aa4  24October2012_C2_raw.atf
In [3]:
for name, c in cells.items():
    atf = c["atf"]
    print(f"== {name}: {atf['version']}, {atf['n_header']} optional header records, {atf['n_cols']} columns")
    for h in atf["header"]:
        print("   ", h[:160])
    print("   columns:", atf["columns"][:5], "...", atf["columns"][-2:])
    dt = np.diff(c["t"])
    print(f"   samples {len(c['t'])}, time {c['t'][0]:.4f} to {c['t'][-1]:.4f} s, "
          f"interval {dt.min()*1e6:.3f} to {dt.max()*1e6:.3f} us, sampling {c['fs']:.0f} Hz")
    print(f"   sweeps {c['v'].shape[1]}; voltage {c['v'].min():.1f} to {c['v'].max():.1f} mV; "
          f"current {c['i'].min():.1f} to {c['i'].max():.1f} pA")
== 1May2012_P3: ATF	1.0, 8 optional header records, 31 columns
    "AcquisitionMode=Episodic Stimulation"
    "Comment="
    "YTop=200,4291.76"
    "YBottom=-200,-4334.72"
    "SyncTimeUnits=12.5"
    "SweepStartTimesMS=0.000,5000.000,10000.000,15000.000,20000.000,25000.000,30000.000,35000.000,40000.000,45000.000,50000.000,55000.000,60000.000,65000.000,70000.
    "SignalsExported=2CC mV,AO #1"
    "Signals="	"2CC mV"	"AO #1"	"2CC mV"	"AO #1"	"2CC mV"	"AO #1"	"2CC mV"	"AO #1"	"2CC mV"	"AO #1"	"2CC mV"	"AO #1"	"2CC mV"	"AO #1"	"2CC mV"	"AO #1"	"2CC mV"	"AO 
   columns: ['Time (s)', 'Trace #1 (mV)', 'Trace #1 (pA)', 'Trace #2 (mV)', 'Trace #2 (pA)'] ... ['Trace #15 (mV)', 'Trace #15 (pA)']
   samples 39999, time 0.0000 to 1.9999 s, interval 50.000 to 50.000 us, sampling 20000 Hz
   sweeps 15; voltage -102.8 to 33.9 mV; current -100.0 to 250.0 pA
== 20Sept2011_P2: ATF	1.0, 8 optional header records, 31 columns
    "AcquisitionMode=Episodic Stimulation"
    "Comment="
    "YTop=200,4291.76"
    "YBottom=-200,-4334.72"
    "SyncTimeUnits=12.5"
    "SweepStartTimesMS=0.000,5000.000,10000.000,15000.000,20000.000,25000.000,30000.000,35000.000,40000.000,45000.000,50000.000,55000.000,60000.000,65000.000,70000.
    "SignalsExported=2CC mV,AO #1"
    "Signals="	"2CC mV"	"AO #1"	"2CC mV"	"AO #1"	"2CC mV"	"AO #1"	"2CC mV"	"AO #1"	"2CC mV"	"AO #1"	"2CC mV"	"AO #1"	"2CC mV"	"AO #1"	"2CC mV"	"AO #1"	"2CC mV"	"AO 
   columns: ['Time (s)', 'Trace #1 (mV)', 'Trace #1 (pA)', 'Trace #2 (mV)', 'Trace #2 (pA)'] ... ['Trace #15 (mV)', 'Trace #15 (pA)']
   samples 39999, time 0.0000 to 1.9999 s, interval 50.000 to 50.000 us, sampling 20000 Hz
   sweeps 15; voltage -99.0 to 29.9 mV; current -100.0 to 250.0 pA
== 24October2012_C2: ATF	1.0, 8 optional header records, 31 columns
    "AcquisitionMode=Episodic Stimulation"
    "Comment="
    "YTop=200,4295.29"
    "YBottom=-200,-4340.95"
    "SyncTimeUnits=12.5"
    "SweepStartTimesMS=0.000,5000.000,10000.000,15000.000,20000.000,25000.000,30000.000,35000.000,40000.000,45000.000,50000.000,55000.000,60000.000,65000.000,70000.
    "SignalsExported=2CC mV,AO #1"
    "Signals="	"2CC mV"	"AO #1"	"2CC mV"	"AO #1"	"2CC mV"	"AO #1"	"2CC mV"	"AO #1"	"2CC mV"	"AO #1"	"2CC mV"	"AO #1"	"2CC mV"	"AO #1"	"2CC mV"	"AO #1"	"2CC mV"	"AO 
   columns: ['Time (s)', 'Trace #1 (mV)', 'Trace #1 (pA)', 'Trace #2 (mV)', 'Trace #2 (pA)'] ... ['Trace #15 (mV)', 'Trace #15 (pA)']
   samples 39999, time 0.0000 to 1.9999 s, interval 50.000 to 50.000 us, sampling 20000 Hz
   sweeps 15; voltage -104.6 to 38.3 mV; current -100.0 to 250.0 pA

The ATF header agrees with the description: line 2 gives the number of further header records (8) and of columns (31), the 11th line names the columns, time is in seconds, voltage in mV and the injected current in pA. Each file holds 15 sweeps of 2 s sampled at 20 kHz (the Methods give 10 to 20 kHz), with sweep starts 5 s apart (the description says 3 to 5 s).

In [4]:
steps = pd.concat({name: hc9.step_table(c) for name, c in cells.items()}, names=["cell"]).reset_index(level=0)
print("Step level, onset and offset of every sweep, read from the current channel:")
print(steps.pivot(index="sweep", columns="cell", values="step_pA").to_string())
print("\nOnsets (s):", sorted(steps.on_s.dropna().round(5).unique()),
      " offsets (s):", sorted(steps.off_s.dropna().round(5).unique()),
      " holding current (pA):", sorted(steps.holding_pA.round(2).unique()))
levels = sorted(steps.step_pA.unique())
print("\nLevels present:", levels)
print("Increment:", np.unique(np.diff(levels)), "pA")
stated = list(range(-100, 251, 30))
print("Description, -100:30:+250 pA, would give", len(stated), "levels:", stated)
Step level, onset and offset of every sweep, read from the current channel:
cell   1May2012_P3  20Sept2011_P2  24October2012_C2
sweep                                              
1           -100.0         -100.0            -100.0
2            -75.0          -75.0             -75.0
3            -50.0          -50.0             -50.0
4            -25.0          -25.0             -25.0
5              0.0            0.0               0.0
6             25.0           25.0              25.0
7             50.0           50.0              50.0
8             75.0           75.0              75.0
9            100.0          100.0             100.0
10           125.0          125.0             125.0
11           150.0          150.0             150.0
12           175.0          175.0             175.0
13           200.0          200.0             200.0
14           225.0          225.0             225.0
15           250.0          250.0             250.0

Onsets (s): [np.float64(0.07505)]  offsets (s): [np.float64(1.07505)]  holding current (pA): [np.int64(0)]

Levels present: [np.float64(-100.0), np.float64(-75.0), np.float64(-50.0), np.float64(-25.0), np.float64(0.0), np.float64(25.0), np.float64(50.0), np.float64(75.0), np.float64(100.0), np.float64(125.0), np.float64(150.0), np.float64(175.0), np.float64(200.0), np.float64(225.0), np.float64(250.0)]
Increment: [25.] pA
Description, -100:30:+250 pA, would give 12 levels: [-100, -70, -40, -10, 20, 50, 80, 110, 140, 170, 200, 230]

The description gives the current steps as -100:30:+250 pA. The current channel shows 15 steps from -100 to +250 pA in 25 pA increments, one of them 0 pA, in all three files; -100:30:+250 would give 12 levels ending at +230 pA. The data are trusted here: the description is wrong on the increment, and this is recorded as a finding. Each step lasts 1 s, from 75 ms to 1.075 s; the description says the injection "does not usually start until about half a second into the recording", which does not hold for these files. Both the -100 pA and the +50 pA steps of the Figure 7 legends are present in every file.

The description names the three cells but not their animals. The superseded description kept in the same bucket (hc-9/old/) gives them:

In [5]:
old = hc9.old_description_text()
lines = [l for l in old.splitlines() if "_raw" in l]
for l in lines:
    print(l)
print("\nFigure 7 legend (paper): A, two O-LM cells, P13 and P14, from CXCR4-EGFP mice; B to D, an O-LM "
      "cell, P15, from the Sst-IRES-Cre;ChR2(H134R)-EYFP mouse.")
1May2012_P3_raw comes from a P13 CXCR4-EGFP mouse
20Sept2011_P2_raw comes P14 CXCR4-EGFP mouse
24October2012_C2_raw comes from a P15 Sst-IRES-Cre; ChR2(H134R)-EYFP mouse
- May2012_P3_raw.atf comes from a P13 CXCR4-EGFP mouse
- 20Sept2011_P2_raw.atf comes P14 CXCR4-EGFP mouse
- 24October2012_C2_raw.atf comes from a P15 Sst-IRES-Cre; ChR2(H134R)-EYFP mouse

Figure 7 legend (paper): A, two O-LM cells, P13 and P14, from CXCR4-EGFP mice; B to D, an O-LM cell, P15, from the Sst-IRES-Cre;ChR2(H134R)-EYFP mouse.

The ages and mouse lines of the old description match the Figure 7 legends: 1May2012_P3 (P13) and 20Sept2011_P2 (P14) are consistent with the two cells of Figure 7A, and 24October2012_C2 (P15, Sst-IRES-Cre;ChR2) with the cell of Figure 7B to D. This is consistent with the data set holding the three Figure 7 cells, as its description says, but the files themselves carry no age, line or panel. The old description also calls the files ABF; they are ATF text files.

Next, a check for recordings deposited more than once: every pair of sweeps, within and across files, is compared.

In [6]:
rows = []
keys = [(n, k) for n, c in cells.items() for k in range(c["v"].shape[1])]
for (n1, k1), (n2, k2) in itertools.combinations(keys, 2):
    a, b = cells[n1]["v"][:, k1], cells[n2]["v"][:, k2]
    rows.append({"a": f"{n1}#{k1+1}", "b": f"{n2}#{k2+1}", "same_cell": n1 == n2,
                 "identical": bool(np.array_equal(a, b)), "max_abs_diff_mV": float(np.abs(a - b).max()),
                 "r": float(np.corrcoef(a, b)[0, 1])})
pairs = pd.DataFrame(rows)
print(f"{len(pairs)} pairs of sweeps; identical pairs: {int(pairs.identical.sum())}")
print("Smallest maximum difference between two sweeps (mV):", round(pairs.max_abs_diff_mV.min(), 3))
print("\nAcross cells, the same step level: correlation of the voltage traces")
same_level = []
for lvl in levels:
    ks = {n: hc9.sweep_for(c, lvl) for n, c in cells.items()}
    for n1, n2 in itertools.combinations(cells, 2):
        a, b = cells[n1]["v"][:, ks[n1]], cells[n2]["v"][:, ks[n2]]
        same_level.append({"level_pA": lvl, "pair": f"{n1} / {n2}", "r": np.corrcoef(a, b)[0, 1],
                           "max_abs_diff_mV": np.abs(a - b).max()})
same_level = pd.DataFrame(same_level)
print(same_level.groupby("pair")[["r", "max_abs_diff_mV"]].agg(["min", "max"]).round(3).to_string())
990 pairs of sweeps; identical pairs: 0
Smallest maximum difference between two sweeps (mV): 2.704

Across cells, the same step level: correlation of the voltage traces
                                      r        max_abs_diff_mV         
                                    min    max             min      max
pair                                                                   
1May2012_P3 / 20Sept2011_P2      -0.006  0.990           3.058   96.112
1May2012_P3 / 24October2012_C2    0.126  0.993           3.296  104.083
20Sept2011_P2 / 24October2012_C2  0.230  0.993           2.704  101.483

No two sweeps are identical, within a file or across files, and the three files differ at every step level, so no recording is deposited twice. Each cell is one recording from one animal, so results below are given per cell.

The voltage must share the clock of the current channel. Two checks: the delay from the step onset in the current channel to the first departure of the voltage from its baseline, and the voltage change in the step window against the step level, with the window as the current channel gives it and shifted in time by 1.1 s (after the step).

In [7]:
rows = []
for name, c in cells.items():
    t, v = c["t"], c["v"]
    on, off = hc9.step_window(c)
    st = hc9.step_table(c)
    pre = t < on
    for k, lvl in enumerate(st.step_pA):
        if lvl == 0:
            continue
        base, sd = v[pre, k].mean(), v[pre, k].std(ddof=1)
        dev = np.flatnonzero((np.abs(v[:, k] - base) > 5 * sd) & (t >= on - 0.02))
        rows.append({"cell": name, "level_pA": lvl, "voltage_departs_ms_after_onset": (t[dev[0]] - on) * 1e3})
lag = pd.DataFrame(rows)
print("Delay from the current step's onset to the voltage leaving baseline by 5 SD (ms):")
print(lag.groupby("cell").voltage_departs_ms_after_onset.describe()[["count", "min", "50%", "max"]].round(2).to_string())

def window_change(c, shift):
    on, off = hc9.step_window(c)
    t, v = c["t"], c["v"]
    a, b = on + shift, min(off + shift, t[-1])
    w = (t >= a + 0.2 * (b - a)) & (t < b)
    pre = t < on
    # median over the window, so spikes on depolarizing steps do not dominate
    return np.median(v[w], axis=0) - v[pre].mean(axis=0)

rows = []
for name, c in cells.items():
    lv = hc9.step_table(c).step_pA.to_numpy()
    neg = lv < 0
    for label, shift in [("step window", 0.0), ("shifted by 1.1 s", 1.1 - 0.075)]:
        dv = window_change(c, shift)
        rho, p = stats.spearmanr(lv[neg], dv[neg])
        rows.append({"cell": name, "window": label, "spearman_rho_neg_steps": rho, "p": p,
                     "dV_at_-100pA_mV": dv[lv == -100][0], "dV_at_+250pA_mV": dv[lv == 250][0]})
print("\nVoltage change (median over the window's last 80%, minus baseline) against step level:")
print(pd.DataFrame(rows).round(4).to_string(index=False))
Delay from the current step's onset to the voltage leaving baseline by 5 SD (ms):
                  count   min   50%   max
cell                                     
1May2012_P3        14.0  0.25  1.52  4.00
20Sept2011_P2      14.0  0.35  0.55  3.35
24October2012_C2   14.0  0.10  0.75  2.55
Voltage change (median over the window's last 80%, minus baseline) against step level:
            cell           window  spearman_rho_neg_steps   p  dV_at_-100pA_mV  dV_at_+250pA_mV
     1May2012_P3      step window                     1.0 0.0         -24.4557          17.9590
     1May2012_P3 shifted by 1.1 s                    -0.2 0.8           1.9359          -3.1591
   20Sept2011_P2      step window                     1.0 0.0         -28.0433          17.4155
   20Sept2011_P2 shifted by 1.1 s                    -0.2 0.8           0.3807          -3.3608
24October2012_C2      step window                     1.0 0.0         -34.2732          24.1759
24October2012_C2 shifted by 1.1 s                    -0.8 0.2           4.9114          -5.0050

The voltage leaves its baseline 0.10 to 4.00 ms after the current step's onset in every sweep (medians 0.55 to 1.52 ms per cell), and in the step window the voltage change follows the hyperpolarizing step level monotonically (Spearman rho 1.0 in each cell), while in the window shifted past the step it does not (rho -0.8 to -0.2, p 0.2 to 0.8). The voltage and current channels share one clock. The three cells' pre-step potentials, about -63 mV, are given in the next section.

fig7a-olm-sag¶

O-LM interneurons recorded in CXCR4-EGFP slices show a sag in their voltage response to hyperpolarizing current pulses (pulses of -100 and +50 pA in Figure 7A), and the ChR2-expressing O-LM cell of Figure 7D has a similar firing pattern.

Criterion (claims file): the -100 pA step, judged under two readings. (a) The relaxation from the peak hyperpolarization during the step to the mean voltage over its last 100 ms is larger than twice the SD of the pre-step baseline. (b) The sag ratio (V_end - V_min) / (V_baseline - V_min) is at least 0.1. Every cell, the Figure 7D cell included, is judged alike. Reproduced when every cell shows a sag, partially when most do. The parameters are in hc9.py. Two details are standard choices: the pre-step baseline is all of the sweep before the step (0 to 75 ms), and V_min is the minimum of a 1 ms running mean during the step (the raw minimum is a variant; the noise checks below show why it is not the main measure).

In [8]:
on_off = {name: hc9.step_window(c) for name, c in cells.items()}
panel = {"1May2012_P3": "7A (P13, CXCR4-EGFP)", "20Sept2011_P2": "7A (P14, CXCR4-EGFP)",
         "24October2012_C2": "7D (P15, Sst-Cre;ChR2)"}
rows = []
for name, c in cells.items():
    on, off = on_off[name]
    k = hc9.sweep_for(c, hc9.HYPER_STEP_PA)
    m = hc9.sag_measures(c["t"], c["v"][:, k], on, off)
    rows.append({"cell": name, "panel (old description)": panel[name], **m})
sag = pd.DataFrame(rows)
sag["a_threshold_mV"] = hc9.SD_FACTOR * sag.sd_base
sag["reading_a"] = sag.sag_mV > sag.a_threshold_mV
sag["reading_b"] = sag.sag_ratio >= hc9.RATIO_THRESHOLD
print("-100 pA step, main measures (V in mV, t_min in ms from step onset):")
print(sag.round(3).to_string(index=False))
print(f"\nReading (a): {int(sag.reading_a.sum())} of {len(sag)} cells show a sag; "
      f"reading (b): {int(sag.reading_b.sum())} of {len(sag)}.")
print("Sag over the reading (a) threshold, per cell:", (sag.sag_mV / sag.a_threshold_mV).round(1).tolist())
print("Sag ratio over the reading (b) threshold, per cell:", (sag.sag_ratio / hc9.RATIO_THRESHOLD).round(2).tolist())
-100 pA step, main measures (V in mV, t_min in ms from step onset):
            cell panel (old description)  v_base  sd_base    v_min  t_min_ms   v_end  sag_mV  sag_over_sd  sag_ratio  a_threshold_mV  reading_a  reading_b
     1May2012_P3    7A (P13, CXCR4-EGFP) -63.575    0.279 -102.735    116.40 -85.196  17.539       62.909      0.448           0.558       True       True
   20Sept2011_P2    7A (P14, CXCR4-EGFP) -63.509    0.182  -98.974     88.95 -91.038   7.936       43.674      0.224           0.363       True       True
24October2012_C2  7D (P15, Sst-Cre;ChR2) -62.974    0.093 -104.546    125.65 -95.880   8.666       93.160      0.208           0.186       True       True

Reading (a): 3 of 3 cells show a sag; reading (b): 3 of 3.
Sag over the reading (a) threshold, per cell: [31.5, 21.8, 46.6]
Sag ratio over the reading (b) threshold, per cell: [4.48, 2.24, 2.08]

How the result changes with the choices the criterion leaves to a standard: V_min from the raw trace or a 1 ms running mean, the end window (50, 100 and 200 ms), the multiple of the baseline SD (2, 3, 5 and 10) and the ratio threshold (0.05 to 0.3).

In [9]:
rows = []
for name, c in cells.items():
    on, off = on_off[name]
    k = hc9.sweep_for(c, hc9.HYPER_STEP_PA)
    for smooth in (0.0, hc9.SMOOTH_S):
        for end in (0.05, 0.1, 0.2):
            m = hc9.sag_measures(c["t"], c["v"][:, k], on, off, end_window=end, smooth=smooth)
            rows.append({"cell": name, "smooth_ms": smooth * 1e3, "end_window_ms": end * 1e3,
                         "sag_mV": m["sag_mV"], "sag_over_sd": m["sag_over_sd"], "sag_ratio": m["sag_ratio"]})
var = pd.DataFrame(rows)
print(var.round(3).to_string(index=False))
print("\nSmallest sag over baseline SD in any variant:", round(var.sag_over_sd.min(), 1),
      "; smallest sag ratio:", round(var.sag_ratio.min(), 3))
for f in (2, 3, 5, 10):
    print(f"  reading (a) with {f} SD: cells with a sag in every variant:",
          int(var.groupby("cell").sag_over_sd.min().gt(f).sum()), "of 3")
for r in (0.05, 0.1, 0.15, 0.2, 0.25, 0.3):
    print(f"  reading (b) with ratio >= {r}: cells with a sag in every variant:",
          int(var.groupby("cell").sag_ratio.min().ge(r).sum()), "of 3")
            cell  smooth_ms  end_window_ms  sag_mV  sag_over_sd  sag_ratio
     1May2012_P3        0.0           50.0  17.874       64.110      0.455
     1May2012_P3        0.0          100.0  17.636       63.257      0.449
     1May2012_P3        0.0          200.0  17.126       61.428      0.436
     1May2012_P3        1.0           50.0  17.777       63.762      0.454
     1May2012_P3        1.0          100.0  17.539       62.909      0.448
     1May2012_P3        1.0          200.0  17.029       61.080      0.435
   20Sept2011_P2        0.0           50.0   7.979       43.907      0.225
   20Sept2011_P2        0.0          100.0   7.992       43.979      0.225
   20Sept2011_P2        0.0          200.0   7.759       42.697      0.218
   20Sept2011_P2        1.0           50.0   7.923       43.602      0.223
   20Sept2011_P2        1.0          100.0   7.936       43.674      0.224
   20Sept2011_P2        1.0          200.0   7.703       42.392      0.217
24October2012_C2        0.0           50.0   8.707       93.599      0.209
24October2012_C2        0.0          100.0   8.765       94.223      0.210
24October2012_C2        0.0          200.0   8.394       90.236      0.201
24October2012_C2        1.0           50.0   8.608       92.535      0.207
24October2012_C2        1.0          100.0   8.666       93.160      0.208
24October2012_C2        1.0          200.0   8.295       89.173      0.200

Smallest sag over baseline SD in any variant: 42.4 ; smallest sag ratio: 0.2
  reading (a) with 2 SD: cells with a sag in every variant: 3 of 3
  reading (a) with 3 SD: cells with a sag in every variant: 3 of 3
  reading (a) with 5 SD: cells with a sag in every variant: 3 of 3
  reading (a) with 10 SD: cells with a sag in every variant: 3 of 3
  reading (b) with ratio >= 0.05: cells with a sag in every variant: 3 of 3
  reading (b) with ratio >= 0.1: cells with a sag in every variant: 3 of 3
  reading (b) with ratio >= 0.15: cells with a sag in every variant: 3 of 3
  reading (b) with ratio >= 0.2: cells with a sag in every variant: 2 of 3
  reading (b) with ratio >= 0.25: cells with a sag in every variant: 1 of 3
  reading (b) with ratio >= 0.3: cells with a sag in every variant: 1 of 3

The same measures on the other hyperpolarizing steps, shown as a variant (the legends name the -100 pA step), and the firing on the +50 pA step, the legends' depolarizing pulse. The firing is described only: the criterion judges the similarity of the Figure 7D cell by the sag.

In [10]:
rows = []
for name, c in cells.items():
    on, off = on_off[name]
    for lvl in (-100, -75, -50, -25):
        k = hc9.sweep_for(c, lvl)
        m = hc9.sag_measures(c["t"], c["v"][:, k], on, off)
        rows.append({"cell": name, "level_pA": lvl, "v_min": m["v_min"], "sag_mV": m["sag_mV"],
                     "sag_over_sd": m["sag_over_sd"], "sag_ratio": m["sag_ratio"]})
print(pd.DataFrame(rows).round(3).to_string(index=False))

rows = []
for name, c in cells.items():
    on, off = on_off[name]
    for lvl in (hc9.DEPOL_STEP_PA, 100.0, 250.0):
        k = hc9.sweep_for(c, lvl)
        sp = hc9.spike_times(c["t"], c["v"][:, k])
        inside = sp[(sp >= on) & (sp < off)]
        isi = np.diff(inside) * 1e3
        rows.append({"cell": name, "level_pA": lvl, "spikes_in_step": len(inside),
                     "rate_Hz": len(inside) / (off - on),
                     "first_spike_ms": (inside[0] - on) * 1e3 if len(inside) else np.nan,
                     "first_isi_ms": isi[0] if len(isi) else np.nan,
                     "last_isi_ms": isi[-1] if len(isi) else np.nan,
                     "spikes_outside_step": len(sp) - len(inside)})
fire = pd.DataFrame(rows)
print("\nFiring on depolarizing steps (spikes are upward crossings of 0 mV):")
print(fire.round(1).to_string(index=False))
            cell  level_pA    v_min  sag_mV  sag_over_sd  sag_ratio
     1May2012_P3      -100 -102.735  17.539       62.909      0.448
     1May2012_P3       -75  -91.979  10.283       76.014      0.360
     1May2012_P3       -50  -84.011   5.646       37.544      0.267
     1May2012_P3       -25  -74.898   2.747       11.928      0.210
   20Sept2011_P2      -100  -98.974   7.936       43.674      0.224
   20Sept2011_P2       -75  -89.160   4.256        7.599      0.171
   20Sept2011_P2       -50  -81.042   3.176        4.006      0.174
   20Sept2011_P2       -25  -72.935   0.674        1.518      0.070
24October2012_C2      -100 -104.546   8.666       93.160      0.208
24October2012_C2       -75  -93.756   4.444       12.007      0.139
24October2012_C2       -50  -85.131   5.334       65.609      0.230
24October2012_C2       -25  -74.939   1.475        9.462      0.117

Firing on depolarizing steps (spikes are upward crossings of 0 mV):
            cell  level_pA  spikes_in_step  rate_Hz  first_spike_ms  first_isi_ms  last_isi_ms  spikes_outside_step
     1May2012_P3      50.0              15     15.0            33.4          56.6         74.7                    0
     1May2012_P3     100.0              27     27.0            18.2          31.0         43.8                    0
     1May2012_P3     250.0              48     48.0             7.9          15.6         22.9                    0
   20Sept2011_P2      50.0              15     15.0            19.2          58.5         93.8                    0
   20Sept2011_P2     100.0              35     35.0             6.8          23.4         31.6                    0
   20Sept2011_P2     250.0              83     83.0             3.0          11.0         10.2                    0
24October2012_C2      50.0              18     18.0            27.6          49.3         61.5                    0
24October2012_C2     100.0              36     36.0            13.7          24.4         29.6                    0
24October2012_C2     250.0              60     60.0             5.2          10.6         21.9                    0

How large a "sag" the measures find where there is none. Two checks of reading (a)'s false-positive level, with V_min from the raw trace and from the 1 ms running mean. First, each cell's 0 pA sweep, measured in the same step window: no current is injected, so any "sag" is noise. Second, a synthetic trace with no sag: a 40 mV exponential hyperpolarization (time constant 20 ms) in the same window, plus Gaussian noise of 0.2 mV SD, 100 draws from a fixed seed. The 0 pA sweep also gives a reference for the -100 pA sag in mV: the same measure with no current injected, and the SD of the voltage over the 1 s step window, which includes slow fluctuations that the 75 ms pre-step baseline does not. For reading (b) the 0 pA "sag" is also divided by the -100 pA hyperpolarization of the same cell (V_base minus V_min): this is the part of the -100 pA sag ratio that noise of the 0 pA size could make. The sag ratio of the 0 pA sweep itself is printed too, but its denominator, the 0 pA sweep's deflection below baseline, is noise as well, so it does not measure a false-positive level of reading (b).

In [11]:
rows = []
for name, c in cells.items():
    on, off = on_off[name]
    k = hc9.sweep_for(c, 0.0)
    for smooth in (0.0, hc9.SMOOTH_S):
        m = hc9.sag_measures(c["t"], c["v"][:, k], on, off, smooth=smooth)
        rows.append({"cell": name, "sweep": k + 1, "smooth_ms": smooth * 1e3, "sd_base": m["sd_base"],
                     "sag_mV": m["sag_mV"], "sag_over_sd": m["sag_over_sd"], "sag_ratio": m["sag_ratio"]})
zero = pd.DataFrame(rows)
print("0 pA sweep, step window as for the -100 pA step:")
print(zero.round(3).to_string(index=False))

rng = np.random.default_rng(20130327)
c0 = cells[hc9.CELLS[0]]
t = c0["t"]
on, off = on_off[hc9.CELLS[0]]
clean = np.where(t >= on, -40.0 * (1 - np.exp(-np.clip(t - on, 0, None) / 0.020)), 0.0) - 63.0
clean[t >= off] = -63.0
rows = []
for draw in range(100):
    v = clean + rng.normal(0, 0.2, len(t))
    for smooth in (0.0, hc9.SMOOTH_S):
        m = hc9.sag_measures(t, v, on, off, smooth=smooth)
        rows.append({"draw": draw, "smooth_ms": smooth * 1e3, "sag_over_sd": m["sag_over_sd"],
                     "sag_ratio": m["sag_ratio"]})
syn = pd.DataFrame(rows)
print("\nSynthetic trace without sag, 100 draws: sag over baseline SD and sag ratio")
print(syn.groupby("smooth_ms")[["sag_over_sd", "sag_ratio"]].describe(percentiles=[0.5])
      .round(3).to_string())
print("Fraction of draws passing reading (a) at 2 SD:",
      syn.groupby("smooth_ms").sag_over_sd.apply(lambda x: float((x > hc9.SD_FACTOR).mean())).to_dict(),
      "; passing reading (b):",
      syn.groupby("smooth_ms").sag_ratio.apply(lambda x: float((x >= hc9.RATIO_THRESHOLD).mean())).to_dict())
rows = []
for name, c in cells.items():
    on, off = on_off[name]
    k = hc9.sweep_for(c, 0.0)
    t, v = c["t"], c["v"][:, k]
    during = (t >= on) & (t < off)
    z = zero[(zero.cell == name) & (zero.smooth_ms == 1.0)].iloc[0]
    real = sag.set_index("cell").loc[name]
    rows.append({"cell": name, "spikes_0pA": len(hc9.spike_times(t, v)),
                 "sd_0pA_in_step_mV": v[during].std(ddof=1), "sag_0pA_mV": z.sag_mV,
                 "sag_0pA_over_sd": z.sag_over_sd, "sag_-100pA_mV": real.sag_mV,
                 "sag_-100pA_over_sag_0pA": real.sag_mV / z.sag_mV,
                 "sag_-100pA_over_2sd_0pA_in_step": real.sag_mV / (2 * v[during].std(ddof=1)),
                 "sag_ratio_0pA_own": z.sag_ratio,
                 "sag_0pA_over_hyper_-100pA": z.sag_mV / (real.v_base - real.v_min),
                 "sag_ratio_-100pA": real.sag_ratio})
ref = pd.DataFrame(rows)
print("\nThe -100 pA sag against the same measure on the 0 pA sweep (V_min from the 1 ms running mean):")
print(ref.round(3).to_string(index=False))
print("\nReading (b) on the 0 pA sweeps. sag_ratio_0pA_own divides the 0 pA 'sag' by the 0 pA sweep's own")
print("deflection below baseline, which is noise too, so it is not a sag ratio of a hyperpolarization;")
print("sag_0pA_over_hyper_-100pA divides it by the -100 pA hyperpolarization, the ratio noise alone would")
print("add to the -100 pA sag ratio. Threshold of reading (b):", hc9.RATIO_THRESHOLD)
print("  0 pA sweeps passing reading (b) on their own deflection:", int((ref.sag_ratio_0pA_own >= hc9.RATIO_THRESHOLD).sum()), "of 3")
print("  noise ratio against the -100 pA hyperpolarization:", ref["sag_0pA_over_hyper_-100pA"].round(3).tolist(),
      "; passing:", int((ref["sag_0pA_over_hyper_-100pA"] >= hc9.RATIO_THRESHOLD).sum()), "of 3")
0 pA sweep, step window as for the -100 pA step:
            cell  sweep  smooth_ms  sd_base  sag_mV  sag_over_sd  sag_ratio
     1May2012_P3      5        0.0    0.286   2.498        8.732      2.885
     1May2012_P3      5        1.0    0.286   2.393        8.364      3.146
   20Sept2011_P2      5        0.0    0.473   1.425        3.012      0.719
   20Sept2011_P2      5        1.0    0.473   1.348        2.848      0.707
24October2012_C2      5        0.0    0.095   2.023       21.258      4.092
24October2012_C2      5        1.0    0.095   1.916       20.136      4.943
Synthetic trace without sag, 100 draws: sag over baseline SD and sag ratio
          sag_over_sd                                    sag_ratio                                   
                count   mean    std    min    50%    max     count   mean    std    min    50%    max
smooth_ms                                                                                            
0.0             100.0  3.984  0.312  3.438  3.960  4.912     100.0  0.020  0.001  0.017  0.019  0.024
1.0             100.0  0.856  0.091  0.692  0.844  1.274     100.0  0.004  0.000  0.003  0.004  0.006
Fraction of draws passing reading (a) at 2 SD: {0.0: 1.0, 1.0: 0.0} ; passing reading (b): {0.0: 0.0, 1.0: 0.0}
The -100 pA sag against the same measure on the 0 pA sweep (V_min from the 1 ms running mean):
            cell  spikes_0pA  sd_0pA_in_step_mV  sag_0pA_mV  sag_0pA_over_sd  sag_-100pA_mV  sag_-100pA_over_sag_0pA  sag_-100pA_over_2sd_0pA_in_step  sag_ratio_0pA_own  sag_0pA_over_hyper_-100pA  sag_ratio_-100pA
     1May2012_P3           0              0.814       2.393            8.364         17.539                    7.330                           10.776              3.146                      0.061             0.448
   20Sept2011_P2           0              0.592       1.348            2.848          7.936                    5.889                            6.706              0.707                      0.038             0.224
24October2012_C2           0              1.027       1.916           20.136          8.666                    4.522                            4.220              4.943                      0.046             0.208

Reading (b) on the 0 pA sweeps. sag_ratio_0pA_own divides the 0 pA 'sag' by the 0 pA sweep's own
deflection below baseline, which is noise too, so it is not a sag ratio of a hyperpolarization;
sag_0pA_over_hyper_-100pA divides it by the -100 pA hyperpolarization, the ratio noise alone would
add to the -100 pA sag ratio. Threshold of reading (b): 0.1
  0 pA sweeps passing reading (b) on their own deflection: 3 of 3
  noise ratio against the -100 pA hyperpolarization: [0.061, 0.038, 0.046] ; passing: 0 of 3
In [12]:
fig, axes = plt.subplots(1, 3, figsize=(13, 3.6), sharey=True)
for ax, (name, c) in zip(axes, cells.items()):
    on, off = on_off[name]
    for lvl, col in ((hc9.HYPER_STEP_PA, "tab:blue"), (hc9.DEPOL_STEP_PA, "tab:red")):
        k = hc9.sweep_for(c, lvl)
        ax.plot(c["t"] * 1e3, c["v"][:, k], color=col, lw=0.6, label=f"{lvl:+.0f} pA")
    row = sag.set_index("cell").loc[name]
    ax.axhline(row.v_min, color="0.5", ls=":", lw=0.8)
    ax.axhline(row.v_end, color="0.5", ls="--", lw=0.8)
    ax.set_title(f"{name}\n{panel[name]}", fontsize=9)
    ax.set_xlabel("time (ms)")
axes[0].set_ylabel("membrane potential (mV)")
axes[0].legend(fontsize=8, loc="upper right")
plt.tight_layout()
plt.show()
yr = [(c["v"][:, hc9.sweep_for(c, l)].min(), c["v"][:, hc9.sweep_for(c, l)].max())
      for c in cells.values() for l in (hc9.HYPER_STEP_PA, hc9.DEPOL_STEP_PA)]
print("Plotted: full 0 to 2000 ms sweeps; voltage range over the plotted traces",
      round(min(a for a, _ in yr), 1), "to", round(max(b for _, b in yr), 1), "mV; nothing falls outside the axes.",
      "Dotted line V_min, dashed line V_end (mean of the last 100 ms of the -100 pA step).")
No description has been provided for this image
Plotted: full 0 to 2000 ms sweeps; voltage range over the plotted traces -104.6 to 36.8 mV; nothing falls outside the axes. Dotted line V_min, dashed line V_end (mean of the last 100 ms of the -100 pA step).

The first 300 ms of the -100 pA step in each cell, the part of the trace where the sag appears, with the pre-step baseline, V_min and V_end marked, for comparison with the traces printed in Figure 7A and 7D of the paper (the figure is not reproduced here).

In [13]:
fig, axes = plt.subplots(1, 3, figsize=(13, 3.4), sharey=True)
ranges = []
for ax, (name, c) in zip(axes, cells.items()):
    on, off = on_off[name]
    k = hc9.sweep_for(c, hc9.HYPER_STEP_PA)
    w = (c["t"] >= on - 0.05) & (c["t"] < on + 0.3)
    x, y = (c["t"][w] - on) * 1e3, c["v"][w, k]
    ranges.append((name, x.min(), x.max(), y.min(), y.max()))
    ax.plot(x, y, color="tab:blue", lw=0.7)
    row = sag.set_index("cell").loc[name]
    ax.axhline(row.v_base, color="0.3", lw=0.8, label="baseline")
    ax.axhline(row.v_min, color="0.5", ls=":", lw=0.8, label="V_min")
    ax.axhline(row.v_end, color="0.5", ls="--", lw=0.8, label="V_end")
    ax.axvline(row.t_min_ms, color="0.7", lw=0.6)
    ax.set_title(f"{name}\n{panel[name]}, -100 pA", fontsize=9)
    ax.set_xlabel("time from step onset (ms)")
axes[0].set_ylabel("membrane potential (mV)")
axes[0].legend(fontsize=8, loc="upper right")
plt.tight_layout()
plt.show()
for name, x0, x1, y0, y1 in ranges:
    print(f"{name}: plotted {x0:.0f} to {x1:.0f} ms from onset, {y0:.1f} to {y1:.1f} mV; nothing outside the axes")
No description has been provided for this image
1May2012_P3: plotted -50 to 300 ms from onset, -102.8 to -62.9 mV; nothing outside the axes
20Sept2011_P2: plotted -50 to 300 ms from onset, -99.0 to -63.1 mV; nothing outside the axes
24October2012_C2: plotted -50 to 300 ms from onset, -104.6 to -62.8 mV; nothing outside the axes
In [14]:
print(f"Total run time of the notebook: {time.time() - started:.0f} s")
Total run time of the notebook: 7 s

Result: reproduced, resting on the measure the claims file assumes. On the -100 pA step all three cells show a sag under both readings, with V_min from a 1 ms running mean. (a) The relaxation from the peak hyperpolarization to the last 100 ms is 17.5, 7.9 and 8.7 mV for 1May2012_P3, 20Sept2011_P2 and 24October2012_C2, 22 to 47 times the threshold of twice the pre-step SD (0.19 to 0.56 mV). (b) The sag ratios are 0.448, 0.224 and 0.208, all above 0.1. Both readings give 3 of 3, and so does every variant shown (raw V_min, end windows of 50 to 200 ms): the smallest sag is 42.4 baseline SDs and the smallest ratio 0.200, so the outcome holds up to a 10 SD threshold or a ratio threshold of 0.15; at 0.2, two of three cells pass in every variant.

The 2 SD threshold of reading (a) does not separate a sag from noise in these recordings. On each cell's 0 pA sweep, with no current injected, the same measure gives 1.3 to 2.4 mV, 2.8 to 20.1 pre-step SDs, so all three 0 pA sweeps pass it; the 75 ms pre-step baseline misses the slow fluctuations of the voltage over the 1 s step (SD 0.59 to 1.03 mV in the 0 pA step window). On a synthetic trace without sag the raw minimum passes reading (a) in 100 of 100 draws (median 3.96 SD) and the 1 ms running mean in none (median 0.84 SD); reading (b) passes in none either way. Against the 0 pA sweep the -100 pA sags are 4.5 to 7.3 times larger, and 4.2 to 10.8 times twice the SD of the 0 pA step window. For reading (b), the 0 pA "sag" is 0.038 to 0.061 of each cell's -100 pA hyperpolarization, below the 0.1 threshold, while the -100 pA sag ratios are 0.208 to 0.448; reading (b) gives 3 of 3. Applied to the 0 pA sweep alone, reading (b) passes in all three cells (0.71 to 4.94), because with no hyperpolarization the ratio divides noise by noise; that number is not a false-positive level of the ratio on a hyperpolarizing step.

The cell the old description places in Figure 7D (24October2012_C2, P15, Sst-Cre;ChR2) has a sag ratio of 0.208, close to the 0.224 of 20Sept2011_P2 and below the 0.448 of 1May2012_P3, the two Figure 7A cells; that the 7D cell shows a sag too is the similarity the criterion judges. At +50 pA the three cells fire 15, 15 and 18 spikes in the 1 s step, which is shown only as a description. The paper defines no sag measure, so the outcome rests on the measure the claims file chose.

Summary¶

Claim Outcome
fig7a-olm-sag reproduced (assumed measure of sag)
The other 25 claims (Figures 1 to 6 and 8 to 12, the light responses of Figure 7, the age analysis) blocked, data missing

How the analysis differs from the paper's. The paper shows the sag only in example traces and defines no measure; the two readings of the claims file are used, and both agree. V_min is taken from a 1 ms running mean, a standard choice. The 2 SD threshold of reading (a) is also passed by each cell's 0 pA sweep, so that reading alone does not separate a sag from noise in these recordings. On the 0 pA sweeps the noise "sag" is 0.038 to 0.061 of the -100 pA hyperpolarization, below reading (b)'s 0.1, and the -100 pA sags are 4.5 to 7.3 times the 0 pA value; on these references the outcome is the same. Which data set cell is in which panel of Figure 7 comes from the superseded description in the bucket's hc-9/old/ folder, not from the paper or the files; the criterion judges every cell alike, so the outcome does not depend on it.

What was found about the data set:

  • The data set holds only the current-step recordings of the three O-LM cells of Figure 7. None of the Cajal-Retzius recordings behind the paper's main findings are in it, nor the light-evoked responses of Figure 7C and 7D.
  • The description gives the current steps as -100:30:+250 pA; the files hold 15 steps from -100 to +250 pA in 25 pA increments.
  • The description says the injection usually starts about half a second into the sweep; in these files it starts at 75 ms and lasts 1 s.
  • The current description does not say which animal, age or mouse line each cell comes from, or which panel of Figure 7 shows it. The superseded description in hc-9/old/ does (P13 and P14 CXCR4-EGFP, P15 Sst-IRES-Cre;ChR2(H134R)-EYFP), and it matches the Figure 7 legends. That old description calls the files ABF; they are ATF text.
  • The zip's MD5 matches checksums.md5; no recording is deposited twice; the voltage and current channels share one clock.