Reproducing Hickman et al. (2026) from dandiset 000774

Paper. Hickman, Hughes, Sahai, Miles and Denman (2026), "Neural population dynamics of direct electrical stimulation of neocortex", Cell Reports 45, 117420, doi:10.1016/j.celrep.2026.117420.

Dataset. DANDI:000774, version 0.260520.1753 (doi:10.48324/dandi.000774/0.260520.1753), streamed from the DANDI API. This is the version the paper cites in its key resources table, and the only published version of the dandiset; the draft has the same ten assets.

Text the claims were written from. The NIH author manuscript (PMC13404931) with the STAR Methods and Document S1 (legends of Figures S1 to S12); Table S1 was not in the text read (claims/000774/10.1016-j.celrep.2026.117420.yaml). Method details the claims file does not give are taken from the paper and cited where used; where the paper is silent, from its code (github.com/denmanlab/estim_populations at c4aaa608), as the claims file does.

What the dandiset holds. Spike times, Kilosort templates, CCF coordinates and regions of the units and the table of electrical pulses, for ten acute recordings. The raw AP band, the LFP band and the behaviour are not in it, as the paper says ("Raw data (LFP band, AP band, and NiDAQ binary files) are not included in the NWB"), so the evoked-potential claims (Figures 1 to 3, S1, S2) and the behaviour claims (Figure 7, S10) are blocked. The three claims that compare units with the evoked-potential boundaries (fig4f, fig4i, fig5i) are tested with the boundaries the paper's code computed (intermediates/EP_distances.pkl at commit c4aaa608), the outside input their criteria name; the units' side of these tests comes from the dataset.

Claim	Tested here
fig1l-saturation-amplitude, figs1b-saturation-polarity	blocked, raw AP band not in the dandiset
fig2g-extent-increases, fig2g-extent-sublinear, fig2-extent-polarity, fig2h-probe-identity, figs2gh-magnitude-linear	blocked, raw AP band
fig3c-volume, fig3e-ventral-compression, fig3fg-anisotropy-strength, fig3h-com-displacement	blocked, raw AP band
fig4f-responders-within-ep, fig4i-proportion-within-ep, fig5i-fs-within-ep	tested, with the code's evoked-potential boundaries (outside input)
fig4d-latency-trough, figs3cd-gmm-window	tested
fig4h-extent-constant, fig4h-distance-decreases, fig4h-ks-distinct, fig4j-drp-increases, figs3e-early-windows	tested
fig5bc-amplitude-clusters, fig5ef-polarity-clusters	tested
fig5-fsrs-proportion, fig5j-fs-closer, figs6b-fsrs-early-firing	tested
fig6c-modulation-cortex, fig6d-hva-subcortical, fig6e-circuit-responders, fig6f-no-propagating-wave	tested
figs9-state-drp, fig6g-dimensionality-collapse	tested
fig7d to fig7j, figs10c, figs10e, no-learning-across-sessions (10 claims)	blocked, behaviour not in the dandiset

Every random step has its own seeded generator. The analysis takes about 13 minutes, most of it the ZETA tests (7 minutes on four cores) and the k-means stability check (2 minutes).

import os
for _v in ["OMP_NUM_THREADS", "OPENBLAS_NUM_THREADS", "MKL_NUM_THREADS"]:
    os.environ.setdefault(_v, "1")
import time
import warnings
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy import stats
from scipy.cluster.hierarchy import linkage, fcluster, cophenet
from scipy.spatial.distance import pdist
from sklearn.metrics import silhouette_score
from sklearn.mixture import GaussianMixture
import statsmodels.formula.api as smf
from statsmodels.stats.anova import anova_lm, AnovaRM
from statsmodels.stats.multicomp import pairwise_tukeyhsd
from statsmodels.stats.multitest import multipletests
import estim

warnings.filterwarnings("ignore", category=RuntimeWarning)
pd.set_option("display.width", 200)
pd.set_option("display.max_columns", 30)
T0 = time.time()

ASSETS = estim.list_assets()
DATA = {}
for a in ASSETS:
    DATA[a["session"]] = estim.load_session(a)
print(f"loaded {len(DATA)} files in {time.time() - T0:.0f} s")

loaded 10 files in 65 s

# distances from the dataset's coordinates, and, as a variant, from the paper's code coordinates
CODE_COORDS = {s: estim.code_coords(s) for s in estim.PAPER_SESSIONS}
for s in estim.PAPER_SESSIONS:
    u, tr, trains, meta = DATA[s]
    contact_xyz = tr[["neg_ap", "neg_dv", "neg_ml"]].iloc[0].to_numpy()
    u["distance"], u["norm_distance"] = estim.signed_distances(u, contact_xyz)
    u["distance_code"], u["norm_distance_code"] = estim.code_distances(
        u, CODE_COORDS[s], tr.contact_negative.iloc[0])
    # the raw CCF distance signed as the code's `distance` column: the code makes it negative from the
    # channel closest to the contact on, as norm_distance (figure4_figure5_figure6.ipynb cell 7,
    # `dists[cross_index:] = -np.abs(dists[cross_index:])`); checked on the code's direct_df.csv below
    u["sdistance"] = np.where(u.norm_distance < 0, -u.distance, u.distance)
    u["sdistance_code"] = np.where(u.norm_distance_code < 0, -u.distance_code, u.distance_code)
    u["keep"] = (u.group == "good") & ~u.unstable
    u["sub_cort"] = u.region.map(estim.classify_region)
    u["region_group"] = u.region.map(estim.map_region_to_group)

Data check

rows = []
for s, (u, tr, trains, meta) in DATA.items():
    probes = "".join(sorted(u.probe.unique()))
    rows.append(dict(session=s, subject=meta["subject"], probes=probes, units=len(u),
                     good=int((u.group == "good").sum()),
                     good_stable=int(((u.group == "good") & ~u.unstable).sum()),
                     has_ccf="region" in u, trials=len(tr),
                     spikes_s=f"{min(t.min() for t in trains if len(t)):.0f} to {max(t.max() for t in trains if len(t)):.0f}",
                     trials_s=f"{tr.start_time.min():.0f} to {tr.start_time.max():.0f}",
                     description=meta["description"][:60]))
files = pd.DataFrame(rows)
print(files.to_string(index=False))
paper = files[files.session.isin(estim.PAPER_SESSIONS)]
print("\nPaper recordings:", len(paper), "| probes:", sum(len(p) for p in paper.probes),
      "| probe A:", sum("A" in p for p in paper.probes), "B:", sum("B" in p for p in paper.probes),
      "C:", sum("C" in p for p in paper.probes))
print("Subjects of the paper recordings:", paper.subject.nunique(), sorted(paper.subject.unique()))
print("jlh48 and jlh49v1 (same subject id):",
      [(s, DATA[s][3]["session_id"]) for s in ["jlh48", "jlh49v1"]])

print("\nUnit groups:", pd.concat([d[0].group for d in DATA.values()]).value_counts().to_dict())
print("Unstable among good units, by recording:",
      {s: int((DATA[s][0].group.eq("good") & DATA[s][0].unstable).sum()) for s in estim.PAPER_SESSIONS})

print("\nTrials per condition (single pulses, 100 us per phase), the counts used (first 75):")
cond = {}
for s in estim.PAPER_SESSIONS:
    tr = DATA[s][1]
    cond[s] = {f"{p[0]}{a}": len(estim.condition_trials(tr, p, a, n=None)) for p in estim.POLARITIES for a in estim.AMPS}
print(pd.DataFrame(cond).T.to_string())
print("\nTrial table values: pulse_duration", sorted(pd.concat([d[1].pulse_duration for d in DATA.values()]).unique()),
      "| amplitude", sorted(pd.concat([d[1].amplitude for d in DATA.values()]).unique()),
      "| polarity", sorted(pd.concat([d[1].polarity for d in DATA.values()]).unique()))

print("\nRuns and contacts of the cathodal single-pulse conditions (all matching trials, and the trials used):")
rows = []
for s in estim.PAPER_SESSIONS:
    tr = DATA[s][1]
    for a in estim.AMPS:
        sub = estim.condition_trials(tr, "cathodal", a, n=None)
        if len(sub):
            used = estim.condition_trials(tr, "cathodal", a)
            rows.append(dict(session=s, amp=a, runs=sorted(sub.run.unique().tolist()), contact_negative=sorted(sub.contact_negative.unique().tolist()),
                             contact_positive=sorted(sub.contact_positive.unique().tolist()), trials=len(sub), used=len(used), used_runs=sorted(used.run.unique().tolist())))
print(pd.DataFrame(rows).to_string(index=False))
print("Cathodal contacts per recording (monopolar pulses):", {s: sorted(DATA[s][1][DATA[s][1].pol == "cathodal"].contact_negative.unique().tolist()) for s in estim.PAPER_SESSIONS})

allu = pd.concat([DATA[s][0] for s in estim.PAPER_SESSIONS])
print("\nUnit CCF ranges (um): AP %.0f to %.0f, DV %.0f to %.0f, ML %.0f to %.0f" % (
    allu.ap.min(), allu.ap.max(), allu.dv.min(), allu.dv.max(), allu.ml.min(), allu.ml.max()))
print("Distance to the first trial's cathodal contact (um): %.0f to %.0f, median %.0f" % (
    allu.distance.min(), allu.distance.max(), allu.distance.median()))
print("Signed distance (um): %.0f to %.0f" % (allu.norm_distance.min(), allu.norm_distance.max()))
print("Waveform duration (ms): %.2f to %.2f" % (allu.wf_duration.min() * 1e3, allu.wf_duration.max() * 1e3))
print("Regions of good units:", allu[allu.group == "good"].sub_cort.value_counts(dropna=False).to_dict())

    session subject probes  units  good  good_stable  has_ccf  trials       spikes_s       trials_s                                                  description
      jlh31   jlh31    ABC    785   409          405     True    2700   1860 to 8502   1992 to 8439                                                   craniotomy
      jlh32   jlh32    ABC    725   403          381     True    1923   2478 to 7534   2569 to 7475                                                   burr holes
      jlh33   jlh33     BC    258   175          175     True    1440   2585 to 6806   2669 to 6777                                                   burr holes
      jlh34   jlh34    ABC    589   263          263     True    2100  5026 to 10925  5133 to 10908                                                   craniotomy
      jlh39   jlh39    ABC    580   209          197     True     750   1559 to 3260   1732 to 4014                                             small craniotomy
      jlh40   jlh40     BC   1353   200          200     True     825 -4071 to -1540 -3951 to -1565                                             small craniotomy
      jlh41   jlh41    ABC   1542   216          207    False     525      7 to 3573    408 to 3113 small craniotomy, lots of noise issues. potentially related 
      jlh48   jlh48    ABC   1162   321          306     True    1200    565 to 3867    772 to 3833 burr holes. goal was to test the brains board but this nwb c
    jlh49v1   jlh48    ABC   1094   459          452     True     300   2024 to 6368   2525 to 3229 burr holes. goal was to test the brains board but this nwb c
jlh53concat   jlh53    ABC    857   391          384     True    1800  -1799 to 4947  -1136 to 4841 burr holes. typical experiment with a very slightly modified

Paper recordings: 9 | probes: 25 | probe A: 7 B: 9 C: 9
Subjects of the paper recordings: 8 ['jlh31', 'jlh32', 'jlh33', 'jlh34', 'jlh39', 'jlh40', 'jlh48', 'jlh53']
jlh48 and jlh49v1 (same subject id): [('jlh48', 'jlh48_2024-11-05_13-08-08'), ('jlh49v1', 'jlh49v1_2024-07-11_16-18-42')]

Unit groups: {'mua': 3453, 'good': 3046, 'nan': 1753, 'noise': 693}
Unstable among good units, by recording: {'jlh31': 4, 'jlh32': 22, 'jlh33': 0, 'jlh34': 0, 'jlh39': 12, 'jlh40': 0, 'jlh48': 15, 'jlh49v1': 7, 'jlh53concat': 7}

Trials per condition (single pulses, 100 us per phase), the counts used (first 75):
              c5  c25  c50  c100  a5  a25  a50  a100   b5  b25  b50  b100
jlh31         75  150  150   150  75  150  150   150  150  150  150   150
jlh32         75   75   75    75  75   75   75    75  150  150  150   150
jlh33         60   60   60    60  60   60    0    60  120  180  180   180
jlh34        150  150  150   150  75  150  150    75   75  150  150   150
jlh39          0    0    0     0   0    0    0     0  150  225  225   150
jlh40         75   75   75    75  75   75   75    75    0    0    0     0
jlh48         75   75   75    75  75   75   75    75  150  150  150   150
jlh49v1        0    0    0     0   0    0    0     0   75   75   75    75
jlh53concat   75   75   75    75  75   75   75    75  150  150  150   150

Trial table values: pulse_duration [np.int32(30), np.int32(50), np.int32(100), np.int32(200)] | amplitude [np.int32(-100), np.int32(-50), np.int32(-25), np.int32(-10), np.int32(-5), np.int32(5), np.int32(10), np.int32(25), np.int32(50), np.int32(100)] | polarity ['bipolar', 'monopolar']

Runs and contacts of the cathodal single-pulse conditions (all matching trials, and the trials used):
    session  amp     runs contact_negative contact_positive  trials  used used_runs
      jlh31    5     [16]             [16]              [0]      75    75      [16]
      jlh31   25 [18, 28]             [16]              [0]     150    75      [18]
      jlh31   50 [20, 30]             [16]              [0]     150    75      [20]
      jlh31  100 [22, 32]             [16]              [0]     150    75      [22]
      jlh32    5     [12]              [3]              [0]      75    75      [12]
      jlh32   25     [14]              [3]              [0]      75    75      [14]
      jlh32   50     [16]              [3]              [0]      75    75      [16]
      jlh32  100     [18]              [3]              [0]      75    75      [18]
      jlh33    5      [8]              [5]              [0]      60    60       [8]
      jlh33   25     [10]              [5]              [0]      60    60      [10]
      jlh33   50     [12]              [5]              [0]      60    60      [12]
      jlh33  100     [13]              [5]              [0]      60    60      [13]
      jlh34    5  [7, 15]             [15]              [0]     150    75       [7]
      jlh34   25  [9, 16]             [15]              [0]     150    75       [9]
      jlh34   50 [11, 19]             [15]              [0]     150    75      [11]
      jlh34  100 [13, 18]             [15]              [0]     150    75      [13]
      jlh40    5      [0]             [16]              [0]      75    75       [0]
      jlh40   25      [1]             [16]              [0]      75    75       [1]
      jlh40   50      [2]             [16]              [0]      75    75       [2]
      jlh40  100      [6]             [16]              [0]      75    75       [6]
      jlh48    5      [8]             [14]              [0]      75    75       [8]
      jlh48   25     [10]             [14]              [0]      75    75      [10]
      jlh48   50     [12]             [14]              [0]      75    75      [12]
      jlh48  100     [14]             [14]              [0]      75    75      [14]
jlh53concat    5      [0]             [16]              [0]      75    75       [0]
jlh53concat   25      [2]             [16]              [0]      75    75       [2]
jlh53concat   50      [4]             [16]              [0]      75    75       [4]
jlh53concat  100      [6]             [16]              [0]      75    75       [6]
Cathodal contacts per recording (monopolar pulses): {'jlh31': [16], 'jlh32': [3], 'jlh33': [5], 'jlh34': [15], 'jlh39': [], 'jlh40': [16], 'jlh48': [14], 'jlh49v1': [], 'jlh53concat': [16]}

Unit CCF ranges (um): AP 6859 to 9635, DV 368 to 2633, ML 1527 to 4376
Distance to the first trial's cathodal contact (um): 86 to 2281, median 664
Signed distance (um): -882 to 1912
Waveform duration (ms): 0.03 to 1.48
Regions of good units: {'Cortical': 1761, nan: 561, 'Subcortical': 508}

# The dataset's coordinates against the paper's code coordinates (outside input, a check only)
rows = []
for s in estim.PAPER_SESSIONS:
    u, tr, trains, meta = DATA[s]
    cc = CODE_COORDS[s]
    off = []
    for c in tr.contact_negative.unique():
        ds = tr.loc[tr.contact_negative == c, ["neg_ap", "neg_dv", "neg_ml"]].iloc[0].to_numpy()
        off.append(np.round(ds - cc["stim"][estim.NN_CONTACTS.index(int(c))]).astype(int).tolist())
    dd = []
    for i, r in u.iterrows():
        arr = cc.get("probe" + r.probe)
        if arr is None:
            continue
        dd.append(np.linalg.norm(arr[min(r.ch, len(arr) - 1)] - np.array([r.ap, r.dv, r.ml])))
    dd = np.array(dd)
    k = u.keep.to_numpy()
    rows.append(dict(session=s, contact_offset_um=off[0], contacts_same_offset=all(o == off[0] for o in off),
                     unit_offset_median_um=round(float(np.median(dd)), 1),
                     units_over_100um=int((dd > 100).sum()),
                     distance_diff_median_um=round(float(np.nanmedian(np.abs(u.distance - u.distance_code)[k])), 1)))
print(pd.DataFrame(rows).to_string(index=False))

    session contact_offset_um  contacts_same_offset  unit_offset_median_um  units_over_100um  distance_diff_median_um
      jlh31  [100, 150, -100]                  True                   24.0                 0                    129.9
      jlh32     [370, 0, -50]                  True                   24.5                 1                    144.0
      jlh33         [0, 0, 0]                  True                   24.0                 1                     18.3
      jlh34         [0, 0, 0]                  True                   24.0                 1                      3.6
      jlh39         [0, 0, 0]                  True                   24.8                 3                      4.8
      jlh40         [0, 0, 0]                  True                   24.0                 5                      5.0
      jlh48         [0, 0, 0]                  True                   25.3                 3                      7.2
    jlh49v1         [0, 0, 0]                  True                   24.0                 3                      9.3
jlh53concat         [0, 0, 0]                  True                   24.0                 7                     10.1

# Clock check: population spikes in 0 to 20 ms after 100 uA pulses (cathodal, or bipolar where
# there is none) against the 20 ms before, with the true onsets and with onsets moved 1 s later
# (half way between pulses, which are 2 s apart). Paired Wilcoxon over trials.
rows = []
for s in estim.PAPER_SESSIONS:
    u, tr, trains, meta = DATA[s]
    pol = "cathodal" if len(estim.condition_trials(tr, "cathodal", 100)) else "bipolar"
    pop = np.sort(np.concatenate([trains[i] for i in np.where(u.keep)[0]]))
    res = {"session": s, "polarity": pol}
    for label, shift in [("true", 0.0), ("shifted", 1.0)]:
        st = estim.condition_trials(tr, pol, 100).start_time.to_numpy() + shift
        post = np.searchsorted(pop, st + 0.02) - np.searchsorted(pop, st)
        pre = np.searchsorted(pop, st) - np.searchsorted(pop, st - 0.02)
        res[f"{label}_post_over_pre"] = round(post.mean() / pre.mean(), 2)
        res[f"{label}_p"] = stats.wilcoxon(post, pre).pvalue if (post != pre).any() else 1.0
    rows.append(res)
print(pd.DataFrame(rows).to_string(index=False))
# Population spike count per trial in 1 ms bins from 8 ms before to 2 ms after the logged onset, 100 uA
edges = np.round(np.arange(-0.008, 0.0021, 0.001), 4)
rows = {}
for s in estim.PAPER_SESSIONS:
    u, tr, trains, meta = DATA[s]
    pol = "cathodal" if len(estim.condition_trials(tr, "cathodal", 100)) else "bipolar"
    pop = np.sort(np.concatenate([trains[i] for i in np.where(u.keep)[0]]))
    rows[f"{s} ({pol})"] = np.round(estim.mean_counts(pop, estim.condition_trials(tr, pol, 100).start_time.to_numpy(), edges), 2)
print("\nPopulation spikes per trial in 1 ms bins starting at", (edges[:-1] * 1000).astype(int).tolist(), "ms:")
print(pd.DataFrame(rows, index=(edges[:-1] * 1000).astype(int)).T.to_string())

    session polarity  true_post_over_pre       true_p  shifted_post_over_pre  shifted_p
      jlh31 cathodal                4.55 5.261854e-14                   0.96   0.214150
      jlh32 cathodal                1.44 7.157224e-09                   1.04   0.442914
      jlh33 cathodal                1.96 3.180128e-10                   0.96   0.587237
      jlh34 cathodal                4.03 5.203523e-14                   1.06   0.487245
      jlh39  bipolar                2.00 1.557926e-10                   1.07   0.294861
      jlh40 cathodal                1.29 2.888356e-05                   0.98   0.483911
      jlh48 cathodal                1.55 8.214772e-11                   1.01   0.920567
    jlh49v1  bipolar                1.78 1.058268e-13                   0.99   0.978502
jlh53concat cathodal                0.58 3.560547e-11                   1.03   0.607904

Population spikes per trial in 1 ms bins starting at [-8, -7, -6, -5, -4, -3, -2, -1, 0, 1] ms:
                          -8    -7    -6    -5    -4    -3    -2    -1     0     1
jlh31 (cathodal)        1.73  1.87  1.44  1.41  1.11  0.99  0.61  1.48  4.97  2.87
jlh32 (cathodal)        1.60  1.75  1.49  1.28  0.80  0.41  0.07  0.45  8.05  4.96
jlh33 (cathodal)        0.98  0.95  0.97  0.72  0.63  0.50  0.12  0.15  0.70  1.32
jlh34 (cathodal)        0.81  1.15  1.09  0.79  0.64  0.36  0.16  0.36  6.43  5.11
jlh39 (bipolar)         0.84  0.69  0.69  0.69  0.44  0.24  0.11  0.41  1.07  0.13
jlh40 (cathodal)        1.12  1.08  1.12  0.95  0.67  0.37  0.12  0.13  0.55  0.17
jlh48 (cathodal)        1.68  1.28  1.32  0.80  0.59  0.35  0.12  0.36  1.47  0.72
jlh49v1 (bipolar)       2.03  2.12  1.71  1.71  1.12  0.37  0.20  1.77  3.73  0.27
jlh53concat (cathodal)  1.31  1.45  1.37  1.29  0.89  0.51  0.12  0.12  1.55  0.35

The checks against the description:

• Ten files; nine are the recordings the paper's code loads, and the tenth (jlh41, "lots of noise issues ... not much analyzed yet here") has no CCF coordinates or regions; it is left out. The nine give the paper's 25 probes (A 7, B 9, C 9). They come from eight subject ids: the file of session jlh49v1 is filed under subject jlh48, which the paper and its code treat as a separate recording and mouse; its age (P381D) is one day older than the jlh48 file recorded four months later.
• The electrodes table holds placeholders (location "none", constant coordinates), so channel CCF coordinates come from the units' own coordinates.
• Conditions are labelled "monopolar" or "bipolar" with signed amplitudes; cathodal and anodal are the sign. Cathodal single pulses are in seven of the nine recordings (jlh39 and jlh49v1 hold only bipolar pulses, jlh40 no bipolar), so the claims on cathodal pulses rest on seven recordings, where the paper's legends say nine. jlh33 has 60 trials per condition, where the Methods say 75; its anodal 50 uA single pulse is missing.
• Each recording has one cathodal contact for monopolar pulses. jlh31 (25 to 100 uA) and jlh34 (all four amplitudes) ran each cathodal condition twice at that contact, 150 trials; the first 75 are the first run, which the analysis uses, and the claims file's "layer 5" reading takes the same trials. The second run is shown as a variant in fig4j-drp-increases. In jlh53concat, a concatenated recording, the run numbers 0 to 7 are used both by monopolar and by bipolar runs.
• Trial onsets and spike times share a clock: population firing in the 20 ms after a 100 uA pulse differs from the 20 ms before it with the true onsets, and not with onsets moved 1 s later (table above). Two files (jlh40, jlh53concat) have negative times; this does not affect the alignment. In every recording the population fires less in the 3 ms before the logged onset, down to 0.07 to 0.20 spikes per trial from 2 to 1 ms before it (0.61 in jlh31), against 0.81 to 2.03 from 8 to 7 ms before (second table). The trials table does not say what happens there; the direct window and the nDRP baseline (7 to 5 ms before onset) lie outside it.
• Unit coordinates agree with the paper's code coordinates at the unit's channel to a median of about 24 um; units beyond the end of the code's channel arrays differ more. The stimulation contact coordinates of jlh31 and jlh32 differ from the code's by 206 and 373 um; the other seven agree exactly. Distances are computed from the dataset's coordinates, with the code's shown as a variant where they matter.

# Direct response table: nDRP per good stable unit and condition (Methods, DRP)
def drp_table(sessions=estim.PAPER_SESSIONS, win=0.002, base=(0.007, 0.005), polarities=estim.POLARITIES,
              contact_choice=None):
    rows = []
    for s in sessions:
        u, tr, trains, meta = DATA[s]
        idx = np.where(u.keep)[0]
        for pol in polarities:
            for a in estim.AMPS:
                contact = contact_choice.get((s, pol)) if contact_choice else None
                c = estim.condition_trials(tr, pol, a, contact=contact)
                if len(c) == 0:
                    continue
                st = c.start_time.to_numpy()
                for i in idx:
                    p, pb, n = estim.spike_prob(trains[i], st, win=win, base=base)
                    if np.isnan(p):
                        continue
                    rows.append((s, i, pol, a, p, pb, n))
    df = pd.DataFrame(rows, columns=["session", "unit", "pol", "amp", "drp", "base", "ndrp"])
    meta_cols = ["distance", "norm_distance", "distance_code", "norm_distance_code", "probe", "region",
                 "sub_cort", "wf_kmeans", "wf_threshold", "sdistance", "sdistance_code"]
    parts = []
    for s, g in df.groupby("session"):
        parts.append(g.join(DATA[s][0][meta_cols], on="unit"))
    return pd.concat(parts, ignore_index=True)

DRP = drp_table()
CAT = DRP[DRP.pol == "cathodal"]
RESP = CAT[CAT.ndrp > 0.15]
print("Unit and condition values:", len(DRP), "| cathodal:", len(CAT))
print(CAT.groupby("amp").agg(units=("ndrp", "size"), responders=("ndrp", lambda x: int((x > 0.15).sum())),
                             proportion=("ndrp", lambda x: (x > 0.15).mean())).round(4).to_string())
print("Responders by recording and amplitude:")
print(RESP.groupby(["session", "amp"]).size().unstack(fill_value=0).to_string())
print("Unit values per condition:", DRP.groupby(["pol", "amp"]).size().to_dict())

Unit and condition values: 26385 | cathodal: 8265
     units  responders  proportion
amp                               
5     2056          29      0.0141
25    2068          34      0.0164
50    2070          45      0.0217
100   2071          45      0.0217
Responders by recording and amplitude:
amp          5    25   50   100
session                        
jlh31          5   11   10   12
jlh32          2    5    9   11
jlh33          0    0    0    2
jlh34          2   10   18   14
jlh40         11    3    1    1
jlh48          5    3    5    2
jlh53concat    4    2    2    3
Unit values per condition: {('anodal', 5): 2082, ('anodal', 25): 2082, ('anodal', 50): 1901, ('anodal', 100): 2067, ('bipolar', 5): 2467, ('bipolar', 25): 2501, ('bipolar', 50): 2514, ('bipolar', 100): 2506, ('cathodal', 5): 2056, ('cathodal', 25): 2068, ('cathodal', 50): 2070, ('cathodal', 100): 2071}

# Check against the code's intermediate table direct_df.csv (outside input, a check only): the
# number of good stable units per condition, and the sign of its raw `distance` column
CDF = estim.code_direct_df().assign(pol=lambda x: x.polarity, amp=lambda x: x.amplitude.abs())
cmp_n = pd.concat([DRP.groupby(["pol", "amp"]).size().rename("this notebook"),
                   CDF.groupby(["pol", "amp"]).size().rename("code direct_df.csv")], axis=1)
print(cmp_n.T.to_string())
print("Equal in", int((cmp_n.iloc[:, 0] == cmp_n.iloc[:, 1]).sum()), "of", len(cmp_n), "conditions;",
      "cathodal equal in", int((cmp_n.loc["cathodal"].iloc[:, 0] == cmp_n.loc["cathodal"].iloc[:, 1]).sum()), "of 4")
print(f"Code's table: distance negative in {(CDF.distance < 0).mean():.3f} of rows; distance <= 0 exactly where "
      f"norm_distance <= 0 in {((CDF.distance <= 0) == (CDF.norm_distance <= 0)).mean():.4f} of rows")

pol                anodal                   bipolar                   cathodal                  
amp                   5     25    50    100     5     25    50    100      5     25    50    100
this notebook        2082  2082  1901  2067    2467  2501  2514  2506     2056  2068  2070  2071
code direct_df.csv   2082  2082  1901  2067    2465  2502  2515  2511     2056  2068  2070  2071
Equal in 8 of 12 conditions; cathodal equal in 4 of 4
Code's table: distance negative in 0.294 of rows; distance <= 0 exactly where norm_distance <= 0 in 1.0000 of rows

fig4d-latency-trough

First-spike latencies after a pulse are bimodal, with an early peak, a pause at a median of 1.95 ms across parameters and a broader peak from about 2.5 ms, which sets the 2 ms direct window.

First spike after each pulse (within 500 ms) for all good units, all trials of the condition's contact (the code's choose_stim_parameter without the 75-trial cut, figure4_figure5_figure6.ipynb cell 20), 0.1 ms bins over 0 to 15 ms, trough the lowest bin with centre in 0.5 to 5 ms (cell 22).

def latencies(pol, a, n=None, stable=False):
    out = []
    for s in estim.PAPER_SESSIONS:
        u, tr, trains, meta = DATA[s]
        c = estim.condition_trials(tr, pol, a, n=n)
        if len(c) == 0:
            continue
        st = c.start_time.to_numpy()
        sel = u.keep if stable else (u.group == "good")
        for i in np.where(sel)[0]:
            out.append(estim.first_latencies(trains[i], st))
    return np.concatenate(out)

EDGES = np.arange(0, 15.1, 0.1)
def trough(lat):
    cnt, e = np.histogram(lat, EDGES)
    ctr = e[:-1] + 0.05
    m = (ctr >= 0.5) & (ctr <= 5)
    j = np.where(m)[0][np.argmin(cnt[m])]
    return ctr[j], cnt[j], cnt[:j].max(), cnt[j + 1:].max(), ctr[j + 1:][np.argmax(cnt[j + 1:])]

LAT = {}
rows = []
for pol in estim.POLARITIES:
    for a in estim.AMPS:
        lat = latencies(pol, a)
        LAT[(pol, a)] = lat
        t, c, before, after, after_at = trough(lat)
        t75 = trough(latencies(pol, a, n=75, stable=True))[0]
        rows.append(dict(pol=pol, amp=a, n_latencies=int(((lat >= 0) & (lat <= 15)).sum()), trough_ms=round(t, 2),
                         trough_count=c, max_before=before, max_after=after, after_peak_ms=round(after_at, 2),
                         trough_first75_stable=round(t75, 2)))
lt = pd.DataFrame(rows)
print(lt.to_string(index=False))
print("\nMedian trough, cathodal (criterion):", lt[lt.pol == "cathodal"].trough_ms.median(),
      "| all 12 conditions:", lt.trough_ms.median(),
      "| first 75 trials, stable units, cathodal:", lt[lt.pol == "cathodal"].trough_first75_stable.median())
print("Conditions with a higher bin before and after the trough:", int(((lt.max_before > lt.trough_count) & (lt.max_after > lt.trough_count)).sum()), "of", len(lt))
fig, ax = plt.subplots(1, 4, figsize=(11, 2.5), sharey=False)
for k, a in enumerate(estim.AMPS):
    ax[k].hist(LAT[("cathodal", a)], EDGES, color="k")
    ax[k].axvline(lt[(lt.pol == "cathodal") & (lt.amp == a)].trough_ms.iloc[0], color="r")
    ax[k].set_title(f"cathodal {a} uA"); ax[k].set_xlabel("first spike latency (ms)")
plt.tight_layout(); plt.show()
print("All latencies plotted lie in 0 to 15 ms by the bins; latencies beyond 15 ms are not shown.")

     pol  amp  n_latencies  trough_ms  trough_count  max_before  max_after  after_peak_ms  trough_first75_stable
cathodal    5         9045       1.95            23         374         81           5.25                   1.95
cathodal   25        12672       1.75            25         323        111          11.85                   1.75
cathodal   50        16182       1.95            51         297        141           7.65                   1.45
cathodal  100        19842       1.95            76         453        165           3.75                   1.95
  anodal    5         8896       1.05            22         397         77           8.05                   1.05
  anodal   25        11317       1.45            17         459        106           8.25                   1.45
  anodal   50        12463       1.45            27         347        107           7.15                   2.25
  anodal  100        14693       1.45            34         376        129           7.25                   1.45
 bipolar    5        19032       1.25            58         927        169           5.95                   1.25
 bipolar   25        21245       1.35            42         974        169           9.05                   1.35
 bipolar   50        25778       1.25            57         842        228           5.85                   1.25
 bipolar  100        28580       1.45            67         632        249           8.85                   1.25

Median trough, cathodal (criterion): 1.95 | all 12 conditions: 1.45 | first 75 trials, stable units, cathodal: 1.85
Conditions with a higher bin before and after the trough: 12 of 12

All latencies plotted lie in 0 to 15 ms by the bins; latencies beyond 15 ms are not shown.

Result: reproduced. The trough of the first-spike latency histogram lies at 1.95, 1.75, 1.95 and 1.95 ms for 5 to 100 uA cathodal pulses, a median of 1.95 ms, with a higher bin before and after it in all 12 conditions. Across all 12 conditions the median is 1.45 ms (anodal and bipolar troughs at 1.05 to 1.45 ms); the criterion takes the cathodal amplitudes. The later cathodal peaks lie at 3.75 to 11.85 ms.

figs3cd-gmm-window

Gaussian mixtures of first-spike latencies give amplitude-specific direct windows (0.68, 1.14, 1.22 and 1.47 ms at 5, 25, 50 and 100 uA cathodal), and the proportion of direct responders does not differ from that of the 2 ms window.

Latencies of fig4d (0 to 15 ms), one-dimensional GMM (sklearn, maximum likelihood, 10 initialisations, seeded), readings (a) 4 components and (b) the number of components (1 to 4) with the lowest BIC. The window is where the posterior of the second-earliest component first equals that of the earliest (Methods, Gaussian mixture models). Responders per recording with each window and the 2 ms window, Wilcoxon signed-rank over recordings, FDR (Benjamini-Hochberg) over amplitudes.

rng_gmm = np.random.default_rng(3)
GRID = np.linspace(0, 15, 15001)
def gmm_window(x, k, seed):
    g = GaussianMixture(k, n_init=10, random_state=seed).fit(x[:, None])
    order = np.argsort(g.means_.ravel())
    if k < 2:
        return np.nan, g
    post = g.predict_proba(GRID[:, None])[:, order]
    m0 = g.means_.ravel()[order][0]
    after = GRID > m0
    cross = np.where(after & (post[:, 1] >= post[:, 0]))[0]
    return (GRID[cross[0]] if len(cross) else np.nan), g

WIN = {}
rows = []
for a in estim.AMPS:
    x = LAT[("cathodal", a)]
    x = x[(x >= 0) & (x <= 15)]
    seed = int(rng_gmm.integers(1e9))
    w4, g4 = gmm_window(x, 4, seed)
    bics = {k: GaussianMixture(k, n_init=10, random_state=seed).fit(x[:, None]).bic(x[:, None]) for k in range(1, 5)}
    kb = min(bics, key=bics.get)
    wb, _ = gmm_window(x, kb, seed)
    WIN[a] = {"a": w4, "b": wb}
    rows.append(dict(amp=a, n=len(x), window_k4_ms=round(w4, 3), means_k4=np.round(np.sort(g4.means_.ravel()), 2).tolist(),
                     bic_k=kb, window_bic_ms=round(wb, 3) if np.isfinite(wb) else np.nan))
print(pd.DataFrame(rows).to_string(index=False))
print("Expected 0.68, 1.14, 1.22, 1.47 ms (tolerance 10%)")

CAT_SESS = sorted(CAT.session.unique())
def prop_responders(win, base):
    out = {}
    for s in CAT_SESS:
        u, tr, trains, meta = DATA[s]
        idx = np.where(u.keep)[0]
        out[s] = {}
        for a in estim.AMPS:
            st = estim.condition_trials(tr, "cathodal", a).start_time.to_numpy()
            w = win[a] if isinstance(win, dict) else win
            v = np.array([estim.spike_prob(trains[i], st, win=w / 1000, base=base(w))[2] for i in idx])
            v = v[np.isfinite(v)]
            out[s][a] = np.mean(v > 0.15)
    return pd.DataFrame(out).T

P2 = prop_responders(2.0, lambda w: (0.007, 0.005))
for reading in ["a", "b"]:
    win = {a: WIN[a][reading] for a in estim.AMPS}
    if not all(np.isfinite(list(win.values()))):
        print(f"reading {reading}: no window at some amplitude, {win}")
        continue
    for bname, base in [("(i) baseline as long as the window", lambda w: (0.005 + w / 1000, 0.005)),
                        ("(ii) 2 ms baseline", lambda w: (0.007, 0.005))]:
        PG = prop_responders(win, base)
        p = [stats.wilcoxon(PG[a], P2[a]).pvalue if (PG[a] != P2[a]).any() else 1.0 for a in estim.AMPS]
        padj = multipletests(p, method="fdr_bh")[1]
        allp = stats.wilcoxon(PG.to_numpy().ravel(), P2.to_numpy().ravel()).pvalue
        print(f"\nGMM reading {reading}, {bname}: mean proportion GMM {PG.mean().round(4).to_dict()} | 2 ms {P2.mean().round(4).to_dict()}")
        print("   Wilcoxon p", np.round(p, 4).tolist(), "FDR", np.round(padj, 4).tolist(), "| pooled over amplitudes p", round(allp, 4))
print("\nRecordings:", len(CAT_SESS))

 amp     n  window_k4_ms                  means_k4  bic_k  window_bic_ms
   5  9045         0.683 [0.24, 4.34, 9.21, 13.17]      4          0.683
  25 12672         1.139 [0.42, 4.57, 8.88, 12.87]      4          1.139
  50 16182         1.221  [0.5, 4.52, 9.03, 13.03]      4          1.221
 100 19842         1.470 [0.68, 4.13, 8.93, 12.86]      4          1.470
Expected 0.68, 1.14, 1.22, 1.47 ms (tolerance 10%)

GMM reading a, (i) baseline as long as the window: mean proportion GMM {5: 0.0149, 25: 0.0145, 50: 0.015, 100: 0.0178} | 2 ms {5: 0.0156, 25: 0.0157, 50: 0.0208, 100: 0.0207}
   Wilcoxon p [1.0, 0.5, 0.0625, 0.125] FDR [1.0, 0.6667, 0.25, 0.25] | pooled over amplitudes p 0.0022

GMM reading a, (ii) 2 ms baseline: mean proportion GMM {5: 0.0149, 25: 0.0145, 50: 0.015, 100: 0.017} | 2 ms {5: 0.0156, 25: 0.0157, 50: 0.0208, 100: 0.0207}
   Wilcoxon p [1.0, 0.5, 0.0625, 0.125] FDR [1.0, 0.6667, 0.25, 0.25] | pooled over amplitudes p 0.0022

GMM reading b, (i) baseline as long as the window: mean proportion GMM {5: 0.0149, 25: 0.0145, 50: 0.015, 100: 0.0178} | 2 ms {5: 0.0156, 25: 0.0157, 50: 0.0208, 100: 0.0207}
   Wilcoxon p [1.0, 0.5, 0.0625, 0.125] FDR [1.0, 0.6667, 0.25, 0.25] | pooled over amplitudes p 0.0022

GMM reading b, (ii) 2 ms baseline: mean proportion GMM {5: 0.0149, 25: 0.0145, 50: 0.015, 100: 0.017} | 2 ms {5: 0.0156, 25: 0.0157, 50: 0.0208, 100: 0.0207}
   Wilcoxon p [1.0, 0.5, 0.0625, 0.125] FDR [1.0, 0.6667, 0.25, 0.25] | pooled over amplitudes p 0.0022

Recordings: 7

Result: reproduced. Four-component mixtures put the direct window at 0.683, 1.139, 1.221 and 1.470 ms for 5 to 100 uA, and BIC also chooses four components at every amplitude. Over the seven recordings with cathodal pulses, the proportion of responders with these windows does not differ from the 2 ms window at any amplitude (smallest FDR-adjusted p 0.25, under both baselines). Pooled over the four amplitudes, the GMM windows give fewer responders (Wilcoxon p = 0.0022), a test the legend does not name.

fig4f-responders-within-ep

Direct responders are sparse and mostly lie within the average boundaries of the evoked potential (33.3%, 74.1%, 88.2% and 100% at 5, 25, 50 and 100 uA), with more outside at 5 uA.

The evoked-potential boundaries need the raw AP band, which the dandiset does not hold. The criterion names the code's EP_distances.pkl as an outside input for the author to judge; it is used here (estim.code_ep_bounds: the code's per-probe boundaries, averaged over probes as in its cell 25), and everything on the units' side comes from the dataset. A cathodal responder is within when its signed distance lies between the mean superficial and deep boundaries at its amplitude (code cell 27). Readings (a) all responders, (b) RS responders (code cell 54), with the code's threshold classes and the Methods' k-means; distances from the dataset's coordinates and, as a variant, the code's.

EPB = estim.code_ep_bounds()["cathodal"]
print("Cathodal evoked-potential boundaries (outside input; mean superficial, mean deep, probes):",
      {a: (round(lo), round(hi), n) for a, (lo, hi, n) in EPB.items()})
def within_ep(d, col, bounds_amp=None):
    lo = d.amp.map(lambda a: EPB[bounds_amp or a][0])
    hi = d.amp.map(lambda a: EPB[bounds_amp or a][1])
    return (d[col] >= lo) & (d[col] <= hi)
rows = []
for col in ["norm_distance", "norm_distance_code"]:
    for reading, cls_col in [("a", None), ("b", "wf_threshold"), ("b", "wf_kmeans")]:
        d = RESP if cls_col is None else RESP[RESP[cls_col] == "rs"]
        w = within_ep(d, col).groupby(d.amp).agg(["mean", "size"])
        rows.append(dict(distances=col, reading=reading, classes=cls_col or "-",
                         **{f"{a} uA": f"{w.loc[a, 'mean']:.3f} (n={int(w.loc[a, 'size'])})" for a in estim.AMPS}))
print(pd.DataFrame(rows).to_string(index=False))
print("By recording, reading (b), threshold classes, dataset coordinates (fraction within, n):")
d = RESP[RESP.wf_threshold == "rs"]
print(within_ep(d, "norm_distance").groupby([d.session, d.amp]).agg(lambda x: f"{x.mean():.2f} ({len(x)})").unstack().to_string())
print("Expected 0.333, 0.741, 0.882, 1.0 (10%), 5 uA the lowest; the code prints 0.483, 0.765, 0.889, 1.0 for all responders")

Cathodal evoked-potential boundaries (outside input; mean superficial, mean deep, probes): {5: (-424, 510, 5), 25: (-725, 605, 31), 50: (-812, 747, 30), 100: (-750, 888, 19)}
         distances reading      classes         5 uA        25 uA        50 uA       100 uA
     norm_distance       a            - 0.552 (n=29) 0.794 (n=34) 0.889 (n=45) 1.000 (n=45)
     norm_distance       b wf_threshold 0.400 (n=15) 0.778 (n=27) 0.882 (n=34) 1.000 (n=30)
     norm_distance       b    wf_kmeans 0.520 (n=25) 0.774 (n=31) 0.900 (n=40) 1.000 (n=37)
norm_distance_code       a            - 0.483 (n=29) 0.765 (n=34) 0.889 (n=45) 1.000 (n=45)
norm_distance_code       b wf_threshold 0.333 (n=15) 0.741 (n=27) 0.882 (n=34) 1.000 (n=30)
norm_distance_code       b    wf_kmeans 0.480 (n=25) 0.742 (n=31) 0.900 (n=40) 1.000 (n=37)
By recording, reading (b), threshold classes, dataset coordinates (fraction within, n):
amp               5          25         50         100
session                                               
jlh31        1.00 (3)  1.00 (11)   1.00 (8)   1.00 (9)
jlh32        1.00 (2)   1.00 (3)   1.00 (5)   1.00 (7)
jlh33             NaN        NaN        NaN   1.00 (1)
jlh34        0.00 (2)   0.67 (9)  0.88 (17)  1.00 (12)
jlh40        0.00 (5)   0.00 (2)   0.00 (1)        NaN
jlh48        0.50 (2)   0.00 (1)   0.50 (2)        NaN
jlh53concat  0.00 (1)   1.00 (1)   1.00 (1)   1.00 (1)
Expected 0.333, 0.741, 0.882, 1.0 (10%), 5 uA the lowest; the code prints 0.483, 0.765, 0.889, 1.0 for all responders

Result: not reproduced, resting on an assumption. With the dataset's coordinates the fractions of RS responders within the evoked potential (reading b, the code's classes) are 0.400, 0.778, 0.882 and 1.000 at 5 to 100 uA: three are within 10% of 0.741, 0.882 and 1.0, the 5 uA fraction is the lowest, and 0.400 is outside 0.333 +- 0.033. Over all responders (reading a) they are 0.552, 0.794, 0.889 and 1.000. With the code's coordinates reading (b) gives 0.333, 0.741, 0.882 and 1.000 and reading (a) 0.483, 0.765, 0.889 and 1.000, the code's printed values; the difference at 5 uA is the contact coordinates of jlh31 and jlh32. Under the Methods' k-means classes RS fractions are 0.520, 0.774, 0.900 and 1.000. The boundaries are the code's (outside input). The 5 uA fraction is outside 0.333 +- 0.033 under every reading run (0.400, 0.552, 0.520), and the criterion has no partial clause.

fig4h-extent-constant

The spatial extent of direct responders stays relatively constant across amplitudes (1180 +- 559, 884 +- 652, 965 +- 451 and 522 +- 313 um), with no significant amplitude effect.

Per recording and cathodal amplitude, extent = |minimum| + maximum of the responders' signed distances (code cell 30); one-way ANOVA across amplitudes of the per-recording extents.

for col in ["norm_distance", "norm_distance_code"]:
    ext = RESP.groupby(["amp", "session"])[col].agg(lambda d: abs(d.min()) + d.max()).reset_index()
    groups = [g[col].to_numpy() for _, g in ext.groupby("amp")]
    F, p = stats.f_oneway(*groups)
    print(f"{col}: mean +- SD per amplitude", {a: f"{g[col].mean():.0f} +- {g[col].std(ddof=0):.0f} (n={len(g)})" for a, g in ext.groupby("amp")},
          f"| one-way ANOVA F = {F:.3f}, p = {p:.3f}")
    if col == "norm_distance":
        print(ext.pivot(index="session", columns="amp", values=col).round(0).to_string())
print("Expected 1180, 884, 965, 522 um (10%); paper F = 1.265, p = 0.312")

norm_distance: mean +- SD per amplitude {5: '1121 +- 626 (n=6)', 25: '825 +- 664 (n=6)', 50: '813 +- 683 (n=6)', 100: '452 +- 252 (n=7)'} | one-way ANOVA F = 1.250, p = 0.317
amp             5       25      50     100
session                                   
jlh31         415.0   366.0   252.0  788.0
jlh32         278.0   284.0   284.0  530.0
jlh33           NaN     NaN     NaN  112.0
jlh34        2131.0  1108.0  1148.0  827.0
jlh40        1186.0  2045.0  2048.0  314.0
jlh48        1411.0  1049.0  1049.0  310.0
jlh53concat  1307.0   100.0   100.0  280.0
norm_distance_code: mean +- SD per amplitude {5: '1179 +- 559 (n=6)', 25: '884 +- 652 (n=6)', 50: '868 +- 671 (n=6)', 100: '522 +- 313 (n=7)'} | one-way ANOVA F = 1.265, p = 0.312
Expected 1180, 884, 965, 522 um (10%); paper F = 1.265, p = 0.312

Result: reproduced. The per-recording extent of responders does not change with amplitude (one-way ANOVA F = 1.250, p = 0.317), with means of 1121, 825, 813 and 452 um; the last two are outside 10% of the paper's 965 and 522 um. With the paper's code coordinates the means are 1179, 884, 868 and 522 um and F = 1.265, p = 0.312: the paper's F and p and three of its four means (1180, 884 and 522 um); 868 um is just outside 10% of the paper's 965 um at 50 uA.

fig4h-distance-decreases

Higher amplitudes bring direct responders closer to the stimulation site (mean distance 426, 300, 137 and 146 um at 5, 25, 50 and 100 uA), driven by 5 against 50 and 5 against 100 uA.

for col in ["norm_distance", "norm_distance_code"]:
    g = RESP.groupby("amp")[col]
    F, p = stats.f_oneway(*[x.to_numpy() for _, x in g])
    print(f"{col}: mean (SD, n)", {a: f"{x.mean():.0f} ({x.std():.0f}, {len(x)})" for a, x in g}, f"| ANOVA F = {F:.2f}, p = {p:.4f}")
    tk = pairwise_tukeyhsd(RESP[col], RESP.amp)
    print(pd.DataFrame(tk.summary().data[1:], columns=tk.summary().data[0]).to_string(index=False))
    print("By recording (mean, n):", RESP.groupby(["session", "amp"])[col].agg(lambda x: f"{x.mean():.0f} ({len(x)})").unstack().to_string())
print("Expected 426, 300, 137, 146 um (10%); paper F = 4.15, p = 0.007; Tukey 5-50 p = 0.015, 5-100 p = 0.019")
print("Range of signed distances of responders: %.0f to %.0f um" % (RESP.norm_distance.min(), RESP.norm_distance.max()))

norm_distance: mean (SD, n) {5: '414 (499, 29)', 25: '264 (442, 34)', 50: '129 (389, 45)', 100: '142 (234, 45)'} | ANOVA F = 3.97, p = 0.0094

 group1  group2  meandiff  p-adj     lower    upper  reject
      5      25 -149.3185 0.4279 -404.4829 105.8458   False
      5      50 -284.1923 0.0134 -524.5725 -43.8120    True
      5     100 -271.6323 0.0199 -512.0125 -31.2520    True
     25      50 -134.8737 0.4236 -364.2540  94.5065   False
     25     100 -122.3138 0.5104 -351.6940 107.0665   False
     50     100   12.5600 0.9987 -200.2524 225.3724   False
By recording (mean, n): amp               5         25        50        100
session                                            
jlh31         173 (5)   98 (11)  -26 (10)   46 (12)
jlh32         137 (2)   102 (5)    42 (9)  128 (11)
jlh33             NaN       NaN       NaN  -224 (2)
jlh34        1065 (2)  538 (10)  387 (18)  322 (14)
jlh40        631 (11)  1023 (3)  1024 (1)   157 (1)
jlh48         -34 (5)  -267 (3)  -407 (5)    26 (2)
jlh53concat   490 (4)  -124 (2)  -124 (2)    54 (3)

norm_distance_code: mean (SD, n) {5: '426 (509, 29)', 25: '300 (443, 34)', 50: '137 (399, 45)', 100: '146 (263, 45)'} | ANOVA F = 4.15, p = 0.0074

 group1  group2  meandiff  p-adj     lower    upper  reject
      5      25 -125.8157 0.5986 -388.1590 136.5275   False
      5      50 -288.8159 0.0148 -535.9591 -41.6727    True
      5     100 -280.0835 0.0194 -527.2267 -32.9403    True
     25      50 -163.0001 0.2794 -398.8339  72.8336   False
     25     100 -154.2678 0.3275 -390.1015  81.5659   False
     50     100    8.7323 0.9996 -210.0674 227.5321   False
By recording (mean, n): amp               5         25        50        100
session                                            
jlh31         285 (5)  254 (11)   16 (10)  101 (12)
jlh32           3 (2)    31 (5)    36 (9)   78 (11)
jlh33             NaN       NaN       NaN  -267 (2)
jlh34        1051 (2)  525 (10)  392 (18)  340 (14)
jlh40        644 (11)  1042 (3)  1043 (1)   169 (1)
jlh48         -32 (5)  -270 (3)  -420 (5)     5 (2)
jlh53concat   478 (4)  -146 (2)  -146 (2)    36 (3)
Expected 426, 300, 137, 146 um (10%); paper F = 4.15, p = 0.007; Tukey 5-50 p = 0.015, 5-100 p = 0.019
Range of signed distances of responders: -868 to 1233 um

Result: partially reproduced. Mean signed distances are 414, 264, 129 and 142 um (ANOVA F = 3.97, p = 0.0094; Tukey 5 against 50 uA p = 0.0134, 5 against 100 uA p = 0.0199, no other pair below 0.05). The tests hold; the 25 uA mean (264 um) is outside 10% of 300 um. With the paper's code coordinates the means are 426, 300, 137 and 146 um and F = 4.15, p = 0.0074, the paper's values: the difference is the stimulation contact coordinates of jlh31 and jlh32 in the dataset.

fig4h-ks-distinct

5 and 100 uA recruit distinct spatial populations of direct responders; the other pairs the paper reports (5-25, 25-50, 25-100 and 50-100 uA) do not differ after correction.

for col in ["norm_distance", "norm_distance_code"]:
    pairs = [(5, 25), (5, 50), (5, 100), (25, 50), (25, 100), (50, 100)]
    p = [stats.ks_2samp(RESP[RESP.amp == x][col], RESP[RESP.amp == y][col]).pvalue for x, y in pairs]
    adj = estim.holm_sidak(p)
    D = [stats.ks_2samp(RESP[RESP.amp == x][col], RESP[RESP.amp == y][col]).statistic for x, y in pairs]
    print(col)
    print(pd.DataFrame(dict(pair=[f"{x}-{y}" for x, y in pairs], D=np.round(D, 3), p=np.round(p, 4), p_holm_sidak=np.round(adj, 4))).to_string(index=False))
print("Paper adjusted p: 5-25 0.566, 5-100 0.016, 25-50 0.566, 25-100 0.363, 50-100 0.083")

norm_distance
  pair     D      p  p_holm_sidak
  5-25 0.268 0.1700        0.5255
  5-50 0.336 0.0277        0.1309
 5-100 0.349 0.0193        0.1105
 25-50 0.166 0.5865        0.8290
25-100 0.198 0.3650        0.7439
50-100 0.133 0.8245        0.8290
norm_distance_code
  pair     D      p  p_holm_sidak
  5-25 0.247 0.2430        0.5662
  5-50 0.372 0.0105        0.0516
 5-100 0.416 0.0028        0.0164
 25-50 0.219 0.2577        0.5662
25-100 0.263 0.1064        0.3625
50-100 0.133 0.8245        0.8245
Paper adjusted p: 5-25 0.566, 5-100 0.016, 25-50 0.566, 25-100 0.363, 50-100 0.083

Result: not reproduced. With the dataset's coordinates, 5 against 100 uA has KS p = 0.0193 and Holm-Sidak adjusted p = 0.1105, so the 5 and 100 uA populations do not differ after correction; the other four pairs do not differ (adjusted p 0.5255 to 0.8290). With the paper's code coordinates the adjusted p are 0.5662 (5-25), 0.0164 (5-100), 0.5662 (25-50), 0.3625 (25-100) and 0.8245 (50-100), which reproduce the claim (the paper prints 0.083 for 50-100).

fig4i-proportion-within-ep

The proportion of direct responders among units within the evoked potential is small and increases with amplitude (0.014, 0.020, 0.031 and 0.030), so fewer than 5% of neurons are directly activated.

Cathodal good stable units with a value, within the code's mean evoked-potential boundaries (outside input, see fig4f-responders-within-ep); responder indicator nDRP > 0.15, pooled over recordings; OLS of the indicator on numeric amplitude (the "linear ANOVA"). Readings of the criterion: (a) the amplitude's own boundaries on the signed distance (norm_distance); (b) the 100 uA boundaries; (c) the amplitude's own boundaries on the raw CCF distance, signed as the code's distance column (negative on the superficial side of the probe's closest point). Dataset and code coordinates.

# (c unsigned): the raw CCF distance without the code's sign, a variant of reading (c)
FIG4I = {"a": ("norm_distance", None), "b": ("norm_distance", 100), "c": ("sdistance", None),
         "c unsigned": ("distance", None)}
for coords in ["", "_code"]:
    for reading, (col, ba) in FIG4I.items():
        d = CAT[within_ep(CAT, col + coords, ba)].assign(resp=lambda x: (x.ndrp > 0.15).astype(float))
        m = smf.ols("resp ~ amp", data=d).fit()
        t = anova_lm(m, typ=2)
        print(f"reading {reading}{' (code coordinates)' if coords else ''}:", {a: f"{x.mean():.4f} +- {x.std():.3f} (n={len(x)})" for a, x in d.groupby("amp").resp},
              f"| F = {t.loc['amp', 'F']:.2f}, p = {t.loc['amp', 'PR(>F)']:.4f}, slope {m.params['amp']:.2e}")
        if reading == "c" and not coords:
            print("   by recording:", d.groupby(["session", "amp"]).resp.mean().round(3).unstack().to_dict("index"))
print("Expected 0.014, 0.020, 0.031, 0.030 (10%), all below 0.05; paper F = 5.41, p = 0.020")

reading a: {5: '0.0124 +- 0.111 (n=1292)', 25: '0.0186 +- 0.135 (n=1449)', 50: '0.0242 +- 0.154 (n=1651)', 100: '0.0254 +- 0.157 (n=1773)'} | F = 5.87, p = 0.0154, slope 1.22e-04
reading b: {5: '0.0119 +- 0.108 (n=1766)', 25: '0.0152 +- 0.123 (n=1771)', 50: '0.0225 +- 0.148 (n=1775)', 100: '0.0254 +- 0.157 (n=1773)'} | F = 10.10, p = 0.0015, slope 1.44e-04
reading c: {5: '0.0133 +- 0.115 (n=678)', 25: '0.0208 +- 0.143 (n=1059)', 50: '0.0307 +- 0.173 (n=1271)', 100: '0.0301 +- 0.171 (n=1427)'} | F = 4.99, p = 0.0256, slope 1.49e-04
   by recording: {'jlh31': {5: 0.015, 25: 0.044, 50: 0.037, 100: 0.034}, 'jlh32': {5: 0.008, 25: 0.011, 50: 0.028, 100: 0.034}, 'jlh33': {5: 0.0, 25: 0.0, 50: 0.0, 100: 0.017}, 'jlh34': {5: 0.0, 25: 0.033, 50: 0.105, 100: 0.077}, 'jlh40': {5: 0.053, 25: 0.0, 50: 0.0, 100: 0.01}, 'jlh48': {5: 0.0, 25: 0.019, 50: 0.027, 100: 0.013}, 'jlh53concat': {5: 0.028, 25: 0.015, 50: 0.009, 100: 0.011}}
reading c unsigned: {5: '0.0129 +- 0.113 (n=778)', 25: '0.0221 +- 0.147 (n=997)', 50: '0.0312 +- 0.174 (n=1248)', 100: '0.0293 +- 0.169 (n=1467)'} | F = 4.55, p = 0.0330, slope 1.39e-04
reading a (code coordinates): {5: '0.0122 +- 0.110 (n=1147)', 25: '0.0197 +- 0.139 (n=1322)', 50: '0.0258 +- 0.159 (n=1549)', 100: '0.0259 +- 0.159 (n=1738)'} | F = 5.34, p = 0.0209, slope 1.23e-04
reading b (code coordinates): {5: '0.0121 +- 0.110 (n=1730)', 25: '0.0156 +- 0.124 (n=1735)', 50: '0.0230 +- 0.150 (n=1739)', 100: '0.0259 +- 0.159 (n=1738)'} | F = 10.08, p = 0.0015, slope 1.47e-04
reading c (code coordinates): {5: '0.0137 +- 0.116 (n=731)', 25: '0.0196 +- 0.139 (n=1073)', 50: '0.0306 +- 0.172 (n=1276)', 100: '0.0299 +- 0.170 (n=1440)'} | F = 5.41, p = 0.0201, slope 1.52e-04
reading c unsigned (code coordinates): {5: '0.0135 +- 0.116 (n=812)', 25: '0.0203 +- 0.141 (n=1033)', 50: '0.0307 +- 0.173 (n=1270)', 100: '0.0295 +- 0.169 (n=1457)'} | F = 5.26, p = 0.0219, slope 1.47e-04
Expected 0.014, 0.020, 0.031, 0.030 (10%), all below 0.05; paper F = 5.41, p = 0.020

Result: reproduced. Under reading (c), the amplitude's own boundaries on the raw signed distance, the proportions of responders among units within the evoked potential are 0.0133, 0.0208, 0.0307 and 0.0301 at 5 to 100 uA, all within 10% of 0.014, 0.020, 0.031 and 0.030 and below 0.05, and rise with amplitude (F = 4.99, p = 0.0256, slope positive). With the code's coordinates reading (c) gives F = 5.41, p = 0.0201, the paper's values. Readings (a) and (b) give the same trend (p = 0.0154 and 0.0015) but values below the paper's at 50 and 100 uA (0.0242 and 0.0254 under a). The boundaries are the code's (outside input). Reading (c) signs the raw distance as the code's distance column does; the code's direct_df.csv has distance <= 0 exactly where norm_distance <= 0 (check above). The unsigned raw distance gives 0.0129, 0.0221, 0.0312 and 0.0293 (F = 4.55, p = 0.0330): 0.0221 is just outside 0.020 +- 0.002, so the outcome rests on the code's sign, which the code applies.

fig4j-drp-increases

Among units responsive at any amplitude (99 units), direct response probability increases with amplitude (0.12, 0.14, 0.19 and 0.23).

PIV = CAT.pivot_table(index=["session", "unit"], columns="amp", values="ndrp")
resp_ids = RESP.set_index(["session", "unit"]).index.unique()
R99 = PIV.loc[resp_ids].dropna()
print("Responders at any cathodal amplitude:", len(resp_ids), "| with a value at all four:", len(R99))
print("Mean +- SD nDRP:", {a: f"{R99[a].mean():.3f} +- {R99[a].std():.3f}" for a in estim.AMPS})
long = R99.stack().rename("ndrp").reset_index()
m = smf.ols("ndrp ~ amp", data=long).fit()
print("(a) OLS on numeric amplitude:", anova_lm(m, typ=2).loc["amp", ["F", "PR(>F)"]].round(5).to_dict(), "slope", round(m.params["amp"], 5))
print("(b) one-way ANOVA, amplitude categorical: F = %.2f, p = %.2e" % stats.f_oneway(*[R99[a] for a in estim.AMPS]))
rm = AnovaRM(long.assign(uid=long.session + "_" + long.unit.astype(str)), "ndrp", "uid", within=["amp"]).fit()
print("(c) repeated-measures ANOVA:", rm.anova_table.round(5).to_dict("records"))
print("By recording, responders:", R99.groupby(level=0).size().to_dict())
print("Expected 99 units (10%), 0.12, 0.14, 0.19, 0.23; paper F = 12.78, p = 0.0004")
# Variant: jlh31 and jlh34 ran most cathodal conditions twice at the same contact; the second run in place of the first
rows = []
for s in ["jlh31", "jlh34"]:
    u, tr, trains, meta = DATA[s]
    for a in estim.AMPS:
        c = estim.condition_trials(tr, "cathodal", a, n=None)
        runs = c.run.unique()
        if len(runs) < 2:
            continue
        st = c[c.run == runs[1]].start_time.to_numpy()
        for i in np.where(u.keep)[0]:
            n = estim.spike_prob(trains[i], st)[2]
            if np.isfinite(n):
                rows.append((s, i, a, n))
SECOND = pd.DataFrame(rows, columns=["session", "unit", "amp", "ndrp"]).set_index(["session", "unit", "amp"]).ndrp
CAT_2 = CAT.set_index(["session", "unit", "amp"]).ndrp.copy()
CAT_2.loc[CAT_2.index.intersection(SECOND.index)] = SECOND
P2V = CAT_2.unstack("amp")
ids2 = P2V.index[(P2V > 0.15).any(axis=1)]
R2 = P2V.loc[ids2].dropna()
long2 = R2.stack().rename("ndrp").reset_index()
t2 = anova_lm(smf.ols("ndrp ~ amp", data=long2).fit(), typ=2)
print(f"\nSecond run for jlh31 and jlh34 ({len(SECOND)} unit values replaced): responders {len(ids2)}, with all four {len(R2)}, mean nDRP",
      {a: round(float(R2[a].mean()), 3) for a in estim.AMPS}, f"| OLS F = {t2.loc['amp', 'F']:.2f}, p = {t2.loc['amp', 'PR(>F)']:.5f}",
      "| by recording", R2.groupby(level=0).size().to_dict())

Responders at any cathodal amplitude: 103 | with a value at all four: 99
Mean +- SD nDRP: {5: '0.116 +- 0.201', 25: '0.142 +- 0.207', 50: '0.194 +- 0.253', 100: '0.233 +- 0.310'}
(a) OLS on numeric amplitude: {'F': 12.77902, 'PR(>F)': 0.00039} slope 0.00124
(b) one-way ANOVA, amplitude categorical: F = 4.44, p = 4.39e-03
(c) repeated-measures ANOVA: [{'F Value': 4.65709, 'Num DF': 3.0, 'Den DF': 294.0, 'Pr > F': 0.00339}]
By recording, responders: {'jlh31': 26, 'jlh32': 16, 'jlh33': 2, 'jlh34': 28, 'jlh40': 13, 'jlh48': 7, 'jlh53concat': 7}
Expected 99 units (10%), 0.12, 0.14, 0.19, 0.23; paper F = 12.78, p = 0.0004

Second run for jlh31 and jlh34 (2163 unit values replaced): responders 93, with all four 91, mean nDRP {5: 0.125, 25: 0.092, 50: 0.157, 100: 0.223} | OLS F = 12.41, p = 0.00048 | by recording {'jlh31': 31, 'jlh32': 16, 'jlh33': 2, 'jlh34': 15, 'jlh40': 13, 'jlh48': 7, 'jlh53concat': 7}

Result: reproduced. 99 units respond at some cathodal amplitude and have values at all four; their mean nDRP is 0.116, 0.142, 0.194 and 0.233, and the code's regression on numeric amplitude gives F = 12.779, p = 0.00039 (paper F = 12.78, p = 0.0004). The number of good stable units per cathodal condition (2056, 2068, 2070, 2071) equals that of the code's direct_df.csv (check above). A one-way ANOVA with amplitude categorical (F = 4.44, p = 0.0044) and a repeated-measures ANOVA (F = 4.66, p = 0.0034) agree. With the second run of jlh31 and jlh34 in place of the first, 91 units respond with values at all four amplitudes, the mean nDRP is 0.125, 0.092, 0.157 and 0.223 and F = 12.41, p = 0.00048: the rise holds, and the 25 uA mean (0.092) would be outside 10% of 0.14.

figs3e-early-windows

Broader windows (5 and 10 ms) give more responsive units than the 2 ms window (0.092 and 0.175 at 100 uA), and their amplitude dependence is stronger.

Proportion of good stable units with nDRP > 0.15 per recording, windows of 2, 5 and 10 ms; baseline readings (a) 2 ms ending 3 ms plus the window before onset, (b) as long as the window, ending 3 ms before onset; the 2 ms comparator is the main DRP. Paired Wilcoxon over recordings, FDR over the eight comparisons.

for reading, base in [("a", lambda w: (0.005 + w / 1000, 0.003 + w / 1000)), ("b", lambda w: (0.003 + w / 1000, 0.003))]:
    PW = {w: prop_responders(float(w), base) for w in (5, 10)}
    print(f"Reading {reading}: mean (SEM) proportion over {len(P2)} recordings")
    tab = pd.DataFrame({f"{w} ms": [f"{d[a].mean():.3f} ({d[a].sem():.3f})" for a in estim.AMPS] for w, d in [(2, P2)] + list(PW.items())}, index=estim.AMPS)
    print(tab.to_string())
    p, lab = [], []
    for w in (5, 10):
        for a in estim.AMPS:
            p.append(stats.wilcoxon(PW[w][a], P2[a]).pvalue if (PW[w][a] != P2[a]).any() else 1.0)
            lab.append(f"{w} ms vs 2 ms, {a} uA")
    adj = multipletests(p, method="fdr_bh")[1]
    print(pd.DataFrame(dict(comparison=lab, p=np.round(p, 4), p_fdr=np.round(adj, 4))).to_string(index=False))
    print("Increase 5 to 100 uA in mean proportion: 2 ms %.3f, 5 ms %.3f, 10 ms %.3f\n" % (
        P2[100].mean() - P2[5].mean(), PW[5][100].mean() - PW[5][5].mean(), PW[10][100].mean() - PW[10][5].mean()))
print("Expected at 100 uA: 5 ms 0.092, 10 ms 0.175 (10%)")

Reading a: mean (SEM) proportion over 7 recordings
              2 ms           5 ms          10 ms
5    0.016 (0.007)  0.017 (0.007)  0.033 (0.007)
25   0.016 (0.005)  0.024 (0.006)  0.047 (0.010)
50   0.021 (0.009)  0.040 (0.013)  0.094 (0.023)
100  0.021 (0.007)  0.058 (0.016)  0.106 (0.026)
           comparison      p  p_fdr
   5 ms vs 2 ms, 5 uA 0.5000 0.5000
  5 ms vs 2 ms, 25 uA 0.0156 0.0208
  5 ms vs 2 ms, 50 uA 0.0625 0.0714
 5 ms vs 2 ms, 100 uA 0.0156 0.0208
  10 ms vs 2 ms, 5 uA 0.0156 0.0208
 10 ms vs 2 ms, 25 uA 0.0156 0.0208
 10 ms vs 2 ms, 50 uA 0.0156 0.0208
10 ms vs 2 ms, 100 uA 0.0156 0.0208
Increase 5 to 100 uA in mean proportion: 2 ms 0.005, 5 ms 0.041, 10 ms 0.073

Reading b: mean (SEM) proportion over 7 recordings
              2 ms           5 ms          10 ms
5    0.016 (0.007)  0.014 (0.006)  0.016 (0.006)
25   0.016 (0.005)  0.019 (0.006)  0.032 (0.009)
50   0.021 (0.009)  0.037 (0.013)  0.068 (0.021)
100  0.021 (0.007)  0.057 (0.016)  0.091 (0.023)
           comparison      p  p_fdr
   5 ms vs 2 ms, 5 uA 0.5000 0.5714
  5 ms vs 2 ms, 25 uA 0.1250 0.1667
  5 ms vs 2 ms, 50 uA 0.0625 0.1000
 5 ms vs 2 ms, 100 uA 0.0156 0.0417
  10 ms vs 2 ms, 5 uA 1.0000 1.0000
 10 ms vs 2 ms, 25 uA 0.0156 0.0417
 10 ms vs 2 ms, 50 uA 0.0312 0.0625
10 ms vs 2 ms, 100 uA 0.0156 0.0417
Increase 5 to 100 uA in mean proportion: 2 ms 0.005, 5 ms 0.043, 10 ms 0.076

Expected at 100 uA: 5 ms 0.092, 10 ms 0.175 (10%)

Result: not reproduced. Broader windows give more responders and a stronger amplitude dependence: from 5 to 100 uA the mean proportion rises by 0.005 with 2 ms, 0.041 with 5 ms and 0.073 with 10 ms (reading a), and at 100 uA both broader windows exceed 2 ms (FDR p = 0.0208 under reading a, 0.0417 under b). The 100 uA proportions are 0.058 and 0.106 (reading a) or 0.057 and 0.091 (b), below 10% of the paper's 0.092 and 0.175. With seven recordings the smallest exact p is 0.0156. Both 100 uA proportions are outside tolerance under both readings, and the criterion has no partial clause.

fig5bc-amplitude-clusters

Clustering the direct responders by response probability across amplitude gives six clusters, broadly dose-dependent (58.6% of responders) or amplitude-selective, responding almost only at one amplitude (19.2% at 5 uA, 15.2% at 25 uA, 7.1% at 100 uA).

The 99 responders of fig4j; Ward linkage, maxclust 6; features (a) the four nDRP values, as the code, (b) each unit's vector divided by its maximum. A cluster is selective when its mean peaks at one amplitude and is below half the peak at every other (the claims file's rule). The fractions are also checked over every assignment of clusters to the four labels (the criterion's "some assignment"), and over those in which a cluster labelled selective for an amplitude peaks at that amplitude.

import itertools
EXP_AMP = {"dose": 0.586, 5: 0.192, 25: 0.152, 100: 0.071}

def cluster_summary(X, k=6):
    Z = linkage(X, method="ward")
    lab = fcluster(Z, k, criterion="maxclust")
    ccc = cophenet(Z, pdist(X))[0]
    sil = silhouette_score(X, lab)
    return Z, lab, ccc, sil

for reading, X in [("a", R99.to_numpy()), ("b", (R99.T / R99.max(axis=1)).T.to_numpy())]:
    Z, lab, ccc, sil = cluster_summary(X)
    rows = []
    for c in np.unique(lab):
        m = R99.to_numpy()[lab == c].mean(0)
        pk = estim.AMPS[int(np.argmax(m))]
        sel = all(m[j] < 0.5 * m.max() for j in range(4) if estim.AMPS[j] != pk)
        rows.append(dict(cluster=c, n=int((lab == c).sum()), frac=round((lab == c).mean(), 3),
                         mean_ndrp=np.round(m, 3).tolist(), peak_uA=pk, selective=sel))
    cl = pd.DataFrame(rows)
    fr = {"dose": cl[~cl.selective].frac.sum()}
    for a in estim.AMPS:
        fr[a] = cl[cl.selective & (cl.peak_uA == a)].frac.sum()
    print(f"Reading {reading}: CCC = {ccc:.3f}, silhouette = {sil:.3f}")
    print(cl.to_string(index=False))
    print("   fractions by the half-peak rule:", {k: round(v, 3) for k, v in fr.items()})
    ok = []
    for assign in itertools.product(["dose", 5, 25, 100], repeat=len(cl)):
        if any(l != "dose" and l != p for l, p in zip(assign, cl.peak_uA)):
            continue
        f = {key: sum(cl.frac[i] for i, l in enumerate(assign) if l == key) for key in EXP_AMP}
        if all(abs(f[key] - EXP_AMP[key]) <= 0.1 * EXP_AMP[key] for key in EXP_AMP):
            ok.append(assign)
    print("   assignments with selective labels at the cluster's peak, matching all four fractions within 10%:", len(ok), ok[:3])
    ok_any = []
    for assign in itertools.product(["dose", 5, 25, 100], repeat=len(cl)):
        f = {key: sum(cl.frac[i] for i, l in enumerate(assign) if l == key) for key in EXP_AMP}
        if all(abs(f[key] - EXP_AMP[key]) <= 0.1 * EXP_AMP[key] for key in EXP_AMP):
            ok_any.append(assign)
    print("   any assignment of clusters to the labels matching within 10%:", len(ok_any), ok_any[:3])
print("Expected CCC 0.702, silhouette 0.352 (10%)")

Reading a: CCC = 0.702, silhouette = 0.352
 cluster  n  frac                     mean_ndrp  peak_uA  selective
       1 15 0.152   [0.27, 0.474, 0.276, 0.088]       25      False
       2 19 0.192   [0.35, 0.034, 0.006, 0.008]        5       True
       3 40 0.404     [0.016, 0.1, 0.12, 0.117]       50      False
       4 13 0.131 [-0.004, 0.034, 0.106, 0.748]      100       True
       5  7 0.071 [-0.004, 0.021, 0.651, 0.343]       50      False
       6  5 0.051   [0.048, 0.339, 0.84, 0.955]      100      False
   fractions by the half-peak rule: {'dose': np.float64(0.678), 5: np.float64(0.192), 25: np.float64(0.0), 50: np.float64(0.0), 100: np.float64(0.131)}
   assignments with selective labels at the cluster's peak, matching all four fractions within 10%: 0 []

   any assignment of clusters to the labels matching within 10%: 2 [(25, 'dose', 'dose', 5, 100, 5), (25, 5, 'dose', 'dose', 100, 'dose')]
Reading b: CCC = 0.840, silhouette = 0.501
 cluster  n  frac                    mean_ndrp  peak_uA  selective
       1 28 0.283 [0.006, 0.034, 0.102, 0.506]      100       True
       2 10 0.101 [-0.008, 0.016, 0.44, 0.059]       50       True
       3 13 0.131  [0.026, 0.16, 0.514, 0.523]      100      False
       4 19 0.192  [0.35, 0.034, 0.006, 0.008]        5       True
       5  9 0.091 [0.453, 0.361, 0.255, 0.127]        5      False
       6 20 0.202  [0.018, 0.346, 0.142, 0.01]       25       True
   fractions by the half-peak rule: {'dose': np.float64(0.222), 5: np.float64(0.192), 25: np.float64(0.202), 50: np.float64(0.101), 100: np.float64(0.283)}
   assignments with selective labels at the cluster's peak, matching all four fractions within 10%: 0 []

   any assignment of clusters to the labels matching within 10%: 0 []
Expected CCC 0.702, silhouette 0.352 (10%)

Result: reproduced. On the code's features (reading a) the 99 responders give CCC = 0.702 and silhouette = 0.352, the paper's values, and clusters of 15, 19, 40, 13, 7 and 5 units. Labelling the clusters 25 uA, 5 uA, dose, dose, 100 uA and dose gives the paper's 15.2%, 19.2%, 58.6% and 7.1%, the assignment the criterion allows. By the half-peak rule the cluster of 7 peaks at 50 uA (mean nDRP 0.651 against 0.343 at 100 uA) and the cluster of 13 is the one selective for 100 uA, so no assignment that follows the clusters' peaks matches; reading (b) gives CCC 0.840 and silhouette 0.501 and no matching assignment.

fig5ef-polarity-clusters

Clustering responders by their response to cathodal, anodal and bipolar pulses shows a polarity bias: 69.1% are polarity-specific, preferring bipolar (27.7%), anodal (23.6%) or cathodal (17.7%) pulses.

As the code (cells 77 to 81): per polarity the maximum nDRP over amplitudes, missing polarities 0, divided by the unit's maximum; Ward, maxclust 6. Populations (a) units with nDRP > 0.15 in any polarity and amplitude, (b) the 99 cathodal responders. A cluster prefers a polarity when its mean is highest there and below half of that for both others (the claims file's rule); the fractions are also checked over every assignment of clusters to the three polarities or none, the criterion's "some assignment", as for fig5bc. Variant: for bipolar pulses, the contact pair whose cathodal contact lies nearest layer 5 units (the Methods' "in or near Layer 5"), in place of the first pair.

EXP_POL = {"bipolar": 0.277, "anodal": 0.236, "cathodal": 0.177}

def polarity_clusters(drp, ids):
    sub = drp.set_index(["session", "unit"]).loc[ids].reset_index()
    feat = sub.groupby(["session", "unit", "pol"]).ndrp.max().unstack(fill_value=0.0)
    feat = feat.reindex(columns=["bipolar", "cathodal", "anodal"], fill_value=0.0)
    feat = feat.div(feat.max(axis=1), axis=0).fillna(0.0)
    X = feat.to_numpy()
    Z, lab, ccc, sil = cluster_summary(X)
    rows = []
    for c in np.unique(lab):
        m = X[lab == c].mean(0)
        top = int(np.argmax(m))
        pref = feat.columns[top] if all(m[j] < 0.5 * m[top] for j in range(3) if j != top) else "none"
        rows.append(dict(cluster=c, n=int((lab == c).sum()), frac=round((lab == c).mean(), 3),
                         mean=np.round(m, 2).tolist(), prefers=pref))
    cl = pd.DataFrame(rows)
    cl["top"] = [feat.columns[int(np.argmax(m))] for m in cl["mean"]]
    fr = {p: round(cl[cl.prefers == p].frac.sum(), 3) for p in EXP_POL}
    return len(X), ccc, sil, cl, fr

def polarity_assignments(cl):
    """Every assignment of the clusters to the three polarities or none (the criterion's 'some
    assignment') whose three fractions are within 10%; marked when each labelled cluster's highest
    mean is at its label."""
    out = []
    for assign in itertools.product(["none"] + list(EXP_POL), repeat=len(cl)):
        f = {p: sum(cl.frac[i] for i, l in enumerate(assign) if l == p) for p in EXP_POL}
        if all(abs(f[p] - EXP_POL[p]) <= 0.1 * EXP_POL[p] for p in EXP_POL):
            out.append(dict(assignment=assign, at_top=all(l == "none" or l == t for l, t in zip(assign, cl.top)),
                            specific=round(sum(f.values()), 3), **{p: round(f[p], 3) for p in EXP_POL}))
    return pd.DataFrame(out)

all_ids = DRP[DRP.ndrp > 0.15].set_index(["session", "unit"]).index.unique()
# layer 5 variant: per recording, the bipolar cathodal contact nearest layer 5 units
L5 = {}
for s in estim.PAPER_SESSIONS:
    u, tr = DATA[s][0], DATA[s][1]
    regs = estim.contact_regions(u, tr)
    bip = tr[(tr.pol == "bipolar") & (tr.pulse_duration == 100) & (tr.pulse_number == 1)]
    cands = [int(c) for c in bip.contact_negative.unique()]
    l5 = [c for c in cands if "layer 5" in regs[c][0]]
    if l5:
        L5[(s, "bipolar")] = l5[0]
print("Bipolar contacts (first used by the code, layer 5 variant):",
      {s: (int(DATA[s][1][DATA[s][1].pol == 'bipolar'].contact_negative.iloc[0]) if (DATA[s][1].pol == 'bipolar').any() else None, L5.get((s, 'bipolar'))) for s in estim.PAPER_SESSIONS})
DRP_L5 = pd.concat([DRP[DRP.pol != "bipolar"], drp_table(polarities=["bipolar"], contact_choice=L5)], ignore_index=True)
for name, drp, ids in [("(a) any polarity", DRP, all_ids), ("(b) 99 cathodal responders", DRP, R99.index),
                       ("(a) any polarity, layer 5 bipolar contact", DRP_L5, DRP_L5[DRP_L5.ndrp > 0.15].set_index(["session", "unit"]).index.unique())]:
    n, ccc, sil, cl, fr = polarity_clusters(drp, ids)
    print(f"\n{name}: n = {n}, CCC = {ccc:.3f}, silhouette = {sil:.3f}, polarity-specific = {sum(fr.values()):.3f}, {fr}")
    print(cl.to_string(index=False))
    asg = polarity_assignments(cl)
    print(f"   assignments matching the three fractions within 10%: {len(asg)}, of which with each label at its cluster's highest mean: {int(asg.at_top.sum()) if len(asg) else 0}")
    if len(asg):
        print(asg.to_string(index=False))
print("\nExpected CCC 0.806, silhouette 0.519, polarity-specific 0.691, bipolar 0.277, anodal 0.236, cathodal 0.177 (10%)")

Bipolar contacts (first used by the code, layer 5 variant): {'jlh31': (16, 16), 'jlh32': (3, None), 'jlh33': (5, 5), 'jlh34': (15, 11), 'jlh39': (15, 15), 'jlh40': (None, None), 'jlh48': (14, 14), 'jlh49v1': (15, None), 'jlh53concat': (16, 16)}

(a) any polarity: n = 220, CCC = 0.773, silhouette = 0.538, polarity-specific = 0.623, {'bipolar': np.float64(0.205), 'anodal': np.float64(0.232), 'cathodal': np.float64(0.186)}
 cluster  n  frac               mean  prefers      top
       1 51 0.232   [0.21, 0.2, 1.0]   anodal   anodal
       2 45 0.205  [1.0, 0.04, 0.03]  bipolar  bipolar
       3 41 0.186  [0.16, 1.0, 0.13] cathodal cathodal
       4 27 0.123 [0.53, 0.98, 0.82]     none cathodal
       5 18 0.082  [0.98, 0.4, 0.79]     none  bipolar
       6 38 0.173  [0.98, 0.76, 0.2]     none  bipolar

   assignments matching the three fractions within 10%: 5, of which with each label at its cluster's highest mean: 2
                                      assignment  at_top  specific  bipolar  anodal  cathodal
(anodal, none, bipolar, none, bipolar, cathodal)   False     0.673    0.268   0.232     0.173
(anodal, none, cathodal, none, bipolar, bipolar)    True     0.673    0.255   0.232     0.186
(anodal, none, cathodal, bipolar, none, bipolar)   False     0.714    0.296   0.232     0.186
(anodal, bipolar, none, none, bipolar, cathodal)   False     0.692    0.287   0.232     0.173
(anodal, bipolar, cathodal, none, bipolar, none)    True     0.705    0.287   0.232     0.186

(b) 99 cathodal responders: n = 99, CCC = 0.711, silhouette = 0.393, polarity-specific = 0.434, {'bipolar': np.float64(0.0), 'anodal': np.float64(0.0), 'cathodal': np.float64(0.434)}
 cluster  n  frac               mean  prefers      top
       1 12 0.121  [0.48, 1.0, 0.15] cathodal cathodal
       2 31 0.313  [0.03, 1.0, 0.18] cathodal cathodal
       3 14 0.141  [0.98, 0.77, 0.1]     none  bipolar
       4 11 0.111 [0.23, 0.75, 0.97]     none   anodal
       5 13 0.131 [0.98, 0.86, 0.53]     none  bipolar
       6 18 0.182 [0.71, 0.98, 0.81]     none cathodal

   assignments matching the three fractions within 10%: 4, of which with each label at its cluster's highest mean: 0
                                        assignment  at_top  specific  bipolar  anodal  cathodal
(bipolar, none, bipolar, anodal, anodal, cathodal)   False     0.686    0.262   0.242     0.182
(bipolar, none, anodal, anodal, bipolar, cathodal)   False     0.686    0.252   0.252     0.182
(anodal, none, bipolar, bipolar, anodal, cathodal)   False     0.686    0.252   0.252     0.182
(anodal, none, bipolar, anodal, bipolar, cathodal)   False     0.686    0.272   0.232     0.182

(a) any polarity, layer 5 bipolar contact: n = 235, CCC = 0.773, silhouette = 0.524, polarity-specific = 0.669, {'bipolar': np.float64(0.26), 'anodal': np.float64(0.213), 'cathodal': np.float64(0.196)}
 cluster  n  frac               mean  prefers      top
       1 61 0.260  [1.0, 0.04, 0.06]  bipolar  bipolar
       2 50 0.213  [0.22, 0.22, 1.0]   anodal   anodal
       3 46 0.196 [0.08, 0.99, 0.29] cathodal cathodal
       4 35 0.149 [0.89, 0.79, 0.14]     none  bipolar
       5 14 0.060  [0.99, 0.26, 0.8]     none  bipolar
       6 29 0.123 [0.88, 0.91, 0.71]     none cathodal

   assignments matching the three fractions within 10%: 1, of which with each label at its cluster's highest mean: 0
                                       assignment  at_top  specific  bipolar  anodal  cathodal
(bipolar, anodal, none, none, cathodal, cathodal)   False     0.656     0.26   0.213     0.183

Expected CCC 0.806, silhouette 0.519, polarity-specific 0.691, bipolar 0.277, anodal 0.236, cathodal 0.177 (10%)

Result: reproduced. Under reading (a), the code's population, 220 units give CCC = 0.773 and silhouette = 0.538 (within 10% of 0.806 and 0.519). Of the assignments of the six clusters to the three polarities or none, five give all three fractions within 10%, and two of them label each cluster with its highest mean polarity; one gives bipolar 0.287, anodal 0.232 and cathodal 0.186 (polarity-specific 0.705), the criterion's "some assignment". It labels as bipolar both the cluster of 45 units (mean 1.0, 0.04, 0.03 for bipolar, cathodal, anodal) and one of 18 units (0.98, 0.40, 0.79) whose anodal response is 0.79 of its bipolar one. By the half-maximum rule only the first is bipolar-specific and the fractions are 0.205, 0.232 and 0.186 (total 0.623). With the bipolar contact nearest layer 5 (235 units) the half-maximum rule gives 0.260, 0.213 and 0.196. Among the 99 cathodal responders (reading b) CCC is 0.711 and silhouette 0.393, outside tolerance.

fig5-fsrs-proportion

Direct responders stay sparse for both FS and RS units across amplitudes, with a main effect of amplitude but no effect of waveform class or interaction.

Responder indicator (nDRP > 0.15) over FS or RS good stable units pooled across recordings, cathodal; two-way ANOVA (type II) of the indicator on amplitude and class. Classes under the Methods' k-means and the code's 0.4 ms threshold. Reading (b) of the criterion, units within the evoked potential, uses the code's evoked-potential boundaries (outside input, see fig4f-responders-within-ep) under the three readings of fig4i-proportion-within-ep. The paper's F = 3.04 with p = 0.028 fits 3 numerator degrees of freedom, amplitude categorical.

for cls_col in ["wf_kmeans", "wf_threshold"]:
    d = CAT[CAT[cls_col].isin(["fs", "rs"])].assign(resp=lambda x: (x.ndrp > 0.15).astype(float), cls=lambda x: x[cls_col])
    print(f"Classes by {cls_col}: units per class", d.groupby("cls").unit.size().div(4).round().to_dict())
    print(d.groupby(["cls", "amp"]).resp.agg(["mean", "std", "size"]).round(4).unstack(0).to_string())
    for amp_term in ["C(amp)", "amp"]:
        t = anova_lm(smf.ols(f"resp ~ {amp_term} * C(cls)", data=d).fit(), typ=2)
        print(f"   amplitude {'categorical' if 'C' in amp_term else 'numeric'}:",
              {k: f"F = {t.loc[k, 'F']:.2f}, p = {t.loc[k, 'PR(>F)']:.3f}" for k in t.index[:-1]})
print("Expected FS 0.032 (5 uA), 0.040 (100 uA); RS 0.013 (5 uA), 0.037 (100 uA) (10%); paper amplitude F = 3.04, p = 0.028; class p = 0.262; interaction p = 0.210")
print("\nReading (b): units within the evoked potential, under the three readings of fig4i-proportion-within-ep (boundaries: outside input)")
for reading, (col, ba) in FIG4I.items():
    din = CAT[within_ep(CAT, col, ba)]
    for cls_col in ["wf_kmeans", "wf_threshold"]:
        d = din[din[cls_col].isin(["fs", "rs"])].assign(resp=lambda x: (x.ndrp > 0.15).astype(float), cls=lambda x: x[cls_col])
        tab = d.groupby(["cls", "amp"]).resp.agg(lambda x: f"{x.mean():.4f} (n={len(x)})").unstack()
        print(f"({reading}) {col}{', 100 uA boundaries' if ba else ''}, {cls_col}:", tab.to_dict("index"))
        for amp_term in ["C(amp)", "amp"]:
            t = anova_lm(smf.ols(f"resp ~ {amp_term} * C(cls)", data=d).fit(), typ=2)
            print(f"   amplitude {'categorical' if 'C' in amp_term else 'numeric'}:",
                  {k: f"F = {t.loc[k, 'F']:.2f}, p = {t.loc[k, 'PR(>F)']:.3f}" for k in t.index[:-1]})

Classes by wf_kmeans: units per class {'fs': 332.0, 'rs': 1637.0}
       mean             std         size      
cls      fs      rs      fs      rs   fs    rs
amp                                           
5    0.0061  0.0153  0.0780  0.1228  328  1633
25   0.0060  0.0189  0.0773  0.1364  334  1636
50   0.0120  0.0244  0.1091  0.1544  333  1639
100  0.0211  0.0225  0.1441  0.1485  331  1641
   amplitude categorical: {'C(amp)': 'F = 1.93, p = 0.122', 'C(cls)': 'F = 4.85, p = 0.028', 'C(amp):C(cls)': 'F = 0.42, p = 0.736'}
   amplitude numeric: {'C(cls)': 'F = 4.84, p = 0.028', 'amp': 'F = 4.33, p = 0.037', 'amp:C(cls)': 'F = 0.71, p = 0.400'}
Classes by wf_threshold: units per class {'fs': 616.0, 'rs': 1352.0}
       mean             std         size      
cls      fs      rs      fs      rs   fs    rs
amp                                           
5    0.0196  0.0111  0.1387  0.1049  613  1348
25   0.0097  0.0200  0.0981  0.1400  619  1351
50   0.0162  0.0251  0.1263  0.1565  618  1354
100  0.0227  0.0221  0.1492  0.1471  616  1356
   amplitude categorical: {'C(amp)': 'F = 1.93, p = 0.123', 'C(cls)': 'F = 0.60, p = 0.437', 'C(amp):C(cls)': 'F = 1.77, p = 0.150'}
   amplitude numeric: {'C(cls)': 'F = 0.60, p = 0.438', 'amp': 'F = 4.33, p = 0.037', 'amp:C(cls)': 'F = 0.09, p = 0.764'}
Expected FS 0.032 (5 uA), 0.040 (100 uA); RS 0.013 (5 uA), 0.037 (100 uA) (10%); paper amplitude F = 3.04, p = 0.028; class p = 0.262; interaction p = 0.210

Reading (b): units within the evoked potential, under the three readings of fig4i-proportion-within-ep (boundaries: outside input)

(a) norm_distance, wf_kmeans: {'fs': {5: '0.0098 (n=205)', 25: '0.0087 (n=230)', 50: '0.0156 (n=256)', 100: '0.0253 (n=277)'}, 'rs': {5: '0.0126 (n=1035)', 25: '0.0207 (n=1160)', 50: '0.0270 (n=1331)', 100: '0.0259 (n=1426)'}}

   amplitude categorical: {'C(amp)': 'F = 2.79, p = 0.039', 'C(cls)': 'F = 1.74, p = 0.187', 'C(amp):C(cls)': 'F = 0.34, p = 0.795'}
   amplitude numeric: {'C(cls)': 'F = 1.75, p = 0.186', 'amp': 'F = 6.18, p = 0.013', 'amp:C(cls)': 'F = 0.20, p = 0.655'}
(a) norm_distance, wf_threshold: {'fs': {5: '0.0233 (n=387)', 25: '0.0115 (n=434)', 50: '0.0209 (n=479)', 100: '0.0278 (n=504)'}, 'rs': {5: '0.0070 (n=853)', 25: '0.0220 (n=956)', 50: '0.0271 (n=1108)', 100: '0.0250 (n=1199)'}}
   amplitude categorical: {'C(amp)': 'F = 2.80, p = 0.039', 'C(cls)': 'F = 0.00, p = 0.983', 'C(amp):C(cls)': 'F = 1.91, p = 0.126'}
   amplitude numeric: {'C(cls)': 'F = 0.00, p = 0.985', 'amp': 'F = 6.20, p = 0.013', 'amp:C(cls)': 'F = 0.13, p = 0.723'}
(b) norm_distance, 100 uA boundaries, wf_kmeans: {'fs': {5: '0.0072 (n=279)', 25: '0.0071 (n=280)', 50: '0.0142 (n=281)', 100: '0.0253 (n=277)'}, 'rs': {5: '0.0120 (n=1419)', 25: '0.0169 (n=1421)', 50: '0.0253 (n=1423)', 100: '0.0259 (n=1426)'}}

   amplitude categorical: {'C(amp)': 'F = 4.32, p = 0.005', 'C(cls)': 'F = 2.18, p = 0.139', 'C(amp):C(cls)': 'F = 0.28, p = 0.836'}
   amplitude numeric: {'C(cls)': 'F = 2.18, p = 0.140', 'amp': 'F = 11.16, p = 0.001', 'amp:C(cls)': 'F = 0.22, p = 0.642'}
(b) norm_distance, 100 uA boundaries, wf_threshold: {'fs': {5: '0.0237 (n=506)', 25: '0.0099 (n=506)', 50: '0.0196 (n=509)', 100: '0.0278 (n=504)'}, 'rs': {5: '0.0059 (n=1192)', 25: '0.0176 (n=1195)', 50: '0.0251 (n=1195)', 100: '0.0250 (n=1199)'}}

   amplitude categorical: {'C(amp)': 'F = 4.33, p = 0.005', 'C(cls)': 'F = 0.26, p = 0.608', 'C(amp):C(cls)': 'F = 2.56, p = 0.053'}
   amplitude numeric: {'C(cls)': 'F = 0.26, p = 0.607', 'amp': 'F = 11.18, p = 0.001', 'amp:C(cls)': 'F = 0.72, p = 0.397'}
(c) sdistance, wf_kmeans: {'fs': {5: '0.0085 (n=117)', 25: '0.0110 (n=181)', 50: '0.0189 (n=212)', 100: '0.0304 (n=230)'}, 'rs': {5: '0.0150 (n=533)', 25: '0.0227 (n=836)', 50: '0.0348 (n=1007)', 100: '0.0307 (n=1140)'}}
   amplitude categorical: {'C(amp)': 'F = 2.58, p = 0.052', 'C(cls)': 'F = 1.74, p = 0.188', 'C(amp):C(cls)': 'F = 0.32, p = 0.810'}
   amplitude numeric: {'C(cls)': 'F = 1.74, p = 0.187', 'amp': 'F = 4.77, p = 0.029', 'amp:C(cls)': 'F = 0.36, p = 0.547'}
(c) sdistance, wf_threshold: {'fs': {5: '0.0275 (n=218)', 25: '0.0123 (n=325)', 50: '0.0265 (n=377)', 100: '0.0324 (n=401)'}, 'rs': {5: '0.0069 (n=432)', 25: '0.0246 (n=692)', 50: '0.0344 (n=842)', 100: '0.0299 (n=969)'}}

   amplitude categorical: {'C(amp)': 'F = 2.60, p = 0.051', 'C(cls)': 'F = 0.05, p = 0.821', 'C(amp):C(cls)': 'F = 1.46, p = 0.223'}
   amplitude numeric: {'C(cls)': 'F = 0.05, p = 0.815', 'amp': 'F = 4.80, p = 0.028', 'amp:C(cls)': 'F = 0.02, p = 0.899'}
(c unsigned) distance, wf_kmeans: {'fs': {5: '0.0076 (n=131)', 25: '0.0115 (n=174)', 50: '0.0190 (n=211)', 100: '0.0299 (n=234)'}, 'rs': {5: '0.0145 (n=619)', 25: '0.0243 (n=781)', 50: '0.0355 (n=985)', 100: '0.0298 (n=1176)'}}
   amplitude categorical: {'C(amp)': 'F = 2.73, p = 0.042', 'C(cls)': 'F = 1.88, p = 0.171', 'C(amp):C(cls)': 'F = 0.38, p = 0.767'}

   amplitude numeric: {'C(cls)': 'F = 1.84, p = 0.175', 'amp': 'F = 4.33, p = 0.037', 'amp:C(cls)': 'F = 0.44, p = 0.507'}
(c unsigned) distance, wf_threshold: {'fs': {5: '0.0282 (n=248)', 25: '0.0129 (n=309)', 50: '0.0269 (n=372)', 100: '0.0316 (n=411)'}, 'rs': {5: '0.0060 (n=502)', 25: '0.0263 (n=646)', 50: '0.0352 (n=824)', 100: '0.0290 (n=999)'}}
   amplitude categorical: {'C(amp)': 'F = 2.74, p = 0.042', 'C(cls)': 'F = 0.01, p = 0.922', 'C(amp):C(cls)': 'F = 1.83, p = 0.139'}
   amplitude numeric: {'C(cls)': 'F = 0.01, p = 0.920', 'amp': 'F = 4.39, p = 0.036', 'amp:C(cls)': 'F = 0.08, p = 0.775'}

Result: not reproduced. No reading gives the four proportions within 10% of FS 0.032 and 0.040 and RS 0.013 and 0.037. Over all FS and RS units (reading a) they are FS 0.0196 and 0.0227 and RS 0.0111 and 0.0221 at 5 and 100 uA under the code's classes. Within the evoked potential (reading b, boundaries from the code) they are closest under reading (c) of fig4i: FS 0.0275 and 0.0324, RS 0.0069 and 0.0299. Responders are sparse in both classes and class and interaction have no effect under the code's classes (class p 0.437 over all units; within the evoked potential 0.983, 0.608 and 0.821 under readings a, b and c of fig4i). The amplitude term, which the paper's F = 3.04 with p = 0.028 implies is categorical, has p = 0.123 over all units and 0.039, 0.005 and 0.051 within the evoked potential under the three readings; numeric, its p is 0.037 or less. Under the Methods' k-means (one run) class has an effect over all units (p = 0.028), and over 20 k-means runs (that run and 19 others; stability check below) class p is below 0.05 in 4.

fig5j-fs-closer

FS direct responders lie closer to the stimulation site than RS responders, and higher amplitude brings both closer.

As the code (cell 57): signed distance of cathodal responders classed FS or RS, ols('value ~ C(amplitude) * C(cell_type)'), type II.

for cls_col in ["wf_kmeans", "wf_threshold"]:
    for col in ["norm_distance", "norm_distance_code"]:
        d = RESP[RESP[cls_col].isin(["fs", "rs"])].assign(cls=lambda x: x[cls_col], value=lambda x: x[col])
        t = anova_lm(smf.ols("value ~ C(amp) * C(cls)", data=d).fit(), typ=2)
        means = d.groupby(["cls", "amp"]).value.agg(lambda x: f"{x.mean():.0f} (n={len(x)})").unstack()
        print(f"{cls_col}, {col}:", {k: f"F = {t.loc[k, 'F']:.2f}, p = {t.loc[k, 'PR(>F)']:.4f}" for k in t.index[:-1]},
              "| mean FS %.0f, RS %.0f" % (d[d.cls == "fs"].value.mean(), d[d.cls == "rs"].value.mean()))
        print(means.to_string())
print("Expected FS 205 (5 uA), RS 567 (5 uA), RS 185 (100 uA) um (10%); paper waveform F = 8.04, p = 0.005; amplitude F = 4.10, p = 0.008")

wf_kmeans, norm_distance: {'C(amp)': 'F = 2.91, p = 0.0366', 'C(cls)': 'F = 2.96, p = 0.0874', 'C(amp):C(cls)': 'F = 0.08, p = 0.9733'} | mean FS 38, RS 244
amp         5           25          50          100
cls                                                
fs    121 (n=2)   106 (n=2)    -8 (n=4)    21 (n=7)
rs   422 (n=25)  287 (n=31)  168 (n=40)  169 (n=37)
wf_kmeans, norm_distance_code: {'C(amp)': 'F = 3.07, p = 0.0299', 'C(cls)': 'F = 2.47, p = 0.1181', 'C(amp):C(cls)': 'F = 0.14, p = 0.9356'} | mean FS 56, RS 255
amp         5           25          50          100
cls                                                
fs    114 (n=2)   207 (n=2)   -30 (n=4)    44 (n=7)
rs   430 (n=25)  320 (n=31)  179 (n=40)  164 (n=37)

wf_threshold, norm_distance: {'C(amp)': 'F = 3.95, p = 0.0097', 'C(cls)': 'F = 7.02, p = 0.0090', 'C(amp):C(cls)': 'F = 0.67, p = 0.5715'} | mean FS 111, RS 267
amp         5           25          50          100
cls                                                
fs   205 (n=12)   214 (n=6)    1 (n=10)   65 (n=14)
rs   556 (n=15)  290 (n=27)  197 (n=34)  183 (n=30)
wf_threshold, norm_distance_code: {'C(amp)': 'F = 4.10, p = 0.0080', 'C(cls)': 'F = 8.04, p = 0.0052', 'C(amp):C(cls)': 'F = 0.56, p = 0.6429'} | mean FS 107, RS 285
amp         5           25          50          100
cls                                                
fs   205 (n=12)   204 (n=6)   -2 (n=10)   60 (n=14)
rs   567 (n=15)  337 (n=27)  208 (n=34)  185 (n=30)
Expected FS 205 (5 uA), RS 567 (5 uA), RS 185 (100 uA) um (10%); paper waveform F = 8.04, p = 0.005; amplitude F = 4.10, p = 0.008

Result: reproduced, resting on an assumption. With the classes of the paper's code (duration threshold), FS responders lie closer than RS responders (mean 111 against 267 um; waveform F = 7.02, p = 0.0090) and amplitude has an effect (F = 3.95, p = 0.0097); FS 205 um and RS 556 and 183 um at 5 and 100 uA are within 10% of 205, 567 and 185 um. With the code's coordinates the paper's F = 8.04 and 4.10 are matched exactly. Under the Methods' k-means classes the waveform effect is not significant (p = 0.0874), and in none of 20 k-means runs (that run and 19 from other seeds and initialisations) is it (p 0.087 to 0.401; stability check below).

fig5i-fs-within-ep

FS responders lie almost entirely within the evoked potential even at low amplitude (75%, 83.3%, 100%, 100%), while only 33% of RS responders do at 5 uA, rising with amplitude.

As fig4f-responders-within-ep, split by waveform class (code cell 54): the fraction of cathodal FS and RS responders whose signed distance lies within the code's mean evoked-potential boundaries (outside input). Classes by the code's threshold and the Methods' k-means; dataset and code coordinates.

rows = []
for col in ["norm_distance", "norm_distance_code"]:
    for cls_col in ["wf_threshold", "wf_kmeans"]:
        for c in ["fs", "rs"]:
            d = RESP[RESP[cls_col] == c]
            w = within_ep(d, col).groupby(d.amp).agg(["mean", "size"])
            rows.append(dict(distances=col, classes=cls_col, cls=c, **{f"{a} uA": f"{w.loc[a, 'mean']:.3f} (n={int(w.loc[a, 'size'])})" for a in estim.AMPS}))
print(pd.DataFrame(rows).to_string(index=False))
print("Expected FS 0.75 (5 uA), 0.833 (25 uA); RS 0.33 (5 uA), 0.741 (25 uA) (10%); FS at or above RS at every amplitude")

         distances      classes cls         5 uA        25 uA        50 uA       100 uA
     norm_distance wf_threshold  fs 0.750 (n=12)  0.833 (n=6) 1.000 (n=10) 1.000 (n=14)
     norm_distance wf_threshold  rs 0.400 (n=15) 0.778 (n=27) 0.882 (n=34) 1.000 (n=30)
     norm_distance    wf_kmeans  fs  1.000 (n=2)  1.000 (n=2)  1.000 (n=4)  1.000 (n=7)
     norm_distance    wf_kmeans  rs 0.520 (n=25) 0.774 (n=31) 0.900 (n=40) 1.000 (n=37)
norm_distance_code wf_threshold  fs 0.750 (n=12)  0.833 (n=6) 1.000 (n=10) 1.000 (n=14)
norm_distance_code wf_threshold  rs 0.333 (n=15) 0.741 (n=27) 0.882 (n=34) 1.000 (n=30)
norm_distance_code    wf_kmeans  fs  1.000 (n=2)  1.000 (n=2)  1.000 (n=4)  1.000 (n=7)
norm_distance_code    wf_kmeans  rs 0.480 (n=25) 0.742 (n=31) 0.900 (n=40) 1.000 (n=37)
Expected FS 0.75 (5 uA), 0.833 (25 uA); RS 0.33 (5 uA), 0.741 (25 uA) (10%); FS at or above RS at every amplitude

Result: not reproduced, resting on an assumption. Under the code's threshold classes and the dataset's coordinates, FS responders lie within the evoked potential at 0.750, 0.833, 1.000 and 1.000 and RS responders at 0.400, 0.778, 0.882 and 1.000 at 5 to 100 uA: FS at 5 and 25 uA and RS at 25 uA are within 10%, FS is at or above RS at every amplitude, and RS at 5 uA (0.400) is outside 0.33 +- 0.033. With the code's coordinates RS gives 0.333 and 0.741, so all four values match. Under the Methods' k-means classes FS responders are 2, 2, 4 and 7 units, all within (1.000), outside 10% of 0.75 and 0.833. The boundaries are the code's (outside input). No class choice gives the four values, and the criterion has no partial clause.

figs6b-fsrs-early-firing

In the 10 ms after a pulse, the firing of FS and RS cortical units differs at every amplitude.

From the S6B legend and the code (cells 62 to 67): per recording, cortical good stable units of each class, first 75 cathodal trials; spikes of the class pooled, counted per trial in bins, averaged over trials and divided by the number of units; the bins from 0 to 10 ms of every recording pooled per class; two-sample KS test per amplitude, Holm-Sidak over the four. Bin widths 1 ms (code) and 2 ms (legend).

for cls_col in ["wf_kmeans", "wf_threshold"]:
    for bw in [0.001, 0.002]:
        edges = np.arange(0, 0.010 + bw / 2, bw)
        p, ns, med = [], [], []
        for a in estim.AMPS:
            samp = {"fs": [], "rs": []}
            for s in CAT_SESS:
                u, tr, trains, meta = DATA[s]
                st = estim.condition_trials(tr, "cathodal", a).start_time.to_numpy()
                for c in ["fs", "rs"]:
                    idx = np.where(u.keep & (u.sub_cort == "Cortical") & (u[cls_col] == c))[0]
                    if len(idx) == 0:
                        continue
                    pop = np.sort(np.concatenate([trains[i] for i in idx]))
                    samp[c].extend(estim.counts_matrix(pop, st, edges).mean(0) / len(idx))
            p.append(stats.ks_2samp(samp["fs"], samp["rs"]).pvalue)
            ns.append((len(samp["fs"]), len(samp["rs"])))
            med.append((round(np.mean(samp["fs"]), 4), round(np.mean(samp["rs"]), 4)))
        adj = estim.holm_sidak(p)
        print(f"{cls_col}, {bw * 1000:.0f} ms bins: p", np.round(p, 4).tolist(), "| Holm-Sidak", np.round(adj, 4).tolist(),
              "| bins (FS, RS)", ns, "| mean spikes per unit and bin (FS, RS)", med)
print("Paper adjusted p: 0.010, 0.007, 0.019, 0.010")

wf_kmeans, 1 ms bins: p [0.0192, 0.0513, 0.0318, 0.0035] | Holm-Sidak [0.0564, 0.0627, 0.0627, 0.0139] | bins (FS, RS) [(70, 70), (70, 70), (70, 70), (70, 70)] | mean spikes per unit and bin (FS, RS) [(np.float64(0.0039), np.float64(0.0036)), (np.float64(0.005), np.float64(0.005)), (np.float64(0.0091), np.float64(0.0086)), (np.float64(0.0141), np.float64(0.0123))]

wf_kmeans, 2 ms bins: p [0.1991, 0.0625, 0.1151, 0.0625] | Holm-Sidak [0.2277, 0.2277, 0.2277, 0.2277] | bins (FS, RS) [(35, 35), (35, 35), (35, 35), (35, 35)] | mean spikes per unit and bin (FS, RS) [(np.float64(0.0078), np.float64(0.0072)), (np.float64(0.0099), np.float64(0.01)), (np.float64(0.0182), np.float64(0.0173)), (np.float64(0.0282), np.float64(0.0245))]

wf_threshold, 1 ms bins: p [0.0019, 0.001, 0.0112, 0.0064] | Holm-Sidak [0.0056, 0.0038, 0.0127, 0.0127] | bins (FS, RS) [(70, 70), (70, 70), (70, 70), (70, 70)] | mean spikes per unit and bin (FS, RS) [(np.float64(0.0051), np.float64(0.0032)), (np.float64(0.0062), np.float64(0.0048)), (np.float64(0.0115), np.float64(0.008)), (np.float64(0.0168), np.float64(0.0115))]

wf_threshold, 2 ms bins: p [0.0153, 0.0153, 0.0319, 0.1151] | Holm-Sidak [0.0598, 0.0598, 0.0629, 0.1151] | bins (FS, RS) [(35, 35), (35, 35), (35, 35), (35, 35)] | mean spikes per unit and bin (FS, RS) [(np.float64(0.0102), np.float64(0.0065)), (np.float64(0.0123), np.float64(0.0095)), (np.float64(0.0231), np.float64(0.016)), (np.float64(0.0336), np.float64(0.0229))]
Paper adjusted p: 0.010, 0.007, 0.019, 0.010

Result: reproduced, resting on an assumption. With the code's classes and 1 ms bins, FS and RS early firing differ at every amplitude (Holm-Sidak p 0.0056, 0.0038, 0.0127, 0.0127). With 2 ms bins (legend) no adjusted p is below 0.05 (0.0598 to 0.1151), and under the Methods' k-means classes only 100 uA differs (1 ms bins, adjusted p 0.0139). Over 20 k-means runs (that run and 19 from other seeds and initialisations) all four amplitudes differ in 9 and one in the other 11 (stability check below).

Stability of the k-means classes (fig5-fsrs-proportion, fig5j-fs-closer, figs6b-fsrs-early-firing)

The Methods' k-means reading above is one run of kmeans2 per recording (seed 0, random points as initial centres). Here it is run from 10 seeds with each of two initialisations (random points and k-means++), and as the run of lowest within-cluster sum of squares over those 20 per recording (the usual way to pick a k-means solution), and the three claims' tests are repeated on each.

Stability. The k-means classes depend on the initialisation: over the 20 runs the recordings hold 296 to 470 good stable FS units (924 under the threshold). fig5j-fs-closer has a waveform p below 0.05 in none of them, figs6b-fsrs-early-firing differs at all four amplitudes in 9, and the class term of fig5-fsrs-proportion has p below 0.05 in 4. The run of lowest within-cluster sum of squares per recording gives 401 FS units, waveform p = 0.1267 and all four amplitudes in figs6b. The code's threshold classes are deterministic; they are the ones the FS/RS outcomes rest on.

KM = {s: [] for s in estim.PAPER_SESSIONS}
for minit in ["points", "++"]:
    for seed in range(10):
        for s in estim.PAPER_SESSIONS:
            u = DATA[s][0]
            lab, w = estim.kmeans_fs_rs(u.wf_duration, u.wf_ptratio, u.wf_slope, seed=seed, minit=minit)
            KM[s].append((f"{minit} {seed}", estim.label_up_axon(lab, u.wf_duration, u.wf_ptratio), w))
runs = [r[0] for r in KM[estim.PAPER_SESSIONS[0]]] + ["lowest WCSS"]
def km_labels(run):
    if run == "lowest WCSS":
        return {s: min(KM[s], key=lambda r: r[2])[1] for s in estim.PAPER_SESSIONS}
    j = runs.index(run)
    return {s: KM[s][j][1] for s in estim.PAPER_SESSIONS}
edges1 = np.arange(0, 0.0105, 0.001)
rows = []
for run in runs:
    lab = km_labels(run)
    cls = lambda df: [lab[s][i] for s, i in zip(df.session, df.unit)]
    nfs = sum(int(((lab[s] == "fs") & DATA[s][0].keep).sum()) for s in estim.PAPER_SESSIONS)
    nrs = sum(int(((lab[s] == "rs") & DATA[s][0].keep).sum()) for s in estim.PAPER_SESSIONS)
    dc = CAT.assign(cls=cls(CAT), resp=lambda x: (x.ndrp > 0.15).astype(float))
    dc = dc[dc.cls.isin(["fs", "rs"])]
    tf = anova_lm(smf.ols("resp ~ C(amp) * C(cls)", data=dc).fit(), typ=2)
    dj = RESP.assign(cls=cls(RESP), value=lambda x: x.norm_distance)
    dj = dj[dj.cls.isin(["fs", "rs"])]
    tj = anova_lm(smf.ols("value ~ C(amp) * C(cls)", data=dj).fit(), typ=2)
    p6 = []
    for a in estim.AMPS:
        samp = {"fs": [], "rs": []}
        for s in CAT_SESS:
            u, tr, trains, meta = DATA[s]
            st = estim.condition_trials(tr, "cathodal", a).start_time.to_numpy()
            for c in ["fs", "rs"]:
                idx = np.where(u.keep & (u.sub_cort == "Cortical") & (lab[s] == c))[0]
                if len(idx):
                    pop = np.sort(np.concatenate([trains[i] for i in idx]))
                    samp[c].extend(estim.counts_matrix(pop, st, edges1).mean(0) / len(idx))
        p6.append(stats.ks_2samp(samp["fs"], samp["rs"]).pvalue)
    adj6 = estim.holm_sidak(p6)
    rows.append(dict(run=run, fs_units=nfs, rs_units=nrs,
                     fsrs_FS5=round(dc[(dc.cls == "fs") & (dc.amp == 5)].resp.mean(), 4), fsrs_RS100=round(dc[(dc.cls == "rs") & (dc.amp == 100)].resp.mean(), 4),
                     fsrs_class_p=round(tf.loc["C(cls)", "PR(>F)"], 3), fsrs_amp_p=round(tf.loc["C(amp)", "PR(>F)"], 3),
                     fig5j_wave_p=round(tj.loc["C(cls)", "PR(>F)"], 4), fig5j_amp_p=round(tj.loc["C(amp)", "PR(>F)"], 4),
                     fig5j_FS_mean=round(dj[dj.cls == "fs"].value.mean()), fig5j_RS_mean=round(dj[dj.cls == "rs"].value.mean()),
                     s6b_adj_max=round(adj6.max(), 4), s6b_sig=int((adj6 < 0.05).sum())))
KMS = pd.DataFrame(rows)
print(KMS.to_string(index=False))
print("Run 'points 0' is the k-means run of the sections above; the other 19 use other seeds and initialisations.")
print(f"\nOver the 20 runs (the run above and 19 others): FS units {KMS.fs_units[:-1].min()} to {KMS.fs_units[:-1].max()} (threshold classes: "
      f"{sum(int(((DATA[s][0].wf_threshold == 'fs') & DATA[s][0].keep).sum()) for s in estim.PAPER_SESSIONS)});"
      f" fig5j waveform p < 0.05 in {int((KMS.fig5j_wave_p[:-1] < 0.05).sum())}; figs6b significant at all four amplitudes in {int((KMS.s6b_sig[:-1] == 4).sum())};"
      f" fig5-fsrs class p < 0.05 in {int((KMS.fsrs_class_p[:-1] < 0.05).sum())}")

        run  fs_units  rs_units  fsrs_FS5  fsrs_RS100  fsrs_class_p  fsrs_amp_p  fig5j_wave_p  fig5j_amp_p  fig5j_FS_mean  fig5j_RS_mean  s6b_adj_max  s6b_sig
   points 0       470      2155    0.0061      0.0225         0.028       0.122        0.0874       0.0366             38            244       0.0627        1
   points 1       405      2220    0.0074      0.0223         0.138       0.124        0.1267       0.0335             53            241       0.0104        4
   points 2       417      2208    0.0071      0.0219         0.147       0.124        0.0874       0.0366             38            244       0.0627        1
   points 3       372      2253    0.0083      0.0230         0.079       0.124        0.2094       0.0271             76            235       0.0127        4
   points 4       421      2204    0.0070      0.0219         0.134       0.124        0.0874       0.0366             38            244       0.0627        1
   points 5       404      2221    0.0074      0.0223         0.142       0.124        0.1267       0.0335             53            241       0.0104        4
   points 6       363      2262    0.0086      0.0229         0.102       0.124        0.2094       0.0271             76            235       0.0127        4
   points 7       466      2159    0.0062      0.0225         0.031       0.122        0.0874       0.0366             38            244       0.0627        1
   points 8       420      2205    0.0070      0.0219         0.138       0.124        0.0874       0.0366             38            244       0.0627        1
   points 9       362      2263    0.0086      0.0229         0.106       0.124        0.2094       0.0271             76            235       0.0127        4
       ++ 0       362      2263    0.0086      0.0229         0.106       0.124        0.2094       0.0271             76            235       0.0127        4
       ++ 1       421      2204    0.0070      0.0219         0.134       0.124        0.0874       0.0366             38            244       0.0627        1
       ++ 2       311      2314    0.0054      0.0229         0.172       0.123        0.2838       0.0303             71            233       0.0318        4
       ++ 3       470      2155    0.0061      0.0225         0.028       0.122        0.0874       0.0366             38            244       0.0627        1
       ++ 4       296      2329    0.0059      0.0233         0.167       0.123        0.4012       0.0280            100            230       0.0037        4
       ++ 5       379      2246    0.0080      0.0226         0.101       0.124        0.1461       0.0291             56            238       0.0925        1
       ++ 6       422      2203    0.0070      0.0231         0.023       0.123        0.1461       0.0291             56            238       0.0925        1
       ++ 7       420      2205    0.0070      0.0219         0.138       0.124        0.0874       0.0366             38            244       0.0627        1
       ++ 8       420      2205    0.0070      0.0219         0.138       0.124        0.0874       0.0366             38            244       0.0627        1
       ++ 9       296      2329    0.0059      0.0233         0.167       0.123        0.4012       0.0280            100            230       0.0037        4
lowest WCSS       401      2224    0.0075      0.0222         0.151       0.124        0.1267       0.0335             53            241       0.0104        4
Run 'points 0' is the k-means run of the sections above; the other 19 use other seeds and initialisations.

Over the 20 runs (the run above and 19 others): FS units 296 to 470 (threshold classes: 924); fig5j waveform p < 0.05 in 0; figs6b significant at all four amplitudes in 9; fig5-fsrs class p < 0.05 in 4

fig6c-modulation-cortex

Spiking modulation after a single pulse increases with amplitude for cortical units (12.60, 13.37, 15.11, 20.21) but stays about constant for subcortical units (12.27, 12.60, 12.08, 13.48), differing between them at 50 and 100 uA.

MI per bin against the median rate 1 s to 0.5 s before onset; magnitude = trapezoid area of |MI| over bins (code, figure7.ipynb cells 28 and 33). Readings (a) 0 to 500 ms, 10 ms bins (code); (b) 0 to 500 ms, 5 ms bins (legend); (c) 500 to 1000 ms, 5 ms bins (Methods). Good stable units, first 75 cathodal trials, cortical or subcortical by the code's classify_region.

READINGS = {"a": (0.0, 0.5, 0.010), "b": (0.0, 0.5, 0.005), "c": (0.5, 1.0, 0.005)}
mi_rows = []
MI50 = {}
for s in CAT_SESS:
    u, tr, trains, meta = DATA[s]
    idx = np.where(u.keep & u.sub_cort.notna())[0]
    for a in estim.AMPS:
        st = estim.condition_trials(tr, "cathodal", a).start_time.to_numpy()
        for i in idx:
            row = dict(session=s, unit=i, amp=a, Region=u.sub_cort[i], region=u.region[i], distance=u.distance[i])
            for k, (t0, t1, bw) in READINGS.items():
                mi = estim.modulation_index(trains[i], st, t0, t1, bw)
                row[f"auc_{k}"] = np.trapezoid(np.abs(mi))
            mi_rows.append(row)
MIDF = pd.DataFrame(mi_rows)
print("Units:", MIDF.groupby(["Region"]).unit.size().div(4).to_dict(), "per amplitude")
for k in READINGS:
    col = f"auc_{k}"
    print(f"\nReading {k}: mean +- SD")
    print(MIDF.groupby(["Region", "amp"])[col].agg(lambda x: f"{x.mean():.2f} +- {x.std():.2f}").unstack().to_string())
    t = anova_lm(smf.ols(f"{col} ~ amp * C(Region)", data=MIDF).fit(), typ=2)
    print("   ANOVA:", {i: f"F = {t.loc[i, 'F']:.2f}, p = {t.loc[i, 'PR(>F)']:.2e}" for i in t.index[:-1]})
    for a in estim.AMPS:
        c = MIDF[(MIDF.amp == a) & (MIDF.Region == "Cortical")][col]
        sc = MIDF[(MIDF.amp == a) & (MIDF.Region == "Subcortical")][col]
        print(f"   {a} uA: Student p = {stats.ttest_ind(c, sc).pvalue:.2e}, Welch p = {stats.ttest_ind(c, sc, equal_var=False).pvalue:.2e}")
    slope = {r: np.polyfit(g.amp, g[col], 1)[0] for r, g in MIDF.groupby("Region")}
    print("   slope on amplitude per region:", {r: round(v, 4) for r, v in slope.items()})
print("By recording, reading a, cortical mean at 5 and 100 uA:",
      MIDF[MIDF.Region == "Cortical"].groupby(["session", "amp"]).auc_a.mean().unstack()[[5, 100]].round(2).to_dict("index"))
print("\nExpected cortical 12.60 (5 uA), 20.21 (100 uA); subcortical 12.27, 13.48 (10%)")

Units: {'Cortical': 1287.0, 'Subcortical': 351.0} per amplitude

Reading a: mean +- SD
amp                    5              25             50              100
Region                                                                  
Cortical     12.60 +- 8.43  13.37 +- 8.41  15.11 +- 8.80  20.21 +- 10.32
Subcortical  12.27 +- 7.21  12.60 +- 8.00  12.08 +- 7.52   13.48 +- 7.70
   ANOVA: {'C(Region)': 'F = 106.27, p = 9.94e-25', 'amp': 'F = 487.36, p = 3.14e-104', 'amp:C(Region)': 'F = 90.88, p = 2.10e-21'}
   5 uA: Student p = 5.01e-01, Welch p = 4.62e-01
   25 uA: Student p = 1.26e-01, Welch p = 1.16e-01
   50 uA: Student p = 4.65e-09, Welch p = 2.39e-10
   100 uA: Student p = 5.50e-29, Welch p = 6.45e-37
   slope on amplitude per region: {'Cortical': np.float64(0.0823), 'Subcortical': np.float64(0.0116)}

Reading b: mean +- SD
amp                     5               25              50              100
Region                                                                     
Cortical     24.72 +- 17.04  26.50 +- 16.90  29.33 +- 18.58  37.04 +- 21.47
Subcortical  27.20 +- 15.95  27.15 +- 15.76  26.58 +- 15.29  27.83 +- 15.45
   ANOVA: {'C(Region)': 'F = 16.57, p = 4.76e-05', 'amp': 'F = 280.85, p = 9.35e-62', 'amp:C(Region)': 'F = 67.53, p = 2.47e-16'}
   5 uA: Student p = 1.42e-02, Welch p = 1.11e-02
   25 uA: Student p = 5.16e-01, Welch p = 4.99e-01
   50 uA: Student p = 1.10e-02, Welch p = 4.62e-03
   100 uA: Student p = 8.75e-14, Welch p = 1.27e-18
   slope on amplitude per region: {'Cortical': np.float64(0.1317), 'Subcortical': np.float64(0.0064)}

Reading c: mean +- SD
amp                     5               25              50              100
Region                                                                     
Cortical     24.64 +- 16.69  25.90 +- 16.70  26.85 +- 16.95  27.49 +- 17.41
Subcortical  27.05 +- 15.72  27.11 +- 15.89  26.46 +- 15.21  25.75 +- 15.29
   ANOVA: {'C(Region)': 'F = 0.55, p = 4.57e-01', 'amp': 'F = 10.72, p = 1.06e-03', 'amp:C(Region)': 'F = 9.40, p = 2.18e-03'}
   5 uA: Student p = 1.52e-02, Welch p = 1.22e-02
   25 uA: Student p = 2.26e-01, Welch p = 2.14e-01
   50 uA: Student p = 6.99e-01, Welch p = 6.81e-01
   100 uA: Student p = 8.86e-02, Welch p = 6.71e-02
   slope on amplitude per region: {'Cortical': np.float64(0.0282), 'Subcortical': np.float64(-0.015)}
By recording, reading a, cortical mean at 5 and 100 uA: {'jlh31': {5: 11.27, 100: 20.67}, 'jlh32': {5: 13.95, 100: 25.0}, 'jlh33': {5: 11.7, 100: 23.9}, 'jlh34': {5: 14.89, 100: 22.52}, 'jlh40': {5: 12.56, 100: 16.42}, 'jlh48': {5: 11.24, 100: 16.51}, 'jlh53concat': {5: 12.62, 100: 15.07}}

Expected cortical 12.60 (5 uA), 20.21 (100 uA); subcortical 12.27, 13.48 (10%)

Result: not reproduced. Under the code's reading (a) the magnitudes are 12.60, 13.37, 15.11 and 20.21 (cortical) and 12.27, 12.60, 12.08 and 13.48 (subcortical), the paper's values, with region F = 106.27, amplitude F = 487.36 and interaction F = 90.88 (all p < 1e-20) and cortical against subcortical Student p = 4.65e-09 at 50 uA and 5.50e-29 at 100 uA. Under the legend's 5 ms bins (b) the effects hold but the 50 uA comparison has p = 1.10e-02. Under the Methods' window of 500 to 1000 ms (c) there is no region effect (p = 0.457), the amplitude term has p = 1.06e-03, above 0.001, and no amplitude comparison has p < 0.001. The criterion asks for every reading and has no partial clause.

fig6d-hva-subcortical

At 100 uA higher visual areas, beyond the reach of direct activation, are more modulated than subcortical units (15.82 against 13.48) although both lie at similar distances from the stimulation site.

Groups by the Figure 6 legend: V1 (primary visual area), HVA (LM, AM, PM), subcortical ("hippocampus (CA1, subiculum, and dentate gyrus) and superior colliculus": Field CA1, Postsubiculum, the dentate gyrus and superior colliculus layers). Variants: the legend's HVA with the code's subcortical list (which has no CA1), and the code's groups (HVA any CCF name with "layer" and not "primary"). Tukey HSD over the three at 100 uA cathodal (magnitude under reading (a)); Mann-Whitney test of distances, HVA against subcortical, over the good stable units of all nine recordings (as the code, cells 40 and 41).

# Magnitude at 100 uA (reading a) for every good stable unit of the recordings with cathodal pulses,
# so that the legend's subcortical regions outside the code's list (CA1) are included
rows = []
for s in CAT_SESS:
    u, tr, trains, meta = DATA[s]
    st = estim.condition_trials(tr, "cathodal", 100).start_time.to_numpy()
    for i in np.where(u.keep)[0]:
        mi = estim.modulation_index(trains[i], st, 0.0, 0.5, 0.010)
        rows.append(dict(session=s, unit=i, region=u.region[i], distance=u.distance[i], auc_a=np.trapezoid(np.abs(mi))))
M100 = pd.DataFrame(rows)
def hva_code(r):
    return ("layer" in r) and ("primary" not in r.lower())
def sub_code(r):
    return estim.classify_region(r) == "Subcortical"
allkeep = pd.concat([DATA[s][0].assign(session=s) for s in estim.PAPER_SESSIONS])
allkeep = allkeep[allkeep.keep]
print("Legend subcortical group, good stable units of the nine recordings by CCF name:",
      allkeep[allkeep.region.map(estim.legend_subcortical)].region.value_counts().to_dict())
print("Field CA1 units among them, which the code's subcortical list leaves out:", int((allkeep.region == "Field CA1").sum()), "\n")
for name, hva, sub in [("legend (HVA LM, AM, PM; subcortical CA1, subiculum, dentate gyrus, superior colliculus)", estim.legend_hva, estim.legend_subcortical),
                       ("legend HVA, the code's subcortical list", estim.legend_hva, sub_code),
                       ("code (HVA any layered name outside V1, the code's subcortical list)", hva_code, sub_code)]:
    d = M100.copy()
    d["group"] = np.where(d.region.str.contains("Primary visual area"), "V1",
                          np.where(d.region.map(hva), "HVA", np.where(d.region.map(sub), "Subcortical", "other")))
    d = d[d.group != "other"]
    print(f"Groups: {name}: n", d.group.value_counts().to_dict(), "| of the subcortical, Field CA1:",
          int((d[d.group == "Subcortical"].region == "Field CA1").sum()), "| of the HVA, superior colliculus or dentate gyrus:",
          int(d[d.group == "HVA"].region.str.contains("Superior colliculus|Dentate gyrus").sum()), "| mean +- SD",
          {g: f"{x.mean():.2f} +- {x.std():.2f}" for g, x in d.groupby("group").auc_a})
    tk = pairwise_tukeyhsd(d.auc_a, d.group)
    print(pd.DataFrame(tk.summary().data[1:], columns=tk.summary().data[0]).to_string(index=False))
    k = allkeep.copy()
    k["group"] = np.where(k.region.map(hva) & ~k.region.str.contains("Primary visual area"), "HVA",
                          np.where(k.region.map(sub), "Subcortical", "other"))
    h, sc = k[k.group == "HVA"].distance, k[k.group == "Subcortical"].distance
    print(f"   distance HVA median {h.median():.0f} um (n={len(h)}), subcortical {sc.median():.0f} um (n={len(sc)}), Mann-Whitney p = {stats.mannwhitneyu(h, sc).pvalue:.3g}\n")
print("Expected HVA 15.82, subcortical 13.48 (10%); paper Tukey p = 0.004; Mann-Whitney p = 0.511")

Legend subcortical group, good stable units of the nine recordings by CCF name: {'Field CA1': 315, 'Postsubiculum': 154, 'Superior colliculus, motor related, intermediate white layer': 68, 'Superior colliculus, motor related, intermediate gray layer': 43, 'Dentate gyrus molecular layer': 36, 'Dentate gyrus granule cell layer': 4}
Field CA1 units among them, which the code's subcortical list leaves out: 315 

Groups: legend (HVA LM, AM, PM; subcortical CA1, subiculum, dentate gyrus, superior colliculus): n {'V1': 1029, 'Subcortical': 460, 'HVA': 120} | of the subcortical, Field CA1: 273 | of the HVA, superior colliculus or dentate gyrus: 0 | mean +- SD {'HVA': '13.91 +- 9.85', 'Subcortical': '13.26 +- 8.36', 'V1': '21.06 +- 10.13'}
     group1      group2  meandiff  p-adj   lower  upper  reject
        HVA Subcortical   -0.6559 0.7843 -2.9723 1.6604   False
        HVA          V1    7.1414 0.0000  4.9615 9.3212    True
Subcortical          V1    7.7973 0.0000  6.5299 9.0647    True
   distance HVA median 1219 um (n=139), subcortical 1053 um (n=620), Mann-Whitney p = 0.54

Groups: legend HVA, the code's subcortical list: n {'V1': 1029, 'Subcortical': 351, 'HVA': 120} | of the subcortical, Field CA1: 0 | of the HVA, superior colliculus or dentate gyrus: 0 | mean +- SD {'HVA': '13.91 +- 9.85', 'Subcortical': '13.48 +- 7.70', 'V1': '21.06 +- 10.13'}

     group1      group2  meandiff  p-adj   lower  upper  reject
        HVA Subcortical   -0.4387  0.902 -2.8180 1.9407   False
        HVA          V1    7.1414  0.000  4.9709 9.3118    True
Subcortical          V1    7.5800  0.000  6.1892 8.9709    True
   distance HVA median 1219 um (n=139), subcortical 1033 um (n=499), Mann-Whitney p = 0.577

Groups: code (HVA any layered name outside V1, the code's subcortical list): n {'V1': 1029, 'HVA': 349, 'Subcortical': 260} | of the subcortical, Field CA1: 0 | of the HVA, superior colliculus or dentate gyrus: 91 | mean +- SD {'HVA': '15.82 +- 9.97', 'Subcortical': '13.68 +- 7.64', 'V1': '21.06 +- 10.13'}
     group1      group2  meandiff  p-adj   lower   upper  reject
        HVA Subcortical   -2.1404 0.0202 -4.0122 -0.2685    True
        HVA          V1    5.2385 0.0000  3.8232  6.6539    True
Subcortical          V1    7.3789 0.0000  5.7930  8.9648    True
   distance HVA median 1264 um (n=428), subcortical 939 um (n=348), Mann-Whitney p = 5.92e-08

Expected HVA 15.82, subcortical 13.48 (10%); paper Tukey p = 0.004; Mann-Whitney p = 0.511

Result: not reproduced. With the legend's groups (HVA: LM, AM and PM, 120 units; subcortical: CA1, subiculum, dentate gyrus and superior colliculus, 460 units in the seven recordings with cathodal pulses, of which 273 are CA1 units that the code's list leaves out) HVA units are modulated at 13.91 and subcortical units at 13.26 at 100 uA (Tukey p = 0.784); the distances do not differ (Mann-Whitney p = 0.54), as in the paper. With the code's subcortical list the values are 13.91 and 13.48 (p = 0.902). With the code's broader HVA group (any layered name outside V1, 349 units, of which 91 are superior colliculus or dentate gyrus units) they are 15.82 and 13.68 and Tukey p = 0.0202, but then the distances differ (p = 5.92e-08).

fig6e-circuit-responders

The proportion of cortical units in the circuit response grows with amplitude (0.14, 0.21, 0.52, 0.79); subcortical units respond less but also increase (0.08, 0.09, 0.17, 0.35).

ZETA test (zetapy 2.7.2 getZeta, 300 ms after onset, seeded per unit and condition), p < 0.05. As the code (figure7.ipynb cells 1 and 3 to 8, jlh_ephys/analysis.py get_zetas), the trials are all those of the run that holds the condition's first trial (tr.run == run, with no other filter; the table below gives their composition), the units good and stable, and the regions the code's map_region_to_group. The tests are run for the cathodal conditions, which give the claim's values and the cathodal model Value ~ Amplitude * C(Region) (type II), whose F values the Results also print. Running all three polarities for the code's full model Value ~ Amplitude * C(Polarity) * C(Region) would take about 20 more minutes; the full model is shown with the paper's cached p-values (outside input, a check only). In jlh53concat, a concatenated recording, the run numbers of the monopolar runs are also used by bipolar runs, so a run there holds 150 trials of both polarities; the condition's own 75 trials are run as a variant.

jobs, keys = [], []
for si, s in enumerate(estim.PAPER_SESSIONS):
    u, tr, trains, meta = DATA[s]
    idx = np.where(u.keep & u.region_group.notna())[0]
    for ai, a in enumerate(estim.AMPS):
        c = estim.condition_trials(tr, "cathodal", a, n=None)
        if len(c) == 0:
            continue
        run = c.run.iloc[0]
        st = tr[tr.run == run].start_time.to_numpy()
        for i in idx:
            jobs.append((trains[i], st, int(1_000_000 * si + 10_000 * ai + i)))
            keys.append((s, i, "cathodal", a, int(run)))
comp = []
for s in estim.PAPER_SESSIONS:
    tr = DATA[s][1]
    for a in estim.AMPS:
        c = estim.condition_trials(tr, "cathodal", a, n=None)
        if len(c):
            r = tr[tr.run == c.run.iloc[0]]
            comp.append(dict(session=s, amp=a, run=int(c.run.iloc[0]), trials=len(r), amplitude=sorted(r.amplitude.unique().tolist()),
                             polarity=sorted(r.polarity.unique()), pulse_duration=sorted(r.pulse_duration.unique().tolist()),
                             pulse_number=sorted(r.pulse_number.unique().tolist()), contact_negative=sorted(r.contact_negative.unique().tolist())))
print("Trials of each run used for ZETA:")
print(pd.DataFrame(comp).to_string(index=False))
t = time.time()
ZP = estim.zeta_pvalues(jobs)
print(f"{len(jobs)} ZETA tests in {time.time() - t:.0f} s; failed: {int(np.isnan(ZP).sum())}")
Z = pd.DataFrame(keys, columns=["session", "unit", "pol", "amp", "run"]).assign(p=ZP)
Z = Z.join(pd.concat([DATA[s][0].assign(session=s).set_index(["session", "unit"], drop=False)[["region_group", "sub_cort", "distance"]]
                      for s in estim.PAPER_SESSIONS]), on=["session", "unit"])
Z["Region"] = np.where(Z.region_group.str.contains("Area"), "cortical", np.where(Z.region_group == "subcortical", "subcortical", None))
Z["sig"] = Z.p < 0.05
props = Z[Z.Region.notna()].groupby(["Region", "amp", "session"]).sig.mean().rename("Value").reset_index()
print(props.groupby(["Region", "amp"]).Value.agg(lambda x: f"{x.mean():.3f} +- {x.std():.3f} (n={len(x)})").unstack().to_string())
print("By recording, cortical:", props[props.Region == "cortical"].pivot(index="session", columns="amp", values="Value").round(2).to_dict("index"))
t = anova_lm(smf.ols("Value ~ amp * C(Region)", data=props).fit(), typ=2)
print("Cathodal model:", {i: f"F = {t.loc[i, 'F']:.2f}, p = {t.loc[i, 'PR(>F)']:.2e}" for i in t.index[:-1]})

# The paper's cached ZETA p-values (outside input, a check only)
CZ = {s: estim.code_zeta(s) for s in estim.PAPER_SESSIONS}
zc = Z.assign(p_code=[CZ[s][r][i] if r in CZ[s] else np.nan for s, r, i in zip(Z.session, Z.run, Z.unit)])
zc = zc[zc.p_code.notna()]
print(f"\nCached p-values found for {len(zc)} of {len(Z)} tests; same decision (p < 0.05) for {((zc.p < 0.05) == (zc.p_code < 0.05)).mean():.3f}; "
      f"responders here {(zc.p < 0.05).mean():.3f}, in the cache {(zc.p_code < 0.05).mean():.3f}; Pearson r of p {stats.pearsonr(zc.p, zc.p_code)[0]:.3f}")
print("Same decision by recording:", zc.groupby("session").apply(lambda g: round(((g.p < 0.05) == (g.p_code < 0.05)).mean(), 3)).to_dict())

# Variant: jlh53concat with the condition's own 75 cathodal trials in place of the whole run
s = "jlh53concat"; si = estim.PAPER_SESSIONS.index(s)
u, tr, trains, meta = DATA[s]
idx = np.where(u.keep & u.region_group.notna())[0]
jobs2, keys2 = [], []
for ai, a in enumerate(estim.AMPS):
    st = estim.condition_trials(tr, "cathodal", a, n=None).start_time.to_numpy()
    for i in idx:
        jobs2.append((trains[i], st, int(1_000_000 * si + 10_000 * ai + i)))
        keys2.append((s, i, a))
ZP2 = pd.Series(estim.zeta_pvalues(jobs2), index=pd.MultiIndex.from_tuples(keys2))
Zv = Z.copy()
m = Zv.session == s
Zv.loc[m, "p"] = [ZP2[(s, i, a)] for i, a in zip(Zv[m].unit, Zv[m].amp)]
Zv["sig"] = Zv.p < 0.05
pv = Zv[Zv.Region.notna()].groupby(["Region", "amp", "session"]).sig.mean().rename("Value").reset_index()
print(f"\njlh53concat with its own trials ({len(jobs2)} tests): proportions",
      pv.groupby(["Region", "amp"]).Value.mean().round(3).unstack().to_dict("index"),
      "| jlh53concat alone", pv[pv.session == s].pivot(index="Region", columns="amp", values="Value").round(2).to_dict("index"))
t = anova_lm(smf.ols("Value ~ amp * C(Region)", data=pv).fit(), typ=2)
print("   cathodal model:", {i: f"F = {t.loc[i, 'F']:.2f}, p = {t.loc[i, 'PR(>F)']:.2e}" for i in t.index[:-1]})
rows = []
for s in estim.PAPER_SESSIONS:
    u, tr = DATA[s][0], DATA[s][1]
    idx = np.where(u.keep & u.region_group.notna())[0]
    reg = np.where(u.region_group.str.contains("Area", na=False), "cortical", np.where(u.region_group == "subcortical", "subcortical", None))
    for pol in estim.POLARITIES:
        for a in estim.AMPS:
            c = estim.condition_trials(tr, pol, a, n=None)
            if len(c) == 0 or int(c.run.iloc[0]) not in CZ[s]:
                continue
            pc = CZ[s][int(c.run.iloc[0])]
            for R in ["cortical", "subcortical"]:
                ii = idx[reg[idx] == R]
                rows.append(dict(session=s, pol=pol, amp=a, Region=R, Value=np.mean(pc[ii] < 0.05)))
cp = pd.DataFrame(rows)
print("Cached p-values, cathodal:", cp[cp.pol == "cathodal"].groupby(["Region", "amp"]).Value.mean().round(2).unstack().to_dict("index"))
t = anova_lm(smf.ols("Value ~ amp * C(pol) * C(Region)", data=cp).fit(), typ=2)
print("Cached p-values, full model:", {i: f"F = {t.loc[i, 'F']:.2f}, p = {t.loc[i, 'PR(>F)']:.2e}" for i in t.index[:-1]})
print("Units with 'posteromedial visual area' names, which the code's grouping drops:",
      int(allkeep.region.str.contains("posteromedial visual area").sum()))
print("Expected cortical 0.14, 0.21, 0.52, 0.79; subcortical 0.08 (5 uA), 0.35 (100 uA) (10%); paper F amplitude 175.18, region 81.50 (full model)")

Trials of each run used for ZETA:
    session  amp  run  trials amplitude             polarity pulse_duration pulse_number contact_negative
      jlh31    5   16      75      [-5]          [monopolar]          [100]          [1]             [16]
      jlh31   25   18      75     [-25]          [monopolar]          [100]          [1]             [16]
      jlh31   50   20      75     [-50]          [monopolar]          [100]          [1]             [16]
      jlh31  100   22      75    [-100]          [monopolar]          [100]          [1]             [16]
      jlh32    5   12      75      [-5]          [monopolar]          [100]          [1]              [3]
      jlh32   25   14      75     [-25]          [monopolar]          [100]          [1]              [3]
      jlh32   50   16      75     [-50]          [monopolar]          [100]          [1]              [3]
      jlh32  100   18      75    [-100]          [monopolar]          [100]          [1]              [3]
      jlh33    5    8      60      [-5]          [monopolar]          [100]          [1]              [5]
      jlh33   25   10      60     [-25]          [monopolar]          [100]          [1]              [5]
      jlh33   50   12      60     [-50]          [monopolar]          [100]          [1]              [5]
      jlh33  100   13      60    [-100]          [monopolar]          [100]          [1]              [5]
      jlh34    5    7      75      [-5]          [monopolar]          [100]          [1]             [15]
      jlh34   25    9      75     [-25]          [monopolar]          [100]          [1]             [15]
      jlh34   50   11      75     [-50]          [monopolar]          [100]          [1]             [15]
      jlh34  100   13      75    [-100]          [monopolar]          [100]          [1]             [15]
      jlh40    5    0      75      [-5]          [monopolar]          [100]          [1]             [16]
      jlh40   25    1      75     [-25]          [monopolar]          [100]          [1]             [16]
      jlh40   50    2      75     [-50]          [monopolar]          [100]          [1]             [16]
      jlh40  100    6      75    [-100]          [monopolar]          [100]          [1]             [16]
      jlh48    5    8      75      [-5]          [monopolar]          [100]          [1]             [14]
      jlh48   25   10      75     [-25]          [monopolar]          [100]          [1]             [14]
      jlh48   50   12      75     [-50]          [monopolar]          [100]          [1]             [14]
      jlh48  100   14      75    [-100]          [monopolar]          [100]          [1]             [14]
jlh53concat    5    0     150      [-5] [bipolar, monopolar]          [100]          [1]             [16]
jlh53concat   25    2     150     [-25] [bipolar, monopolar]          [100]          [1]             [16]
jlh53concat   50    4     150     [-50] [bipolar, monopolar]          [100]          [1]             [16]
jlh53concat  100    6     150    [-100] [bipolar, monopolar]          [100]          [1]             [16]

7980 ZETA tests in 518 s; failed: 0
amp                           5                     25                    50                    100
Region                                                                                             
cortical     0.125 +- 0.051 (n=7)  0.209 +- 0.087 (n=7)  0.520 +- 0.255 (n=7)  0.785 +- 0.215 (n=7)
subcortical  0.089 +- 0.015 (n=7)  0.094 +- 0.036 (n=7)  0.164 +- 0.061 (n=7)  0.360 +- 0.197 (n=7)
By recording, cortical: {'jlh31': {5: 0.18, 25: 0.25, 50: 0.56, 100: 0.88}, 'jlh32': {5: 0.06, 25: 0.36, 50: 0.86, 100: 0.97}, 'jlh33': {5: 0.12, 25: 0.12, 50: 0.74, 100: 0.96}, 'jlh34': {5: 0.08, 25: 0.17, 50: 0.65, 100: 0.93}, 'jlh40': {5: 0.16, 25: 0.27, 50: 0.37, 100: 0.39}, 'jlh48': {5: 0.19, 25: 0.16, 50: 0.28, 100: 0.72}, 'jlh53concat': {5: 0.09, 25: 0.14, 50: 0.16, 100: 0.64}}
Cathodal model: {'C(Region)': 'F = 36.59, p = 1.61e-07', 'amp': 'F = 89.48, p = 6.84e-13', 'amp:C(Region)': 'F = 15.26, p = 2.72e-04'}

Cached p-values found for 7980 of 7980 tests; same decision (p < 0.05) for 0.970; responders here 0.325, in the cache 0.326; Pearson r of p 0.992
Same decision by recording: {'jlh31': 0.978, 'jlh32': 0.968, 'jlh33': 0.968, 'jlh34': 0.971, 'jlh40': 0.965, 'jlh48': 0.964, 'jlh53concat': 0.97}

jlh53concat with its own trials (1476 tests): proportions {'cortical': {5: 0.121, 25: 0.203, 50: 0.514, 100: 0.756}, 'subcortical': {5: 0.095, 25: 0.09, 50: 0.141, 100: 0.339}} | jlh53concat alone {'cortical': {5: 0.06, 25: 0.09, 50: 0.12, 100: 0.44}, 'subcortical': {5: 0.1, 25: 0.05, 50: 0.05, 100: 0.39}}
   cathodal model: {'C(Region)': 'F = 32.74, p = 5.27e-07', 'amp': 'F = 71.96, p = 2.21e-11', 'amp:C(Region)': 'F = 13.82, p = 4.93e-04'}
Cached p-values, cathodal: {'cortical': {5: 0.14, 25: 0.21, 50: 0.52, 100: 0.79}, 'subcortical': {5: 0.08, 25: 0.09, 50: 0.17, 100: 0.35}}
Cached p-values, full model: {'C(pol)': 'F = 5.13, p = 6.89e-03', 'C(Region)': 'F = 85.99, p = 1.10e-16', 'C(pol):C(Region)': 'F = 1.07, p = 3.45e-01', 'amp': 'F = 172.99, p = 2.40e-27', 'amp:C(pol)': 'F = 3.93, p = 2.16e-02', 'amp:C(Region)': 'F = 46.27, p = 1.89e-10', 'amp:C(pol):C(Region)': 'F = 0.09, p = 9.16e-01'}
Units with 'posteromedial visual area' names, which the code's grouping drops: 16
Expected cortical 0.14, 0.21, 0.52, 0.79; subcortical 0.08 (5 uA), 0.35 (100 uA) (10%); paper F amplitude 175.18, region 81.50 (full model)

Result: partially reproduced. With zetapy's getZeta the proportion of cortical circuit responders is 0.125, 0.209, 0.520 and 0.785 and of subcortical ones 0.089, 0.094, 0.164 and 0.360 at 5 to 100 uA, rising with amplitude and above one half for cortex at 50 and 100 uA; the cathodal model gives amplitude F = 89.48 (p = 6.84e-13) and region F = 36.59 (p = 1.61e-07). Two values are just outside 10% (cortical 5 uA 0.125 against 0.14, subcortical 5 uA 0.089 against 0.08). The recomputed p-values agree with the paper's cached ones (Pearson r = 0.992, same decision for 0.970 of 7980 tests), and the cached ones give the paper's 0.14, 0.21, 0.52, 0.79 and 0.08 to 0.35 on the dataset's units, with amplitude F = 172.99 and region F = 85.99 in the code's full model. The runs used hold only the condition's trials, except in jlh53concat, where a run number covers 150 trials of monopolar and bipolar pulses; with its own 75 cathodal trials the proportions are 0.121, 0.203, 0.514 and 0.756 (cortical) and 0.095 and 0.339 (subcortical at 5 and 100 uA), the same two values outside 10%.

fig6f-no-propagating-wave

Response timing shows no propagating wave: time to peak inhibition and to peak rebound excitation do not correlate with distance from the stimulation site.

As the code (figure7.ipynb cells 20 to 23): cortical good stable units that are ZETA responders to 50 uA cathodal pulses; MI in 5 ms bins over 0 to 500 ms (first 75 trials), Gaussian smoothing (sigma 2.5 ms), the most prominent peak after 20 ms (prominence 0.2), time by the local centre of mass. The code's 4 ms minimum width and 6 ms centre-of-mass half-width are converted to whole 5 ms bins (its _to_bins, rounding 0.8 and 1.2 bins to 1), so it applies no width filter (only above 1 bin) and takes the centre of mass over +-1 bin, and a unit without a prominent peak falls back to its largest bin after 20 ms. The criterion lists the width of at least 4 ms; applied, it changes the outcome, so the width filter and dropping units without a peak are run as variants, alone and together. A centre of mass over +-6 ms and over +-1 bin of 5 ms take the same bins. The distance is the code's: Euclidean distance to the first trial's cathodal contact, signed negative where the unit is dorsal to it (the code takes the anodal contact for jlh32; both shown). jlh39 and jlh49v1, which the code leaves out, have no cathodal pulses, so the two readings are the same.

from scipy.ndimage import gaussian_filter1d
from scipy.signal import find_peaks
BIN = 0.005
def pick_peaks(resp, sign, width=False, drop=False):
    """The code's pick_peak_times: its 4 ms width and 6 ms centre-of-mass half-width round to 1 bin of
    5 ms (_to_bins), so the width filter is skipped (applied only when above 1 bin) and the centre of
    mass spans +-1 bin; units without a prominent peak fall back to the largest bin after 20 ms.
    width: apply the criterion's width of at least 4 ms (0.8 bins); drop: leave out units without
    such a peak in place of the fallback."""
    sig = resp if sign > 0 else -resp
    out = np.full(len(resp), np.nan)
    minlat = int(np.floor(0.02 / BIN))
    for k, y0 in enumerate(sig):
        y = y0.copy()
        y[:minlat] = -np.inf
        pk, props = find_peaks(y, prominence=0.2, width=0.004 / BIN) if width else find_peaks(y, prominence=0.2)
        if len(pk) == 0:
            if drop:
                continue
            idx = int(np.nanargmax(y[minlat:]) + minlat)
        else:
            idx = int(pk[np.argmax(props["prominences"])])
        w = np.maximum(y0, 0)
        lo, hi = max(0, idx - 1), min(len(w), idx + 2)
        out[k] = (np.sum(np.arange(lo, hi) * w[lo:hi]) / w[lo:hi].sum() if w[lo:hi].sum() > 0 else idx) * BIN
    return out

RULES = {"code's rule": dict(), "width >= 4 ms, fallback kept": dict(width=True),
         "no width filter, units without a peak dropped": dict(drop=True),
         "criterion as written (width >= 4 ms, units without a peak dropped)": dict(width=True, drop=True)}
z50 = Z[(Z.pol == "cathodal") & (Z.amp == 50) & Z.sig]
for variant in ["code (jlh32 anodal contact)", "cathodal contact for all"]:
    mis, dists, sess = [], [], []
    for s in CAT_SESS:
        u, tr, trains, meta = DATA[s]
        st = estim.condition_trials(tr, "cathodal", 50).start_time.to_numpy()
        f0 = tr.iloc[0]
        stim = np.array([f0.neg_ap, f0.neg_dv, f0.neg_ml])
        if s == "jlh32" and variant.startswith("code"):
            f = estim.open_asset([a for a in ASSETS if a["session"] == s][0]["asset_id"])
            order = np.argsort(f["intervals/trials/start_time"][:], kind="stable")
            stim = f["intervals/trials/contact_positive_coords"][:][order][0]
            f.close()
        ids = z50[z50.session == s].unit.to_numpy()
        ids = [i for i in ids if "area" in u.region[i].lower()]
        for i in ids:
            mis.append(estim.modulation_index(trains[i], st, 0, 0.5, BIN))
            xyz = np.array([u.ap[i], u.dv[i], u.ml[i]])
            dists.append(np.linalg.norm(xyz - stim) * (-1 if xyz[1] < stim[1] else 1))
            sess.append(s)
    sm = gaussian_filter1d(np.array(mis), sigma=0.0025 / BIN, axis=1, mode="nearest")
    d = np.array(dists)
    print(f"{variant}: units {len(d)} by recording {pd.Series(sess).value_counts().to_dict()}; signed distance {d.min():.0f} to {d.max():.0f} um")
    for rule, kw in RULES.items():
        t_exc, t_inh = pick_peaks(sm, 1, **kw), pick_peaks(sm, -1, **kw)
        for name, tt in [("excitatory", t_exc), ("inhibitory", t_inh)]:
            ok = np.isfinite(tt) & np.isfinite(d)
            r, p = stats.pearsonr(d[ok], tt[ok])
            print(f"   {rule}, time to peak {name}: median {np.median(tt[ok]) * 1000:.0f} ms, Pearson r = {r:.3f}, p = {p:.3f}, n = {ok.sum()} ({len(d) - ok.sum()} dropped)")
print("Paper: Pearson r = -0.03 (p = 0.51) and 0.09 (p = 0.64)")

code (jlh32 anodal contact): units 614 by recording {'jlh31': 174, 'jlh32': 150, 'jlh34': 113, 'jlh33': 61, 'jlh40': 53, 'jlh48': 39, 'jlh53concat': 24}; signed distance -1218 to 1396 um
   code's rule, time to peak excitatory: median 205 ms, Pearson r = -0.001, p = 0.987, n = 614 (0 dropped)
   code's rule, time to peak inhibitory: median 70 ms, Pearson r = 0.047, p = 0.244, n = 614 (0 dropped)
   width >= 4 ms, fallback kept, time to peak excitatory: median 215 ms, Pearson r = 0.131, p = 0.001, n = 614 (0 dropped)
   width >= 4 ms, fallback kept, time to peak inhibitory: median 115 ms, Pearson r = -0.115, p = 0.004, n = 614 (0 dropped)
   no width filter, units without a peak dropped, time to peak excitatory: median 205 ms, Pearson r = 0.002, p = 0.969, n = 595 (19 dropped)
   no width filter, units without a peak dropped, time to peak inhibitory: median 70 ms, Pearson r = 0.050, p = 0.219, n = 595 (19 dropped)

   criterion as written (width >= 4 ms, units without a peak dropped), time to peak excitatory: median 220 ms, Pearson r = 0.153, p = 0.000, n = 595 (19 dropped)
   criterion as written (width >= 4 ms, units without a peak dropped), time to peak inhibitory: median 125 ms, Pearson r = -0.121, p = 0.004, n = 580 (34 dropped)

cathodal contact for all: units 614 by recording {'jlh31': 174, 'jlh32': 150, 'jlh34': 113, 'jlh33': 61, 'jlh40': 53, 'jlh48': 39, 'jlh53concat': 24}; signed distance -1218 to 1396 um
   code's rule, time to peak excitatory: median 205 ms, Pearson r = 0.004, p = 0.918, n = 614 (0 dropped)
   code's rule, time to peak inhibitory: median 70 ms, Pearson r = 0.047, p = 0.240, n = 614 (0 dropped)
   width >= 4 ms, fallback kept, time to peak excitatory: median 215 ms, Pearson r = 0.137, p = 0.001, n = 614 (0 dropped)
   width >= 4 ms, fallback kept, time to peak inhibitory: median 115 ms, Pearson r = -0.118, p = 0.004, n = 614 (0 dropped)
   no width filter, units without a peak dropped, time to peak excitatory: median 205 ms, Pearson r = 0.008, p = 0.853, n = 595 (19 dropped)
   no width filter, units without a peak dropped, time to peak inhibitory: median 70 ms, Pearson r = 0.051, p = 0.211, n = 595 (19 dropped)

   criterion as written (width >= 4 ms, units without a peak dropped), time to peak excitatory: median 220 ms, Pearson r = 0.161, p = 0.000, n = 595 (19 dropped)
   criterion as written (width >= 4 ms, units without a peak dropped), time to peak inhibitory: median 125 ms, Pearson r = -0.121, p = 0.003, n = 580 (34 dropped)
Paper: Pearson r = -0.03 (p = 0.51) and 0.09 (p = 0.64)

Result: reproduced, resting on an assumption. With the code's peak rule, among 614 cortical circuit responders at 50 uA, time to peak excitation (Pearson r = -0.001, p = 0.987) and to peak inhibition (r = 0.047, p = 0.244) do not correlate with signed distance, with either contact for jlh32; dropping the 19 units without a prominent peak gives r = 0.002 and 0.050 (p = 0.969 and 0.219). With the 4 ms width filter applied, as the criterion lists it, both correlate: r = 0.131 (p = 0.001) and -0.115 (p = 0.004), and r = 0.153 and -0.121 (p < 0.001 and 0.004) when units without such a peak are dropped. The outcome rests on the code's rule, which skips the width filter.

figs9-state-drp

Direct response probability is slightly but significantly higher on high than on low cortical-state trials at every cathodal amplitude.

Trials (first 75 cathodal) split at the median of the population spike count in the 500 ms before onset (high: above the median); DRP or nDRP per unit within each half; paired Wilcoxon over units, FDR over amplitudes. Readings: population (a) all good stable units, (b) cortical ones; measure (i) DRP, (ii) nDRP; units (x) all good stable units, (y) the 99 direct responders.

rows = []
for s in CAT_SESS:
    u, tr, trains, meta = DATA[s]
    idx = np.where(u.keep)[0]
    for popname, pidx in [("a", idx), ("b", np.where(u.keep & (u.sub_cort == "Cortical"))[0])]:
        pop = np.sort(np.concatenate([trains[i] for i in pidx]))
        for a in estim.AMPS:
            st = estim.condition_trials(tr, "cathodal", a).start_time.to_numpy()
            cnt = np.searchsorted(pop, st) - np.searchsorted(pop, st - 0.5)
            hi = cnt > np.median(cnt)
            for i in idx:
                for state, m in [("high", hi), ("low", ~hi)]:
                    p, pb, n = estim.spike_prob(trains[i], st[m])
                    rows.append((s, i, popname, a, state, p, n))
ST = pd.DataFrame(rows, columns=["session", "unit", "pop", "amp", "state", "drp", "ndrp"])
resp_set = set(map(tuple, R99.index))
for popname in ["a", "b"]:
    for meas in ["drp", "ndrp"]:
        for uset in ["x", "y"]:
            d = ST[ST["pop"] == popname]
            if uset == "y":
                d = d[[k in resp_set for k in zip(d.session, d.unit)]]
            w = d.pivot_table(index=["session", "unit", "amp"], columns="state", values=meas).dropna().reset_index()
            res, ps = {}, []
            for a in estim.AMPS:
                x = w[w.amp == a]
                ps.append(stats.wilcoxon(x.high, x.low).pvalue if (x.high != x.low).any() else 1.0)
                res[a] = f"low {x.low.mean():.4f}, high {x.high.mean():.4f}"
            adj = multipletests(ps, method="fdr_bh")[1]
            print(f"population {popname}, {meas}, units {uset} (n={w.groupby('amp').size().iloc[0]}):", res, "| FDR p", np.round(adj, 4).tolist(),
                  "| high above low at", [a for a in estim.AMPS if w[w.amp == a].high.mean() > w[w.amp == a].low.mean()])
print("Expected DRP low 0.051, high 0.056 at 5 uA; low 0.104, high 0.107 at 100 uA (10%); paper p < 0.001, < 0.001, 0.010, 0.030")

population a, drp, units x (n=2056): {5: 'low 0.0120, high 0.0124', 25: 'low 0.0126, high 0.0143', 50: 'low 0.0150, high 0.0159', 100: 'low 0.0159, high 0.0165'} | FDR p [0.0, 0.0, 0.0003, 0.0016] | high above low at [5, 25, 50, 100]
population a, drp, units y (n=99): {5: 'low 0.1328, high 0.1202', 25: 'low 0.1492, high 0.1542', 50: 'low 0.2023, high 0.2128', 100: 'low 0.2441, high 0.2463'} | FDR p [0.3008, 0.3601, 0.3008, 0.9463] | high above low at [25, 50, 100]
population a, ndrp, units x (n=2056): {5: 'low 0.0053, high 0.0043', 25: 'low 0.0054, high 0.0049', 50: 'low 0.0072, high 0.0057', 100: 'low 0.0081, high 0.0065'} | FDR p [0.0336, 0.0336, 0.0, 0.0] | high above low at []
population a, ndrp, units y (n=99): {5: 'low 0.1253, high 0.1068', 25: 'low 0.1411, high 0.1422', 50: 'low 0.1893, high 0.1982', 100: 'low 0.2329, high 0.2325'} | FDR p [0.0382, 0.7443, 0.7443, 0.7443] | high above low at [25, 50]
population b, drp, units x (n=2054): {5: 'low 0.0118, high 0.0127', 25: 'low 0.0131, high 0.0138', 50: 'low 0.0152, high 0.0157', 100: 'low 0.0160, high 0.0163'} | FDR p [0.0, 0.0, 0.0003, 0.0156] | high above low at [5, 25, 50, 100]
population b, drp, units y (n=99): {5: 'low 0.1267, high 0.1264', 25: 'low 0.1531, high 0.1501', 50: 'low 0.2088, high 0.2061', 100: 'low 0.2489, high 0.2410'} | FDR p [0.7679, 0.7679, 0.7679, 0.7679] | high above low at []
population b, ndrp, units x (n=2054): {5: 'low 0.0049, high 0.0047', 25: 'low 0.0058, high 0.0046', 50: 'low 0.0072, high 0.0057', 100: 'low 0.0079, high 0.0067'} | FDR p [0.1843, 0.0042, 0.0002, 0.0002] | high above low at []
population b, ndrp, units y (n=99): {5: 'low 0.1184, high 0.1139', 25: 'low 0.1465, high 0.1367', 50: 'low 0.1966, high 0.1907', 100: 'low 0.2372, high 0.2277'} | FDR p [0.7557, 0.7557, 0.7557, 0.4155] | high above low at []
Expected DRP low 0.051, high 0.056 at 5 uA; low 0.104, high 0.107 at 100 uA (10%); paper p < 0.001, < 0.001, 0.010, 0.030

Result: partially reproduced. With all good stable units, the whole population and DRP (readings a, i, x), DRP is higher on high-state trials at every amplitude (0.0124 against 0.0120 at 5 uA, 0.0165 against 0.0159 at 100 uA; FDR p 0.0016 or less); the cortical population gives the same (FDR p up to 0.0156). The values are about a quarter of the paper's at 5 uA (0.051 and 0.056) and a sixth at 100 uA (0.104 and 0.107). Among the 99 direct responders no difference is significant (FDR p 0.3008 to 0.9463), and with nDRP over all units the low state is higher at every amplitude (FDR p 0.0336 or less).

fig6g-dimensionality-collapse

After a pulse the population's dimensionality drops (participation ratio and effective dimensionality fall), with a brief rise in effective dimensionality within 20 ms at 5, 25 and 50 uA.

Per recording, cortical good stable units (CCF name with "area"; the code's "layer" filter as a variant), first 75 cathodal trials, spike counts in bins from -200 to +500 ms; per bin the participation ratio of the neuron covariance across trials and the number of components for 95% of the variance; each recording's trace divided by its median over pre-stimulus bins, then averaged over recordings. Epochs as the code: middle 12 to 150 ms. Bins 2 ms (Methods, code) and 5 ms (legend).

DIM = {}
for unit_filter in ["area", "layer"]:
    for bw in [0.002, 0.005]:
        edges = np.round(np.arange(-0.2, 0.5 + bw / 2, bw), 6)
        ctr = edges[:-1] + bw / 2
        res = {}
        for a in estim.AMPS:
            prs, eds = [], []
            for s in CAT_SESS:
                u, tr, trains, meta = DATA[s]
                idx = np.where(u.keep & u.region.str.contains(unit_filter))[0]
                st = estim.condition_trials(tr, "cathodal", a).start_time.to_numpy()
                C = np.stack([estim.counts_matrix(trains[i], st, edges) for i in idx], axis=2)  # trials x bins x units
                pr, ed = np.array([estim.participation_ratio(C[:, b, :]) for b in range(C.shape[1])]).T
                pre = ctr < 0
                prs.append(pr / np.nanmedian(pr[pre]))
                eds.append(ed / np.nanmedian(ed[pre]))
            res[a] = (np.nanmean(prs, 0), np.nanmean(eds, 0))
        DIM[(unit_filter, bw)] = (ctr, res)
        mid = (ctr >= 0.012) & (ctr <= 0.150)
        early = (ctr >= 0) & (ctr <= 0.020)
        print(f"units '{unit_filter}', {bw * 1000:.0f} ms bins:")
        for a in estim.AMPS:
            pr, ed = res[a]
            print(f"   {a} uA: middle-epoch PR {np.nanmean(pr[mid]):.3f}, ED {np.nanmean(ed[mid]):.3f} | 0-20 ms ED mean {np.nanmean(ed[early]):.3f}, max {np.nanmax(ed[early]):.3f}"
                  f" | late 170-300 ms PR {np.nanmean(pr[(ctr >= 0.17) & (ctr <= 0.3)]):.3f}")
ctr, res = DIM[("area", 0.002)]
pre = ctr < -0.004
print("\nPre-stimulus spread of the normalised traces (2 ms bins, 'area' units, -200 to -4 ms: firing falls in the 4 ms before the logged onset,"
      " see the data check), the reference for the 0 to 20 ms rise:")
print("   ED in the two bins before onset:", {a: np.round(res[a][1][(ctr < 0) & (ctr > -0.004)], 3).tolist() for a in estim.AMPS})
for a in estim.AMPS:
    pr, ed = res[a]
    print(f"   {a} uA: ED pre-stimulus min {np.nanmin(ed[pre]):.3f}, max {np.nanmax(ed[pre]):.3f}, 99th percentile {np.nanpercentile(ed[pre], 99):.3f};"
          f" bins of 0 to 20 ms above the pre-stimulus max: {int((ed[(ctr >= 0) & (ctr <= 0.020)] > np.nanmax(ed[pre])).sum())}"
          f" | PR pre-stimulus min {np.nanmin(pr[pre]):.3f}, max {np.nanmax(pr[pre]):.3f}")
fig, ax = plt.subplots(1, 2, figsize=(9, 2.5))
for a in estim.AMPS:
    ax[0].plot(ctr * 1000, res[a][0], label=f"{a} uA"); ax[1].plot(ctr * 1000, res[a][1])
ax[0].set_ylabel("normalised PR"); ax[1].set_ylabel("normalised ED"); ax[0].legend(fontsize=6)
for x in ax: x.set_xlabel("time from onset (ms)"); x.axhline(1, color="k", lw=0.5)
plt.tight_layout(); plt.show()
print("Plotted ranges: PR %.2f to %.2f, ED %.2f to %.2f" % (min(np.nanmin(res[a][0]) for a in estim.AMPS), max(np.nanmax(res[a][0]) for a in estim.AMPS),
                                                       min(np.nanmin(res[a][1]) for a in estim.AMPS), max(np.nanmax(res[a][1]) for a in estim.AMPS)))

units 'area', 2 ms bins:
   5 uA: middle-epoch PR 0.982, ED 0.986 | 0-20 ms ED mean 0.929, max 1.024 | late 170-300 ms PR 1.033
   25 uA: middle-epoch PR 0.951, ED 0.961 | 0-20 ms ED mean 0.994, max 1.087 | late 170-300 ms PR 1.013
   50 uA: middle-epoch PR 0.765, ED 0.778 | 0-20 ms ED mean 0.960, max 1.094 | late 170-300 ms PR 0.990
   100 uA: middle-epoch PR 0.459, ED 0.476 | 0-20 ms ED mean 0.884, max 1.043 | late 170-300 ms PR 0.979

units 'area', 5 ms bins:
   5 uA: middle-epoch PR 0.966, ED 0.973 | 0-20 ms ED mean 0.954, max 0.977 | late 170-300 ms PR 1.019
   25 uA: middle-epoch PR 0.959, ED 0.974 | 0-20 ms ED mean 1.010, max 1.038 | late 170-300 ms PR 0.998
   50 uA: middle-epoch PR 0.781, ED 0.813 | 0-20 ms ED mean 0.976, max 1.021 | late 170-300 ms PR 0.946
   100 uA: middle-epoch PR 0.478, ED 0.515 | 0-20 ms ED mean 0.894, max 0.984 | late 170-300 ms PR 0.872

units 'layer', 2 ms bins:
   5 uA: middle-epoch PR 0.983, ED 0.988 | 0-20 ms ED mean 0.931, max 1.015 | late 170-300 ms PR 1.030
   25 uA: middle-epoch PR 0.951, ED 0.963 | 0-20 ms ED mean 1.000, max 1.107 | late 170-300 ms PR 1.010
   50 uA: middle-epoch PR 0.773, ED 0.787 | 0-20 ms ED mean 0.954, max 1.079 | late 170-300 ms PR 0.991
   100 uA: middle-epoch PR 0.484, ED 0.505 | 0-20 ms ED mean 0.886, max 1.048 | late 170-300 ms PR 0.974

units 'layer', 5 ms bins:
   5 uA: middle-epoch PR 0.967, ED 0.972 | 0-20 ms ED mean 0.958, max 0.972 | late 170-300 ms PR 1.014
   25 uA: middle-epoch PR 0.959, ED 0.976 | 0-20 ms ED mean 1.016, max 1.045 | late 170-300 ms PR 0.996
   50 uA: middle-epoch PR 0.789, ED 0.826 | 0-20 ms ED mean 0.979, max 1.025 | late 170-300 ms PR 0.950
   100 uA: middle-epoch PR 0.504, ED 0.541 | 0-20 ms ED mean 0.894, max 0.986 | late 170-300 ms PR 0.870

Pre-stimulus spread of the normalised traces (2 ms bins, 'area' units, -200 to -4 ms: firing falls in the 4 ms before the logged onset, see the data check), the reference for the 0 to 20 ms rise:
   ED in the two bins before onset: {5: [0.805, 0.683], 25: [0.823, 0.574], 50: [0.801, 0.475], 100: [0.658, 0.333]}
   5 uA: ED pre-stimulus min 0.877, max 1.079, 99th percentile 1.078; bins of 0 to 20 ms above the pre-stimulus max: 0 | PR pre-stimulus min 0.856, max 1.087
   25 uA: ED pre-stimulus min 0.923, max 1.056, 99th percentile 1.047; bins of 0 to 20 ms above the pre-stimulus max: 5 | PR pre-stimulus min 0.911, max 1.093
   50 uA: ED pre-stimulus min 0.930, max 1.091, 99th percentile 1.074; bins of 0 to 20 ms above the pre-stimulus max: 3 | PR pre-stimulus min 0.917, max 1.098
   100 uA: ED pre-stimulus min 0.910, max 1.054, 99th percentile 1.050; bins of 0 to 20 ms above the pre-stimulus max: 0 | PR pre-stimulus min 0.913, max 1.091

Plotted ranges: PR 0.30 to 1.13, ED 0.33 to 1.13

Result: reproduced. With 2 ms bins, the normalised participation ratio in 12 to 150 ms is 0.982, 0.951, 0.765 and 0.459 and the effective dimensionality 0.986, 0.961, 0.778 and 0.476 at 5 to 100 uA, all below 1. Within 0 to 20 ms the largest bin of effective dimensionality is 1.024, 1.087 and 1.094 at 5, 25 and 50 uA (reading b; also 1.043 at 100 uA), while its mean over the window is below 1 (reading a). With 5 ms bins the largest bin exceeds 1 only at 25 and 50 uA. Against the pre-stimulus spread (-200 to -4 ms; normalised effective dimensionality up to 1.079, 1.056, 1.091 and 1.054), the largest bin in 0 to 20 ms exceeds the pre-stimulus maximum at 25 uA (5 bins) and 50 uA (3 bins), not at 5 uA (1.024 against 1.079) or 100 uA. The criterion does not compare with the pre-stimulus spread.

Summary

Claim	Outcome
fig4d-latency-trough	reproduced
figs3cd-gmm-window	reproduced
fig4f-responders-within-ep	not reproduced (assumption: the code's FS/RS classes; outside input: evoked-potential boundaries)
fig4h-extent-constant	reproduced
fig4h-distance-decreases	partially reproduced
fig4h-ks-distinct	not reproduced
fig4i-proportion-within-ep	reproduced (outside input: evoked-potential boundaries)
fig4j-drp-increases	reproduced
figs3e-early-windows	not reproduced
fig5bc-amplitude-clusters	reproduced
fig5ef-polarity-clusters	reproduced
fig5-fsrs-proportion	not reproduced
fig5j-fs-closer	reproduced (assumption: the code's FS/RS classes)
fig5i-fs-within-ep	not reproduced (assumption: the code's FS/RS classes; outside input: evoked-potential boundaries)
figs6b-fsrs-early-firing	reproduced (assumption: the code's FS/RS classes)
fig6c-modulation-cortex	not reproduced
fig6d-hva-subcortical	not reproduced
fig6e-circuit-responders	partially reproduced
fig6f-no-propagating-wave	reproduced (assumption: the code's peak rule, without the width filter)
figs9-state-drp	partially reproduced
fig6g-dimensionality-collapse	reproduced
fig1l, figs1b, fig2g (2), fig2-extent-polarity, fig2h, figs2gh, fig3c, fig3e, fig3fg, fig3h (11 claims)	blocked: the raw AP band is not in the dandiset
fig7d, fig7e, fig7g, fig7h (2), figs10c, fig7i, fig7j, no-learning-across-sessions, figs10e (10 claims)	blocked: the behaviour is not in the dandiset

A failure to reproduce here may be a failure of this attempt, not of the paper: see the reasoning and the differences from the paper's method.

How the analysis differs from the paper's. Units, trials, the stability filter and the direct response follow the paper and its code; the unit counts per condition match the code's intermediate table direct_df.csv for the cathodal and anodal conditions, and differ by 1 to 5 units for bipolar. Distances use the dataset's coordinates; the paper's code uses its own, which differ for the stimulation contacts of jlh31 and jlh32, and the distance claims are also shown with the code's coordinates, which give the paper's statistics and most of its values (not the 50 uA extent, 868 against 965 um, or the 50 against 100 uA KS p, 0.8245 against 0.083). The signed distance needs channel coordinates that the dandiset does not give; they are fitted from the units' coordinates. The FS and RS classes follow the paper's code (a 0.4 ms duration threshold), where the Methods describe k-means; both are shown; the k-means classes change with the initialisation. The evoked-potential boundaries of fig4f, fig4i and fig5i are the code's, an outside input. ZETA uses zetapy 2.7.2, whose getZeta the code calls, and only the cathodal conditions are tested. The claims on cathodal pulses rest on seven recordings: two of the nine hold no cathodal pulses.

What was found about the dataset. The dandiset holds the spikes, templates, CCF coordinates and regions of the nine recordings the paper analyses, and a tenth (jlh41) without coordinates that the paper does not use. The raw AP band, LFP and behaviour are not deposited, so 21 of the 42 claims cannot be tested, and three more only with the boundaries the paper's code computed. The file of session jlh49v1 is filed under subject jlh48. The stimulation contact coordinates of jlh31 and jlh32 differ from those the paper's code uses by 206 and 373 um, which changes the outcome of the distance claims. The electrodes table holds placeholder locations and coordinates. jlh33 has 60 trials per condition, where the Methods state 75. Neither FS/RS labels nor the stability flags the paper uses are in the files; both are recomputed here. In jlh53concat the run numbers of monopolar and bipolar runs coincide. Population firing falls in the 3 ms before the logged onset of the pulses in every recording.