Reproducing Pasquereau & Turner (2013) from dandiset 001636¶
Paper: Pasquereau B, Turner RS (2013). Primary motor cortex of the parkinsonian monkey: altered neuronal responses to muscle stretch. Frontiers in Systems Neuroscience 7:98. doi:10.3389/fnsys.2013.00098. The claims file was written from the full text (open access, PubMed Central PMC3840326) and its supplementary data sheet, whose table gives the neuronal results of each monkey separately.
Dandiset: DANDI:001636, version draft, as asked
for this attempt. The draft was last modified on 2026-10-21 at 20:16:41 UTC (the first cell prints
the modification time read at run time); a draft can change, and when it was read the archive
listed 4 uploads in progress. Between 23:00 and 23:10 UTC on the same day every asset of the draft
was uploaded again; this run read the draft after that (the first cell prints the time), and every
number recorded from the earlier state (the attempt record and the conclusions) is unchanged. This is a repeat attempt: the earlier one,
001636-pasquereau-2013-01, used the published version 0.260518.1819. Since then the files were
reconverted to fix what that attempt found, and the fixes are not yet published, so the draft is
tested here. The paper (2013) cites no DANDI version.
How the draft differs from 0.260518.1819, in what the analysis uses (shown in the data check): it
holds 439 files where the published version held 447. Source files recorded on the same day and
linked through a neuron recorded across them are merged into one session, and a unit is observed
only in the trials of its own source files (its obs_intervals), so each unit's perturbations are
counted within those trials. The merge adds trials of second files that the published version did
not hold (for 12 monkey V units), and removes the four files that repeated another file's trials
and the two unit copies they held. The analog streams now share the clock of the trials, the
angle's reference frame says positive is flexion, the monkey L series are named WristAngle and
WristTorque, unit_name gives the source file and unit (v2701.1+v2702.2), and
unit_also_in_session_id is gone. The monkey L perturbation type is not_available where it was
none. Every file was replaced; the angle, torque command and triceps EMG values compared in one
file (v4527) while preparing are unchanged (not run here).
What can be tested. The files of monkey V (Venus, elbow) mark the torque perturbations, with the joint angle, the torque command and triceps EMG. The files of monkey L (Leu, wrist) mark none, so every claim is tested on monkey V alone, against the paper's monkey V values where it gives them (supplementary table) and against its pooled values otherwise.
| Claim | Tested |
|---|---|
cell-counts: 112 and 159 monkey V cells, 49 and 54 PTNs, 31 and 48 CSNs, 3 dual |
yes |
monkey-l-perturbations: monkey L torque responses |
blocked (data missing) |
baseline-rate: lower resting rate after MPTP in PTNs, not CSNs |
yes |
ptn-faster-than-csn: PTNs fire faster than CSNs before MPTP |
yes |
antidromic-latency: CSN antidromic latency longer after MPTP, PTN unchanged |
yes |
false-positive-rate: the detection threshold's false-positive rate near 0.01 |
run, blocked (methods underspecified) |
torque-related-fractions: monkey V fractions of torque-related cells |
run, blocked (methods underspecified) |
prevalence-unchanged: MPTP does not change the fraction of torque-related cells |
run, blocked (methods underspecified) |
ptn-more-responsive: PTNs respond more often than CSNs |
run, blocked (methods underspecified) |
one-direction-half: 44% of torque-related cells respond to one direction only |
run, blocked (methods underspecified) |
excitatory-majority: most responses are increases; PTN balance unchanged |
run, blocked (methods underspecified) |
magnitude-unchanged: response magnitudes unchanged by MPTP |
run, blocked (methods underspecified) |
venus-response-magnitudes: monkey V response magnitudes |
run, blocked (methods underspecified) |
latency-shorter: response latencies shorter after MPTP |
run, blocked (methods underspecified) |
latency-not-training: latencies do not follow days of training |
run, blocked (methods underspecified) |
fwhm-shorter-ptn: PTN increases briefer after MPTP |
run, blocked (methods underspecified) |
dsi-reduced: fewer directional responses after MPTP |
run, blocked (methods underspecified) |
venus-dsi: monkey V mean DSI |
run, blocked (methods underspecified) |
dsi-not-baseline: resting rate does not differ between directional categories |
run, blocked (methods underspecified) |
kinematics-slower: torque-induced movements slower and smaller after MPTP |
yes |
kinematic-latency-unchanged: movement latency unchanged |
yes |
emg-m2-larger: larger long-latency (m2) triceps response after MPTP |
yes |
emg-duration-longer: longer EMG responses after MPTP |
yes |
emg-onset-unchanged: EMG onset unchanged |
yes |
emg-m1-weak: m1 much smaller than m2 |
yes |
velocity-matched-magnitude: magnitudes unchanged in velocity-matched sessions |
run, blocked (methods underspecified) |
session-correlations: no session-level correlation of PTN latency with velocity or m2 |
run, blocked (methods underspecified) |
trial-correlations: no trial-level correlation |
blocked (methods underspecified) |
The detection's false-positive rate (false-positive-rate) is 0.23 here with the detection fixed
in the claims file (0.46 with the t-test variant), far from the 0.01 the paper's threshold implies,
so the procedure differs from the paper's and false-positive-rate is blocked with
methods_underspecified. The claims file blocks the 15 claims that rest on the selection of
responses for the same reason when that rate falls outside 0.0033 to 0.03. Their analyses are run all the same, and their numbers are reported
for reference; none of them decides an outcome.
Code read. The paper links no analysis code. The dandiset links its conversion code,
catalystneuro/turner-lab-to-nwb. Its main
branch is still at commit e80a192917490b50db6871180b6d1d9a24772e35, read for the earlier attempt;
the code that built the draft is in an open pull request,
#33 (branch fix_errors, head
5f9bf6320a0160e4e56301d5409fb9242c1d8cf4), whose description was read: it merges linked source
files into sessions, reads the second file of each neuron recorded across files, drops repeated
trials, and lays the analog streams by trial. It holds no analysis of the torque responses, so
every method choice here comes from the paper's Methods or is fixed in the claims file.
The claims, with their criteria, are in claims/001636/10.3389-fnsys.2013.00098.yaml. The
loading and the measures are in stretch.py, whose docstring says which parameters come from the
paper and which were chosen here.
import time
import warnings
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import statsmodels.api as sm
from scipy import stats
from statsmodels.formula.api import ols
import stretch
pd.set_option("display.width", 200)
pd.set_option("display.max_columns", 40)
pd.set_option("display.max_rows", 200)
warnings.filterwarnings("ignore", category=RuntimeWarning)
started = time.time()
VERSION = stretch.version_info()
print(f"dandiset {stretch.DANDISET_ID}, version {VERSION['version']}, last modified {VERSION['modified']}, "
f"{VERSION['assets']} assets")
S = stretch.load_all(stretch.asset_urls())
print(f"loaded {len(S)} files of dandiset {stretch.DANDISET_ID} version {stretch.DANDISET_VERSION} "
f"in {time.time() - started:.0f} s")
dandiset 001636, version draft, last modified 2026-10-21 23:10:32.357921+00:00, 439 assets
loaded 439 files of dandiset 001636 version draft in 52 s
Data check¶
Each check of the metadata against the values is stated under its output.
sessions = pd.DataFrame([{
"si": si, "session": s["session_id"], "code": s["session_id"].split("++")[1],
"subject": s["subject"], "post": s["post"], "date": s["date"], "n_trials": len(s["trials"]),
"n_flex": int((s["pert_direction"] == "flexion").sum()),
"n_ext": int((s["pert_direction"] == "extension").sum()), "n_units": len(s["units"]),
"analog_dt_ms": s.get("analog_dt_ms", np.nan), "has_analog": "windows" in s or "analog_dt_ms" in s,
"angle_all_nan": s.get("angle_all_nan", True), "torque_all_nan": s.get("torque_all_nan", True),
"emg_channels": len(s.get("emg_channels", [])), "emg_same_clock": s.get("emg_same_clock", np.nan),
} for si, s in enumerate(S)])
units = pd.DataFrame([{
"si": si, "ui": ui, "session": s["session_id"], "subject": s["subject"], "post": s["post"],
"is_post_mptp": u["is_post_mptp"], "projection": u["type"], "ptype": stretch.short_type(u["type"]),
"unit_name": u["unit_name"], "also_in": u["also_in"], "date": s["date"],
"n_flex": int((u["pert_direction"] == "flexion").sum()),
"n_ext": int((u["pert_direction"] == "extension").sum()),
"n_flex_session": int((s["pert_direction"] == "flexion").sum()),
"n_ext_session": int((s["pert_direction"] == "extension").sum()),
"source_files": len(u["unit_name"].split("+")), "outside_obs": u.get("outside_obs", 0),
"antidromic_latency_ms": u["antidromic_latency_ms"], "n_spikes": len(u["spikes"]),
"isi_below_2ms": u["isi_below_2ms"], "outside_trials": u["outside_trials"], "sorted": u["sorted"],
} for si, s in enumerate(S) for ui, u in enumerate(s["units"])])
units["state"] = np.where(units.post, "post", "pre")
units["studied"] = (units.n_flex >= stretch.MIN_REPS) & (units.n_ext >= stretch.MIN_REPS)
print("Files and units by subject and MPTP state (from the session id):")
print(pd.concat([sessions.groupby(["subject", "post"]).size().rename("files"),
units.groupby(["subject", "post"]).size().rename("units")], axis=1).to_string())
print("\nUnits whose is_post_mptp differs from the state in the session id:",
int((units.post != units.is_post_mptp).sum()))
print("\nRecording dates by subject:")
print(sessions.groupby("subject").date.agg(["min", "max"]).to_string())
print("\nProjection types as stored, by subject and state:")
print(pd.crosstab(units.projection, [units.subject, units.state]).to_string())
print("\nperturbation type labels by subject:", {k: dict(v) for k, v in
pd.concat([s["trials"].assign(subject=s["subject"]) for s in S]).groupby("subject").torque.value_counts().unstack(fill_value=0).iterrows()})
print("Units recorded across several source files (unit_name lists them):",
units.groupby("subject").source_files.apply(lambda x: int((x > 1).sum())).to_dict(),
"| spikes outside the unit's obs_intervals:", int(units.outside_obs.sum()))
print("\nTorque perturbations per file (flexion / extension), by subject and state:")
print(sessions.groupby(["subject", "post"])[["n_flex", "n_ext"]].describe()
.loc[:, (slice(None), ["min", "50%", "max"])].to_string())
print("\nUnits studied (at least 5 perturbations of each direction within the unit's obs_intervals), by subject, state and type:")
print(pd.crosstab([units.subject, units.state], [units.studied, units.ptype]).to_string())
units["studied_session"] = (units.n_flex_session >= stretch.MIN_REPS) & (units.n_ext_session >= stretch.MIN_REPS)
print("Units whose selection changes if all perturbations of the session are counted:",
int((units.studied != units.studied_session).sum()))
Files and units by subject and MPTP state (from the session id):
files units
subject post
Leu False 98 138
True 50 72
Venus False 154 194
True 137 162
Units whose is_post_mptp differs from the state in the session id: 0
Recording dates by subject:
min max
subject
Leu 1992-08-06 1993-12-01
Venus 1999-04-06 2000-04-03
Projection types as stored, by subject and state:
subject Leu Venus
state post pre post pre
projection
corticostriatal_neuron 10 25 49 51
no_response 46 84 51 63
not_tested 0 0 3 3
pyramidal_tract_and_corticostriatal_neuron 0 0 1 2
pyramidal_tract_and_subthalamic_neuron 0 1 0 0
pyramidal_tract_and_thalamic_neuron 0 0 1 0
pyramidal_tract_neuron 10 19 55 75
subthalamic_and_corticostriatal_neuron 0 1 0 0
subthalamic_and_thalamic_neuron 0 1 0 0
subthalamic_projection_neuron 6 7 0 0
thalamic_projection_neuron 0 0 2 0
perturbation type labels by subject: {'Leu': {'extension': np.int64(0), 'flexion': np.int64(0), 'none': np.int64(0), 'not_available': np.int64(5952)}, 'Venus': {'extension': np.int64(2792), 'flexion': np.int64(2758), 'none': np.int64(2906), 'not_available': np.int64(0)}}
Units recorded across several source files (unit_name lists them): {'Leu': 1, 'Venus': 12} | spikes outside the unit's obs_intervals: 0
Torque perturbations per file (flexion / extension), by subject and state:
n_flex n_ext
min 50% max min 50% max
subject post
Leu False 0.0 0.0 0.0 0.0 0.0 0.0
True 0.0 0.0 0.0 0.0 0.0 0.0
Venus False 0.0 6.0 19.0 1.0 6.0 17.0
True 1.0 13.0 28.0 3.0 13.0 28.0
Units studied (at least 5 perturbations of each direction within the unit's obs_intervals), by subject, state and type:
studied False True
ptype CSN NA PTN CSN NA PTN dual
subject state
Leu post 10 52 10 0 0 0 0
pre 25 94 19 0 0 0 0
Venus post 1 3 1 48 54 54 1
pre 19 36 24 32 30 51 2
Units whose selection changes if all perturbations of the session are counted: 0
The draft holds 439 files with 566 units: 210 of monkey L (Leu) and 356 of monkey V (Venus). The
MPTP state in the session id agrees with each unit's is_post_mptp. Leu files date from August
1992 to December 1993, as the draft's description now says (earlier finding: the description said
1993 only; fixed).
Every Leu trial has the perturbation type not_available (it was none), so no monkey L trial
can be identified as perturbed. The monkey V units with at least 5 perturbations of each direction
within their own trials number 115 before and 157 after MPTP, with 51 and 54 PTNs and 32 and 48
CSNs. In the published version these counts were 112 and 159, 49 and 54, 31 and 48, the paper's
supplementary table. Counting all perturbations of the session for each unit selects the same
units. Thirteen units (12 of monkey V) are recorded across two source files; their spikes all fall
inside their obs_intervals. The other projection types (thalamic and subthalamic projections)
are counted as non-identified (NA) cells, as the paper has no such category, and the dual PTN and
CSN units are kept apart from both.
# The version used by the earlier attempt, units only: what the selection of cells studied gives
# there, and which units differ.
old = pd.DataFrame(stretch.load_all_brief(stretch.asset_urls(stretch.EARLIER_VERSION)))
old["studied"] = (old.n_flex >= stretch.MIN_REPS) & (old.n_ext >= stretch.MIN_REPS)
old["ptype"] = old.type.map(stretch.short_type)
old["state"] = np.where(old.post, "post", "pre")
print(f"Version {stretch.EARLIER_VERSION}: {old.path.nunique()} files, {len(old)} units; this version: "
f"{len(sessions)} files, {len(units)} units")
print(f"Files of {stretch.EARLIER_VERSION} not in this version:", sorted(set(old.path) - set(s["path"] for s in S)))
print(f"\nMonkey V cells studied in {stretch.EARLIER_VERSION}, by type and state:")
print(pd.crosstab(old[old.studied & (old.subject == "Venus")].ptype, old[old.studied & (old.subject == "Venus")].state)
.reindex(columns=["pre", "post"]).to_string())
# match each old unit to the unit of this version that lists it in its unit_name
part_of = {p: n for n in units.unit_name for p in n.split("+")}
old["now"] = old.unit_name.map(part_of)
now = units.set_index("unit_name")
diff = []
for r in old[old.subject == "Venus"].itertuples():
n = now.loc[r.now] if isinstance(r.now, str) else None
row = {"unit_then": r.unit_name, "type_then": r.ptype, "studied_then": r.studied,
"flex_ext_then": f"{r.n_flex}/{r.n_ext}", "unit_now": r.now if n is not None else "(absent)",
"type_now": n.ptype if n is not None else "", "studied_now": bool(n.studied) if n is not None else False,
"flex_ext_now": f"{n.n_flex}/{n.n_ext}" if n is not None else ""}
if (row["type_then"], row["studied_then"], row["flex_ext_then"]) != (row["type_now"], row["studied_now"], row["flex_ext_now"]):
diff.append(row)
print("\nMonkey V units whose type, perturbation counts or selection differ between the versions:")
print(pd.DataFrame(diff).to_string(index=False))
Version 0.260518.1819: 447 files, 568 units; this version: 439 files, 566 units Files of 0.260518.1819 not in this version: ['sub-Leu/sub-Leu_ses-Leu++l16301b++PostMPTP++Depth23000um++19931129_behavior+ecephys.nwb', 'sub-Venus/sub-Venus_ses-Venus++v2702++PreMPTP++Depth19650um++19990607_behavior+ecephys.nwb', 'sub-Venus/sub-Venus_ses-Venus++v3205++PreMPTP++Depth19400um++19990630_behavior+ecephys.nwb', 'sub-Venus/sub-Venus_ses-Venus++v3607b++PreMPTP++Depth19650um++19990729_behavior+ecephys.nwb', 'sub-Venus/sub-Venus_ses-Venus++v4527b++PostMPTP++Depth21600um++20000126_behavior+ecephys.nwb', 'sub-Venus/sub-Venus_ses-Venus++v4810b++PostMPTP++Depth18900um++20000204_behavior+ecephys.nwb', 'sub-Venus/sub-Venus_ses-Venus++v5604b++PostMPTP++Depth17960um++20000229_behavior+ecephys.nwb', 'sub-Venus/sub-Venus_ses-Venus++v5812++PostMPTP++Depth18300um++20000331_behavior+ecephys.nwb'] Monkey V cells studied in 0.260518.1819, by type and state: state pre post ptype CSN 31 48 NA 30 56 PTN 49 54 dual 2 1 Monkey V units whose type, perturbation counts or selection differ between the versions: unit_then type_then studied_then flex_ext_then unit_now type_now studied_now flex_ext_now v1501.1 CSN True 7/6 v1501.1+v1501b.1 CSN True 11/16 v1601.1 PTN True 6/7 v1601.1+v1601b.1 PTN True 11/15 v1602.1 PTN True 6/7 v1602.1+v1602b.1 PTN True 9/10 v1702.1 PTN True 6/7 v1702.1+v1702b.1 PTN True 12/17 v1703.1 CSN True 7/6 v1703.1+v1703b.1 CSN True 12/10 v1902.1 PTN False 10/4 v1902.1+v1902b.1 PTN True 17/11 v2701.1 CSN False 7/4 v2701.1+v2702.2 CSN True 16/10 v3202.1 CSN True 5/8 v3202.1+v3202b.1 CSN True 10/16 v3204.1 PTN True 8/8 v3204.1+v3205.1 PTN True 16/13 v3607.1 PTN False 10/1 v3607.1+v3607b.1 PTN True 19/7 v4527.2 NA True 16/16 v4527.2 PTN True 16/16 v4527b.2 PTN True 9/9 (absent) False v4603.1 PTN True 13/13 v4603.1+v4603b.1 PTN True 27/26 v4606.2 PTN True 12/12 v4606.2+v4607.1 PTN True 24/25 v5812.2 NA True 13/13 (absent) False
Against version 0.260518.1819, used by the earlier attempt: that version had 447 files and 568 units, and its selection gives the supplementary table's counts exactly (112 and 159 cells, 49 and 54 PTNs, 31 and 48 CSNs). The eight files it held that this version does not were second files of a neuron or repeats of another file's trials; their units are now in the sessions they continue. Twelve monkey V units gain the trials of their second file, which brings three of them (v1902.1, v2701.1 and v3607.1, two PTNs and a CSN, all before MPTP) to at least 5 perturbations of each direction. After MPTP, the two unit copies of repeated files (v4527b.2, a PTN, and v5812.2, labelled no_response) are gone, and v4527.2, labelled not_tested there, is now a PTN. These changes account for every difference in the counts of cells studied.
hold = []
for s in S:
if "hold_excursion_deg" not in s["trials"]:
continue
tr = s["trials"]
hold.append(pd.DataFrame({"subject": s["subject"], "post": s["post"], "perturbed": tr.torque.isin(["flexion", "extension"]),
"excursion": tr.hold_excursion_deg, "movement": tr.movement,
"active_change": tr.active_angle_change, "end_position": tr.end_position}))
hold = pd.concat(hold, ignore_index=True)
print("Largest angle excursion from its median, 1.0 to 2.5 s after the center target appears (deg):")
print(hold.groupby(["subject", "post", "perturbed"]).excursion.describe(percentiles=[.1, .5, .9]).round(2).to_string())
print("\nFraction of trials with an excursion above 2 deg:")
print(hold.assign(big=hold.excursion > 2).groupby(["subject", "post", "perturbed"]).big.mean().round(3).to_string())
leu = sessions[sessions.subject == "Leu"]
print(f"\nLeu files with no analog data: {int((~leu.has_analog).sum())} of {len(leu)} "
f"({int((~leu[leu.post].has_analog).sum())} of {int(leu.post.sum())} post-MPTP)")
print("Leu files with analog data whose torque command is all NaN, by state:")
print(leu[leu.has_analog].groupby("post").torque_all_nan.agg(["sum", "count"]).to_string())
print("Analog sample interval (ms) by subject:")
print(sessions.groupby("subject").analog_dt_ms.describe()[["count", "min", "50%", "max"]].round(4).to_string())
Largest angle excursion from its median, 1.0 to 2.5 s after the center target appears (deg):
count mean std min 10% 50% 90% max
subject post perturbed
Leu False False 2997.0 1.06 3.35 0.00 0.02 0.08 2.66 45.24
True False 2091.0 0.17 1.37 0.00 0.02 0.06 0.16 39.62
Venus False False 1003.0 0.19 0.63 0.01 0.01 0.03 0.47 17.27
True 1947.0 2.83 1.28 0.35 1.29 2.79 4.42 19.56
True False 1807.0 0.89 1.00 0.01 0.07 0.57 2.18 12.34
True 3603.0 3.32 1.45 0.11 1.72 3.19 4.79 31.75
Fraction of trials with an excursion above 2 deg:
subject post perturbed
Leu False False 0.148
True False 0.011
Venus False False 0.007
True 0.671
True False 0.107
True 0.865
Leu files with no analog data: 15 of 148 (15 of 50 post-MPTP)
Leu files with analog data whose torque command is all NaN, by state:
sum count
post
False 98 98
True 1 35
Analog sample interval (ms) by subject:
count min 50% max
subject
Leu 133.0 10.0 10.0 10.0
Venus 291.0 1.0 1.0 1.0
In monkey V the joint moves by more than 2 deg in the hold window in 67% (pre) and 86.5% (post) of
perturbed trials and in 1% and 11% of the others. In monkey L it does so in 15% of trials before
MPTP and 1% after, which does not fit perturbations on two thirds of the trials as the Methods
describe; these windows were set from the monkey V onsets, so the monkey L perturbations could fall
outside them. Either way, the monkey L torque trials cannot be identified from the files, and
monkey-l-perturbations is blocked. The Leu torque command is NaN in every pre-MPTP file (98) and
in one post-MPTP file, and 15 post-MPTP Leu files have no analog data. Earlier finding (no monkey L
perturbations): still present in the data; the draft's description and the torque series now say
that the Leu source files hold no perturbation times. The Leu analog streams are sampled at
100 Hz, where the paper says 200 Hz for monkey L; monkey V is sampled at 1 kHz (the paper digitized
at 500 Hz; the conversion notes say the lab interpolated to 1 kHz).
# Every excursion of the torque command (more than 20 units from rest, lasting at least 10 ms),
# in every file with a torque command; a hold-period impulse starts in the hold window.
exc_rows = []
for s in S:
if "command_excursions" not in s:
continue
e = s["command_excursions"]
e = e[e.duration_ms >= 10]
h = e[e.in_hold]
exc_rows.append({"subject": s["subject"], "state": "post" if s["post"] else "pre", "code": s["session_id"].split("++")[1],
"trials": len(s["trials"]), "excursions": len(e),
"labelled_torque_trials": int(s["trials"].torque.isin(["flexion", "extension"]).sum()),
"trials_hold_excursion": h.trial.nunique(),
"trials_hold_40_70ms": h[(h.duration_ms >= 40) & (h.duration_ms <= 70)].trial.nunique(),
"hold_median_ms": h.duration_ms.median(), "hold_median_abs_peak": h.peak.abs().median(),
"all_median_from_center_s": e.from_center_s.median(),
"all_median_from_lateral_s": e.from_lateral_s.median()})
exc = pd.DataFrame(exc_rows)
print("Files with a torque command, by subject and state (trials with an excursion starting in the hold window,"
" 1.0 to 2.5 s after the center target and before the lateral target):")
print(exc.groupby(["subject", "state"]).agg(files=("code", "size"), trials=("trials", "sum"),
labelled_torque_trials=("labelled_torque_trials", "sum"), trials_hold_excursion=("trials_hold_excursion", "sum"),
trials_hold_40_70ms=("trials_hold_40_70ms", "sum")).to_string())
print("\nEvery monkey L file with a torque command:")
print(exc[exc.subject == "Leu"].drop(columns="subject").round(2).to_string(index=False))
Files with a torque command, by subject and state (trials with an excursion starting in the hold window, 1.0 to 2.5 s after the center target and before the lateral target):
files trials labelled_torque_trials trials_hold_excursion trials_hold_40_70ms
subject state
Leu post 34 2048 0 1 1
Venus post 137 5506 3603 3517 2571
pre 154 2950 1947 1882 1882
Every monkey L file with a torque command:
state code trials excursions labelled_torque_trials trials_hold_excursion trials_hold_40_70ms hold_median_ms hold_median_abs_peak all_median_from_center_s all_median_from_lateral_s
post l14001 61 348 0 0 0 NaN NaN 3.72 0.53
post l14002 60 354 0 0 0 NaN NaN 3.68 0.54
post l14101 60 352 0 0 0 NaN NaN 3.50 0.41
post l14102 60 348 0 0 0 NaN NaN 3.67 0.47
post l14103 60 351 0 0 0 NaN NaN 3.67 0.48
post l14201 60 354 0 0 0 NaN NaN 3.66 0.48
post l14501 85 379 0 0 0 NaN NaN 3.74 0.49
post l14601 60 350 0 0 0 NaN NaN 3.63 0.39
post l14602 61 348 0 0 0 NaN NaN 3.52 0.41
post l14701 62 297 0 0 0 NaN NaN 3.94 0.83
post l14801 88 348 0 0 0 NaN NaN 3.32 0.35
post l14802 18 88 0 0 0 NaN NaN 3.83 0.72
post l14803 60 212 0 0 0 NaN NaN 3.04 0.30
post l14804 60 237 0 0 0 NaN NaN 3.30 0.32
post l15001 62 226 0 0 0 NaN NaN 3.85 0.73
post l15601 60 407 0 0 0 NaN NaN 3.45 0.34
post l15602 37 257 0 0 0 NaN NaN 3.54 0.42
post l15701 60 410 0 0 0 NaN NaN 3.56 0.39
post l15702 62 410 0 0 0 NaN NaN 3.38 0.38
post l15703 60 396 0 0 0 NaN NaN 3.47 0.32
post l15801 60 391 0 0 0 NaN NaN 3.40 0.38
post l15901 61 397 0 0 0 NaN NaN 3.58 0.37
post l16001 60 381 0 0 0 NaN NaN 3.47 0.37
post l16002 60 392 0 0 0 NaN NaN 3.41 0.38
post l16101 60 385 0 0 0 NaN NaN 3.26 0.32
post l16102 43 283 0 0 0 NaN NaN 3.36 0.35
post l16103 47 295 0 1 1 70.0 41.12 3.20 0.34
post l16201 60 387 0 0 0 NaN NaN 3.39 0.31
post l16301 120 767 0 0 0 NaN NaN 3.29 0.36
post l16302 60 379 0 0 0 NaN NaN 3.12 0.38
post l16401 60 402 0 0 0 NaN NaN 3.54 0.36
post l16402 60 377 0 0 0 NaN NaN 3.23 0.34
post l16403 61 380 0 0 0 NaN NaN 3.28 0.34
post l16501 40 265 0 0 0 NaN NaN 3.35 0.34
The torque command shows the same. In monkey V a command excursion starts in the hold window in 1882 pre-MPTP and 3517 post-MPTP trials, against 1947 and 3603 labelled perturbations; before MPTP every one of them lasts 40 to 70 ms. In monkey L the command is NaN in every pre-MPTP file, and 34 post-MPTP files hold one (35 in the published version, where l16301b was a file of its own). Of their 2048 trials, one has an excursion starting in the hold window (70 ms, 41 units). Their 88 to 767 excursions per file start a median 0.30 to 0.83 s after the lateral target appears, during the active movement. No monkey L file holds a hold-period torque impulse.
venus = [s for s in S if s["subject"] == "Venus" and "windows" in s]
AX = np.arange(stretch.ANALOG_WINDOW_MS[0], stretch.ANALOG_WINDOW_MS[1] + 1)
def command_onset(w):
# first ms where the torque command leaves its pre-onset level by more than 6 SD (and 10 units)
base = np.nanmedian(w[:, :400], axis=1, keepdims=True)
sd = np.nanstd(w[:, :400], axis=1, keepdims=True)
dev = np.abs(w - base) > np.maximum(10, 6 * sd)
return np.where(dev.any(1), AX[np.argmax(dev, axis=1)], np.nan)
realigned = np.concatenate([command_onset(s["windows"]["torque_cmd"]) for s in venus])
stored = np.concatenate([command_onset(s["windows"]["torque_cmd_stored"]) for s in venus])
lag = np.concatenate([s["analog_lag_ms"][s["pert_trial"]] for s in venus])
ok = np.isfinite(stored) & np.isfinite(realigned)
print(f"Monkey V perturbations: {len(realigned)}; command step found in {int(np.isfinite(realigned).sum())}")
no_step = [s["session_id"].split("++")[1] for s in venus if not np.isfinite(command_onset(s["windows"]["torque_cmd"])).any()]
print("Monkey V sessions with perturbations but no command step in any of them:", no_step)
print("Torque command step relative to the labelled onset, realigned clock (ms), percentiles 1/50/99:",
np.nanpercentile(realigned, [1, 50, 99]))
print("Same on the stored analog clock (ms), percentiles 1/50/99:", np.nanpercentile(stored, [1, 50, 99]))
print(f"Stored analog clock behind the trials table at each perturbation: {lag.min():.0f} to {lag.max():.0f} ms; "
f"correlation of the stored-clock step time with this lag {np.corrcoef(stored[ok], lag[ok])[0, 1]:.3f}")
per_trial = np.concatenate([np.diff(s["analog_lag_ms"]) for s in S if "analog_lag_ms" in s])
print(f"Lag added per trial: median {np.median(per_trial):.1f} ms (range {per_trial.min():.1f} to {per_trial.max():.1f})")
print("EMG series share the angle's timestamps in every monkey V file:",
bool(sessions[sessions.subject == "Venus"].emg_same_clock.astype(bool).all()))
print("Angle series and their reference_frame, by subject (files):")
print(pd.Series([(s["subject"], s["joint"], s["angle_reference_frame"]) for s in S]).value_counts().to_string())
print("\nStored angle change over active movements (lateral target to reward), median deg:")
print(hold.groupby(["subject", "movement"]).active_change.median().round(1).to_string())
print("derived_end_position, median deg:")
print(hold.groupby(["subject", "movement"]).end_position.median().round(1).to_string())
chg = {d: np.concatenate([s["windows"]["angle"][s["pert_direction"] == d][:, AX == 100][:, 0]
- s["windows"]["angle"][s["pert_direction"] == d][:, (AX >= -50) & (AX < 0)].mean(1)
for s in venus]) for d in ["flexion", "extension"]}
for d, v in chg.items():
print(f"Angle change 100 ms after a {d} torque: median {np.nanmedian(v):.2f} deg, positive in "
f"{np.nanmean(v > 0):.1%} of {np.isfinite(v).sum()} perturbations")
Monkey V perturbations: 5550; command step found in 5499
Monkey V sessions with perturbations but no command step in any of them: ['v1401', 'v1402', 'v1501']
Torque command step relative to the labelled onset, realigned clock (ms), percentiles 1/50/99: [4. 4. 6.]
Same on the stored analog clock (ms), percentiles 1/50/99: [4. 4. 6.]
Stored analog clock behind the trials table at each perturbation: -0 to 0 ms; correlation of the stored-clock step time with this lag 0.012
Lag added per trial: median 0.0 ms (range -0.0 to 0.0)
EMG series share the angle's timestamps in every monkey V file: True
Angle series and their reference_frame, by subject (files):
(Venus, Elbow, Anatomical flexion/extension axis, 0 degrees = neutral position, positive = flexion) 291
(Leu, Wrist, Anatomical flexion/extension axis, 0 degrees = neutral position, positive = flexion) 133
(Leu, None, None) 15
Stored angle change over active movements (lateral target to reward), median deg:
subject movement
Leu extension -16.3
flexion 17.6
Venus extension -20.9
flexion 18.4
derived_end_position, median deg:
subject movement
Leu extension -18.6
flexion 20.1
Venus extension -20.1
flexion 19.1
Angle change 100 ms after a flexion torque: median 3.42 deg, positive in 99.8% of 2758 perturbations Angle change 100 ms after a extension torque: median -2.50 deg, positive in 1.0% of 2792 perturbations
Earlier finding (analog clock): fixed. In the published version the analog streams were laid end
to end by their sample count and fell behind the trials table by about 5 ms per trial, so the
torque command step appeared up to 108 ms before the labelled onset. In the draft each analog
segment starts at its trial's start: the lag stretch.realign would correct is 0 ms at every
perturbation, and on the stored clock the command steps 4 to 6 ms (1st to 99th percentile) after
the labelled onset, the same as after realignment. The EMG series share the angle's timestamps.
Three sessions (v1401, v1402, v1501) have labelled perturbations but no step in the torque command,
as the draft's torque description says; their 51 perturbations are kept (the command is not used
to select trials).
Earlier finding (reference frame): fixed. The angle's reference_frame now says positive is
flexion in every file with an angle, and the values agree: active flexion movements raise the
angle by about 18 deg in both monkeys and extension movements lower it, and flexion torques raise
it in 99.8% of perturbations. The monkey L series are now named WristAngle and WristTorque. The
torque labels agree with the EMG (the triceps, an extensor stretched by a flexion torque, responds
to flexion torques; see the EMG sections).
# Angle scale: active movements measured here (same filter as the torque-evoked movements, around
# derived_movement_onset_time) against the files' derived columns and target amplitudes.
act = pd.concat([s["trials"].assign(subject=s["subject"]) for s in S if "active_peak_velocity" in s["trials"]])
act = act[act.active_peak_velocity.notna() & act.derived_peak_velocity.notna()]
act["ratio_velocity"] = act.active_peak_velocity / act.derived_peak_velocity
act["ratio_amplitude"] = act.active_amplitude / act.derived_amplitude.abs()
print(act.groupby("subject").agg(movements=("ratio_velocity", "size"),
target_amplitude_median=("target_amplitude", "median"),
derived_amplitude_median=("derived_amplitude", lambda x: x.abs().median()),
amplitude_here_median=("active_amplitude", "median"),
derived_peak_velocity_median=("derived_peak_velocity", "median"),
peak_velocity_here_median=("active_peak_velocity", "median"),
ratio_velocity_median=("ratio_velocity", "median"),
ratio_amplitude_median=("ratio_amplitude", "median")).round(2).to_string())
for subj, g in act.groupby("subject"):
print(f"{subj}: Spearman rho of peak velocity here against derived_peak_velocity "
f"{stats.spearmanr(g.active_peak_velocity, g.derived_peak_velocity)[0]:.3f}; "
f"velocity ratio within 10% of 1 in {((g.ratio_velocity - 1).abs() <= 0.1).mean():.1%}")
movements target_amplitude_median derived_amplitude_median amplitude_here_median derived_peak_velocity_median peak_velocity_here_median ratio_velocity_median ratio_amplitude_median subject Leu 4794 20.0 19.18 19.86 100.53 101.40 1.00 1.01 Venus 7566 20.0 19.46 20.42 100.90 100.35 0.99 1.02 Leu: Spearman rho of peak velocity here against derived_peak_velocity 0.989; velocity ratio within 10% of 1 in 93.5% Venus: Spearman rho of peak velocity here against derived_peak_velocity 0.994; velocity ratio within 10% of 1 in 99.4%
The angle's scale agrees with the files' own kinematic columns. Active movements measured here,
with the filter used for the torque-evoked movements, have median peak velocities and amplitudes
within 2% of derived_peak_velocity and derived_movement_amplitude (median ratios 0.99 to 1.02;
Spearman rho 0.994 in monkey V and 0.989 in monkey L for peak velocity), and their median amplitude
(20.4 deg in monkey V) matches the 20 deg targets. The angle is consistent with degrees, as stated,
and the filter keeps the peak velocity of these movements.
# Spikes against the trials table: does the torque direction become decodable from single-cell
# spike counts only after the labelled onset? 20 ms windows, Mann-Whitney flexion vs extension.
edges = np.arange(-200, 201, 20)
studied_v = units[units.studied & (units.subject == "Venus")]
sig = np.zeros(len(edges) - 1)
for r in studied_v.itertuples():
s, u = S[r.si], S[r.si]["units"][r.ui]
counts = np.array([np.histogram(x, edges)[0] for x in u["pert"]])
f, e = counts[u["pert_direction"] == "flexion"], counts[u["pert_direction"] == "extension"]
for j in range(len(edges) - 1):
if np.ptp(np.r_[f[:, j], e[:, j]]) > 0:
sig[j] += stats.mannwhitneyu(f[:, j], e[:, j]).pvalue < 0.01
frac = sig / len(studied_v)
print(f"Fraction of the {len(studied_v)} monkey V cells studied whose spike counts differ between "
"flexion and extension trials (Mann-Whitney p < 0.01), by window (ms from labelled onset):")
print(pd.Series(frac.round(3), index=[f"{a} to {b}" for a, b in zip(edges[:-1], edges[1:])]).to_string())
Fraction of the 272 monkey V cells studied whose spike counts differ between flexion and extension trials (Mann-Whitney p < 0.01), by window (ms from labelled onset): -200 to -180 0.004 -180 to -160 0.000 -160 to -140 0.000 -140 to -120 0.000 -120 to -100 0.011 -100 to -80 0.000 -80 to -60 0.004 -60 to -40 0.000 -40 to -20 0.000 -20 to 0 0.004 0 to 20 0.011 20 to 40 0.070 40 to 60 0.092 60 to 80 0.044 80 to 100 0.062 100 to 120 0.048 120 to 140 0.037 140 to 160 0.037 160 to 180 0.051 180 to 200 0.022
Before the labelled onset, the fraction of cells whose counts differ between directions stays near the 1% expected by chance (at most 1.1%); it rises to 7.0% in the window from 20 to 40 ms and peaks at 9.2% from 40 to 60 ms. The spike times therefore share the trials table's clock to within the 20 ms windows, consistent with the 20 to 60 ms response latencies of the paper.
# Files that repeat another file's recording: same subject and date, and the trial durations of one
# are a contiguous block of the other's.
dur = {si: np.round((s["trials"].stop - s["trials"].start).values, 3) for si, s in enumerate(S)}
repeats = []
for (subj, date), g in sessions.groupby(["subject", "date"]):
for a in g.si:
for b in g.si:
if a != b and len(dur[b]) <= len(dur[a]) and (len(dur[b]) < len(dur[a]) or a < b):
n = len(dur[b])
if any(np.array_equal(dur[a][i:i + n], dur[b]) for i in range(len(dur[a]) - n + 1)):
repeats.append({"file": sessions.code[a], "repeated_in": sessions.code[b],
"trials": len(dur[a]), "repeated_trials": n,
"units_file": [u["unit_name"] for u in S[a]["units"]],
"units_repeat": [u["unit_name"] for u in S[b]["units"]],
"types_file": [stretch.short_type(u["type"]) for u in S[a]["units"]],
"types_repeat": [stretch.short_type(u["type"]) for u in S[b]["units"]]})
repeats = pd.DataFrame(repeats, columns=["file", "repeated_in", "trials", "repeated_trials", "units_file",
"units_repeat", "types_file", "types_repeat"])
repeated_files = set(repeats.repeated_in)
print("Files whose trials repeat a block of another file's trials:", "none" if repeats.empty else "")
print(repeats.to_string(index=False))
# Within a merged session: does a later source file repeat trials of an earlier one?
within = []
for si, s in enumerate(S):
tr = s["trials"]
if tr.source_file.nunique() < 2:
continue
key = list(zip(np.round((tr.stop - tr.start).values, 3), tr.torque, tr.movement))
files = tr.source_file.values
first_file = {}
for k, f in zip(key, files):
first_file.setdefault(k, f)
within.append({"session": sessions.code[si], "source_files": list(pd.unique(files)), "trials": len(tr),
"trials_repeating_an_earlier_file": sum(first_file[k] != f for k, f in zip(key, files))})
print("\nMerged sessions (several source files) and trials of a later file that repeat a trial of an earlier one "
"(same duration, torque and movement):")
print(pd.DataFrame(within).to_string(index=False))
def trial_spikes(s, u):
# spike times relative to the start of their trial, one array per trial
tr = s["trials"]
k = np.searchsorted(tr.start.values, u["spikes"], side="right") - 1
return [np.round(u["spikes"][k == i] - tr.start.values[i], 4) for i in range(len(tr))]
same_unit = []
for r in repeats.itertuples():
a = S[int(sessions[sessions.code == r.file].si.iloc[0])]
b = S[int(sessions[sessions.code == r.repeated_in].si.iloc[0])]
da, db = dur[int(sessions[sessions.code == r.file].si.iloc[0])], dur[int(sessions[sessions.code == r.repeated_in].si.iloc[0])]
off = [i for i in range(len(da) - len(db) + 1) if np.array_equal(da[i:i + len(db)], db)][0]
for ub in b["units"]:
tb = trial_spikes(b, ub)
for ua in a["units"]:
ta = trial_spikes(a, ua)[off:off + len(db)]
n_same = sum(np.array_equal(x, y) for x, y in zip(ta, tb))
same_unit.append({"file": r.file, "unit": ua["unit_name"], "type": ua["type"],
"repeat": r.repeated_in, "repeat_unit": ub["unit_name"], "repeat_type": ub["type"],
"identical_trials": f"{n_same}/{len(db)}",
"spikes_in_block": (sum(map(len, ta)), sum(map(len, tb)))})
same_unit = pd.DataFrame(same_unit, columns=["file", "unit", "type", "repeat", "repeat_unit", "repeat_type",
"identical_trials", "spikes_in_block"])
print("\nUnits of each repeated recording compared trial by trial over the repeated block:")
print(same_unit.to_string(index=False))
dup_units = same_unit[same_unit.identical_trials.str.split("/", regex=False).str[0].astype(int) > 0]
# the copy in the repeated file is the one dropped in the variant that counts each unit once
dup_keys = {(r.repeat, r.repeat_unit) for r in dup_units.itertuples()}
units["duplicate_copy"] = [(sess.split("++")[1], un) in dup_keys for sess, un in zip(units.session, units.unit_name)]
print("Unit copies in repeated files (dropped in the variant that counts each unit once):", sorted(dup_keys))
print("\nUnits with unit_also_in_session_id set:", int((units.also_in != "").sum()),
"| files with the unit_also_in_session_id column:", sum(s["has_also_in_column"] for s in S), "of", len(S))
# Units recorded across several source files: their firing rate over the trials of each file, as a
# check that the parts are one neuron (spike times are only inside the unit's obs_intervals).
rows = []
for r in units[units.source_files > 1].itertuples():
s, u = S[r.si], S[r.si]["units"][r.ui]
tr = s["trials"][u["observed_trials"]]
rates = {}
for f, g in tr.groupby("source_file", sort=False):
n = sum(int(((u["spikes"] >= a) & (u["spikes"] <= b)).sum()) for a, b in zip(g.start, g.stop))
rates[f] = round(n / float((g.stop - g.start).sum()), 2)
rows.append({"unit": r.unit_name, "subject": r.subject, "type": r.ptype, "studied": r.studied,
"rate_by_file (spikes/s)": rates})
merged_units = pd.DataFrame(rows)
print(merged_units.to_string(index=False))
print("\nSpike times sorted in every unit:", bool(units["sorted"].all()),
"| spikes outside the trials:", int(units.outside_trials.sum()))
print("Within-trial ISIs below 2 ms, fraction per unit: median "
f"{units.isi_below_2ms.median():.4f}, units above 1%: {int((units.isi_below_2ms > 0.01).sum())} of {len(units)} "
f"({int((units[units.studied & (units.subject == 'Venus')].isi_below_2ms > 0.01).sum())} of the cells studied)")
Files whose trials repeat a block of another file's trials: none
Empty DataFrame
Columns: [file, repeated_in, trials, repeated_trials, units_file, units_repeat, types_file, types_repeat]
Index: []
Merged sessions (several source files) and trials of a later file that repeat a trial of an earlier one (same duration, torque and movement):
session source_files trials trials_repeating_an_earlier_file
l16301 [l16301, l16301b] 120 1
v1501 [v1501, v1501b] 40 0
v1601 [v1601, v1601b] 40 0
v1602 [v1602, v1602b] 30 0
v1702 [v1702, v1702b] 40 0
v1703 [v1703, v1703b] 40 0
v1902 [v1902, v1902b] 40 0
v2701 [v2701, v2702] 38 0
v3202 [v3202, v3202b] 40 0
v3204 [v3204, v3205] 40 0
v3607 [v3607, v3607b] 40 0
v4603 [v4603, v4603b] 80 0
v4606 [v4606, v4607] 76 0
Units of each repeated recording compared trial by trial over the repeated block:
Empty DataFrame
Columns: [file, unit, type, repeat, repeat_unit, repeat_type, identical_trials, spikes_in_block]
Index: []
Unit copies in repeated files (dropped in the variant that counts each unit once): []
Units with unit_also_in_session_id set: 0 | files with the unit_also_in_session_id column: 0 of 439
unit subject type studied rate_by_file (spikes/s)
l16301.1+l16301b.1 Leu CSN False {'l16301': 10.91, 'l16301b': 6.42}
v1501.1+v1501b.1 Venus CSN True {'v1501': 0.52, 'v1501b': 0.93}
v1601.1+v1601b.1 Venus PTN True {'v1601': 16.08, 'v1601b': 18.54}
v1602.1+v1602b.1 Venus PTN True {'v1602': 7.75, 'v1602b': 7.76}
v1702.1+v1702b.1 Venus PTN True {'v1702': 27.86, 'v1702b': 28.23}
v1703.1+v1703b.1 Venus CSN True {'v1703': 6.6, 'v1703b': 7.29}
v1902.1+v1902b.1 Venus PTN True {'v1902': 41.66, 'v1902b': 60.66}
v2701.1+v2702.2 Venus CSN True {'v2701': 3.54, 'v2702': 5.73}
v3202.1+v3202b.1 Venus CSN True {'v3202': 2.81, 'v3202b': 9.12}
v3204.1+v3205.1 Venus PTN True {'v3204': 10.85, 'v3205': 8.87}
v3607.1+v3607b.1 Venus PTN True {'v3607': 19.33, 'v3607b': 17.01}
v4603.1+v4603b.1 Venus PTN True {'v4603': 17.87, 'v4603b': 17.64}
v4606.2+v4607.1 Venus PTN True {'v4606': 10.61, 'v4607': 6.03}
Spike times sorted in every unit: True | spikes outside the trials: 0
Within-trial ISIs below 2 ms, fraction per unit: median 0.0006, units above 1%: 68 of 566 (36 of the cells studied)
Earlier finding (repeated files and unit copies): fixed. No file repeats a block of another file's trials. The four repeated files of the published version (v4527b, v4810b, v5604b, v5812) are folded into the sessions they repeat, their repeated trials left out, and the two unit copies they held (v4527b.2 and v5812.2) are stored once, as v4527.2 and v5811.2, both PTNs; the files' notes say which projection label was kept and why. Within the 13 merged sessions, trials of a later file match no trial of an earlier one in monkey V; in l16301 one trial of 120 matches on duration and movement alone, which a coincidence explains as well as a repeat. The variants below that drop unit copies are kept from the earlier attempt; there is nothing left to drop, so they equal the main results.
Earlier finding (unit_also_in_session_id): fixed. The column is in none of the files. Units
recorded across two source files are now one unit, named after both (v2701.1+v2702.2) and
observed in the trials of both. Treating each as one cell is the dataset's statement, taken here as
an assumption; their rates over the trials of each file differ by less than a factor of 1.5 in 8 of
the 13, and by up to 3.2 (v3202.1+v3202b.1, 2.81 and 9.12 spikes/s). A variant of cell-counts
counts each from its first file only.
Every unit's spike times are sorted and fall within the trials. The paper verified isolation by a refractory period over 2 ms; a few units have more than 1% of their ISIs below 2 ms (68 of 566, 36 of the cells studied), and they are kept, as the dandiset holds the units the paper selected.
emg_rows = []
for s in venus:
for ch in stretch.EMG_CHANNELS[:3]:
if ch not in s["windows"]:
state = "absent"
else:
w = s["windows"][ch]
if np.isnan(w).any():
state = "NaN"
elif np.nanstd(w[:, AX < 0]) < 1:
state = "flat"
else:
state = "used"
emg_rows.append({"code": s["session_id"].split("++")[1], "post": s["post"], "channel": ch, "state": state})
emg_avail = pd.DataFrame(emg_rows)
emg_avail["repeat"] = emg_avail.code.isin(repeated_files)
print("Monkey V triceps channels per file (repeated files excluded), by state:")
print(pd.crosstab([emg_avail[~emg_avail.repeat].channel], [emg_avail[~emg_avail.repeat].post, emg_avail[~emg_avail.repeat].state]).to_string())
vals = np.concatenate([s["windows"][ch].ravel() for s in venus for ch in stretch.EMG_CHANNELS[:3] if ch in s["windows"]])
print(f"\nTriceps EMG values (A/D units, as stored in a.u.): {np.nanmin(vals):.0f} to {np.nanmax(vals):.0f}, "
f"median {np.nanmedian(vals):.0f}")
ang = np.concatenate([s["windows"]["angle"].ravel() for s in venus])
print(f"Monkey V angle around perturbations: {np.nanmin(ang):.1f} to {np.nanmax(ang):.1f} deg")
tq = np.concatenate([s["windows"]["torque_cmd"].ravel() for s in venus])
print(f"Monkey V torque command (arbitrary units): {np.nanmin(tq):.0f} to {np.nanmax(tq):.0f}")
iv = np.concatenate([(s["trials"].onset - s["trials"].center_on).dropna().values for s in venus])
# size of the torque command step (sign-corrected, 10 to 40 ms after onset, less the 100 ms before)
command = []
for s in venus:
w = s["windows"]["torque_cmd"]
sign = np.where(s["pert_direction"] == "flexion", 1.0, -1.0)[:, None]
rel = (w - np.nanmedian(w[:, (AX >= -100) & (AX < 0)], axis=1, keepdims=True)) * sign
m = np.nanmean(rel, axis=0)
command.append({"session": s["session_id"], "code": s["session_id"].split("++")[1], "post": s["post"],
"date": s["date"], "step": float(np.nanmedian(rel[:, (AX >= 10) & (AX < 40)])),
"width_ms": int((m > m.max() / 2).sum()), "after_pulse": float(m[(AX >= 70) & (AX < 150)].mean())})
command = pd.DataFrame(command)
command["level"] = np.round(command.step / command.step[command.step > 10].min()).astype(int)
print("\nTorque command step (units above the pre-onset level) by state; level = step in multiples of the smallest:")
print(pd.crosstab(command.level, command.post).to_string())
print(command.groupby(["post", "level"])[["step", "width_ms", "after_pulse"]].median().round(1).to_string())
# per-trial step, against the level of its session
trial_steps = []
for s, lv in zip(venus, command.level):
w = s["windows"]["torque_cmd"]
sign = np.where(s["pert_direction"] == "flexion", 1.0, -1.0)
rel = (w - np.nanmedian(w[:, (AX >= -100) & (AX < 0)], axis=1, keepdims=True)) * sign[:, None]
step = np.nanmedian(rel[:, (AX >= 10) & (AX < 40)], axis=1)
trial_steps.append(pd.DataFrame({"post": s["post"], "session_level": lv,
"trial_level": np.round(step / command.step[command.step > 10].min())}))
trial_steps = pd.concat(trial_steps, ignore_index=True)
print("Per-trial step in multiples of the smallest level, against its session's level (perturbations):")
print(pd.crosstab([trial_steps.post, trial_steps.session_level], trial_steps.trial_level).to_string())
print("Dates of the level-2 sessions by state:", command[command.level == 2].groupby("post").date.agg(["min", "max", "count"]).to_string())
print("Torque onset after the center target appears (s), percentiles 0/5/50/95/100:",
np.percentile(iv, [0, 5, 50, 95, 100]).round(3))
Monkey V triceps channels per file (repeated files excluded), by state: post False True state NaN absent flat used NaN absent flat used channel EMGTricepsLateralis 6 112 0 36 0 97 4 36 EMGTricepsLongus1 2 0 9 143 0 4 3 130 EMGTricepsLongus2 0 0 26 128 13 4 4 116 Triceps EMG values (A/D units, as stored in a.u.): 1705 to 4294, median 2072
Monkey V angle around perturbations: -28.9 to 29.8 deg Monkey V torque command (arbitrary units): 1846 to 2233
Torque command step (units above the pre-onset level) by state; level = step in multiples of the smallest:
post False True
level
0 3 0
1 105 0
2 43 12
4 3 125
step width_ms after_pulse
post level
False 0 0.0 47.0 -0.4
1 45.7 50.0 -0.3
2 91.3 50.0 -0.0
4 182.7 50.0 -0.2
True 2 91.2 50.0 -30.3
4 182.5 50.0 -36.8
Per-trial step in multiples of the smallest level, against its session's level (perturbations):
trial_level -0.0 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0
post session_level
False 0 51 0 0 0 0 0 0 0 0
1 0 1359 0 0 0 0 0 0 0
2 0 13 483 0 13 0 0 0 0
4 0 0 0 0 28 0 0 0 0
True 2 0 16 241 13 0 0 0 0 0
4 0 34 181 446 2002 425 197 45 3
Dates of the level-2 sessions by state: min max count
post
False 1999-04-06 1999-06-30 43
True 2000-01-13 2000-02-16 12
Torque onset after the center target appears (s), percentiles 0/5/50/95/100: [ 1.257 1.696 1.696 1.732 11.841]
EMG units are A/D values stored as arbitrary units, so EMG is compared only within the dandiset. Earlier finding (sampling and EMG description): partly fixed. The EMG description now gives the paper's processing (rectified and low-pass filtered at 100 Hz) where it said sample-hold integrated over 10 ms; the values are the same as in the published version (checked on v4527 while preparing, not run here). The A/D units and the 100 Hz monkey L sampling remain, and the draft's description now states both.
Earlier finding (torque command sizes): partly fixed. The data are unchanged: the command's 50 ms pulse has three session sizes, about 46, 91 and 183 units above rest; 105 of 154 pre-MPTP sessions have the smallest and 125 of 137 post-MPTP sessions the largest, with the middle size in 43 and 12 sessions, and post-MPTP pulses are followed by an opposite swing of about 30 to 37 units. Within post-MPTP sessions the step also varies from trial to trial (1 to 8 times the smallest level in the sessions of the largest size). The torque description now says all this, and that no conversion to newton-metres is available, so whether the torque or only the command's scaling changed still cannot be told. The paper says every impulse was 0.1 Nm for 50 ms. The pre and post perturbations are not shown to be the same; the kinematic and EMG claims add a variant restricted to sessions of the middle size, which both states have.
Earlier finding (torque timing): still present, now described. The torque onset follows the appearance of the center target by 1.696 to 1.732 s in 90% of monkey V perturbations (5th to 95th percentile). The Methods put it 1 to 2 s, uniformly at random, after the monkey captured the center target; the draft's description now states the difference. This does not enter any analysis here.
Responses of the cells studied¶
Every monkey V cell studied is analysed once, below, and the claims draw on these tables. For each cell: the single-trial SDFs (causal gamma kernel, k = 1.5, theta = 20 ms), its baseline rate (mean SDF over the 500 ms before onset, all perturbations), and for each direction the response found by the main detection (threshold on the trial-averaged SDF) and by the variant (t-test across trials), at the true onset and at a pseudo onset 600 ms earlier.
cells = units[units.studied & (units.subject == "Venus")].copy().reset_index(drop=True)
resp_rows, base, mean_sdf = [], [], {}
for c in cells.itertuples():
s, u = S[c.si], S[c.si]["units"][c.ui]
X = stretch.sdfs(u["pert"])
base.append(X[:, stretch._mask(*stretch.BASELINE_MS)].mean())
for d in ("flexion", "extension"):
Xd = X[u["pert_direction"] == d]
mean_sdf[(c.Index, d)] = Xd.mean(0)
for method, shift in (("mean", 0), ("trials", 0), ("mean", stretch.PSEUDO_ONSET_MS),
("trials", stretch.PSEUDO_ONSET_MS)):
r = stretch.detect(Xd, shift_ms=shift, method=method)
r.update(cell=c.Index, direction=d, method=method, pseudo=shift != 0)
resp_rows.append(r)
cells["baseline_rate"] = base
R = pd.DataFrame(resp_rows).merge(cells[["state", "ptype", "si", "date", "duplicate_copy"]],
left_on="cell", right_index=True)
main = R[(R.method == "mean") & ~R.pseudo]
by_cell = main.pivot(index="cell", columns="direction")
cells["torque_related"] = by_cell["responds"].any(axis=1).values
cells["n_directions"] = by_cell["responds"].sum(axis=1).values
cells["latency"] = by_cell["latency"].min(axis=1).values
dsis = [stretch.dsi(main[(main.cell == i) & (main.direction == "flexion")].iloc[0].to_dict(),
main[(main.cell == i) & (main.direction == "extension")].iloc[0].to_dict())
for i in cells.index]
cells["dsi"] = [d for d, _ in dsis]
cells["category"] = [c for _, c in dsis]
resp = main[main.responds].copy()
print(f"{len(cells)} cells studied, {len(main)} cell-directions; responses found: {len(resp)} "
f"in {int(cells.torque_related.sum())} cells")
flat = main[main.baseline_sd == 0]
print("Cell-directions without a usable baseline (SD of the baseline SDF 0; they can never respond):",
len(flat), "| by type and state:", flat.groupby(["ptype", "state"]).size().to_dict())
print("Trials per cell and direction: ", main.n_trials.min(), "to", main.n_trials.max(),
"| baseline rate range", round(cells.baseline_rate.min(), 2), "to", round(cells.baseline_rate.max(), 2), "spikes/s")
print("Response latencies (ms):", resp.latency.min(), "to", resp.latency.max(),
"| offsets:", resp.offset.min(), "to", resp.offset.max())
272 cells studied, 544 cell-directions; responses found: 204 in 138 cells
Cell-directions without a usable baseline (SD of the baseline SDF 0; they can never respond): 14 | by type and state: {('CSN', 'post'): 3, ('CSN', 'pre'): 7, ('NA', 'post'): 3, ('NA', 'pre'): 1}
Trials per cell and direction: 5 to 28 | baseline rate range 0.0 to 48.96 spikes/s
Response latencies (ms): 20.0 to 60.0 | offsets: 24.0 to 600.0
# Per-cell results of the main detection, for recomputing any group statistic: baseline rate, and per
# direction the baseline SD of the mean SDF, the threshold half-width, and for a response its
# latency, sign and largest change.
t_main_ = stats.t.ppf(1 - stretch.P_THRESHOLD / 2, int(np.diff(stretch.BASELINE_MS)[0]) - 1)
per_cell = cells[["unit_name", "state", "ptype", "baseline_rate", "torque_related", "latency", "dsi", "category"]].copy()
for d, k in (("flexion", "F"), ("extension", "E")):
x = main[main.direction == d].set_index("cell")
per_cell[f"{k}_n"] = x.n_trials
per_cell[f"{k}_thr"] = t_main_ * x.baseline_sd
per_cell[f"{k}_lat"] = x.latency
per_cell[f"{k}_sign"] = x.sign
per_cell[f"{k}_max"] = x.max_change
print(per_cell.round(2).to_string())
unit_name state ptype baseline_rate torque_related latency dsi category F_n F_thr F_lat F_sign F_max E_n E_thr E_lat E_sign E_max 0 v0502.1 pre PTN 9.37 True 54.0 0.06 unidirectional 5 11.27 NaN NaN NaN 11 6.40 54.0 1.0 10.65 1 v0601.1 pre CSN 3.72 False NaN 0.41 none 6 13.76 NaN NaN NaN 6 5.10 NaN NaN NaN 2 v1002.1 pre PTN 25.35 True 43.0 0.23 non-directional 6 9.97 45.0 -1.0 16.72 8 10.30 43.0 -1.0 21.75 3 v1003.1 pre CSN 3.76 False NaN 1.27 none 8 7.57 NaN NaN NaN 6 6.39 NaN NaN NaN 4 v1004.1 pre CSN 4.87 False NaN 0.04 none 5 10.55 NaN NaN NaN 8 8.43 NaN NaN NaN 5 v1005.1 pre NA 2.53 True 42.0 0.00 unidirectional 5 5.16 NaN NaN NaN 6 9.32 42.0 1.0 21.69 6 v1006.1 pre PTN 13.60 False NaN 0.55 none 6 16.69 NaN NaN NaN 7 10.52 NaN NaN NaN 7 v1302.1 pre PTN 25.82 True 60.0 0.32 unidirectional 7 14.12 NaN NaN NaN 7 14.38 60.0 -1.0 21.13 8 v1303.1 pre PTN 23.86 True 37.0 0.15 unidirectional 9 10.23 37.0 -1.0 24.96 5 16.17 NaN NaN NaN 9 v1303.2 pre NA 13.93 False NaN 0.30 none 9 8.57 NaN NaN NaN 5 20.34 NaN NaN NaN 10 v1304.1 pre PTN 16.63 False NaN 0.06 none 6 10.26 NaN NaN NaN 8 8.81 NaN NaN NaN 11 v1401.1 pre CSN 0.55 False NaN 1.11 none 5 6.54 NaN NaN NaN 6 0.00 NaN NaN NaN 12 v1402.1 pre CSN 10.88 False NaN 0.04 none 7 18.02 NaN NaN NaN 6 19.25 NaN NaN NaN 13 v1403.1 pre PTN 21.27 True 28.0 0.09 non-directional 5 12.42 39.0 -1.0 17.58 5 9.54 28.0 -1.0 10.50 14 v1501.1+v1501b.1 pre CSN 0.22 False NaN 1.00 none 11 0.00 NaN NaN NaN 16 1.98 NaN NaN NaN 15 v1601.1+v1601b.1 pre PTN 14.80 False NaN 0.01 none 11 6.51 NaN NaN NaN 15 7.51 NaN NaN NaN 16 v1602.1+v1602b.1 pre PTN 8.16 False NaN 0.43 none 9 7.74 NaN NaN NaN 10 9.29 NaN NaN NaN 17 v1702.1+v1702b.1 pre PTN 25.31 True 25.0 0.38 non-directional 12 10.19 38.0 1.0 28.22 17 10.55 25.0 1.0 43.88 18 v1703.1+v1703b.1 pre CSN 4.21 False NaN 1.17 none 12 3.64 NaN NaN NaN 10 6.92 NaN NaN NaN 19 v1802.1 pre NA 23.54 True 36.0 1.43 unidirectional 5 25.15 36.0 1.0 42.50 8 21.95 NaN NaN NaN 20 v1802.2 pre PTN 17.52 True 58.0 0.35 unidirectional 5 12.69 NaN NaN NaN 8 10.65 58.0 1.0 27.79 21 v1901.1 pre PTN 48.96 True 46.0 0.26 unidirectional 6 26.45 NaN NaN NaN 7 28.26 46.0 -1.0 32.54 22 v1902.1+v1902b.1 pre PTN 43.80 True 55.0 0.02 unidirectional 17 18.47 55.0 1.0 56.72 11 37.42 NaN NaN NaN 23 v2301.1 pre NA 10.40 True 20.0 0.37 non-directional 7 14.27 34.0 1.0 75.83 8 20.92 20.0 1.0 120.43 24 v2302.1 pre PTN 13.92 True 36.0 0.75 unidirectional 9 12.77 NaN NaN NaN 5 12.88 36.0 1.0 74.84 25 v2401.1 pre NA 4.98 True 27.0 0.23 non-directional 7 7.47 27.0 1.0 85.22 8 14.71 33.0 1.0 111.00 26 v2502.1 pre NA 6.31 True 38.0 0.34 non-directional 12 12.23 44.0 1.0 61.70 12 7.47 38.0 1.0 93.78 27 v2701.1+v2702.2 pre CSN 3.38 False NaN 1.58 none 16 5.66 NaN NaN NaN 10 8.89 NaN NaN NaN 28 v2702.1 pre NA 0.00 False NaN 1.00 none 9 0.00 NaN NaN NaN 6 0.00 NaN NaN NaN 29 v2703.1 pre CSN 4.70 False NaN 1.26 none 5 10.47 NaN NaN NaN 8 9.81 NaN NaN NaN 30 v2704.1 pre PTN 12.01 True 55.0 1.39 unidirectional 10 8.57 55.0 -1.0 11.68 5 13.20 NaN NaN NaN 31 v2802.1 pre NA 17.43 True 25.0 0.14 non-directional 5 15.15 32.0 1.0 33.35 6 9.65 25.0 1.0 9.99 32 v2804.1 pre NA 0.29 True 22.0 0.95 unidirectional 6 2.73 NaN NaN NaN 8 2.05 22.0 1.0 89.79 33 v2805.1 pre PTN 32.86 True 32.0 0.13 non-directional 9 12.57 34.0 1.0 24.83 5 13.26 32.0 1.0 28.51 34 v2809.1 pre PTN 12.36 False NaN 1.77 none 7 9.39 NaN NaN NaN 6 13.39 NaN NaN NaN 35 v2810.1 pre NA 41.61 True 32.0 0.23 non-directional 5 47.06 32.0 1.0 211.86 6 42.40 40.0 1.0 163.84 36 v2810.2 pre CSN 0.39 True 57.0 0.01 non-directional 5 3.26 59.0 1.0 18.56 6 2.73 57.0 1.0 18.80 37 v2811.1 pre PTN 29.58 True 37.0 1.80 bidirectional 7 13.09 37.0 -1.0 25.05 8 12.49 46.0 -1.0 19.32 38 v2902.1 pre CSN 0.12 False NaN 1.00 none 6 2.68 NaN NaN NaN 6 0.00 NaN NaN NaN 39 v2904.1 pre CSN 3.19 False NaN 0.45 none 7 10.05 NaN NaN NaN 5 7.52 NaN NaN NaN 40 v2906.1 pre NA 4.72 True 21.0 0.28 non-directional 7 14.39 21.0 1.0 162.31 6 9.42 26.0 1.0 117.61 41 v2907.1 pre dual 14.30 True 24.0 0.09 non-directional 6 13.15 24.0 1.0 14.93 7 9.82 50.0 1.0 28.97 42 v3001.1 pre CSN 3.47 False NaN 1.35 none 5 12.44 NaN NaN NaN 6 5.91 NaN NaN NaN 43 v3002.1 pre NA 2.48 True 56.0 0.88 unidirectional 5 10.19 NaN NaN NaN 9 5.33 56.0 1.0 9.53 44 v3002.2 pre NA 6.07 True 36.0 0.04 non-directional 5 11.47 39.0 1.0 68.83 9 11.16 36.0 1.0 71.35 45 v3004.1 pre CSN 1.31 False NaN 0.62 none 7 6.13 NaN NaN NaN 7 3.11 NaN NaN NaN 46 v3006.1 pre PTN 17.25 True 40.0 1.89 unidirectional 7 11.60 40.0 -1.0 11.80 5 14.20 NaN NaN NaN 47 v3007.1 pre PTN 30.29 False NaN 1.74 none 5 24.34 NaN NaN NaN 7 15.76 NaN NaN NaN 48 v3007.2 pre PTN 12.96 True 54.0 1.85 unidirectional 5 10.28 NaN NaN NaN 7 9.54 54.0 1.0 10.87 49 v3103.1 pre PTN 6.44 True 37.0 0.80 bidirectional 6 9.22 51.0 1.0 52.20 8 5.73 37.0 1.0 10.40 50 v3103.2 pre PTN 10.49 False NaN 1.25 none 6 9.70 NaN NaN NaN 8 10.73 NaN NaN NaN 51 v3105.1 pre PTN 18.03 True 46.0 0.68 unidirectional 5 11.40 NaN NaN NaN 7 9.13 46.0 1.0 38.15 52 v3107.1 pre PTN 25.78 True 35.0 0.35 non-directional 6 14.63 35.0 1.0 69.82 5 15.34 48.0 1.0 45.06 53 v3108.1 pre PTN 40.97 True 41.0 1.61 bidirectional 6 12.77 41.0 -1.0 31.45 8 20.48 45.0 -1.0 31.32 54 v3108.2 pre NA 18.52 True 30.0 0.58 bidirectional 6 13.76 30.0 1.0 20.67 8 14.93 41.0 1.0 48.84 55 v3202.1+v3202b.1 pre CSN 7.07 False NaN 2.00 none 10 9.97 NaN NaN NaN 16 6.96 NaN NaN NaN 56 v3203.1 pre NA 10.60 False NaN 1.84 none 5 18.91 NaN NaN NaN 8 14.30 NaN NaN NaN 57 v3203.2 pre NA 4.23 True 54.0 0.17 unidirectional 5 8.15 54.0 1.0 26.32 8 13.19 NaN NaN NaN 58 v3204.1+v3205.1 pre PTN 4.98 True 54.0 1.12 unidirectional 16 6.37 NaN NaN NaN 13 4.90 54.0 1.0 36.98 59 v3205.2 pre CSN 3.84 True 50.0 1.14 unidirectional 8 10.43 50.0 1.0 23.93 5 12.07 NaN NaN NaN 60 v3206.1 pre PTN 15.67 True 53.0 0.50 unidirectional 6 12.94 53.0 1.0 18.29 9 13.87 NaN NaN NaN 61 v3209.1 pre PTN 4.57 True 48.0 0.57 bidirectional 5 5.93 48.0 1.0 17.64 8 10.12 60.0 1.0 40.82 62 v3209.2 pre PTN 20.52 True 53.0 1.37 unidirectional 5 14.46 53.0 1.0 49.37 8 13.21 NaN NaN NaN 63 v3210.1 pre PTN 3.80 False NaN 0.73 none 9 6.78 NaN NaN NaN 8 5.63 NaN NaN NaN 64 v3212.1 pre PTN 9.19 True 37.0 0.74 unidirectional 8 8.28 37.0 1.0 80.03 8 10.84 NaN NaN NaN 65 v3214.1 pre CSN 0.17 False NaN 1.00 none 5 3.28 NaN NaN NaN 7 0.00 NaN NaN NaN 66 v3216.1 pre dual 15.67 False NaN 0.44 none 8 9.47 NaN NaN NaN 5 15.10 NaN NaN NaN 67 v3217.1 pre CSN 0.69 True 23.0 0.20 unidirectional 9 2.89 23.0 1.0 4.18 6 3.55 NaN NaN NaN 68 v3218.1 pre CSN 6.13 False NaN 0.68 none 7 6.21 NaN NaN NaN 6 14.67 NaN NaN NaN 69 v3219.1 pre CSN 20.73 False NaN 1.55 none 7 19.15 NaN NaN NaN 6 27.85 NaN NaN NaN 70 v3220.1 pre PTN 16.94 True 34.0 1.32 unidirectional 6 10.11 34.0 -1.0 15.54 6 11.66 NaN NaN NaN 71 v3220.2 pre NA 19.54 False NaN 0.18 none 6 16.98 NaN NaN NaN 6 10.99 NaN NaN NaN 72 v3221.1 pre CSN 3.02 False NaN 1.90 none 8 6.73 NaN NaN NaN 9 6.51 NaN NaN NaN 73 v3223.1 pre CSN 5.58 False NaN 0.14 none 5 10.15 NaN NaN NaN 5 11.18 NaN NaN NaN 74 v3225.1 pre CSN 2.83 False NaN 0.43 none 6 8.80 NaN NaN NaN 7 7.35 NaN NaN NaN 75 v3227.1 pre NA 0.34 True 32.0 0.69 unidirectional 9 3.56 32.0 1.0 6.11 8 0.00 NaN NaN NaN 76 v3302.1 pre NA 8.46 True 38.0 0.22 non-directional 7 13.49 38.0 1.0 37.14 6 13.61 56.0 1.0 47.88 77 v3302.2 pre NA 0.73 False NaN 0.11 none 7 3.54 NaN NaN NaN 6 4.15 NaN NaN NaN 78 v3303.1 pre NA 21.97 True 58.0 0.14 unidirectional 5 24.91 58.0 1.0 29.05 5 34.98 NaN NaN NaN 79 v3402.1 pre PTN 27.53 True 43.0 0.39 non-directional 10 29.90 59.0 1.0 67.58 6 22.11 43.0 1.0 110.39 80 v3403.1 pre PTN 2.56 True 28.0 0.90 unidirectional 5 9.35 28.0 1.0 24.33 7 8.00 NaN NaN NaN 81 v3403.2 pre NA 24.92 False NaN 0.28 none 5 18.39 NaN NaN NaN 7 10.20 NaN NaN NaN 82 v3404.1 pre NA 16.16 False NaN 0.40 none 7 11.49 NaN NaN NaN 6 11.23 NaN NaN NaN 83 v3405.1 pre PTN 13.80 False NaN 0.13 none 7 9.74 NaN NaN NaN 6 11.05 NaN NaN NaN 84 v3501.1 pre CSN 0.30 False NaN 1.03 none 8 2.93 NaN NaN NaN 9 1.14 NaN NaN NaN 85 v3503.1 pre PTN 21.54 True 43.0 0.68 bidirectional 6 10.50 51.0 1.0 21.94 8 9.96 43.0 1.0 69.03 86 v3505.1 pre CSN 6.18 False NaN 0.32 none 7 8.12 NaN NaN NaN 10 13.20 NaN NaN NaN 87 v3601.1 pre CSN 0.00 False NaN NaN none 8 0.00 NaN NaN NaN 5 0.00 NaN NaN NaN 88 v3601.2 pre PTN 7.28 True 41.0 1.34 unidirectional 8 10.15 41.0 1.0 13.04 5 13.52 NaN NaN NaN 89 v3602.1 pre NA 2.39 True 37.0 0.00 non-directional 8 8.94 39.0 1.0 37.46 7 2.33 37.0 1.0 37.45 90 v3602.2 pre PTN 14.74 False NaN 0.03 none 8 10.64 NaN NaN NaN 7 9.91 NaN NaN NaN 91 v3603.1 pre CSN 3.09 False NaN 0.02 none 7 7.82 NaN NaN NaN 5 10.01 NaN NaN NaN 92 v3603.2 pre NA 28.78 False NaN 1.99 none 7 23.95 NaN NaN NaN 5 14.64 NaN NaN NaN 93 v3607.1+v3607b.1 pre PTN 17.28 True 41.0 0.82 bidirectional 19 9.02 57.0 1.0 11.56 7 10.73 41.0 1.0 62.71 94 v3607b.2 pre CSN 0.98 False NaN 1.86 none 9 1.82 NaN NaN NaN 6 6.46 NaN NaN NaN 95 v3609.1 pre PTN 18.19 False NaN 0.09 none 7 8.44 NaN NaN NaN 6 10.30 NaN NaN NaN 96 v3702.1 pre PTN 10.32 True 46.0 0.20 unidirectional 10 7.78 46.0 1.0 27.58 6 13.59 NaN NaN NaN 97 v3702.2 pre NA 20.94 True 35.0 0.35 unidirectional 10 22.96 35.0 1.0 75.61 6 19.66 NaN NaN NaN 98 v3704.1 pre PTN 25.41 True 38.0 0.54 unidirectional 9 9.31 38.0 1.0 43.91 6 15.35 NaN NaN NaN 99 v3704.2 pre PTN 20.64 True 58.0 1.83 unidirectional 9 9.18 58.0 1.0 12.64 6 11.39 NaN NaN NaN 100 v3802.1 pre CSN 2.75 False NaN 0.38 none 6 14.95 NaN NaN NaN 7 2.36 NaN NaN NaN 101 v3802.2 pre NA 1.01 True 52.0 0.16 unidirectional 6 4.05 52.0 1.0 12.18 7 4.96 NaN NaN NaN 102 v3901.1 pre CSN 0.79 True 33.0 0.92 unidirectional 5 8.15 33.0 1.0 30.13 10 0.00 NaN NaN NaN 103 v3903.1 pre PTN 16.84 False NaN 0.20 none 6 15.24 NaN NaN NaN 5 17.62 NaN NaN NaN 104 v3903.2 pre PTN 26.46 True 59.0 0.42 unidirectional 6 14.81 59.0 1.0 19.23 5 15.06 NaN NaN NaN 105 v4003.1 pre NA 0.47 True 27.0 0.40 non-directional 6 3.52 35.0 1.0 35.53 8 3.13 27.0 1.0 58.88 106 v4003.2 pre CSN 1.41 False NaN 0.42 none 6 4.09 NaN NaN NaN 8 4.81 NaN NaN NaN 107 v4004.1 pre PTN 21.03 True 45.0 1.43 unidirectional 6 14.61 NaN NaN NaN 8 11.35 45.0 1.0 38.13 108 v4004.2 pre PTN 17.83 False NaN 2.00 none 6 11.10 NaN NaN NaN 8 12.12 NaN NaN NaN 109 v4005.1 pre PTN 12.96 True 36.0 0.53 bidirectional 8 11.35 36.0 1.0 63.24 5 15.71 46.0 1.0 29.45 110 v4101.1 pre PTN 15.77 False NaN 1.59 none 7 14.34 NaN NaN NaN 6 10.15 NaN NaN NaN 111 v4203.1 pre CSN 1.18 True 43.0 1.04 unidirectional 8 4.52 NaN NaN NaN 8 4.45 43.0 1.0 22.95 112 v4204.1 pre PTN 19.66 True 52.0 0.27 unidirectional 9 11.58 NaN NaN NaN 6 11.60 52.0 -1.0 18.96 113 v4205.1 pre NA 36.97 False NaN 1.66 none 5 26.90 NaN NaN NaN 6 26.97 NaN NaN NaN 114 v4207.1 pre NA 21.44 False NaN 0.46 none 10 11.87 NaN NaN NaN 5 19.95 NaN NaN NaN 115 v4301.1 post NA 4.18 False NaN 1.36 none 7 9.33 NaN NaN NaN 6 8.19 NaN NaN NaN 116 v4301.2 post NA 8.21 True 30.0 0.17 non-directional 7 22.13 32.0 1.0 94.71 6 14.35 30.0 1.0 114.74 117 v4302.1 post NA 27.89 True 50.0 0.26 unidirectional 8 35.08 50.0 1.0 93.57 7 38.11 NaN NaN NaN 118 v4401.1 post PTN 20.64 True 23.0 0.55 reciprocal 10 6.53 23.0 -1.0 6.82 10 11.74 45.0 1.0 19.69 119 v4401.2 post PTN 2.16 False NaN 0.77 none 10 3.49 NaN NaN NaN 10 4.05 NaN NaN NaN 120 v4402.1 post CSN 6.62 False NaN 0.00 none 14 7.04 NaN NaN NaN 12 6.18 NaN NaN NaN 121 v4403.1 post NA 1.84 False NaN 1.18 none 6 5.02 NaN NaN NaN 6 10.16 NaN NaN NaN 122 v4406.1 post PTN 16.67 True 54.0 0.06 unidirectional 14 9.40 54.0 -1.0 16.95 13 8.21 NaN NaN NaN 123 v4407.1 post PTN 18.61 True 51.0 1.39 reciprocal 13 7.98 51.0 -1.0 9.79 13 5.20 54.0 1.0 25.09 124 v4408.1 post CSN 1.37 True 46.0 0.68 unidirectional 7 4.67 46.0 1.0 8.51 7 4.73 NaN NaN NaN 125 v4410.1 post PTN 4.80 True 32.0 0.47 unidirectional 13 3.27 32.0 1.0 22.10 13 7.24 NaN NaN NaN 126 v4411.1 post PTN 6.42 True 24.0 0.13 non-directional 13 6.21 29.0 1.0 48.03 13 4.85 24.0 1.0 55.19 127 v4412.1 post CSN 0.15 True 55.0 1.09 unidirectional 13 1.26 NaN NaN NaN 13 1.26 55.0 1.0 1.70 128 v4501.1 post CSN 0.98 False NaN 0.68 none 9 3.57 NaN NaN NaN 10 5.13 NaN NaN NaN 129 v4502.1 post CSN 3.10 True 33.0 0.62 unidirectional 11 16.15 NaN NaN NaN 11 3.32 33.0 1.0 3.81 130 v4503.1 post NA 6.97 True 24.0 0.29 non-directional 12 10.40 24.0 1.0 16.00 13 9.11 36.0 1.0 22.46 131 v4504.1 post CSN 3.21 False NaN 0.64 none 13 9.93 NaN NaN NaN 13 8.24 NaN NaN NaN 132 v4505.1 post CSN 2.58 False NaN 1.49 none 6 8.45 NaN NaN NaN 6 6.06 NaN NaN NaN 133 v4506.1 post NA 24.06 True 23.0 0.45 unidirectional 13 13.81 23.0 1.0 30.70 14 17.78 NaN NaN NaN 134 v4507.1 post NA 4.39 True 29.0 0.63 unidirectional 8 17.00 NaN NaN NaN 8 7.86 29.0 1.0 45.04 135 v4508.1 post NA 0.74 False NaN 0.01 none 13 3.04 NaN NaN NaN 14 2.84 NaN NaN NaN 136 v4508.2 post NA 5.09 False NaN 0.43 none 13 10.59 NaN NaN NaN 14 7.64 NaN NaN NaN 137 v4510.1 post NA 12.99 False NaN 1.91 none 12 8.82 NaN NaN NaN 13 8.71 NaN NaN NaN 138 v4511.1 post NA 4.04 True 35.0 0.30 non-directional 10 7.37 38.0 1.0 26.15 10 8.28 35.0 1.0 37.52 139 v4511.2 post NA 13.83 False NaN 1.94 none 10 10.90 NaN NaN NaN 10 12.91 NaN NaN NaN 140 v4512.1 post CSN 1.11 True 54.0 1.03 unidirectional 12 4.11 54.0 1.0 21.11 12 2.83 NaN NaN NaN 141 v4513.1 post CSN 2.50 False NaN 0.39 none 12 5.54 NaN NaN NaN 12 3.23 NaN NaN NaN 142 v4513.2 post CSN 6.19 False NaN 1.88 none 12 8.94 NaN NaN NaN 12 6.68 NaN NaN NaN 143 v4514.1 post NA 6.08 False NaN 0.22 none 13 6.05 NaN NaN NaN 12 6.72 NaN NaN NaN 144 v4515.1 post NA 0.09 False NaN 0.77 none 12 0.00 NaN NaN NaN 12 1.44 NaN NaN NaN 145 v4516.1 post CSN 1.15 False NaN 0.08 none 6 4.09 NaN NaN NaN 8 5.60 NaN NaN NaN 146 v4517.1 post PTN 12.36 False NaN 0.46 none 13 8.16 NaN NaN NaN 12 7.75 NaN NaN NaN 147 v4518.1 post CSN 0.76 False NaN 0.58 none 11 3.17 NaN NaN NaN 10 2.11 NaN NaN NaN 148 v4519.1 post CSN 3.73 False NaN 0.64 none 6 8.03 NaN NaN NaN 6 4.82 NaN NaN NaN 149 v4520.1 post PTN 45.20 True 33.0 0.28 unidirectional 13 14.15 NaN NaN NaN 14 8.52 33.0 -1.0 37.22 150 v4521.1 post NA 6.64 False NaN 0.07 none 13 6.45 NaN NaN NaN 14 7.23 NaN NaN NaN 151 v4522.1 post CSN 13.36 False NaN 0.11 none 8 14.95 NaN NaN NaN 8 10.44 NaN NaN NaN 152 v4523.1 post NA 15.18 True 44.0 0.05 non-directional 16 9.38 44.0 1.0 23.95 17 6.75 49.0 1.0 25.19 153 v4524.1 post CSN 1.37 False NaN 1.35 none 26 1.65 NaN NaN NaN 27 3.25 NaN NaN NaN 154 v4525.1 post NA 14.08 True 23.0 0.80 unidirectional 13 11.79 23.0 1.0 33.45 14 9.59 NaN NaN NaN 155 v4526.1 post CSN 1.17 True 30.0 0.32 non-directional 20 2.73 34.0 1.0 5.87 20 2.81 30.0 1.0 6.03 156 v4527.1 post PTN 13.85 True 23.0 1.65 reciprocal 16 7.22 23.0 1.0 19.67 16 7.94 36.0 -1.0 12.72 157 v4527.2 post PTN 14.30 True 30.0 0.04 unidirectional 16 10.51 30.0 1.0 48.30 16 8.42 NaN NaN NaN 158 v4528.1 post CSN 0.29 True 39.0 0.21 unidirectional 23 1.05 39.0 1.0 1.48 23 1.12 NaN NaN NaN 159 v4529.1 post CSN 1.88 False NaN 0.53 none 18 4.49 NaN NaN NaN 17 4.19 NaN NaN NaN 160 v4530.1 post NA 14.82 True 26.0 0.22 non-directional 9 9.82 28.0 1.0 88.19 10 9.50 26.0 1.0 68.92 161 v4531.1 post PTN 14.38 True 44.0 1.74 unidirectional 13 6.68 NaN NaN NaN 13 7.83 44.0 -1.0 14.96 162 v4532.1 post NA 12.33 False NaN 0.29 none 10 16.01 NaN NaN NaN 10 7.61 NaN NaN NaN 163 v4533.1 post PTN 32.51 True 22.0 0.44 non-directional 11 7.75 22.0 1.0 70.04 12 10.39 46.0 1.0 39.25 164 v4534.1 post NA 5.72 False NaN 0.14 none 14 6.20 NaN NaN NaN 13 5.42 NaN NaN NaN 165 v4535.1 post CSN 0.81 False NaN 0.10 none 14 2.59 NaN NaN NaN 13 3.56 NaN NaN NaN 166 v4536.1 post NA 41.01 True 37.0 1.43 unidirectional 9 9.89 37.0 1.0 40.96 9 10.45 NaN NaN NaN 167 v4537.1 post PTN 30.49 True 27.0 0.43 non-directional 13 6.56 28.0 1.0 24.41 13 6.64 27.0 1.0 42.58 168 v4601.1 post CSN 1.75 False NaN 0.23 none 12 5.82 NaN NaN NaN 11 4.06 NaN NaN NaN 169 v4602.1 post NA 10.53 True 59.0 0.35 unidirectional 18 6.71 59.0 1.0 12.46 18 8.58 NaN NaN NaN 170 v4602.2 post NA 9.22 False NaN 0.50 none 18 9.01 NaN NaN NaN 18 7.03 NaN NaN NaN 171 v4603.1+v4603b.1 post PTN 15.87 True 36.0 0.70 bidirectional 27 4.88 54.0 1.0 12.44 26 5.87 36.0 1.0 41.53 172 v4604.1 post PTN 22.40 True 39.0 0.38 non-directional 14 9.54 54.0 1.0 28.84 13 13.21 39.0 1.0 46.41 173 v4604.2 post PTN 34.66 False NaN 0.16 none 14 18.29 NaN NaN NaN 13 20.40 NaN NaN NaN 174 v4605.1 post PTN 15.23 True 40.0 0.68 unidirectional 13 8.59 40.0 1.0 30.99 13 9.14 NaN NaN NaN 175 v4605.2 post NA 37.16 True 48.0 1.91 unidirectional 13 15.72 NaN NaN NaN 13 15.74 48.0 1.0 29.09 176 v4606.1 post CSN 0.17 False NaN 1.00 none 12 0.00 NaN NaN NaN 12 2.70 NaN NaN NaN 177 v4606.2+v4607.1 post PTN 8.15 True 60.0 0.28 unidirectional 24 4.80 NaN NaN NaN 25 4.80 60.0 -1.0 5.66 178 v4701.1 post NA 1.08 False NaN 1.00 none 7 0.01 NaN NaN NaN 7 5.21 NaN NaN NaN 179 v4702.1 post PTN 10.46 False NaN 0.53 none 9 7.99 NaN NaN NaN 9 8.36 NaN NaN NaN 180 v4702.2 post NA 1.07 False NaN 0.43 none 9 3.37 NaN NaN NaN 9 5.10 NaN NaN NaN 181 v4703.1 post PTN 27.70 False NaN 1.56 none 12 9.16 NaN NaN NaN 12 12.96 NaN NaN NaN 182 v4704.1 post CSN 4.38 False NaN 0.05 none 10 6.86 NaN NaN NaN 10 7.04 NaN NaN NaN 183 v4801.1 post CSN 0.13 False NaN 0.32 none 26 1.31 NaN NaN NaN 27 0.00 NaN NaN NaN 184 v4802.1 post NA 3.70 False NaN 0.14 none 14 6.36 NaN NaN NaN 13 4.80 NaN NaN NaN 185 v4803.1 post CSN 7.59 False NaN 1.46 none 9 8.73 NaN NaN NaN 9 13.59 NaN NaN NaN 186 v4804.1 post CSN 3.05 False NaN 0.05 none 13 5.99 NaN NaN NaN 14 7.10 NaN NaN NaN 187 v4805.1 post CSN 2.48 False NaN 0.30 none 7 5.06 NaN NaN NaN 7 6.95 NaN NaN NaN 188 v4806.1 post CSN 2.07 False NaN 0.44 none 23 3.10 NaN NaN NaN 24 2.70 NaN NaN NaN 189 v4807.1 post NA 2.59 True 22.0 0.26 non-directional 13 6.95 22.0 1.0 33.00 13 3.22 23.0 1.0 24.34 190 v4808.1 post NA 2.40 True 30.0 0.23 unidirectional 13 4.44 30.0 1.0 6.52 13 5.77 NaN NaN NaN 191 v4809.1 post NA 0.76 True 28.0 0.03 unidirectional 13 3.39 28.0 1.0 33.94 13 1.85 NaN NaN NaN 192 v4809.2 post PTN 27.86 True 35.0 0.27 non-directional 13 7.10 35.0 -1.0 9.02 13 7.95 37.0 -1.0 13.19 193 v4810.1 post PTN 7.43 False NaN 0.03 none 27 4.56 NaN NaN NaN 27 3.63 NaN NaN NaN 194 v4810b.2 post NA 1.90 False NaN 0.68 none 13 3.20 NaN NaN NaN 14 4.05 NaN NaN NaN 195 v4812.1 post PTN 9.55 True 55.0 0.10 unidirectional 13 6.43 55.0 -1.0 7.55 13 6.31 NaN NaN NaN 196 v4813.1 post PTN 21.20 True 45.0 0.07 unidirectional 9 11.34 45.0 -1.0 17.31 8 10.56 NaN NaN NaN 197 v4813.2 post PTN 10.84 True 35.0 0.06 unidirectional 9 6.45 35.0 -1.0 10.59 8 8.64 NaN NaN NaN 198 v4814.1 post PTN 12.66 True 51.0 0.07 non-directional 27 5.26 53.0 -1.0 8.79 27 4.82 51.0 -1.0 8.15 199 v4815.1 post PTN 4.43 True 24.0 0.86 bidirectional 28 4.11 41.0 1.0 5.54 26 2.50 24.0 1.0 40.18 200 v4901.1 post PTN 28.13 True 26.0 0.38 unidirectional 14 8.54 26.0 1.0 14.39 13 8.36 NaN NaN NaN 201 v4902.1 post CSN 5.12 False NaN 0.18 none 28 4.34 NaN NaN NaN 26 6.16 NaN NaN NaN 202 v4903.1 post NA 7.87 True 22.0 0.35 non-directional 5 13.57 22.0 1.0 114.25 5 22.26 25.0 1.0 74.11 203 v5001.1 post NA 0.11 True 26.0 0.21 unidirectional 13 0.00 NaN NaN NaN 13 1.56 26.0 1.0 29.54 204 v5002.1 post PTN 5.16 True 35.0 0.04 unidirectional 13 4.12 35.0 1.0 37.55 13 5.90 NaN NaN NaN 205 v5003.1 post PTN 11.76 False NaN 0.24 none 13 6.34 NaN NaN NaN 14 5.13 NaN NaN NaN 206 v5004.1 post PTN 19.00 True 49.0 0.12 non-directional 11 8.99 54.0 -1.0 10.86 10 7.87 49.0 -1.0 9.52 207 v5005.1 post PTN 1.76 True 38.0 0.67 bidirectional 13 4.08 60.0 1.0 5.46 14 3.50 38.0 1.0 16.39 208 v5006.1 post CSN 2.31 False NaN 0.18 none 13 5.79 NaN NaN NaN 14 5.53 NaN NaN NaN 209 v5101.1 post CSN 0.50 False NaN 0.30 none 16 2.58 NaN NaN NaN 16 3.21 NaN NaN NaN 210 v5102.1 post PTN 19.82 True 44.0 0.11 non-directional 13 8.67 51.0 -1.0 13.71 13 8.08 44.0 -1.0 15.35 211 v5103.1 post CSN 0.93 False NaN 0.42 none 13 3.77 NaN NaN NaN 13 3.31 NaN NaN NaN 212 v5104.1 post PTN 18.62 True 27.0 0.20 non-directional 14 8.13 27.0 -1.0 14.16 13 10.86 38.0 -1.0 16.38 213 v5105.1 post PTN 17.01 True 27.0 0.55 bidirectional 13 7.54 27.0 1.0 10.21 13 7.51 35.0 1.0 22.65 214 v5106.1 post PTN 22.15 True 31.0 0.31 non-directional 20 6.72 32.0 -1.0 14.47 20 6.06 31.0 -1.0 20.84 215 v5107.1 post CSN 22.13 True 49.0 0.29 non-directional 24 8.62 50.0 -1.0 15.76 24 7.58 49.0 -1.0 22.24 216 v5108.1 post dual 11.84 True 41.0 0.24 unidirectional 14 7.19 41.0 -1.0 12.49 13 7.80 NaN NaN NaN 217 v5109.1 post NA 7.43 True 25.0 0.52 unidirectional 13 7.44 NaN NaN NaN 13 6.92 25.0 1.0 17.14 218 v5110.1 post PTN 7.92 False NaN 1.78 none 13 9.08 NaN NaN NaN 13 6.47 NaN NaN NaN 219 v5111.1 post PTN 2.50 False NaN 0.29 none 13 5.56 NaN NaN NaN 13 4.40 NaN NaN NaN 220 v5111.2 post PTN 5.45 False NaN 1.61 none 13 5.65 NaN NaN NaN 13 6.97 NaN NaN NaN 221 v5112.1 post PTN 5.25 False NaN 1.28 none 10 6.17 NaN NaN NaN 10 9.31 NaN NaN NaN 222 v5113.1 post PTN 2.80 False NaN 1.73 none 13 5.01 NaN NaN NaN 12 4.42 NaN NaN NaN 223 v5114.1 post PTN 17.75 True 23.0 0.25 non-directional 13 6.35 23.0 1.0 12.62 13 7.78 35.0 1.0 16.84 224 v5115.1 post PTN 0.15 True 29.0 0.77 bidirectional 13 1.30 32.0 1.0 44.39 13 1.26 29.0 1.0 10.38 225 v5115.2 post PTN 6.17 False NaN 0.32 none 13 6.65 NaN NaN NaN 13 5.51 NaN NaN NaN 226 v5116.1 post PTN 3.42 False NaN 0.15 none 14 3.75 NaN NaN NaN 13 5.38 NaN NaN NaN 227 v5201.1 post NA 4.27 True 39.0 0.36 non-directional 6 6.60 39.0 1.0 17.50 7 15.01 55.0 1.0 27.47 228 v5202.1 post NA 7.65 True 32.0 0.28 non-directional 14 14.67 32.0 1.0 67.74 13 10.63 33.0 1.0 94.03 229 v5203.1 post CSN 7.95 False NaN 0.35 none 19 9.36 NaN NaN NaN 18 8.57 NaN NaN NaN 230 v5301.1 post NA 0.05 True 26.0 0.71 bidirectional 13 0.11 30.0 1.0 3.69 13 0.79 26.0 1.0 12.80 231 v5301.2 post CSN 2.46 False NaN 1.69 none 13 6.13 NaN NaN NaN 13 4.60 NaN NaN NaN 232 v5302.1 post PTN 6.86 False NaN 0.43 none 8 8.71 NaN NaN NaN 8 9.53 NaN NaN NaN 233 v5303.1 post CSN 0.73 False NaN 0.39 none 13 2.25 NaN NaN NaN 12 2.69 NaN NaN NaN 234 v5304.1 post CSN 0.43 False NaN 0.17 none 13 3.37 NaN NaN NaN 13 1.65 NaN NaN NaN 235 v5305.1 post NA 8.35 True 48.0 0.06 non-directional 12 7.08 48.0 1.0 15.81 12 8.25 51.0 1.0 14.84 236 v5401.1 post NA 13.29 True 20.0 0.00 non-directional 9 18.10 23.0 1.0 93.87 8 16.22 20.0 1.0 94.14 237 v5402.1 post PTN 39.55 True 28.0 0.11 non-directional 13 16.52 29.0 1.0 28.31 13 10.05 28.0 1.0 31.93 238 v5403.1 post CSN 10.96 False NaN 1.44 none 14 8.66 NaN NaN NaN 13 5.66 NaN NaN NaN 239 v5404.1 post CSN 6.32 False NaN 1.96 none 14 5.80 NaN NaN NaN 13 5.36 NaN NaN NaN 240 v5501.1 post NA 0.00 False NaN 0.15 none 13 0.00 NaN NaN NaN 14 0.00 NaN NaN NaN 241 v5502.1 post NA 4.52 True 25.0 0.50 bidirectional 13 7.92 25.0 1.0 41.59 14 7.95 27.0 1.0 20.64 242 v5503.1 post CSN 1.64 True 32.0 1.32 unidirectional 13 3.41 32.0 1.0 5.48 13 4.02 NaN NaN NaN 243 v5601.1 post NA 2.47 True 49.0 0.45 unidirectional 8 8.96 NaN NaN NaN 7 5.71 49.0 1.0 6.70 244 v5602.1 post CSN 0.59 False NaN 0.08 none 13 1.65 NaN NaN NaN 14 2.89 NaN NaN NaN 245 v5603.1 post PTN 12.46 True 30.0 0.35 non-directional 13 8.45 39.0 1.0 19.70 14 5.19 30.0 1.0 30.33 246 v5604.1 post PTN 4.31 False NaN 0.15 none 26 3.90 NaN NaN NaN 28 4.19 NaN NaN NaN 247 v5604b.2 post CSN 0.00 False NaN NaN none 13 0.00 NaN NaN NaN 14 0.00 NaN NaN NaN 248 v5605.1 post CSN 0.62 False NaN 1.12 none 8 2.04 NaN NaN NaN 8 4.69 NaN NaN NaN 249 v5606.1 post PTN 3.29 True 32.0 0.32 unidirectional 20 4.21 32.0 1.0 15.59 20 4.49 NaN NaN NaN 250 v5607.1 post CSN 0.21 False NaN 0.63 none 26 0.92 NaN NaN NaN 27 1.04 NaN NaN NaN 251 v5608.1 post PTN 8.56 True 42.0 1.67 reciprocal 13 6.34 42.0 1.0 9.44 13 6.07 50.0 -1.0 6.34 252 v5801.1 post NA 0.00 True 30.0 0.35 unidirectional 7 0.00 NaN NaN NaN 7 0.00 30.0 1.0 14.44 253 v5801.2 post NA 6.91 True 25.0 0.24 non-directional 7 17.23 25.0 1.0 104.25 7 16.36 27.0 1.0 79.37 254 v5802.1 post NA 0.56 False NaN 1.27 none 7 3.50 NaN NaN NaN 7 2.34 NaN NaN NaN 255 v5802.2 post NA 8.35 False NaN 0.44 none 7 19.42 NaN NaN NaN 7 6.85 NaN NaN NaN 256 v5803.1 post NA 1.86 False NaN 0.22 none 13 9.03 NaN NaN NaN 13 3.22 NaN NaN NaN 257 v5804.1 post CSN 8.84 False NaN 0.46 none 13 8.27 NaN NaN NaN 14 7.09 NaN NaN NaN 258 v5805.1 post NA 1.55 False NaN 1.84 none 13 3.73 NaN NaN NaN 14 4.51 NaN NaN NaN 259 v5805.2 post NA 8.18 False NaN 0.28 none 13 10.73 NaN NaN NaN 14 6.21 NaN NaN NaN 260 v5806.1 post CSN 8.05 False NaN 0.02 none 18 4.50 NaN NaN NaN 17 6.21 NaN NaN NaN 261 v5806.2 post NA 0.64 False NaN 1.40 none 18 3.07 NaN NaN NaN 17 0.96 NaN NaN NaN 262 v5807.1 post CSN 2.09 False NaN 1.47 none 9 6.18 NaN NaN NaN 9 3.53 NaN NaN NaN 263 v5807.2 post PTN 18.19 True 38.0 0.49 unidirectional 9 10.54 NaN NaN NaN 9 7.94 38.0 1.0 19.72 264 v5808.1 post NA 7.82 True 20.0 0.20 non-directional 14 10.47 21.0 1.0 95.79 13 8.78 20.0 1.0 119.50 265 v5809.1 post CSN 5.22 False NaN 1.75 none 13 7.40 NaN NaN NaN 13 4.73 NaN NaN NaN 266 v5810.1 post NA 33.46 True 22.0 0.32 non-directional 14 18.16 22.0 1.0 150.03 13 22.79 25.0 1.0 102.12 267 v5811.1 post PTN 13.51 True 47.0 0.04 unidirectional 26 5.32 47.0 -1.0 10.23 26 5.83 NaN NaN NaN 268 v5811.2 post PTN 13.50 False NaN 1.74 none 26 6.19 NaN NaN NaN 26 7.11 NaN NaN NaN 269 v5813.1 post NA 11.76 False NaN 0.38 none 7 10.21 NaN NaN NaN 6 11.41 NaN NaN NaN 270 v5814.1 post CSN 0.45 False NaN 0.27 none 13 1.90 NaN NaN NaN 14 2.38 NaN NaN NaN 271 v5815.1 post NA 1.94 False NaN 0.06 none 7 3.60 NaN NaN NaN 7 5.08 NaN NaN NaN
cell-counts¶
In monkey V, 112 neurons were studied before MPTP and 159 after, of which 49 and 54 were PTNs and 31 and 48 CSNs; three neurons in all were activated from both the peduncle and the putamen.
counts = pd.crosstab(cells.ptype, cells.state).reindex(columns=["pre", "post"])
counts.loc["all"] = counts.sum()
print("Monkey V cells studied by type and state:")
print(counts.to_string())
dual_all = units[units.ptype == "dual"]
print(f"\nUnits activated from both the peduncle and the putamen, all files: {len(dual_all)} "
f"({dual_all.subject.value_counts().to_dict()}); among the cells studied: {int((cells.ptype == 'dual').sum())}")
print("Other dual projections (with thalamus or STN):", int(units.projection.str.contains('_and_').sum() - len(dual_all)))
print("\nVariant, the two unit copies in repeated files dropped:")
print(pd.crosstab(cells[~cells.duplicate_copy].ptype, cells[~cells.duplicate_copy].state)
.reindex(columns=["pre", "post"]).to_string())
# Units recorded across several source files are one cell in this version. Variant: each unit
# counted from the perturbations of its first source file only (the source_file of its first observed
# trial; a unit of a repeated "b" file has its trials stored under the earlier file's name).
first_only = []
for r in units.itertuples():
u = S[r.si]["units"][r.ui]
src = S[r.si]["trials"].source_file.values
f0 = src[np.flatnonzero(u["observed_trials"])[0]]
own = src[u["pert_trial"]] == f0
d = u["pert_direction"][own]
first_only.append(((d == "flexion").sum() >= stretch.MIN_REPS) & ((d == "extension").sum() >= stretch.MIN_REPS))
fo = units.assign(studied_first=first_only)
fo = fo[fo.studied_first & (fo.subject == "Venus")]
print("\nVariant, units across several files counted from their first file only:")
print(pd.crosstab(fo.ptype, fo.state).reindex(columns=["pre", "post"]).to_string())
Monkey V cells studied by type and state:
state pre post
ptype
CSN 32 48
NA 30 54
PTN 51 54
dual 2 1
all 115 157
Units activated from both the peduncle and the putamen, all files: 3 ({'Venus': 3}); among the cells studied: 3
Other dual projections (with thalamus or STN): 4
Variant, the two unit copies in repeated files dropped:
state pre post
ptype
CSN 32 48
NA 30 54
PTN 51 54
dual 2 1
Variant, units across several files counted from their first file only: state pre post ptype CSN 31 48 NA 30 54 PTN 49 54 dual 2 1
Result: reproduced. The monkey V cells studied number 115 before and 157 after MPTP, with 51 and 54 PTNs and 32 and 48 CSNs, all within tolerance of the supplementary table's 112 and 159, 49 and 54, 31 and 48; the files hold 3 units activated from both the peduncle and the putamen, all in monkey V and all among the cells studied. Counting each unit from its first source file only (variant) gives 112 and 157, with 49 and 54 PTNs and 31 and 48 CSNs. The published version gave the table's counts exactly (112 and 159), counting the two unit copies the draft no longer holds and without the trials of the second files.
monkey-l-perturbations¶
In monkey L, 57 of 115 neurons before MPTP and 35 of 73 after responded to torque perturbations of the wrist.
Blocked: the monkey L files mark no torque perturbation (data check above; the type is
not_available in every trial), so neither the trials nor their onsets can be found. The files
hold 138 and 72 monkey L units, against the paper's 115 and 73 cells studied.
baseline-rate¶
Before torque onset, the firing rate of M1 neurons was lower after MPTP, because PTN rates fell (Mann-Whitney p < 0.001) while CSN rates did not change (p = 0.5). In monkey V, 12.6 and 8.9 spikes/s for all cells, 20.3 and 14.8 for PTNs, 3.5 and 3.4 for CSNs.
def mwu(a, b):
return stats.mannwhitneyu(a, b, alternative="two-sided").pvalue
expected_rate = {"M1": (12.6, 8.9, 1.98, 1.33), "PTN": (20.3, 14.8, 2.46, 2.72), "CSN": (3.5, 3.4, 1.47, 1.01)}
rows = []
for g, (e_pre, e_post, t_pre, t_post) in expected_rate.items():
sub = cells if g == "M1" else cells[cells.ptype == g]
pre, post = sub[sub.state == "pre"].baseline_rate, sub[sub.state == "post"].baseline_rate
rows.append({"group": g, "n_pre": len(pre), "pre": pre.mean(), "sd_pre": pre.std(), "expected_pre": e_pre,
"in_tol_pre": abs(pre.mean() - e_pre) <= t_pre, "n_post": len(post), "post": post.mean(),
"sd_post": post.std(), "expected_post": e_post, "in_tol_post": abs(post.mean() - e_post) <= t_post,
"change_%": 100 * (post.mean() / pre.mean() - 1), "mwu_p": mwu(pre, post)})
rate_table = pd.DataFrame(rows)
print(rate_table.round(4).to_string(index=False))
nodup = cells[~cells.duplicate_copy]
p = nodup[nodup.ptype == "PTN"]
print(f"\nVariant without the unit copies, PTN: {p[p.state == 'pre'].baseline_rate.mean():.2f} and "
f"{p[p.state == 'post'].baseline_rate.mean():.2f} spikes/s, p = {mwu(p[p.state == 'pre'].baseline_rate, p[p.state == 'post'].baseline_rate):.4f}")
group n_pre pre sd_pre expected_pre in_tol_pre n_post post sd_post expected_post in_tol_post change_% mwu_p M1 115 12.6618 10.8980 12.6 True 157 8.8463 9.4447 8.9 True -30.1340 0.0022 PTN 51 18.4881 9.7849 20.3 True 54 14.3682 10.1503 14.8 True -22.2842 0.0201 CSN 32 3.4845 4.0101 3.5 True 48 3.3647 4.1454 3.4 True -3.4374 0.6514 Variant without the unit copies, PTN: 18.49 and 14.37 spikes/s, p = 0.0201
Result: reproduced. All six means are within tolerance (M1 12.66 and 8.85, PTN 18.49 and 14.37, CSN 3.48 and 3.36 spikes/s). The PTN rate is lower after MPTP (Mann-Whitney p = 0.020) and the CSN rate does not differ (p = 0.65); the M1 rate is lower (p = 0.0022). The statement's p < 0.001 for PTNs, from both monkeys, is not reached in monkey V alone.
ptn-faster-than-csn¶
Before MPTP, PTNs fired much faster than CSNs at rest (Mann-Whitney p < 0.001).
pre = cells[cells.state == "pre"]
ptn_pre, csn_pre = pre[pre.ptype == "PTN"].baseline_rate, pre[pre.ptype == "CSN"].baseline_rate
print(f"PTN {ptn_pre.mean():.2f} spikes/s (range {ptn_pre.min():.2f} to {ptn_pre.max():.2f}, n = {len(ptn_pre)}); "
f"CSN {csn_pre.mean():.2f} (range {csn_pre.min():.2f} to {csn_pre.max():.2f}, n = {len(csn_pre)}); "
f"Mann-Whitney p = {mwu(ptn_pre, csn_pre):.2e}")
print("Paper (both monkeys): text 16.3 and 3.6, Table 1 18.3 and 3.2 spikes/s")
PTN 18.49 spikes/s (range 2.56 to 48.96, n = 51); CSN 3.48 (range 0.00 to 20.73, n = 32); Mann-Whitney p = 6.33e-12 Paper (both monkeys): text 16.3 and 3.6, Table 1 18.3 and 3.2 spikes/s
Result: reproduced. Before MPTP, monkey V PTNs fire at 18.49 spikes/s and CSNs at 3.48 (Mann-Whitney p = 6.3e-12).
antidromic-latency¶
PTN antidromic latencies did not change after MPTP (1.9 ms, p = 0.9), while CSN latencies were longer (5.1 against 6.1 ms, p < 0.05/3).
def antidromic(df, label):
rows = []
for g, e in (("PTN", (1.9, 1.9, 0.22, 0.2)), ("CSN", (5.1, 6.1, 0.61, 0.89))):
sub = df[df.ptype == g]
a, b = sub[sub.state == "pre"].antidromic_latency_ms.dropna(), sub[sub.state == "post"].antidromic_latency_ms.dropna()
rows.append({"cells": label, "group": g, "n_pre": len(a), "pre": a.mean(), "sd_pre": a.std(),
"in_tol_pre": abs(a.mean() - e[0]) <= e[2], "n_post": len(b), "post": b.mean(),
"sd_post": b.std(), "in_tol_post": abs(b.mean() - e[1]) <= e[3], "mwu_p": mwu(a, b)})
return rows
leu_units = units[(units.subject == "Leu") & units.ptype.isin(["PTN", "CSN"])]
ad = pd.DataFrame(antidromic(cells, "monkey V studied")
+ antidromic(pd.concat([cells, leu_units]), "variant: plus all monkey L units"))
print(ad.round(4).to_string(index=False))
print("Antidromic latency range (ms), cells studied:", cells.antidromic_latency_ms.min(), "to",
cells.antidromic_latency_ms.max(), "| NaN for", int(cells[cells.ptype.isin(['PTN', 'CSN'])].antidromic_latency_ms.isna().sum()), "PTNs and CSNs")
alpha = 0.05 / 3
for r in ad.itertuples():
if r.group == "CSN":
verdict = f"longer after MPTP and below 0.05/3 = {alpha:.4f}: {bool(r.post > r.pre and r.mwu_p < alpha)}"
else:
verdict = f"p >= 0.05 (no difference): {bool(r.mwu_p >= 0.05)}"
print(f"{r.cells}, {r.group}: Mann-Whitney p = {r.mwu_p:.4f}; {verdict}")
cells group n_pre pre sd_pre in_tol_pre n_post post sd_post in_tol_post mwu_p
monkey V studied PTN 49 1.9214 0.8305 True 48 1.8594 0.8072 True 0.6808
monkey V studied CSN 32 4.9109 1.7512 True 38 6.8474 3.7243 True 0.0190
variant: plus all monkey L units PTN 65 1.8038 0.8009 True 55 1.8318 0.7730 True 0.7841
variant: plus all monkey L units CSN 56 5.3848 2.3268 True 47 7.0149 3.4964 False 0.0072
Antidromic latency range (ms), cells studied: 0.65 to 20.5 | NaN for 18 PTNs and CSNs
monkey V studied, PTN: Mann-Whitney p = 0.6808; p >= 0.05 (no difference): True
monkey V studied, CSN: Mann-Whitney p = 0.0190; longer after MPTP and below 0.05/3 = 0.0167: False
variant: plus all monkey L units, PTN: Mann-Whitney p = 0.7841; p >= 0.05 (no difference): True
variant: plus all monkey L units, CSN: Mann-Whitney p = 0.0072; longer after MPTP and below 0.05/3 = 0.0167: True
Result: partially reproduced. The four means are within tolerance (PTN 1.92 and 1.86 ms, CSN 4.91 and 6.85 ms) and the PTN latency does not differ (p = 0.68). The CSN latency is longer after MPTP with p = 0.019, just above the paper's threshold of 0.05/3 = 0.0167, so the CSN test fails (printed above). With all monkey L PTNs and CSNs added (variant), p = 0.0072 for CSNs and 0.78 for PTNs, and the post-MPTP CSN mean (7.01 ms) is outside tolerance.
false-positive-rate¶
The detection threshold, p < 0.01 divided by 5 independent comparisons, holds the chance of a false response over the 200 ms epoch near 0.01.
fp = R[R.pseudo].groupby("method").agg(cell_directions=("responds", "size"),
any_crossing_0_200=("any_crossing_0_200", "mean"),
response_20_60=("responds", "mean"))
real = R[~R.pseudo].groupby("method").agg(any_crossing_0_200=("any_crossing_0_200", "mean"),
response_20_60=("responds", "mean"))
print("At the pseudo onset (600 ms before the torque), fraction of cell-directions:")
print(fp.round(4).to_string())
print("\nAt the true onset, for reference:")
print(real.round(4).to_string())
print("\nPseudo onset, main method, by state:")
print(R[R.pseudo & (R.method == "mean")].groupby("state")[["any_crossing_0_200", "responds"]].mean().round(4).to_string())
print(f"\nPseudo-onset over true-onset rate of responses (main): {fp.loc['mean', 'response_20_60'] / real.loc['mean', 'response_20_60']:.3f}")
At the pseudo onset (600 ms before the torque), fraction of cell-directions:
cell_directions any_crossing_0_200 response_20_60
method
mean 544 0.2279 0.0478
trials 544 0.4632 0.0901
At the true onset, for reference:
any_crossing_0_200 response_20_60
method
mean 0.6636 0.3750
trials 0.8640 0.3566
Pseudo onset, main method, by state:
any_crossing_0_200 responds
state
post 0.2102 0.0478
pre 0.2522 0.0478
Pseudo-onset over true-onset rate of responses (main): 0.127
# Variant null: unperturbed trials, at the time a torque usually comes (the session's median delay
# from center target to torque onset), so the epoch matches the true one but has no torque.
null_rows = []
for c in cells.itertuples():
s, u = S[c.si], S[c.si]["units"][c.ui]
tr = s["trials"]
delay = np.nanmedian(tr.onset - tr.center_on)
t0 = (tr.center_on + delay)[(tr.torque == "none") & u["observed_trials"]].dropna().values
if len(t0) < stretch.MIN_REPS:
continue
rel = [np.round((u["spikes"][(u["spikes"] >= t - 1.8) & (u["spikes"] < t + 0.7)] - t) * 1000) for t in t0]
for method in ("mean", "trials"):
r = stretch.detect(stretch.sdfs(rel), method=method)
null_rows.append({"cell": c.Index, "method": method, "n_trials": len(t0), "state": c.state,
"any_crossing_0_200": r["any_crossing_0_200"], "responds": r["responds"]})
null = pd.DataFrame(null_rows)
print(f"Unperturbed trials at the usual torque time: {null.cell.nunique()} cells with at least 5 such trials "
f"({null.n_trials.min()} to {null.n_trials.max()} trials)")
print(null.groupby("method")[["any_crossing_0_200", "responds"]].mean().round(4).to_string())
Unperturbed trials at the usual torque time: 255 cells with at least 5 such trials (5 to 27 trials)
any_crossing_0_200 responds
method
mean 0.1922 0.0471
trials 0.4431 0.0706
# Where the false crossings come from. (1) The pseudo-onset rate by baseline rate and by number of
# trials. (2) Stationary Poisson spike trains at each cell-direction's baseline rate, with its number
# of trials, through the same detection: no response exists, so every crossing is false.
pm = R[R.pseudo & (R.method == "mean")].copy()
pm["rate_group"] = pd.qcut(pm.cell.map(cells.baseline_rate), 4)
pm["trials_group"] = pd.qcut(pm.n_trials, 3, duplicates="drop")
print("Pseudo onset, main detection, by baseline rate quartile (spikes/s):")
print(pm.groupby("rate_group", observed=True)[["any_crossing_0_200", "responds"]].agg(["mean", "size"]).round(3).to_string())
print("\nBy number of trials per direction:")
print(pm.groupby("trials_group", observed=True)[["any_crossing_0_200", "responds"]].agg(["mean", "size"]).round(3).to_string())
rng_fp = np.random.default_rng(1636)
span = int((stretch.SPIKE_WINDOW_S[1] - stretch.SPIKE_WINDOW_S[0]) * 1000)
sim = []
for r in main.itertuples():
rate = cells.baseline_rate[r.cell]
trains = [np.flatnonzero(rng_fp.random(span) < rate / 1000) + int(stretch.SPIKE_WINDOW_S[0] * 1000)
for _ in range(r.n_trials)]
X = stretch.sdfs(trains)
for method in ("mean", "trials"):
q = stretch.detect(X, method=method)
sim.append({"method": method, "any_crossing_0_200": q["any_crossing_0_200"], "responds": q["responds"]})
sim = pd.DataFrame(sim)
print("\nStationary Poisson trains at each cell-direction's baseline rate (seeded), same detection:")
print(sim.groupby("method")[["any_crossing_0_200", "responds"]].mean().round(4).to_string())
# The same simulation with other seeds, main detection only, to show its spread.
spread = []
for seed in range(1637, 1645):
g = np.random.default_rng(seed)
hits = []
for r in main.itertuples():
rate = cells.baseline_rate[r.cell]
trains = [np.flatnonzero(g.random(span) < rate / 1000) + int(stretch.SPIKE_WINDOW_S[0] * 1000)
for _ in range(r.n_trials)]
hits.append(stretch.detect(stretch.sdfs(trains))["any_crossing_0_200"])
spread.append(np.mean(hits))
print(f"Main detection, seeds 1637 to 1644: any crossing {np.min(spread):.4f} to {np.max(spread):.4f} "
f"(mean {np.mean(spread):.4f}) over {len(main)} cell-directions")
Pseudo onset, main detection, by baseline rate quartile (spikes/s):
any_crossing_0_200 responds
mean size mean size
rate_group
(-0.001, 2.466] 0.176 136 0.044 136
(2.466, 7.02] 0.287 136 0.059 136
(7.02, 15.94] 0.206 136 0.059 136
(15.94, 48.965] 0.243 136 0.029 136
By number of trials per direction:
any_crossing_0_200 responds
mean size mean size
trials_group
(4.999, 8.0] 0.227 233 0.047 233
(8.0, 13.0] 0.213 211 0.038 211
(13.0, 28.0] 0.260 100 0.070 100
Stationary Poisson trains at each cell-direction's baseline rate (seeded), same detection:
any_crossing_0_200 responds
method
mean 0.2096 0.0515
trials 0.4835 0.0993
Main detection, seeds 1637 to 1644: any crossing 0.1581 to 0.1930 (mean 0.1735) over 544 cell-directions
# How the result changes with the threshold. The main detection crosses when the largest change over
# 0 to 200 ms exceeds t times the baseline SD of the mean SDF; the t that gives a pseudo-onset rate
# of 0.01 is read off the pseudo-onset distribution of that ratio and applied at the true onset.
z_pseudo = pm.max_abs_z.fillna(-np.inf).values
t_main = stats.t.ppf(1 - stretch.P_THRESHOLD / 2, int(np.diff(stretch.BASELINE_MS)[0]) - 1)
t_cal = float(np.quantile(z_pseudo, 0.99))
print(f"t of the main detection {t_main:.2f}; t giving a pseudo-onset rate of 0.01: {t_cal:.2f}")
cal_rows = []
for (cell, d), m in mean_sdf.items():
for shift in (0, stretch.PSEUDO_ONSET_MS):
q = stretch.detect(m[None, :], shift_ms=shift, tcrit=t_cal)
cal_rows.append({"cell": cell, "direction": d, "pseudo": shift != 0, "responds": q["responds"],
"any_crossing_0_200": q["any_crossing_0_200"]})
cal = pd.DataFrame(cal_rows)
print(cal.groupby("pseudo")[["any_crossing_0_200", "responds"]].mean().round(4).to_string())
cal_cell = cal[~cal.pseudo].groupby("cell").responds.any()
first_overall = main.assign(early=main.first_crossing < stretch.ONSET_RANGE_MS[0])
print(f"\nDetected responses (main) with an earlier crossing at 0 to {stretch.ONSET_RANGE_MS[0] - 1} ms: "
f"{int((first_overall.responds & first_overall.early).sum())} of {int(main.responds.sum())}")
fo_cell = first_overall.assign(r=first_overall.responds & ~first_overall.early).groupby("cell").r.any()
rows = []
for label, rel in (("main", cells.torque_related), ("calibrated threshold", cal_cell.reindex(cells.index)),
("first crossing overall in 20 to 60 ms", fo_cell.reindex(cells.index))):
for g in ("M1", "PTN", "CSN"):
sub = cells.index if g == "M1" else cells.index[cells.ptype == g]
for st in ("pre", "post"):
idx = [k for k in sub if cells.state[k] == st]
rows.append({"detection": label, "group": g, "state": st, "percent": 100 * rel[idx].mean(),
"expected": {('M1', 'pre'): 75, ('M1', 'post'): 68, ('PTN', 'pre'): 55, ('PTN', 'post'): 68, ('CSN', 'pre'): 45, ('CSN', 'post'): 21}[(g, st)]})
print("\nTorque-related percent by detection (paper's monkey V values for reference):")
print(pd.DataFrame(rows).pivot(index=["group", "state"], columns="detection", values="percent").round(1)
.join(pd.DataFrame(rows).groupby(["group", "state"]).expected.first()).to_string())
# Example mean SDFs with their thresholds, as Figure 3 draws them: the cell-directions at the
# 10th, 50th and 90th percentile of baseline rate among detected responses and among non-responses.
pick = []
for flag in (True, False):
g = main[main.responds == flag].sort_values("baseline")
pick += [g.iloc[int(q * (len(g) - 1))] for q in (0.1, 0.5, 0.9)]
fig, axes = plt.subplots(2, 3, figsize=(10, 5), sharex=True)
for ax, r in zip(axes.ravel(), pick):
m = mean_sdf[(r.cell, r.direction)]
thr = t_main * r.baseline_sd
ax.plot(stretch.GRID, m, color="k", lw=0.8)
for y in (r.baseline - thr, r.baseline + thr):
ax.axhline(y, color="0.5", lw=0.8, ls="--")
ax.axvspan(stretch.PSEUDO_ONSET_MS, stretch.PSEUDO_ONSET_MS + 200, color="0.92")
ax.axvspan(0, 200, color="0.85")
ax.set_xlim(stretch.GRID[0], stretch.GRID[-1])
ax.set_title(f"{cells.ptype[r.cell]} {r.state} {r.direction}, {r.baseline:.1f} sp/s, "
f"{'response' if r.responds else 'none'}", fontsize=8)
axes[1, 0].set_xlabel("ms from torque onset")
axes[0, 0].set_ylabel("mean SDF (spikes/s)")
plt.tight_layout()
plt.show()
for r in pick:
m = mean_sdf[(r.cell, r.direction)]
print(f"cell {r.cell} {r.direction}: SDF {m.min():.1f} to {m.max():.1f} spikes/s, threshold "
f"{r.baseline - t_main * r.baseline_sd:.1f} to {r.baseline + t_main * r.baseline_sd:.1f}; "
f"pseudo-onset crossing {bool(pm[(pm.cell == r.cell) & (pm.direction == r.direction)].any_crossing_0_200.iloc[0])}")
print("Shaded: 200 ms after the pseudo onset (light) and the torque (dark). Every SDF is within its axes.")
t of the main detection 3.11; t giving a pseudo-onset rate of 0.01: 6.67
any_crossing_0_200 responds
pseudo
False 0.3401 0.2132
True 0.0110 0.0018
Detected responses (main) with an earlier crossing at 0 to 19 ms: 5 of 204
Torque-related percent by detection (paper's monkey V values for reference):
calibrated threshold first crossing overall in 20 to 60 ms main expected
group state
CSN post 0.0 16.7 16.7 21
pre 3.1 15.6 15.6 45
M1 post 32.5 47.1 47.8 68
pre 30.4 53.0 54.8 75
PTN post 42.6 66.7 68.5 68
pre 33.3 72.5 72.5 55
cell 191 flexion: SDF 0.0 to 35.1 spikes/s, threshold -2.2 to 4.5; pseudo-onset crossing False cell 267 flexion: SDF 2.7 to 19.0 spikes/s, threshold 7.6 to 18.3; pseudo-onset crossing True cell 200 flexion: SDF 12.9 to 43.6 spikes/s, threshold 20.6 to 37.7; pseudo-onset crossing True cell 234 extension: SDF 0.0 to 3.4 spikes/s, threshold -1.3 to 2.0; pseudo-onset crossing False cell 239 extension: SDF 0.8 to 11.7 spikes/s, threshold -0.3 to 10.4; pseudo-onset crossing False cell 112 flexion: SDF 4.4 to 28.9 spikes/s, threshold 7.7 to 30.8; pseudo-onset crossing False Shaded: 200 ms after the pseudo onset (light) and the torque (dark). Every SDF is within its axes.
# Do the spike windows fall inside the trials? Spikes are stored only inside trials, so a window that
# starts before its trial's start is partly empty by construction. For every perturbation of the cells
# studied: time from the trial start to the torque onset, and the share of perturbations whose
# spike window (-1800 ms), pseudo baseline with the kernel's reach (-1500 ms), true baseline with the
# kernel's reach (-900 ms), and SDF end (+600 ms, against the trial stop) fall outside the trial.
lead, tail = [], []
for c in cells.itertuples():
s, u = S[c.si], S[c.si]["units"][c.ui]
tr = s["trials"].loc[u["pert_trial"]]
lead.append((tr.onset - tr.start).values * 1000)
tail.append((tr.stop - tr.onset).values * 1000)
lead, tail = np.concatenate(lead), np.concatenate(tail)
print(f"{len(lead)} cell-perturbations; onset minus trial start (ms), percentiles 0/1/5/50/100:",
np.percentile(lead, [0, 1, 5, 50, 100]).round(0))
print("Trial stop minus onset (ms), percentiles 0/1/50:", np.percentile(tail, [0, 1, 50]).round(0))
for label, ms in (("spike window, -1800 ms", 1800), ("pseudo baseline with kernel, -1500 ms", 1500),
("true baseline with kernel, -900 ms", 900)):
print(f" {label}: starts before the trial start in {np.mean(lead < ms):.2%}")
print(f" SDF end, +600 ms: after the trial stop in {np.mean(tail < 600):.2%}")
5857 cell-perturbations; onset minus trial start (ms), percentiles 0/1/5/50/100: [ 2968. 3275. 3404. 4049. 14653.]
Trial stop minus onset (ms), percentiles 0/1/50: [1516. 3891. 4777.] spike window, -1800 ms: starts before the trial start in 0.00% pseudo baseline with kernel, -1500 ms: starts before the trial start in 0.00% true baseline with kernel, -900 ms: starts before the trial start in 0.00% SDF end, +600 ms: after the trial stop in 0.00%
Result: blocked (methods underspecified). With the detection fixed in the claims file, a threshold crossing in the 200 ms after a pseudo onset 600 ms before the torque occurs in 0.23 of the 544 cell-directions, against the 0.01 the paper's threshold implies; the t-test variant gives 0.46. Unperturbed trials at the usual torque time give 0.19 and 0.44. Stationary Poisson spike trains at each cell-direction's baseline rate give 0.21 and 0.48 with the seed of the earlier attempt, and 0.16 to 0.19 (main detection) with eight other seeds, so most of the excess comes from the detection procedure, not from activity in the hold period. The rate does not fall with baseline rate or number of trials (0.18 to 0.29 by baseline rate quartile, 0.21 to 0.26 by number of trials). A response with onset between 20 and 60 ms is found in 0.048 of the cell-directions at the pseudo onset and 0.047 on unperturbed trials, against 0.375 at the true onset, a ratio of 0.13. Fourteen cell-directions (10 of them CSNs) have a baseline SD of exactly 0 and can never respond. Near-silent cell-directions with a threshold of about 0 can respond, and add to the false positives: v5801.1 has a baseline rate of 0.00 spikes/s and an extension "response" at 30 ms. The spike windows lie inside their trials: every torque onset comes at least 2968 ms after its trial's start and 1516 ms before its stop, so no window (from -1800 to +600 ms) is cut by a trial boundary. Both rates are about 23 and 46 times the paper's, so the procedure here does not give the rate the paper states for its own: either it differs from the paper's in a way the Methods do not settle, or the paper's stated rate does not hold, and this analysis cannot tell which. The claims that rest on the selection of responses are blocked with the same reason. A threshold calibrated to a pseudo-onset rate of 0.01 (t = 6.67 against 3.11) finds 30.4% and 32.5% of M1 cells torque-related, further from the paper's 75% and 68%; taking the first crossing anywhere in 0 to 200 ms as the onset removes 5 of the 204 responses. The analyses of the blocked claims are run below all the same, and their numbers are reported.
torque-related-fractions¶
In monkey V, 84 of 112 M1 neurons (75%) were torque-related before MPTP and 108 of 159 (68%) after; PTNs 27 of 49 (55%) and 37 of 54 (68%); CSNs 14 of 31 (45%) and 10 of 48 (21%).
Result: blocked (methods underspecified), as the selection of responses is (see false-positive-rate). The detection here finds 54.8% of monkey V cells
torque-related before MPTP and 47.8% after (paper 75% and 68%), PTNs 72.5% and 68.5% (55% and 68%),
CSNs 15.6% and 16.7% (45% and 21%); only the post-MPTP PTN value is within tolerance. The t-test
variant gives 41.7% and 60.5%.
prevalence-unchanged¶
MPTP did not change the fraction of torque-related neurons, in M1, in PTNs or in CSNs (chi-square, p > 0.05/3).
def chi2_yates(a_k, a_n, b_k, b_n):
table = np.array([[a_k, a_n - a_k], [b_k, b_n - b_k]])
chi2, p, _, _ = stats.chi2_contingency(table, correction=True)
return chi2, p
rows = []
for label, col in (("main", "torque_related"), ("t-test variant", "torque_related_trials")):
for g in ("M1", "PTN", "CSN"):
sub = cells if g == "M1" else cells[cells.ptype == g]
a, b = sub[sub.state == "pre"][col], sub[sub.state == "post"][col]
chi2, p = chi2_yates(a.sum(), len(a), b.sum(), len(b))
rows.append({"detection": label, "group": g, "pre_%": 100 * a.mean(), "post_%": 100 * b.mean(),
"change_points": 100 * (b.mean() - a.mean()), "chi2": chi2, "p": p, "p_below_0.05/3": p < 0.05 / 3})
print(pd.DataFrame(rows).round(4).to_string(index=False))
print("\nThe paper's monkey V counts, through the same test:",
{g: round(chi2_yates(*k)[1], 4) for g, k in {"M1": (84, 112, 108, 159), "PTN": (27, 49, 37, 54), "CSN": (14, 31, 10, 48)}.items()})
detection group pre_% post_% change_points chi2 p p_below_0.05/3
main M1 54.7826 47.7707 -7.0119 1.0403 0.3078 False
main PTN 72.5490 68.5185 -4.0305 0.0569 0.8115 False
main CSN 15.6250 16.6667 1.0417 0.0000 1.0000 False
t-test variant M1 41.7391 60.5096 18.7704 8.6420 0.0033 True
t-test variant PTN 49.0196 74.0741 25.0545 5.9595 0.0146 True
t-test variant CSN 31.2500 37.5000 6.2500 0.1122 0.7377 False
The paper's monkey V counts, through the same test: {'M1': np.float64(0.26), 'PTN': np.float64(0.2307), 'CSN': np.float64(0.0408)}
Result: blocked (methods underspecified), as the selection of responses is (see false-positive-rate). By the paper's test the fractions do not change in M1 (p = 0.31),
PTNs (p = 0.81) or CSNs (p = 1.0) with the main detection; the t-test variant finds an increase in
M1 (p = 0.0033) and PTNs (p = 0.015), both below 0.05/3.
ptn-more-responsive¶
PTNs responded to torque more often than CSNs (68% against 31%, chi-square 32.3, p < 0.001).
for label, col in (("main", "torque_related"), ("t-test variant", "torque_related_trials")):
p_, c_ = cells[cells.ptype == "PTN"][col], cells[cells.ptype == "CSN"][col]
chi2, p = chi2_yates(p_.sum(), len(p_), c_.sum(), len(c_))
print(f"{label}: PTN {int(p_.sum())}/{len(p_)} ({100 * p_.mean():.1f}%), CSN {int(c_.sum())}/{len(c_)} "
f"({100 * c_.mean():.1f}%), chi-square {chi2:.1f}, p = {p:.2e}; expected 62.1% and 30.4%")
main: PTN 74/105 (70.5%), CSN 13/80 (16.2%), chi-square 51.4, p = 7.38e-13; expected 62.1% and 30.4% t-test variant: PTN 65/105 (61.9%), CSN 28/80 (35.0%), chi-square 12.1, p = 5.06e-04; expected 62.1% and 30.4%
Result: blocked (methods underspecified), as the selection of responses is (see false-positive-rate). PTNs respond more often than CSNs (70.5% against 16.2%,
chi-square 51.4, p = 7.4e-13); both fractions are outside tolerance, the PTN fraction above the
expected 62.1 +/- 6.2% and the CSN fraction below the expected 30.4 +/- 3.0%. The t-test variant
gives 61.9% and 35.0% (p = 5.1e-4).
one-direction-half¶
About half of the torque-related cells (124 of 284, 44%) responded to only one direction.
tr_cells = cells[cells.torque_related]
one = (tr_cells.n_directions == 1)
print(f"Torque-related cells responding in one direction only: {int(one.sum())}/{len(tr_cells)} "
f"({100 * one.mean():.1f}%); by state: {tr_cells.groupby('state').apply(lambda d: round(100 * (d.n_directions == 1).mean(), 1)).to_dict()}")
vt = cells[cells.torque_related_trials]
n_dir_t = var["responds"].sum(axis=1).values[cells.torque_related_trials.values]
print(f"t-test variant: {int((n_dir_t == 1).sum())}/{len(vt)} ({100 * (n_dir_t == 1).mean():.1f}%)")
Torque-related cells responding in one direction only: 72/138 (52.2%); by state: {'post': 46.7, 'pre': 58.7}
t-test variant: 92/143 (64.3%)
Result: blocked (methods underspecified), as the selection of responses is (see false-positive-rate). 52.2% of torque-related cells (72 of 138) respond to one direction
only (58.7% before MPTP, 46.7% after), against the expected 44 +/- 4.4%; the t-test variant gives
64.3%.
excitatory-majority¶
Most torque responses were increases in firing (86% and 88% of M1 responses; 70% and 75% for PTNs), and MPTP did not change the PTN balance (chi-square 0.1, p = 0.7).
rows = []
for g in ("M1", "PTN", "CSN"):
sub = resp if g == "M1" else resp[resp.ptype == g]
for st in ("pre", "post"):
x = sub[sub.state == st]
rows.append({"group": g, "state": st, "increases": int((x.sign > 0).sum()), "responses": len(x),
"percent": 100 * (x.sign > 0).mean(),
"flexion": f"{int(((x.sign > 0) & (x.direction == 'flexion')).sum())}/{int((x.direction == 'flexion').sum())}",
"extension": f"{int(((x.sign > 0) & (x.direction == 'extension')).sum())}/{int((x.direction == 'extension').sum())}"})
inc = pd.DataFrame(rows)
print(inc.round(1).to_string(index=False))
for g in ("PTN", "CSN"):
a, b = inc[(inc.group == g) & (inc.state == "pre")].iloc[0], inc[(inc.group == g) & (inc.state == "post")].iloc[0]
chi2, p = chi2_yates(a.increases, a.responses, b.increases, b.responses)
print(f"{g} increases pre against post: chi-square {chi2:.2f}, p = {p:.3f}")
print("All responses, both states:", f"{100 * (resp.sign > 0).mean():.1f}% increases of {len(resp)}")
chi2, p = chi2_yates(23, 24, 19, 22)
print(f"Reference, the paper's own CSN counts in Table 1 (23/24 against 19/22): chi-square {chi2:.2f}, p = {p:.2f} "
"(Table 1 marks this p < 0.01)")
group state increases responses percent flexion extension M1 pre 74 89 83.1 41/49 33/40 M1 post 88 115 76.5 48/63 40/52 PTN pre 35 50 70.0 19/27 16/23 PTN post 35 59 59.3 20/33 15/26 CSN pre 6 6 100.0 4/4 2/2 CSN post 8 10 80.0 5/6 3/4 PTN increases pre against post: chi-square 0.92, p = 0.338 CSN increases pre against post: chi-square 0.15, p = 0.696 All responses, both states: 79.4% increases of 204 Reference, the paper's own CSN counts in Table 1 (23/24 against 19/22): chi-square 0.38, p = 0.54 (Table 1 marks this p < 0.01)
Result: blocked (methods underspecified), as the selection of responses is (see false-positive-rate). Increases make up 83.1% and 76.5% of M1 responses and 70.0% and
59.3% of PTN responses (expected 86, 88, 70 and 75%); the PTN balance does not change (chi-square
0.92, p = 0.34). The paper's Table 1 marks its CSN change p < 0.01, but its own counts give p =
0.54 with the test it names.
magnitude-unchanged and venus-response-magnitudes¶
The magnitude of torque responses did not differ before and after MPTP, for M1 cells or for PTNs, by any of three measures (Kolmogorov-Smirnov and Mann-Whitney, all p > 0.15).
In monkey V the mean response magnitudes were 42.4 and 37.9 spikes/s (M1 flexion), 50.9 and 42.6 (M1 extension), 33.4 and 26.1 (PTN flexion), 38.8 and 29.3 (PTN extension).
measures = ["mean_change", "max_change", "auc"]
rows = []
for g in ("M1", "PTN"):
sub = resp if g == "M1" else resp[resp.ptype == g]
for d in ("both", "flexion", "extension"):
x = sub if d == "both" else sub[sub.direction == d]
a, b = x[x.state == "pre"], x[x.state == "post"]
for m in measures:
rows.append({"group": g, "direction": d, "measure": m, "n_pre": len(a), "n_post": len(b),
"mean_pre": a[m].mean(), "mean_post": b[m].mean(),
"ks_p": stats.ks_2samp(a[m], b[m]).pvalue, "mwu_p": mwu(a[m], b[m])})
mag = pd.DataFrame(rows)
print(mag.round(4).to_string(index=False))
main_tests = mag[mag.direction == "both"]
print(f"\nMain (directions pooled): smallest p {main_tests[['ks_p', 'mwu_p']].min().min():.4f}; "
f"tests with p < 0.05: {int((main_tests[['ks_p', 'mwu_p']] < 0.05).sum().sum())} of 12")
for m in measures:
print(f"{m} range: {resp[m].min():.2f} to {resp[m].max():.2f}")
group direction measure n_pre n_post mean_pre mean_post ks_p mwu_p M1 both mean_change 89 115 27.3362 20.3560 0.0002 0.0002 M1 both max_change 89 115 43.5538 32.1179 0.0010 0.0015 M1 both auc 89 115 4.3516 1.9838 0.0254 0.0105 M1 flexion mean_change 49 63 26.3511 20.3303 0.0063 0.0066 M1 flexion max_change 49 63 40.2113 31.2246 0.0111 0.0202 M1 flexion auc 49 63 4.6622 1.8001 0.0732 0.0481 M1 extension mean_change 40 52 28.5429 20.3871 0.0118 0.0097 M1 extension max_change 40 52 47.6484 33.2000 0.0577 0.0221 M1 extension auc 40 52 3.9712 2.2063 0.1637 0.0989 PTN both mean_change 50 59 23.6217 15.0402 0.0003 0.0000 PTN both max_change 50 59 33.8026 21.6409 0.0028 0.0007 PTN both auc 50 59 3.4824 1.1372 0.0113 0.0451 PTN flexion mean_change 27 33 22.5249 14.2037 0.0063 0.0010 PTN flexion max_change 27 33 31.7392 19.9486 0.0209 0.0041 PTN flexion auc 27 33 4.4932 0.9628 0.0209 0.0465 PTN extension mean_change 23 26 24.9092 16.1019 0.0127 0.0075 PTN extension max_change 23 26 36.2250 23.7888 0.1357 0.0557 PTN extension auc 23 26 2.2958 1.3585 0.5238 0.4646 Main (directions pooled): smallest p 0.0000; tests with p < 0.05: 12 of 12 mean_change range: 1.33 to 130.45 max_change range: 1.48 to 211.86 auc range: 0.01 to 25.46
Result: blocked (methods underspecified), as the selection of responses is (see false-positive-rate). Here the responses are smaller after MPTP: all 12 tests with
directions pooled have p < 0.05 (for example M1 maximum change 43.6 against 32.1 spikes/s,
Kolmogorov-Smirnov p = 0.0010; PTN mean change 23.6 against 15.0, Mann-Whitney p < 0.0001).
expected_mag = {("M1", "flexion"): (42.4, 37.9), ("M1", "extension"): (50.9, 42.6),
("PTN", "flexion"): (33.4, 26.1), ("PTN", "extension"): (38.8, 29.3)}
rows = []
for reading, col in (("a max change", "max_change"), ("b mean rate", "mean_rate"), ("c mean change", "mean_change")):
within = 0
for (g, d), (e_pre, e_post) in expected_mag.items():
sub = resp if g == "M1" else resp[resp.ptype == g]
x = sub[sub.direction == d]
for st, e in (("pre", e_pre), ("post", e_post)):
v = x[x.state == st][col].mean()
within += abs(v - e) <= 0.1 * e
rows.append({"reading": reading, "group": g, "direction": d, "state": st, "value": v, "expected": e,
"in_tol": abs(v - e) <= 0.1 * e, "n": int((x.state == st).sum())})
mags = pd.DataFrame(rows)
print(mags.round(2).to_string(index=False))
print("\nWithin tolerance per reading (of 8):", mags.groupby("reading").in_tol.sum().to_dict())
reading group direction state value expected in_tol n
a max change M1 flexion pre 40.21 42.4 True 49
a max change M1 flexion post 31.22 37.9 False 63
a max change M1 extension pre 47.65 50.9 True 40
a max change M1 extension post 33.20 42.6 False 52
a max change PTN flexion pre 31.74 33.4 True 27
a max change PTN flexion post 19.95 26.1 False 33
a max change PTN extension pre 36.22 38.8 True 23
a max change PTN extension post 23.79 29.3 False 26
b mean rate M1 flexion pre 37.24 42.4 False 49
b mean rate M1 flexion post 28.45 37.9 False 63
b mean rate M1 extension pre 37.40 50.9 False 40
b mean rate M1 extension post 28.06 42.6 False 52
b mean rate PTN flexion pre 33.99 33.4 True 27
b mean rate PTN flexion post 21.75 26.1 False 33
b mean rate PTN extension pre 34.45 38.8 False 23
b mean rate PTN extension post 24.58 29.3 False 26
c mean change M1 flexion pre 26.35 42.4 False 49
c mean change M1 flexion post 20.33 37.9 False 63
c mean change M1 extension pre 28.54 50.9 False 40
c mean change M1 extension post 20.39 42.6 False 52
c mean change PTN flexion pre 22.52 33.4 False 27
c mean change PTN flexion post 14.20 26.1 False 33
c mean change PTN extension pre 24.91 38.8 False 23
c mean change PTN extension post 16.10 29.3 False 26
Within tolerance per reading (of 8): {'a max change': 4, 'b mean rate': 1, 'c mean change': 0}
Result: blocked (methods underspecified), as the selection of responses is (see false-positive-rate). The maximum change matches the four pre-MPTP values and
none of the post-MPTP ones (4 of 8 within tolerance); the other readings match 1 and 0 of 8.
latency-shorter¶
Response latencies were shorter after MPTP for M1 cells and for PTNs, for flexion and extension (Kolmogorov-Smirnov, all p < 0.05). Monkey V means 39.2 and 34.6 ms (M1 flexion), 40.4 and 33.4 (M1 extension), 43.9 and 35.5 (PTN flexion), 44.4 and 36.7 (PTN extension).
expected_lat = {("M1", "flexion"): (39.2, 34.56), ("M1", "extension"): (40.4, 33.4),
("PTN", "flexion"): (43.9, 35.5), ("PTN", "extension"): (44.4, 36.7)}
def latency_table(rsp):
rows = []
for (g, d), (e_pre, e_post) in expected_lat.items():
sub = rsp if g == "M1" else rsp[rsp.ptype == g]
x = sub[sub.direction == d]
a, b = x[x.state == "pre"].latency, x[x.state == "post"].latency
rows.append({"group": g, "direction": d, "n_pre": len(a), "pre": a.mean(), "expected_pre": e_pre,
"in_tol_pre": abs(a.mean() - e_pre) <= 0.1 * e_pre, "n_post": len(b), "post": b.mean(),
"expected_post": e_post, "in_tol_post": abs(b.mean() - e_post) <= 0.1 * e_post,
"ks_p": stats.ks_2samp(a, b).pvalue, "mwu_p": mwu(a, b)})
return pd.DataFrame(rows)
lat = latency_table(resp)
print(lat.round(4).to_string(index=False))
print("\nPooled means (paper, both monkeys: M1 41 and 35 ms, PTN 46 and 38 ms):")
for g in ("M1", "PTN"):
sub = resp if g == "M1" else resp[resp.ptype == g]
cs = cells[cells.torque_related] if g == "M1" else cells[cells.torque_related & (cells.ptype == g)]
print(f" {g}: responses {sub[sub.state == 'pre'].latency.mean():.1f} and {sub[sub.state == 'post'].latency.mean():.1f} ms; "
f"earliest per cell {cs[cs.state == 'pre'].latency.mean():.1f} and {cs[cs.state == 'post'].latency.mean():.1f} ms "
f"(KS p = {stats.ks_2samp(cs[cs.state == 'pre'].latency, cs[cs.state == 'post'].latency).pvalue:.4f})")
resp_t = R[(R.method == "trials") & ~R.pseudo & R.responds]
print("\nt-test variant:")
print(latency_table(resp_t)[["group", "direction", "n_pre", "pre", "n_post", "post", "ks_p"]].round(4).to_string(index=False))
fig, axes = plt.subplots(1, 2, figsize=(9, 3.2), sharey=True)
for ax, g in zip(axes, ("M1", "PTN")):
sub = resp if g == "M1" else resp[resp.ptype == g]
for st, col in (("pre", "k"), ("post", "0.6")):
for d, ls in (("flexion", "-"), ("extension", "--")):
v = np.sort(sub[(sub.state == st) & (sub.direction == d)].latency)
ax.step(v, np.arange(1, len(v) + 1) / len(v), color=col, ls=ls, where="post", label=f"{st} {d}")
ax.set_title(g)
ax.set_xlabel("response onset (ms)")
axes[0].set_ylabel("cumulative fraction")
axes[1].legend(fontsize=7)
plt.tight_layout()
plt.show()
print("Plotted latencies span", resp.latency.min(), "to", resp.latency.max(), "ms; none outside the axes.")
group direction n_pre pre expected_pre in_tol_pre n_post post expected_post in_tol_post ks_p mwu_p M1 flexion 49 41.4694 39.2 True 63 36.2063 34.56 True 0.0045 0.0086 M1 extension 40 42.2750 40.4 True 52 36.2115 33.40 True 0.0072 0.0110 PTN flexion 27 44.6296 43.9 True 33 38.1818 35.50 True 0.0075 0.0184 PTN extension 23 45.3043 44.4 True 26 38.5000 36.70 True 0.0127 0.0161 Pooled means (paper, both monkeys: M1 41 and 35 ms, PTN 46 and 38 ms): M1: responses 41.8 and 36.2 ms; earliest per cell 40.9 and 35.0 ms (KS p = 0.0033) PTN: responses 44.9 and 38.3 ms; earliest per cell 44.0 and 35.9 ms (KS p = 0.0042) t-test variant:
group direction n_pre pre n_post post ks_p M1 flexion 31 43.7097 65 40.2615 0.3066 M1 extension 32 42.3750 66 41.6212 0.7331 PTN flexion 19 46.8947 26 41.7308 0.1702 PTN extension 14 45.4286 27 44.2963 0.4625
Plotted latencies span 20.0 to 60.0 ms; none outside the axes.
Result: blocked (methods underspecified), as the selection of responses is (see false-positive-rate). With the main detection, latencies are shorter after MPTP in all four
comparisons (Kolmogorov-Smirnov p = 0.0045, 0.0072, 0.0075 and 0.0127) and all eight means are
within tolerance (for example M1 flexion 41.5 and 36.2 ms, PTN extension 45.3 and 38.5 ms). The
t-test variant finds no significant shortening (p 0.17 to 0.73).
latency-not-training¶
Latencies did not correlate with the days of training (monkey V Spearman |rho| = 0.02, p = 0.85).
tr_cells = cells[cells.torque_related].copy()
first_day = sessions[sessions.subject == "Venus"].date.min()
tr_cells["day"] = (tr_cells.date - first_day).dt.days
rho, p = stats.spearmanr(tr_cells.day, tr_cells.latency)
print(f"All torque-related cells (n = {len(tr_cells)}): Spearman rho = {rho:.3f}, p = {p:.4f}; days {tr_cells.day.min()} to {tr_cells.day.max()}")
for st in ("pre", "post"):
x = tr_cells[tr_cells.state == st]
rho, p = stats.spearmanr(x.day, x.latency)
print(f" {st}: n = {len(x)}, rho = {rho:.3f}, p = {p:.4f}")
All torque-related cells (n = 138): Spearman rho = -0.249, p = 0.0033; days 0 to 360 pre: n = 63, rho = 0.084, p = 0.5106 post: n = 75, rho = -0.147, p = 0.2083
Result: blocked (methods underspecified), as the selection of responses is (see false-positive-rate). Latency correlates with the day of recording over both states
(rho = -0.25, p = 0.0033, n = 138), as post-MPTP sessions come later; within each state it does not
(p = 0.51 and 0.21).
fwhm-shorter-ptn¶
Torque-evoked increases in PTN firing were briefer after MPTP (FWHM 73 to 50 ms over both monkeys; Mann-Whitney p < 0.05 for each direction). Monkey V: 97 and 52 ms (flexion), 57 and 41 ms (extension).
expected_fwhm = {"flexion": (97, 52), "extension": (57, 41)}
ptn = resp[resp.ptype == "PTN"]
rows = []
for subset, x in (("increases", ptn[ptn.sign > 0]), ("all PTN responses", ptn)):
for d, (e_pre, e_post) in expected_fwhm.items():
y = x[x.direction == d]
a, b = y[y.state == "pre"].fwhm, y[y.state == "post"].fwhm
rows.append({"responses": subset, "direction": d, "n_pre": len(a), "pre": a.mean(),
"in_tol_pre": abs(a.mean() - e_pre) <= 0.1 * e_pre, "n_post": len(b), "post": b.mean(),
"in_tol_post": abs(b.mean() - e_post) <= 0.1 * e_post, "mwu_p": mwu(a, b)})
print(pd.DataFrame(rows).round(4).to_string(index=False))
inc = ptn[ptn.sign > 0]
print(f"\nPTN increases, directions pooled: {inc[inc.state == 'pre'].fwhm.mean():.1f} and {inc[inc.state == 'post'].fwhm.mean():.1f} ms "
f"(paper 73 and 50), p = {mwu(inc[inc.state == 'pre'].fwhm, inc[inc.state == 'post'].fwhm):.4f}; FWHM range {ptn.fwhm.min():.0f} to {ptn.fwhm.max():.0f} ms")
responses direction n_pre pre in_tol_pre n_post post in_tol_post mwu_p
increases flexion 19 91.5789 True 20 40.0500 False 0.0273
increases extension 16 68.8125 False 15 41.7333 True 0.3731
all PTN responses flexion 27 95.5556 True 33 50.8485 True 0.0108
all PTN responses extension 23 65.1304 False 26 60.7308 False 0.7637
PTN increases, directions pooled: 81.2 and 40.8 ms (paper 73 and 50), p = 0.0339; FWHM range 17 to 382 ms
Result: blocked (methods underspecified), as the selection of responses is (see false-positive-rate). PTN increases are briefer after MPTP for flexion (91.6 against 40.1
ms, p = 0.027) but not for extension (68.8 against 41.7 ms, p = 0.37); two of the four means are
within tolerance.
dsi-reduced, venus-dsi and dsi-not-baseline¶
Directional selectivity fell after MPTP; reciprocal and unidirectional responses became less common (chi-square 22.2, p < 0.001).
In monkey V the mean DSI was 0.72 and 0.66 (M1), 0.73 and 0.75 (PTN), 0.81 and 0.83 (CSN).
Pre-torque rates did not differ between directional categories (one-way ANOVA, F = 0.77, p = 0.46).
cats = ["reciprocal", "unidirectional", "bidirectional", "non-directional"]
tr_cells = cells[cells.torque_related]
tab = pd.crosstab(tr_cells.category, tr_cells.state).reindex(index=cats, columns=["pre", "post"]).fillna(0)
print(tab.to_string())
chi2, p, dof, _ = stats.chi2_contingency(tab.values)
high = tr_cells.category.isin(["reciprocal", "unidirectional"])
fr = high.groupby(tr_cells.state).mean()
print(f"\nReciprocal or unidirectional: pre {100 * fr['pre']:.1f}%, post {100 * fr['post']:.1f}% "
f"(change {100 * (fr['post'] - fr['pre']):.1f} points); chi-square (2 x 4) {chi2:.2f}, dof {dof}, p = {p:.4f}")
for g in ("PTN", "CSN"):
x = tr_cells[tr_cells.ptype == g]
t = pd.crosstab(x.category, x.state).reindex(index=cats, columns=["pre", "post"]).fillna(0)
t = t[t.sum(axis=1) > 0]
c2, pp, _, _ = stats.chi2_contingency(t.values)
h = x.category.isin(["reciprocal", "unidirectional"]).groupby(x.state).mean()
print(f" {g}: reciprocal or unidirectional pre {100 * h.get('pre', np.nan):.1f}%, post {100 * h.get('post', np.nan):.1f}%; chi-square p = {pp:.3f}")
# Variant: categories from the DSI value as the paper words them: reciprocal DSI > 1, unidirectional
# a response in one direction only (DSI <= 1), bidirectional 0.5 < DSI <= 1, non-directional DSI <= 0.5.
def by_value(r):
if r.dsi > 1:
return "reciprocal"
if r.n_directions == 1:
return "unidirectional"
return "bidirectional" if r.dsi > 0.5 else "non-directional"
tr_v = tr_cells.assign(category_value=[by_value(r) for r in tr_cells.itertuples()])
tab_v = pd.crosstab(tr_v.category_value, tr_v.state).reindex(index=cats, columns=["pre", "post"]).fillna(0)
chi2_v, p_v, _, _ = stats.chi2_contingency(tab_v.values)
fr_v = tr_v.category_value.isin(["reciprocal", "unidirectional"]).groupby(tr_v.state).mean()
print("\nVariant, categories by DSI value:")
print(tab_v.to_string())
print(f"Reciprocal or unidirectional: pre {100 * fr_v['pre']:.1f}%, post {100 * fr_v['post']:.1f}%; chi-square {chi2_v:.2f}, p = {p_v:.4f}")
print("Smallest expected count in the main 2 x 4 table:", round(stats.chi2_contingency(tab.values)[3].min(), 1))
state pre post category reciprocal 0 4 unidirectional 37 35 bidirectional 8 7 non-directional 18 29 Reciprocal or unidirectional: pre 58.7%, post 52.0% (change -6.7 points); chi-square (2 x 4) 5.70, dof 3, p = 0.1274 PTN: reciprocal or unidirectional pre 64.9%, post 51.4%; chi-square p = 0.029 CSN: reciprocal or unidirectional pre 80.0%, post 75.0%; chi-square p = 1.000 Variant, categories by DSI value: state pre post category_value reciprocal 14 9 unidirectional 25 29 bidirectional 6 8 non-directional 18 29 Reciprocal or unidirectional: pre 61.9%, post 50.7%; chi-square 3.22, p = 0.3583 Smallest expected count in the main 2 x 4 table: 1.8
Result: blocked (methods underspecified), as the selection of responses is (see false-positive-rate). Reciprocal and unidirectional cells fall from 58.7% to 52.0% of
torque-related cells, by 6.7 points, without significance (chi-square 5.70, p = 0.13). With
categories set from the DSI value as the paper words them (reciprocal above 1), the fraction falls
from 61.9% to 50.7% (p = 0.36). The main table has no reciprocal cell before MPTP; its smallest
expected count is 1.8, so the chi-square p is approximate. The DSI's P and NP are the largest
changes over 0 to 200 ms, whether or not a response was detected there, a choice made here.
expected_dsi = {("M1", "pre"): (0.72, 0.087), ("M1", "post"): (0.66, 0.077), ("PTN", "pre"): (0.73, 0.115),
("PTN", "post"): (0.75, 0.132), ("CSN", "pre"): (0.81, 0.214), ("CSN", "post"): (0.83, 0.253)}
rows = []
for (g, st), (e, tol) in expected_dsi.items():
x = tr_cells if g == "M1" else tr_cells[tr_cells.ptype == g]
v = x[x.state == st].dsi
rows.append({"group": g, "state": st, "n": len(v), "mean_dsi": v.mean(), "sd": v.std(), "expected": e,
"in_tol": abs(v.mean() - e) <= tol})
print(pd.DataFrame(rows).round(3).to_string(index=False))
print("DSI range:", round(tr_cells.dsi.min(), 3), "to", round(tr_cells.dsi.max(), 3))
group state n mean_dsi sd expected in_tol M1 pre 63 0.620 0.533 0.72 False M1 post 75 0.463 0.437 0.66 False PTN pre 37 0.759 0.575 0.73 True PTN post 37 0.460 0.467 0.75 False CSN pre 5 0.662 0.516 0.81 True CSN post 8 0.695 0.415 0.83 True DSI range: 0.0 to 1.905
Result: blocked (methods underspecified), as the selection of responses is (see false-positive-rate). Three of the six mean DSIs are within tolerance (PTN pre and both CSN
values); M1 is 0.62 and 0.46 and PTN post 0.46, against 0.72, 0.66 and 0.75.
groups = [tr_cells[tr_cells.category == c].baseline_rate for c in cats]
f4, p4 = stats.f_oneway(*groups)
merged = tr_cells.category.replace({"reciprocal": "highly directional", "unidirectional": "highly directional"})
f3, p3 = stats.f_oneway(*[tr_cells[merged == c].baseline_rate for c in merged.unique()])
print("Baseline rate by category (mean, n):", {c: (round(g.mean(), 2), len(g)) for c, g in zip(cats, groups)})
print(f"One-way ANOVA, four categories: F = {f4:.2f}, p = {p4:.4f}; three groups: F = {f3:.2f}, p = {p3:.4f}")
Baseline rate by category (mean, n): {'reciprocal': (np.float64(15.42), 4), 'unidirectional': (np.float64(13.17), 72), 'bidirectional': (np.float64(13.04), 15), 'non-directional': (np.float64(15.36), 47)}
One-way ANOVA, four categories: F = 0.40, p = 0.7497; three groups: F = 0.54, p = 0.5855
Result: blocked (methods underspecified), as the selection of responses is (see false-positive-rate). Baseline rates do not differ between the four categories (F = 0.40,
p = 0.75) or the three groups (F = 0.54, p = 0.59).
Kinematics and EMG¶
Session measures use monkey V files with analog data. No file repeats another's recording in this version, so none is left out (in the published version four were). Kinematics are measured on every perturbation; a trial whose joint was already moving faster than 5 deg/s in the 50 ms before onset is counted and kept.
# Per-trial torque command step, in multiples of the smallest session level (as in the data check)
level_unit = command.step[command.step > 10].min()
trial_level = {}
for s in venus:
w = s["windows"]["torque_cmd"]
sign = np.where(s["pert_direction"] == "flexion", 1.0, -1.0)
rel = (w - np.nanmedian(w[:, (AX >= -100) & (AX < 0)], axis=1, keepdims=True)) * sign[:, None]
trial_level[s["session_id"].split("++")[1]] = np.round(np.nanmedian(rel[:, (AX >= 10) & (AX < 40)], axis=1) / level_unit)
kin_rows = []
for s in venus:
if s["session_id"].split("++")[1] in repeated_files:
continue
tl = trial_level[s["session_id"].split("++")[1]]
for k, d in enumerate(s["pert_direction"]):
r = stretch.kinematics(s["windows"]["angle"][k], d, s["analog_step_ms"])
if r["valid"]:
kin_rows.append({"session": s["session_id"], "post": s["post"], "direction": d,
"trial_level": tl[k], **r})
kin = pd.DataFrame(kin_rows)
kin["state"] = np.where(kin.post, "post", "pre")
print(f"{len(kin)} perturbations in {kin.session.nunique()} sessions; moving before onset: {int(kin.moving_before.sum())}; "
f"no velocity crossing within 150 ms: {int(kin.latency_ms.isna().sum())}")
print("Ranges: latency", kin.latency_ms.min(), "to", kin.latency_ms.max(), "ms; peak velocity",
round(kin.peak_velocity.min(), 1), "to", round(kin.peak_velocity.max(), 1), "deg/s; amplitude",
round(kin.amplitude.min(), 2), "to", round(kin.amplitude.max(), 2), "deg")
sess_kin = kin.groupby(["session", "state", "direction"])[["latency_ms", "peak_velocity", "amplitude"]].mean().reset_index()
sess_kin["level"] = sess_kin.session.map(command.set_index("session").level)
print("\nSession means, mean over sessions:")
print(sess_kin.groupby(["direction", "state"])[["latency_ms", "peak_velocity", "amplitude"]].agg(["mean", "count"]).round(2).to_string())
5550 perturbations in 291 sessions; moving before onset: 503; no velocity crossing within 150 ms: 4
Ranges: latency 0.0 to 40.0 ms; peak velocity -5.0 to 572.4 deg/s; amplitude -1.68 to 13.88 deg
Session means, mean over sessions:
latency_ms peak_velocity amplitude
mean count mean count mean count
direction state
extension post 5.82 137 58.08 137 2.86 137
pre 10.86 154 46.10 154 2.41 154
flexion post 3.05 137 75.37 137 3.78 137
pre 9.65 153 58.11 153 3.39 153
kinematics-slower¶
After MPTP, torque-induced displacements were slower and smaller; mean peak velocity fell by 16.2% and amplitude by 6.6% (ANOVA, F > 515 and F > 49).
def anova(df, measure):
model = ols(f"{measure} ~ C(state) * C(direction)", data=df).fit()
return sm.stats.anova_lm(model, typ=2)
for m in ("peak_velocity", "amplitude"):
print(f"{m}:")
print(anova(sess_kin, m).round(4).to_string())
print()
rows = []
for reading, x in (("a both directions", sess_kin), ("b flexion", sess_kin[sess_kin.direction == "flexion"]),
("extension, for reference", sess_kin[sess_kin.direction == "extension"])):
g = x.groupby("state")[["peak_velocity", "amplitude"]].mean()
rows.append({"reading": reading, "velocity_change_%": 100 * (g.loc["post", "peak_velocity"] / g.loc["pre", "peak_velocity"] - 1),
"amplitude_change_%": 100 * (g.loc["post", "amplitude"] / g.loc["pre", "amplitude"] - 1)})
kin_change = pd.DataFrame(rows)
kin_change["velocity_in_tol"] = (kin_change["velocity_change_%"] + 16.2).abs() <= 1.62
kin_change["amplitude_in_tol"] = (kin_change["amplitude_change_%"] + 6.6).abs() <= 0.66
print(kin_change.round(2).to_string(index=False))
lvl2 = sess_kin[sess_kin.level == 2]
g = lvl2.groupby("state")[["peak_velocity", "amplitude"]].mean()
print(f"\nVariant, sessions with the middle command size only ({lvl2[lvl2.state == 'pre'].session.nunique()} pre, "
f"{lvl2[lvl2.state == 'post'].session.nunique()} post): peak velocity {g.loc['pre', 'peak_velocity']:.1f} and "
f"{g.loc['post', 'peak_velocity']:.1f} deg/s ({100 * (g.loc['post', 'peak_velocity'] / g.loc['pre', 'peak_velocity'] - 1):.1f}%), "
f"amplitude {g.loc['pre', 'amplitude']:.2f} and {g.loc['post', 'amplitude']:.2f} deg "
f"({100 * (g.loc['post', 'amplitude'] / g.loc['pre', 'amplitude'] - 1):.1f}%)")
for m in ("peak_velocity", "amplitude"):
print(f" {m}: MPTP F = {anova(lvl2, m).loc['C(state)', 'F']:.2f}, p = {anova(lvl2, m).loc['C(state)', 'PR(>F)']:.4f}")
print("\nBy command size (session means, both directions):")
print(sess_kin.groupby(["level", "state"])[["peak_velocity", "amplitude", "latency_ms"]].agg(["mean", "count"]).round(1).to_string())
fig, axes = plt.subplots(1, 3, figsize=(10, 3))
for ax, m in zip(axes, ("latency_ms", "peak_velocity", "amplitude")):
for i, (d, st) in enumerate([(d, st) for d in ("flexion", "extension") for st in ("pre", "post")]):
v = sess_kin[(sess_kin.direction == d) & (sess_kin.state == st)][m]
ax.scatter(np.full(len(v), i) + np.random.default_rng(0).uniform(-0.15, 0.15, len(v)), v, s=6,
color="k" if st == "pre" else "0.6")
ax.set_xticks(range(4), ["flex pre", "flex post", "ext pre", "ext post"], fontsize=7)
ax.set_title(m)
plt.tight_layout()
plt.show()
print("All session means are plotted (ranges printed above).")
print(f"\nAbsolute scale: torque-induced amplitude, mean of session means {sess_kin.amplitude.mean():.2f} deg "
f"(paper: about 10 deg); peak velocity {sess_kin.peak_velocity.mean():.1f} deg/s, highest session mean "
f"{sess_kin.peak_velocity.max():.1f} (paper's matched sessions 130 to 200 deg/s). Active movements in the "
f"same files, same filter: median peak velocity {act[act.subject == 'Venus'].active_peak_velocity.median():.1f} "
f"deg/s and amplitude {act[act.subject == 'Venus'].active_amplitude.median():.1f} deg.")
peak_velocity:
sum_sq df F PR(>F)
C(state) 30935.6834 1.0 98.5694 0.0000
C(direction) 30522.8518 1.0 97.2540 0.0000
C(state):C(direction) 1007.2373 1.0 3.2093 0.0737
Residual 181089.5280 577.0 NaN NaN
amplitude:
sum_sq df F PR(>F)
C(state) 25.6013 1.0 24.1037 0.0000
C(direction) 132.3503 1.0 124.6082 0.0000
C(state):C(direction) 0.1651 1.0 0.1554 0.6936
Residual 612.8498 577.0 NaN NaN
reading velocity_change_% amplitude_change_% velocity_in_tol amplitude_in_tol
a both directions 28.11 14.56 False False
b flexion 29.70 11.39 False False
extension, for reference 26.00 18.86 False False
Variant, sessions with the middle command size only (43 pre, 12 post): peak velocity 42.4 and 83.9 deg/s (98.2%), amplitude 2.19 and 4.58 deg (109.0%)
peak_velocity: MPTP F = 41.28, p = 0.0000
amplitude: MPTP F = 44.31, p = 0.0000
By command size (session means, both directions):
peak_velocity amplitude latency_ms
mean count mean count mean count
level state
0 pre 95.7 6 5.1 6 7.8 6
1 pre 54.4 210 3.1 210 11.6 210
2 post 83.9 24 4.6 24 7.6 24
pre 42.4 85 2.2 85 7.7 85
4 post 65.1 250 3.2 250 4.1 250
pre 66.4 6 3.8 6 3.7 6
All session means are plotted (ranges printed above). Absolute scale: torque-induced amplitude, mean of session means 3.10 deg (paper: about 10 deg); peak velocity 59.0 deg/s, highest session mean 145.3 (paper's matched sessions 130 to 200 deg/s). Active movements in the same files, same filter: median peak velocity 100.3 deg/s and amplitude 20.4 deg.
# Variant: pre against post at the same per-trial command step (the session level hides trial-to-trial
# variation after MPTP). Single perturbations, by direction, levels with at least 20 trials in each state.
print("Perturbations by per-trial command level and state:")
print(pd.crosstab(kin.trial_level, kin.state).reindex(columns=["pre", "post"]).fillna(0).astype(int).to_string())
rows = []
for (lv, d), g in kin.groupby(["trial_level", "direction"]):
a, b = g[g.state == "pre"], g[g.state == "post"]
if len(a) < 20 or len(b) < 20:
continue
for m in ("peak_velocity", "amplitude"):
rows.append({"level": lv, "direction": d, "measure": m, "n_pre": len(a), "n_post": len(b),
"sessions_pre": a.session.nunique(), "sessions_post": b.session.nunique(),
"median_pre": a[m].median(), "median_post": b[m].median(),
"change_%": 100 * (b[m].mean() / a[m].mean() - 1), "mwu_p": mwu(a[m], b[m])})
print("\nPre against post at the same per-trial level:")
print(pd.DataFrame(rows).round(4).to_string(index=False))
Perturbations by per-trial command level and state: state pre post trial_level -0.0 51 0 1.0 1372 50 2.0 483 422 3.0 0 459 4.0 41 2002 5.0 0 425 6.0 0 197 7.0 0 45 8.0 0 3 Pre against post at the same per-trial level: level direction measure n_pre n_post sessions_pre sessions_post median_pre median_post change_% mwu_p 1.0 extension peak_velocity 695 36 106 24 43.5740 13.7173 -42.6296 0.0000 1.0 extension amplitude 695 36 106 24 2.2571 0.5205 -52.7935 0.0000 2.0 extension peak_velocity 246 271 42 87 27.7904 29.0602 24.5710 0.3682 2.0 extension amplitude 246 271 42 87 1.3725 1.2495 27.4090 0.1966 2.0 flexion peak_velocity 237 151 41 28 31.4938 112.4136 82.3357 0.0000 2.0 flexion amplitude 237 151 41 28 1.6899 5.6993 84.2642 0.0000 4.0 extension peak_velocity 24 999 4 124 59.6048 59.3326 -2.0402 0.4631 4.0 extension amplitude 24 999 4 124 3.2160 2.9499 -11.2136 0.0008
Result: not reproduced. In monkey V, torque-induced movements are faster and larger after MPTP: session-mean peak velocity rises by 28.1% and amplitude by 14.6% over both directions (29.7% and 11.4% for flexion), with MPTP effects F = 98.6 and 24.1 (p < 0.0001). Restricted to the sessions with the same torque command size (variant), the increase is larger (98.2% and 109.0%). The mean of the session-mean peak velocities is 59.0 deg/s, and only one response comes from a session in the paper's 130 to 200 deg/s window of matched sessions. The angle's scale matches the files' own kinematic columns (data check), so the small torque-induced movements (mean amplitude 3.10 deg, where the paper gives a mean of 10 deg) are in the data and do not come from the measure here. Compared trial by trial at the same command step (variant), the result is mixed: at twice the smallest step, flexion movements are faster and larger after MPTP (median peak velocity 31.5 against 112.4 deg/s, p < 0.0001), extension movements do not differ (p = 0.37 and 0.20), at the smallest step extension movements are slower and smaller after MPTP (43.6 against 13.7 deg/s, 36 post-MPTP trials), and at four times the smallest step extension peak velocity does not differ (p = 0.46, 24 pre-MPTP trials from 4 sessions). The torque command differs between the states, so the change cannot be attributed to MPTP.
kinematic-latency-unchanged¶
The latency of torque-induced displacements did not change after MPTP (F < 0.2, p > 0.05).
print(anova(sess_kin, "latency_ms").round(4).to_string())
print(sess_kin.groupby(["direction", "state"]).latency_ms.mean().round(2).to_string())
a = anova(lvl2, "latency_ms")
print(f"Variant, middle command size only: MPTP F = {a.loc['C(state)', 'F']:.2f}, p = {a.loc['C(state)', 'PR(>F)']:.4f}; "
f"means {lvl2.groupby('state').latency_ms.mean().round(2).to_dict()}")
print(f"\nPerturbations with a latency of 0 ms: {int((kin.latency_ms == 0).sum())} of {len(kin)}")
still = kin[~kin.moving_before]
for label, df, col in (("not moving before onset", still, "latency_ms"),
("causal 10 ms average of the unfiltered velocity", kin, "latency_causal_ms"),
("causal, not moving before onset", still, "latency_causal_ms")):
sk = df.groupby(["session", "state", "direction"])[col].mean().reset_index()
sk["level"] = sk.session.map(command.set_index("session").level)
a = anova(sk.rename(columns={col: "lat"}), "lat")
a2 = anova(sk[sk.level == 2].rename(columns={col: "lat"}), "lat")
print(f"Variant, {label}: means {sk.groupby(['direction', 'state'])[col].mean().round(2).to_dict()}; "
f"MPTP F = {a.loc['C(state)', 'F']:.1f}, p = {a.loc['C(state)', 'PR(>F)']:.2g}; "
f"middle command size only F = {a2.loc['C(state)', 'F']:.2f}, p = {a2.loc['C(state)', 'PR(>F)']:.2g}")
sum_sq df F PR(>F)
C(state) 4909.0313 1.0 1036.6564 0.0
C(direction) 549.9918 1.0 116.1436 0.0
C(state):C(direction) 87.5602 1.0 18.4904 0.0
Residual 2732.3528 577.0 NaN NaN
direction state
extension post 5.82
pre 10.86
flexion post 3.05
pre 9.65
Variant, middle command size only: MPTP F = 0.05, p = 0.8292; means {'post': 7.62, 'pre': 7.68}
Perturbations with a latency of 0 ms: 241 of 5550
Variant, not moving before onset: means {('extension', 'post'): 5.57, ('extension', 'pre'): 10.86, ('flexion', 'post'): 3.15, ('flexion', 'pre'): 9.65}; MPTP F = 1042.7, p = 2e-131; middle command size only F = 0.00, p = 0.99
Variant, causal 10 ms average of the unfiltered velocity: means {('extension', 'post'): 15.31, ('extension', 'pre'): 18.47, ('flexion', 'post'): 12.48, ('flexion', 'pre'): 17.48}; MPTP F = 791.4, p = 2.7e-110; middle command size only F = 1.35, p = 0.25
Variant, causal, not moving before onset: means {('extension', 'post'): 15.36, ('extension', 'pre'): 18.48, ('flexion', 'post'): 13.27, ('flexion', 'pre'): 17.48}; MPTP F = 761.7, p = 1.6e-107; middle command size only F = 5.20, p = 0.025
# The zero-phase filter on synthetic movements: the angle is flat, then moves from a known time at
# a constant velocity (a ramp), or follows a 50 ms smooth pulse. The latency each measure reports
# against the true start shows how far the forwards-backwards FIR moves an onset earlier.
ax_ = np.arange(stretch.ANALOG_WINDOW_MS[0], stretch.ANALOG_WINDOW_MS[1] + 1)
rows = []
for shape in ("ramp", "smooth pulse"):
for start in (5, 10, 20, 40):
for vel in (20, 60, 200): # deg/s
t = np.clip(ax_ - start, 0, None) / 1000
if shape == "ramp":
ang = vel * t
else: # velocity rises and falls over 50 ms (raised cosine), then holds
v = np.where((t > 0) & (t < 0.05), vel * (1 - np.cos(2 * np.pi * t / 0.05)) / 2, 0.0)
ang = np.cumsum(v) / 1000
r = stretch.kinematics(ang, "flexion", 1)
rows.append({"shape": shape, "true_start_ms": start, "velocity": vel, "latency_filtered": r["latency_ms"],
"latency_causal": r["latency_causal_ms"], "moving_before": r["moving_before"]})
syn = pd.DataFrame(rows)
print(syn.to_string(index=False))
print(f"\nFiltered latency minus true start: {(syn.latency_filtered - syn.true_start_ms).min():.0f} to "
f"{(syn.latency_filtered - syn.true_start_ms).max():.0f} ms; causal: {(syn.latency_causal - syn.true_start_ms).min():.0f} "
f"to {(syn.latency_causal - syn.true_start_ms).max():.0f} ms; flagged moving before onset: {int(syn.moving_before.sum())} of {len(syn)}")
# Latency at the same per-trial command step
rows = []
for (lv, d), g in kin.groupby(["trial_level", "direction"]):
a, b = g[g.state == "pre"], g[g.state == "post"]
if len(a) < 20 or len(b) < 20:
continue
for col in ("latency_ms", "latency_causal_ms"):
x, y = a[col].dropna(), b[col].dropna()
rows.append({"level": lv, "direction": d, "measure": col, "n_pre": len(x), "n_post": len(y),
"pre": x.mean(), "post": y.mean(), "mwu_p": mwu(x, y)})
print("\nPerturbations at the same per-trial command step (levels with at least 20 trials in each state):")
print(pd.DataFrame(rows).round(4).to_string(index=False))
shape true_start_ms velocity latency_filtered latency_causal moving_before
ramp 5 20 0.0 8.0 True
ramp 5 60 0.0 6.0 True
ramp 5 200 0.0 6.0 True
ramp 10 20 4.0 13.0 False
ramp 10 60 0.0 11.0 True
ramp 10 200 0.0 11.0 True
ramp 20 20 14.0 23.0 False
ramp 20 60 9.0 21.0 False
ramp 20 200 6.0 21.0 True
ramp 40 20 34.0 43.0 False
ramp 40 60 29.0 41.0 False
ramp 40 200 26.0 41.0 False
smooth pulse 5 20 10.0 18.0 False
smooth pulse 5 60 4.0 14.0 False
smooth pulse 5 200 1.0 11.0 True
smooth pulse 10 20 15.0 23.0 False
smooth pulse 10 60 9.0 19.0 False
smooth pulse 10 200 6.0 16.0 True
smooth pulse 20 20 25.0 33.0 False
smooth pulse 20 60 19.0 29.0 False
smooth pulse 20 200 16.0 26.0 True
smooth pulse 40 20 45.0 53.0 False
smooth pulse 40 60 39.0 49.0 False
smooth pulse 40 200 36.0 46.0 False
Filtered latency minus true start: -14 to 5 ms; causal: 1 to 13 ms; flagged moving before onset: 9 of 24
Perturbations at the same per-trial command step (levels with at least 20 trials in each state): level direction measure n_pre n_post pre post mwu_p 1.0 extension latency_ms 695 33 12.0561 16.3636 0.0030 1.0 extension latency_causal_ms 695 34 19.5108 21.3529 0.6778 2.0 extension latency_ms 246 271 7.9919 9.0000 0.0023 2.0 extension latency_causal_ms 246 271 16.0488 16.7454 0.0000 2.0 flexion latency_ms 237 151 7.1730 7.1656 0.9610 2.0 flexion latency_causal_ms 237 151 15.3165 15.7285 0.0001 4.0 extension latency_ms 24 999 3.7917 4.6857 0.0028 4.0 extension latency_causal_ms 24 999 14.5833 14.7968 0.0380
Result: not reproduced. The movement onset latency is shorter after MPTP (session means 9.65 and 3.05 ms for flexion, 10.86 and 5.82 ms for extension; MPTP F = 1037, p < 0.0001), with an MPTP by direction interaction (F = 18.5). In the sessions with the same command size (variant) it does not differ (7.68 and 7.62 ms, p = 0.83). The zero-phase filter can move an onset earlier, and 241 of 5550 perturbations have a latency of 0 ms. With a causal 10 ms average of the unfiltered velocity, which cannot move it earlier, the session means are 12.5 to 18.5 ms and still shorter after MPTP (F = 791); leaving out the trials already moving before onset gives F = 1043 (filtered) and 762 (causal). At the middle command size the causal variant gives p = 0.25, and p = 0.025 without the trials moving before onset. On synthetic movements the filtered latency is 14 ms early to 5 ms late against the true start, and the causal one 1 to 13 ms late; 9 of 24 synthetic movements that start after the onset are flagged as moving before it. Compared trial by trial at the same command step (variant), latencies are not shorter after MPTP: at twice the smallest step, flexion 7.17 and 7.17 ms (p = 0.96) and extension 7.99 and 9.00 ms (p = 0.0023), and the causal measure is longer after MPTP in both directions (p <= 0.0001, differences under 1 ms). The latency difference goes with the torque command, so it cannot be attributed to MPTP; the outcome is the criterion's comparison.
EMG¶
Session-average responses of the three triceps channels to flexion torques, which stretch the triceps, with the extension torques shown for the contrast.
emg_rows, traces = [], {}
for s in venus:
code_ = s["session_id"].split("++")[1]
if code_ in repeated_files:
continue
for ch in stretch.EMG_CHANNELS[:3]:
if ch not in s["windows"]:
continue
w = s["windows"][ch]
if np.isnan(w).any() or np.nanstd(w[:, AX < 0]) < 1:
continue
for d in ("flexion", "extension"):
r = stretch.emg_response(w[s["pert_direction"] == d], s["analog_step_ms"])
emg_rows.append({"session": code_, "post": s["post"], "channel": ch, "direction": d,
**{k: r[k] for k in ("m1", "m2", "onset_ms", "duration_ms")}})
traces[(code_, ch, d)] = r["trace"]
emg = pd.DataFrame(emg_rows)
emg["state"] = np.where(emg.post, "post", "pre")
emg["level"] = emg.session.map(command.set_index("code").level)
print("Channel-direction-sessions without perturbations of that direction (dropped):", int(emg.m2.isna().sum()))
emg = emg.dropna(subset=["m2"])
print("Channel-sessions used:")
print(pd.crosstab([emg.channel, emg.direction], emg.state).to_string())
print("\nMedian m2 (A/D units) by channel, direction and state:")
print(emg.groupby(["channel", "direction", "state"]).m2.median().unstack().round(2).to_string())
print("\nRanges: m1", round(emg.m1.min(), 1), "to", round(emg.m1.max(), 1), "; m2", round(emg.m2.min(), 1), "to",
round(emg.m2.max(), 1), "; responses found in", int(emg.onset_ms.notna().sum()), "of", len(emg))
fig, axes = plt.subplots(1, 3, figsize=(11, 3), sharex=True)
for ax, ch in zip(axes, stretch.EMG_CHANNELS[:3]):
for st, col in (("pre", "k"), ("post", "0.6")):
keys = [(r.session, ch, "flexion") for r in emg[(emg.channel == ch) & (emg.state == st) & (emg.direction == "flexion")].itertuples()]
if keys:
m = np.mean([traces[k] for k in keys], axis=0)
ax.plot(AX, m, color=col, label=f"{st} (n = {len(keys)})")
ax.set_xlim(-50, 200)
ax.axvspan(41, 80, color="0.9")
ax.set_title(f"{ch}, flexion torque")
ax.set_xlabel("ms from torque onset")
ax.legend(fontsize=7)
axes[0].set_ylabel("EMG less baseline (A/D units)")
plt.tight_layout()
plt.show()
print("Grand averages are plotted from -50 to 200 ms (shaded: m2 window); y range is set by the data.")
Channel-direction-sessions without perturbations of that direction (dropped): 2
Channel-sessions used:
state post pre
channel direction
EMGTricepsLateralis extension 36 36
flexion 36 36
EMGTricepsLongus1 extension 130 143
flexion 130 142
EMGTricepsLongus2 extension 116 128
flexion 116 127
Median m2 (A/D units) by channel, direction and state:
state post pre
channel direction
EMGTricepsLateralis extension 5.41 0.59
flexion 75.97 5.48
EMGTricepsLongus1 extension -2.50 -0.52
flexion 164.04 3.12
EMGTricepsLongus2 extension -3.04 -0.21
flexion 297.14 3.68
Ranges: m1 -76.7 to 271.7 ; m2 -96.8 to 507.1 ; responses found in 1039 of 1176
Grand averages are plotted from -50 to 200 ms (shaded: m2 window); y range is set by the data.
emg-m2-larger¶
The long-latency (m2, 41 to 80 ms) EMG response of the stretched triceps was larger after MPTP (Kolmogorov-Smirnov, p < 0.05).
rows = []
flex = emg[emg.direction == "flexion"]
for ch in stretch.EMG_CHANNELS[:3]:
x = flex[flex.channel == ch]
a, b = x[x.state == "pre"].m2, x[x.state == "post"].m2
rows.append({"channel": ch, "n_pre": len(a), "median_pre": a.median(), "mean_pre": a.mean(), "n_post": len(b),
"median_post": b.median(), "mean_post": b.mean(), "ks_p": stats.ks_2samp(a, b).pvalue,
"larger_and_p<0.05": (b.median() > a.median()) and stats.ks_2samp(a, b).pvalue < 0.05})
m2_table = pd.DataFrame(rows)
print(m2_table.round({c: 2 for c in ("median_pre", "mean_pre", "median_post", "mean_post")})
.to_string(index=False, formatters={"ks_p": "{:.2e}".format}))
print("\nVariant, middle command size only:")
for ch in stretch.EMG_CHANNELS[:3]:
x = flex[(flex.channel == ch) & (flex.level == 2)]
a, b = x[x.state == "pre"].m2, x[x.state == "post"].m2
if len(a) and len(b):
print(f" {ch}: n {len(a)} and {len(b)}, median m2 {a.median():.2f} and {b.median():.2f}, KS p = {stats.ks_2samp(a, b).pvalue:.4f}")
print("\nMedian m2 by command size and state (flexion, all triceps channels):")
print(flex.groupby(["level", "state"]).m2.agg(["median", "count"]).round(2).to_string())
channel n_pre median_pre mean_pre n_post median_post mean_post ks_p larger_and_p<0.05
EMGTricepsLongus1 142 3.12 14.94 130 164.04 156.94 1.09e-56 True
EMGTricepsLongus2 127 3.68 28.91 116 297.14 285.02 6.01e-50 True
EMGTricepsLateralis 36 5.48 9.79 36 75.97 76.07 4.65e-15 True
Variant, middle command size only:
EMGTricepsLongus1: n 40 and 6, median m2 47.56 and 1.33, KS p = 0.0132
EMGTricepsLongus2: n 36 and 4, median m2 49.48 and 123.76, KS p = 0.3200
EMGTricepsLateralis: n 16 and 1, median m2 19.65 and 0.52, KS p = 0.1176
Median m2 by command size and state (flexion, all triceps channels):
median count
level state
0 pre 3.68 9
1 pre 0.86 199
2 post 1.81 11
pre 41.85 92
4 post 194.38 271
pre 17.09 5
# Variant: pre against post at the same per-trial command step. Session averages are rebuilt from the
# flexion perturbations of one per-trial level only (at least 5 in the session), so pre and post
# traces come from the same command step.
lvl_rows = []
for s in venus:
code_ = s["session_id"].split("++")[1]
if code_ in repeated_files:
continue
tl = trial_level[code_]
for ch in stretch.EMG_CHANNELS[:3]:
if ch not in s["windows"]:
continue
w = s["windows"][ch]
if np.isnan(w).any() or np.nanstd(w[:, AX < 0]) < 1:
continue
for lv in np.unique(tl[np.isfinite(tl)]):
k = (s["pert_direction"] == "flexion") & (tl == lv)
if k.sum() < stretch.MIN_REPS:
continue
r = stretch.emg_response(w[k], s["analog_step_ms"])
lvl_rows.append({"session": code_, "state": "post" if s["post"] else "pre", "channel": ch,
"trial_level": int(lv), "n_trials": int(k.sum()),
**{q: r[q] for q in ("m1", "m2", "onset_ms", "duration_ms")}})
emg_lvl = pd.DataFrame(lvl_rows)
print("Channel-sessions with at least 5 flexion perturbations of one per-trial level, by level and state:")
print(pd.crosstab(emg_lvl.trial_level, emg_lvl.state).reindex(columns=["pre", "post"]).fillna(0).astype(int).to_string())
rows = []
for lv, g in emg_lvl.groupby("trial_level"):
a, b = g[g.state == "pre"], g[g.state == "post"]
if len(a) < 5 or len(b) < 5:
continue
for ch in stretch.EMG_CHANNELS[:3] + ("all triceps",):
x, y = (a, b) if ch == "all triceps" else (a[a.channel == ch], b[b.channel == ch])
if len(x) < 3 or len(y) < 3:
continue
xd, yd = x.dropna(subset=["duration_ms"]), y.dropna(subset=["duration_ms"])
rows.append({"level": lv, "channel": ch, "n_pre": len(x), "n_post": len(y),
"m2_pre": x.m2.median(), "m2_post": y.m2.median(), "m2_ks_p": stats.ks_2samp(x.m2, y.m2).pvalue,
"dur_pre": xd.duration_ms.mean(), "dur_post": yd.duration_ms.mean(),
"dur_ks_p": stats.ks_2samp(xd.duration_ms, yd.duration_ms).pvalue if len(xd) > 1 and len(yd) > 1 else np.nan,
"onset_pre": xd.onset_ms.mean(), "onset_post": yd.onset_ms.mean(),
"onset_ks_p": stats.ks_2samp(xd.onset_ms, yd.onset_ms).pvalue if len(xd) > 1 and len(yd) > 1 else np.nan})
print("\nPre against post at the same per-trial level (levels with at least 5 channel-sessions in each state;"
" m2 medians in A/D units, duration and onset means in ms):")
print(pd.DataFrame(rows).round(3).to_string(index=False))
Channel-sessions with at least 5 flexion perturbations of one per-trial level, by level and state:
state pre post
trial_level
0 9 0
1 156 2
2 66 11
3 0 21
4 4 205
5 0 31
6 0 12
Pre against post at the same per-trial level (levels with at least 5 channel-sessions in each state; m2 medians in A/D units, duration and onset means in ms):
level channel n_pre n_post m2_pre m2_post m2_ks_p dur_pre dur_post dur_ks_p onset_pre onset_post onset_ks_p
2 EMGTricepsLongus1 28 6 50.462 1.329 0.011 38.464 19.000 0.104 28.214 53.000 0.121
2 EMGTricepsLongus2 26 4 52.206 113.003 0.363 30.154 29.750 0.977 22.769 51.250 0.267
2 all triceps 66 11 44.919 1.810 0.004 33.455 22.182 0.434 28.455 48.636 0.087
Result: reproduced. For flexion torques, the m2 response is larger after MPTP in all three triceps channels (medians 3.1 against 164.0, 3.7 against 297.1 and 5.5 against 76.0 A/D units; Kolmogorov-Smirnov p < 1e-14 for each). The increase follows the torque command: in sessions of the middle command size (variant), m2 is not larger after MPTP on triceps longus 1 (median 47.6 against 1.3, p = 0.013, in the opposite direction). By command size and state the median m2 is 0.86 at the smallest size (pre only), 41.85 (pre) and 1.81 (post) at the middle, and 17.09 (pre, n = 5) and 194.4 (post) at the largest, so before MPTP it does not rise with command size (the largest size has only 5 channel-sessions), while after MPTP it rises from 1.81 at the middle size to 194.4 at the largest. Rebuilt from the perturbations of one per-trial command step (variant), at twice the smallest step the median m2 over the triceps channels is 44.9 before and 1.8 after MPTP (Kolmogorov-Smirnov p = 0.004, 66 and 11 channel-sessions), though on triceps longus 2 it is larger after MPTP (52.2 against 113.0, n = 26 and 4, p = 0.36); no other step has 5 channel-sessions in both states. As the torque command is about four times larger in most post-MPTP sessions, the increase cannot be attributed to MPTP; the outcome is the criterion's comparison.
emg-duration-longer and emg-onset-unchanged¶
The total duration of EMG responses to torque increased after MPTP, from 57 +/- 6 to 79 +/- 3 ms (Kolmogorov-Smirnov, p < 0.001). The latency of the earliest EMG responses did not change.
fr = flex.dropna(subset=["duration_ms"])
a, b = fr[fr.state == "pre"], fr[fr.state == "post"]
print(f"Responses found: pre {len(a)} of {int((flex.state == 'pre').sum())}, post {len(b)} of {int((flex.state == 'post').sum())} channel-sessions")
print(f"Duration: pre {a.duration_ms.mean():.1f} ms (SEM {a.duration_ms.sem():.1f}, median {a.duration_ms.median():.0f}), "
f"post {b.duration_ms.mean():.1f} ms (SEM {b.duration_ms.sem():.1f}, median {b.duration_ms.median():.0f}); "
f"KS p = {stats.ks_2samp(a.duration_ms, b.duration_ms).pvalue:.4g}; range {fr.duration_ms.min():.0f} to {fr.duration_ms.max():.0f} ms")
print(f"Onset: pre {a.onset_ms.mean():.1f} ms, post {b.onset_ms.mean():.1f} ms; KS p = {stats.ks_2samp(a.onset_ms, b.onset_ms).pvalue:.4g}; "
f"range {fr.onset_ms.min():.0f} to {fr.onset_ms.max():.0f} ms")
l2 = fr[fr.level == 2]
print(f"Variant, middle command size only: duration {l2[l2.state == 'pre'].duration_ms.mean():.1f} and "
f"{l2[l2.state == 'post'].duration_ms.mean():.1f} ms (n {int((l2.state == 'pre').sum())} and {int((l2.state == 'post').sum())}); "
f"onset {l2[l2.state == 'pre'].onset_ms.mean():.1f} and {l2[l2.state == 'post'].onset_ms.mean():.1f} ms; "
f"KS p duration {stats.ks_2samp(l2[l2.state == 'pre'].duration_ms, l2[l2.state == 'post'].duration_ms).pvalue:.4f}, "
f"onset {stats.ks_2samp(l2[l2.state == 'pre'].onset_ms, l2[l2.state == 'post'].onset_ms).pvalue:.4f}")
print(f"Durations of 10 ms or less: pre {(a.duration_ms <= 10).mean():.1%}, post {(b.duration_ms <= 10).mean():.1%}; "
f"median duration pre {a.duration_ms.median():.0f} ms, post {b.duration_ms.median():.0f} ms")
to_end = fr.onset_ms + fr.duration_ms > stretch.ANALOG_WINDOW_MS[1]
print(f"Main (2 SD): responses lasting to the end of the window: {int(to_end.sum())} of {len(fr)}")
print("\nBy channel:")
print(fr.groupby(["channel", "state"])[["duration_ms", "onset_ms"]].agg(["mean", "count"]).round(1).to_string())
print("\nThreshold variants (pooled channels):")
for sd in (1.5, 3.0):
rr = []
post_of = dict(zip(sessions.code, sessions.post))
for k, tr in traces.items():
if k[2] != "flexion":
continue
st = "post" if post_of[k[0]] else "pre"
b0 = (AX >= -200) & (AX < 0)
thr = sd * tr[b0].std(ddof=1)
above = (tr > thr) & (AX >= 10)
cand = np.flatnonzero(above & (AX <= 150))
if len(cand):
rest = np.flatnonzero(~above[cand[0]:])
rr.append({"state": st, "duration": int(rest[0]) if len(rest) else len(AX) - cand[0],
"to_end": not len(rest), "onset": AX[cand[0]]})
rr = pd.DataFrame(rr)
pa, pb = rr[rr.state == "pre"], rr[rr.state == "post"]
print(f" {sd} SD: duration {pa.duration.mean():.1f} and {pb.duration.mean():.1f} ms (KS p = {stats.ks_2samp(pa.duration, pb.duration).pvalue:.4g}); "
f"onset {pa.onset.mean():.1f} and {pb.onset.mean():.1f} ms (KS p = {stats.ks_2samp(pa.onset, pb.onset).pvalue:.4g}); "
f"responses lasting to the end of the window (500 ms) {int(rr.to_end.sum())}")
Responses found: pre 279 of 305, post 281 of 282 channel-sessions
Duration: pre 24.6 ms (SEM 3.8, median 9), post 49.4 ms (SEM 2.4, median 59); KS p = 4.119e-36; range 1 to 491 ms
Onset: pre 40.3 ms, post 24.0 ms; KS p = 9.114e-21; range 10 to 145 ms
Variant, middle command size only: duration 29.2 and 27.9 ms (n 92 and 10); onset 28.8 and 45.5 ms; KS p duration 0.9372, onset 0.1828
Durations of 10 ms or less: pre 63.4%, post 27.0%; median duration pre 9 ms, post 59 ms
Main (2 SD): responses lasting to the end of the window: 5 of 560
By channel:
duration_ms onset_ms
mean count mean count
channel state
EMGTricepsLateralis post 68.6 36 24.6 36
pre 15.5 33 45.0 33
EMGTricepsLongus1 post 46.8 130 23.4 130
pre 28.1 125 40.4 125
EMGTricepsLongus2 post 46.4 115 24.4 115
pre 23.5 121 38.9 121
Threshold variants (pooled channels):
1.5 SD: duration 20.8 and 50.1 ms (KS p = 3.615e-35); onset 32.8 and 22.2 ms (KS p = 2.556e-12); responses lasting to the end of the window (500 ms) 4
3.0 SD: duration 23.4 and 50.1 ms (KS p = 3.902e-42); onset 54.4 and 25.8 ms (KS p = 8.885e-38); responses lasting to the end of the window (500 ms) 3
Result: partially reproduced. EMG responses last longer after MPTP (24.6 against 49.4 ms, Kolmogorov-Smirnov p = 4.1e-36), but both means are outside tolerance (expected 57 and 79 ms); thresholds of 1.5 and 3 SD give the same direction (20.8 against 50.1 and 23.4 against 50.1 ms; responses that last to the end of the window are now kept in these variants, as in the main one, where 5 of 560 do). At the middle command size (variant) the durations do not differ (29.2 and 27.9 ms, p = 0.94). Before MPTP 63.4% of the durations are 10 ms or less (median 9 ms, against 59 ms after), so many pre-MPTP responses are brief crossings of a small signal. At the same per-trial command step (twice the smallest, variant) the durations do not differ (33.5 and 22.2 ms, p = 0.43).
Result: not reproduced. The onset of the EMG response is earlier after MPTP (40.3 against 24.0 ms, Kolmogorov-Smirnov p = 9.1e-21), with the same direction at thresholds of 1.5 and 3 SD. At the middle command size (variant) the onsets do not differ (28.8 and 45.5 ms, p = 0.18), the post-MPTP onset being the later one there; at the same per-trial command step it is also later after MPTP (28.5 and 48.6 ms, p = 0.087).
emg-m1-weak¶
The short-latency m1 component (15 to 40 ms) was poorly formed and generally much smaller than m2.
for st in ("pre", "post"):
x = flex[flex.state == st]
ratio = x.m1 / x.m2
print(f"{st}: m2 > m1 in {100 * (x.m2 > x.m1).mean():.1f}% of {len(x)} channel-sessions; median m1 {x.m1.median():.2f}, "
f"median m2 {x.m2.median():.2f}, median m1/m2 {ratio.median():.2f} (over channel-sessions with m2 > 0: "
f"{(x.m1 / x.m2)[x.m2 > 0].median():.2f})")
print(f" m2 <= 0 in {int((x.m2 <= 0).sum())} of {len(x)}; among channel-sessions with a response found "
f"({int(x.onset_ms.notna().sum())}), m2 > m1 in {100 * (x.m2 > x.m1)[x.onset_ms.notna()].mean():.1f}%")
pre: m2 > m1 in 74.4% of 305 channel-sessions; median m1 0.46, median m2 4.23, median m1/m2 0.10 (over channel-sessions with m2 > 0: 0.09)
m2 <= 0 in 66 of 305; among channel-sessions with a response found (279), m2 > m1 in 76.7%
post: m2 > m1 in 96.5% of 282 channel-sessions; median m1 41.54, median m2 191.13, median m1/m2 0.23 (over channel-sessions with m2 > 0: 0.23)
m2 <= 0 in 8 of 282; among channel-sessions with a response found (281), m2 > m1 in 96.8%
Result: partially reproduced. After MPTP m2 exceeds m1 in 96.5% of 282 channel-sessions with a median m1 / m2 of 0.23 (0.23 over channel-sessions with m2 > 0). Before MPTP the median ratio is 0.10 (0.09 with m2 > 0), but m2 exceeds m1 in 74.4% of 305, just below the 75% of the criterion. Before MPTP m2 is zero or negative in 66 of the 305; over the 279 with a response found, m2 exceeds m1 in 76.7%.
velocity-matched-magnitude¶
In sessions with equivalent torque-induced velocities (130 to 200 deg/s), MPTP did not change the magnitude of responses (Mann-Whitney, p = 0.19 and 0.46; 24 and 33 flexion, 22 and 31 extension responses).
sv = sess_kin.set_index(["session", "direction"]).peak_velocity
resp_v = resp.copy()
resp_v["session"] = [S[si]["session_id"] for si in resp_v.si]
resp_v["session_velocity"] = [sv.get((s, d), np.nan) for s, d in zip(resp_v.session, resp_v.direction)]
print("Responses without a session velocity (repeated files):", int(resp_v.session_velocity.isna().sum()))
matched = resp_v[(resp_v.session_velocity >= 130) & (resp_v.session_velocity <= 200)]
rows = []
for d in ("flexion", "extension"):
x = matched[matched.direction == d]
a, b = x[x.state == "pre"], x[x.state == "post"]
row = {"direction": d, "n_pre": len(a), "n_post": len(b),
"velocity_pre": a.session_velocity.mean(), "velocity_post": b.session_velocity.mean(),
"velocity_p": mwu(a.session_velocity, b.session_velocity)}
for m in ("max_change", "mean_change", "auc"):
row[f"{m}_pre"], row[f"{m}_post"], row[f"{m}_p"] = a[m].mean(), b[m].mean(), mwu(a[m], b[m])
rows.append(row)
print(pd.DataFrame(rows).round(4).T.to_string())
print("Session mean peak velocity range:", round(sess_kin.peak_velocity.min(), 1), "to", round(sess_kin.peak_velocity.max(), 1), "deg/s")
Responses without a session velocity (repeated files): 0
0 1
direction flexion extension
n_pre 0 0
n_post 1 0
velocity_pre NaN NaN
velocity_post 133.1665 NaN
velocity_p NaN NaN
max_change_pre NaN NaN
max_change_post 93.5684 NaN
max_change_p NaN NaN
mean_change_pre NaN NaN
mean_change_post 56.2528 NaN
mean_change_p NaN NaN
auc_pre NaN NaN
auc_post 10.8005 NaN
auc_p NaN NaN
Session mean peak velocity range: 8.5 to 145.3 deg/s
Result: blocked (methods underspecified), as the selection of responses is (see false-positive-rate). It is also blocked by the velocity scale: only one
response (post-MPTP) comes from a session whose mean peak velocity is between 130 and 200 deg/s
(session means reach 145.3 deg/s here), against 22 to 33 per group in the paper.
session-correlations¶
Across sessions, the mean PTN latency did not correlate with peak velocity (|rho| < 0.18) or with the m2 EMG response (|rho| < 0.26; Spearman, threshold 0.05/4).
ptn_cells = cells[cells.torque_related & (cells.ptype == "PTN")].copy()
ptn_cells["code"] = [S[si]["session_id"].split("++")[1] for si in ptn_cells.si]
ptn_cells["session"] = [S[si]["session_id"] for si in ptn_cells.si]
per_session = ptn_cells.groupby(["session", "code"]).latency.mean().reset_index()
vel = sess_kin.groupby("session").peak_velocity.mean()
m2 = flex[flex.channel == "EMGTricepsLongus1"].set_index("session").m2
per_session["velocity"] = per_session.session.map(vel)
per_session["m2"] = per_session.code.map(m2)
print(f"Sessions with a torque-related PTN: {len(per_session)}; with velocity {int(per_session.velocity.notna().sum())}; "
f"with triceps longus 1 m2 {int(per_session.m2.notna().sum())}")
for col in ("velocity", "m2"):
x = per_session.dropna(subset=[col])
rho, p = stats.spearmanr(x.latency, x[col])
print(f"PTN latency against {col}: Spearman rho = {rho:.3f}, p = {p:.4f} (n = {len(x)})")
for st, sub in x.groupby(x.session.str.contains("PostMPTP")):
r2, p2 = stats.spearmanr(sub.latency, sub[col])
print(f" {'post' if st else 'pre'} only: rho = {r2:.3f}, p = {p2:.4f} (n = {len(sub)})")
Sessions with a torque-related PTN: 70; with velocity 70; with triceps longus 1 m2 67 PTN latency against velocity: Spearman rho = -0.157, p = 0.1936 (n = 70) pre only: rho = 0.099, p = 0.5723 (n = 35) post only: rho = 0.028, p = 0.8724 (n = 35) PTN latency against m2: Spearman rho = -0.277, p = 0.0234 (n = 67) pre only: rho = 0.131, p = 0.4601 (n = 34) post only: rho = 0.163, p = 0.3650 (n = 33)
Result: blocked (methods underspecified), as the selection of responses is (see false-positive-rate). Over 70 sessions, mean PTN latency does not reach the paper's
threshold with velocity (rho = -0.16, p = 0.19) or with the triceps longus 1 m2 (rho = -0.28,
p = 0.023, n = 67, above 0.0125); within each state neither correlates (p 0.37 to 0.87).
trial-correlations¶
Trial by trial, PTN response latencies did not correlate with the peak velocity or the m2 EMG response (Spearman |rho| < 0.34, p > 0.05).
Blocked: the Methods define response onset on the trial-averaged SDF and give no way of measuring a latency on a single trial, which this test needs.
Summary¶
A failure to reproduce here may be a failure of this attempt, not of the paper: see the reasoning and the differences from the paper's method.
| Claim | Outcome |
|---|---|
cell-counts |
reproduced |
monkey-l-perturbations |
blocked: data missing |
baseline-rate |
reproduced |
ptn-faster-than-csn |
reproduced |
antidromic-latency |
partially reproduced |
false-positive-rate |
blocked: methods underspecified |
torque-related-fractions |
blocked: methods underspecified |
prevalence-unchanged |
blocked: methods underspecified |
ptn-more-responsive |
blocked: methods underspecified |
one-direction-half |
blocked: methods underspecified |
excitatory-majority |
blocked: methods underspecified |
magnitude-unchanged |
blocked: methods underspecified |
venus-response-magnitudes |
blocked: methods underspecified |
latency-shorter |
blocked: methods underspecified |
latency-not-training |
blocked: methods underspecified |
fwhm-shorter-ptn |
blocked: methods underspecified |
dsi-reduced |
blocked: methods underspecified |
venus-dsi |
blocked: methods underspecified |
dsi-not-baseline |
blocked: methods underspecified |
kinematics-slower |
not reproduced |
kinematic-latency-unchanged |
not reproduced |
emg-m2-larger |
reproduced |
emg-duration-longer |
partially reproduced |
emg-onset-unchanged |
not reproduced |
emg-m1-weak |
partially reproduced |
velocity-matched-magnitude |
blocked: methods underspecified |
session-correlations |
blocked: methods underspecified |
trial-correlations |
blocked: methods underspecified |
How the analysis differs from the paper's. Only monkey V is analysed, as the monkey L files mark no torque perturbation; the paper pools both monkeys except in its supplementary table. The response detection is one reading of the Methods (a threshold on the trial-averaged SDF set from the baseline's SD), with a t-test across trials as a variant. Its false-positive rate over the 200 ms epoch is 0.23, not the 0.01 the paper's threshold implies, and stationary Poisson spike trains give 0.16 to 0.21, most of it, so the claims that rest on the selection of responses are blocked, though their analyses ran. A threshold calibrated to a false-positive rate of 0.01 finds fewer torque-related cells, further from the paper's. With this detection the latency shortening after MPTP appears with the paper's monkey V means, but the response magnitudes fall after MPTP, which the paper did not find. The kinematic and EMG measures follow the Methods where they are given; the response windows, the EMG threshold and the peak velocity window are choices made here. The torque command differs in size between the states, so the kinematic and EMG differences found here cannot be attributed to MPTP, whatever their outcome on the criteria; compared at the same per-trial command step, movements are faster after MPTP for some steps and directions and slower for others. The m2 response can be compared at one step only (twice the smallest, 66 against 11 channel-sessions): its median over the triceps channels is smaller after MPTP (44.9 against 1.8), though on triceps longus 2 it is larger (52.2 against 113.0, 4 post-MPTP channel-sessions, p = 0.36), so a smaller m2 after MPTP rests on that one step and 11 post-MPTP channel-sessions. In synthetic tests the zero-phase filter, which the claims file fixes, moves a movement onset up to 14 ms earlier.
Changes after the first verification. The header now marks false-positive-rate and the 15
claims blocked by it as run and blocked. Added: the antidromic tests against their thresholds, a
per-cell table of the detection, the position of the spike windows within the trials, the movement
and EMG measures compared at the same per-trial command step, and the zero-phase filter on
synthetic movements. The EMG threshold variants now keep responses lasting to the end of the window.
No outcome changed.
Changes after the second verification. Text only: the m2 comparison at the same command step, the m2 by command size within each state, the false-positive rate against the paper's (about 23 and 46 times), and the cell-directions that can never respond are stated as the outputs show them. No outcome changed and none is proposed.
What changed from the earlier attempt (version 0.260518.1819). Every outcome is the same. The cell counts are no longer exactly the paper's (115 and 157 cells against 112 and 159), because the draft adds trials of second source files, so three pre-MPTP units now reach 5 perturbations per direction, and stores the two unit copies once; counting each unit from its first file gives the paper's pre-MPTP counts exactly. The session correlation of PTN latency with m2 no longer reaches the paper's threshold (p = 0.023, where it was 0.011). The other numbers moved by small amounts.
What was found about the dandiset, against the earlier attempt's findings.
- Monkey L perturbations: still present. The monkey L files mark no torque perturbation (type
not_available), the Leu torque command is NaN in every pre-MPTP file, and in the post-MPTP files it moves during the active movements, with one excursion in a hold window in 2048 trials. The draft's description now says the Leu files hold no perturbation times. - Analog clock: fixed. The analog streams share the trials' clock; the torque command steps 4 to 6 ms after every labelled onset on the stored clock.
- Torque command sizes: partly fixed. The pulse is still about four times larger in most post-MPTP sessions than in most pre-MPTP sessions, with an opposite swing after MPTP and trial-to-trial variation; the torque description now says so, and that its scale is not known. The kinematic and EMG results differ between sessions of different command sizes.
- Angle reference frame and series names: fixed. Positive is flexion, as the values show, and the monkey L series are named for the wrist.
- Repeated files and unit copies: fixed. No file repeats another's trials, and each unit is stored once.
unit_also_in_session_id: fixed. The column is gone; units recorded across two files are one unit observed in both files' trials.- Dandiset description: fixed. It now says intracarotid MPTP, wrist movements for monkey L, and Leu recordings from August 1992 to December 1993, as the files do.
- Sampling and EMG description: partly fixed. The EMG description now gives the paper's processing; the monkey L analog streams are still at 100 Hz (the paper says 200 Hz) and the EMG in A/D units, both now stated in the description.
- Torque timing: still present, now stated in the description. The torque onset comes 1.70 to 1.73 s after the center target appears in 90% of perturbations, where the Methods describe an onset 1 to 2 s after center capture, at random.
- New: three sessions (v1401, v1402, v1501) have labelled perturbations but no pulse in the torque command, as the torque description says.
- New: the code that built the draft is in an open pull request of the conversion repository, not on its main branch, which the dandiset links.
- The angle's scale agrees with the files' own movement amplitudes and peak velocities; the torque-induced movements (about 3 deg) are smaller than the paper's mean of 10 deg.
print(f"Notebook run time: {(time.time() - started) / 60:.1f} min")
Notebook run time: 1.9 min