Reproducing Amjad et al. (2025) from dandiset 002013¶
Paper. Amjad U, Mahajan S, Choi J, Shrivastav R, Murray R, Somich A, Coyne O, Schwerdt HN (2025). Microinvasive probes for monitoring electrical and chemical neural activity in nonhuman primates. ACS Chemical Neuroscience 16:2237-2247. doi:10.1021/acschemneuro.5c00071. This attempt used the full text, which is open access (CC-BY 4.0, PMC12183687).
Dandiset. DANDI:002013, version 0.261012.0921. The paper cites no DANDI version (it was published before the dandiset; its data statement points to Zenodo, doi:10.5281/zenodo.14826579). 0.261012.0921 is the only published version, so it is pinned. The draft lists the same three assets with the same sizes.
What the dandiset holds. Three sessions: Monkey T on 09/13/2024 (shape task, raw and trial-aligned FSCV from site c8dg), Monkey T on 09/26/2024 (direction task, raw ephys from c3bs and c3a and their bipolar LFP), and Monkey P on 08/13/2018 (raw ephys and sorted spikes from one site, with no trials or events). These are the example sessions of Figures 3 to 5. The second FSCV channel of Figure 3A (c8ds) is not in it.
Code read. The paper links two code deposits: extractFSCVspikes (GitHub schwerdtlab/extractFSCVspikes, tag v0.2.0-alpha, commit ea215c3; FSCV artifact removal for spike recording, used upstream of the data here) and Zenodo doi:10.5281/zenodo.10397773 (ISI_Plot.m, spikes_PCA.m, percent_recovery.mlx; these plot Figure 5B and 5C and take Plexon unit 1 as the sorted unit). Neither holds the dopamine PCA or the beta-power code. The dandiset links its conversion code (GitHub catalystneuro/schwerdt-lab-to-nwb, commit 6edf729), which shows how the trials, events, trial-aligned FSCV and LFP tables were built.
Outside input. One claim (zenodo-link) compares the dandiset with the paper's own Zenodo files:
four lists of trial numbers and the beta power of one example trial. The data check also reads Monkey P's
Neuralynx events file from the same deposit (session126_events.mat), for the time of its first event,
which places the spike times on the raw signal's clock, and to show its event codes (fig5a-spike-reward).
No other claim uses outside data.
| Claim | Tested here |
|---|---|
insertion-yield: 21.4% and 33% of sensors functional |
Blocked: data missing (no list of implanted sensors) |
chronic-duration: 483, 496 and 390 days post-implant |
Blocked: metadata insufficient (no implant dates) |
fscv-waveform: -0.4 to 1.3 V, 400 V/s, 10 Hz |
Tested |
ephys-acquisition: 30 and 32 kHz, 0.1 to 7500 Hz |
Tested |
shape-task-structure: 25% forced trials, blocks of 50 to 300 trials |
Tested |
zenodo-link: dandiset agrees with the paper's Zenodo files |
Tested, with outside data |
fig3a-dual-channel: dopamine from c8dg and c8ds |
Blocked: data missing (no c8ds) |
fig3a-redox-peaks: redox current near 0.6 and -0.2 V after the cue |
Blocked: metadata insufficient (raw FSCV not on the trial clock) |
fig3b-beta-bursts: beta bursts change around the cue |
Tested |
fig4a-dopamine-reward: more dopamine after large-reward cues |
Tested |
fig4b-beta-reward: more beta suppression after large-reward cues |
Tested |
fig4c-beta-side: more beta suppression for contralateral cues |
Tested |
fig5a-spike-reward: firing differs with reward size |
Blocked: data missing (no events for Monkey P) |
fig5bc-unit-isolation: a well-isolated sorted unit |
Tested |
import json
import re
import time
import warnings
import h5py
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from scipy import signal, stats
import microprobes as mp
warnings.filterwarnings("ignore", category=RuntimeWarning)
pd.set_option("display.width", 200)
pd.set_option("display.max_columns", 30)
LARGE, SMALL = "#eb6834", "#2a78d6" # categorical slots 2 and 1 of the reference palette
started = time.time()
meta = mp.dandiset_metadata()
files = {key: mp.open_file(key) for key in mp.PATHS}
print(f"dandiset {mp.DANDISET_ID} version {mp.DANDISET_VERSION}: {meta['assetsSummary']['numberOfFiles']} files, "
f"{meta['assetsSummary']['numberOfSubjects']} subjects, {meta['assetsSummary']['numberOfBytes'] / 1e9:.2f} GB")
dandiset 002013 version 0.261012.0921: 3 files, 2 subjects, 4.69 GB
Data check¶
Each check of the metadata against the values is stated under its output.
rows = []
for key, f in files.items():
el = f["general/extracellular_ephys/electrodes"]
rows.append({"file": mp.PATHS[key], "subject": mp.text(f, "general/subject/subject_id"),
"sex": mp.text(f, "general/subject/sex"), "age": mp.text(f, "general/subject/age"),
"session_start": mp.text(f, "session_start_time")[:19],
"sites": ", ".join(x.decode() for x in el["location"][:]),
"related_publication": ", ".join(x.decode() for x in f["general/related_publications"][:]),
"trials": f["intervals/trials/id"].shape[0] if "intervals/trials" in f else 0,
"events": "processing/events" in f})
info = pd.DataFrame(rows)
print(info.to_string(index=False))
print("\nContents per file:")
for key, f in files.items():
print(f" {key}: acquisition {sorted(f['acquisition'])}, processing "
f"{ {m: sorted(f['processing'][m]) for m in f['processing']} }, stimulus {sorted(f['stimulus/presentation'])}, "
f"units {'units' in f}")
print("\nDandiset relatedResource:")
for r in meta["relatedResource"]:
print(" ", r["relation"], r.get("url"))
print("\nPaper: Figure 3A sites c8dg and c8ds, 483 days post-implant; Figure 3B and 4B,C site c3bs-c3a, 496 days; "
"Figure 5 site c33, 390 days.")
print("Session dates differ by", (pd.Timestamp("2024-09-26") - pd.Timestamp("2024-09-13")).days,
"days; the paper's days post-implant differ by", 496 - 483)
for key, f in files.items():
print(f"{key}: subject fields { {k: mp.text(f, 'general/subject/' + k)[:60] for k in f['general/subject']} }")
print(f" devices { {k: f['general/devices'][k].attrs.get('description') for k in f['general/devices']} }; "
f"electrode columns {sorted(f['general/extracellular_ephys/electrodes'])}")
texts = {"dandiset metadata": json.dumps(meta, ensure_ascii=False)}
for key, f in files.items():
def visit(name, obj):
for a_, val in obj.attrs.items():
texts[f"{key}: {name} ({a_})"] = str(val)
if isinstance(obj, h5py.Dataset) and obj.dtype.kind in "OS" and obj.size < 1000:
texts[f"{key}: {name}"] = str(obj[()])
f["general"].visititems(visit)
date_like = r"\d{4}-\d{2}-\d{2}|\d{4}/\d{2}/\d{2}|\d{1,2}/\d{1,2}/\d{4}"
print("Every mention of 'implant' in the dandiset metadata and in the general section of the files, with context:")
for where, t in texts.items():
for m_ in re.finditer("(?i)implant", t):
ctx = t[max(0, m_.start() - 60):m_.end() + 60]
print(f" {where}: {ctx!r}; date in context: {bool(re.search(date_like, ctx))}")
print("\nFiles of the paper's Zenodo deposit:", mp.zenodo_files())
file subject sex age session_start sites related_publication trials events
sub-Monkey-T/sub-Monkey-T_ses-09132024_ecephys.nwb Monkey T F P11Y 2024-09-13T09:03:06 c8dg https://doi.org/10.1101/2025.01.30.635139 1468 True
sub-Monkey-T/sub-Monkey-T_ses-09262024_ecephys.nwb Monkey T F P11Y 2024-09-26T09:01:38 c3bs, c3a https://doi.org/10.1101/2025.01.30.635139 1501 True
sub-Monkey-P/sub-Monkey-P_ses-08132018_ecephys.nwb Monkey P F P8Y 2018-08-13T10:20:08 cl3 https://doi.org/10.1101/2025.01.30.635139 0 False
Contents per file:
T-0913: acquisition ['fscv_response_series'], processing {'events': ['events_table'], 'fscv': ['trial_aligned_fscv_table']}, stimulus ['fscv_excitation_series'], units False
T-0926: acquisition ['electrical_series'], processing {'ecephys': ['FilteredEphys'], 'events': ['events_table']}, stimulus [], units False
P-0813: acquisition ['electrical_series'], processing {'ecephys': ['thresholded_units']}, stimulus [], units True
Dandiset relatedResource:
dcite:IsDescribedBy https://doi.org/10.1021/acschemneuro.5c00071
dcite:IsDerivedFrom https://doi.org/10.5281/zenodo.14826579
dcite:IsCompiledBy https://github.com/catalystneuro/schwerdt-lab-to-nwb
Paper: Figure 3A sites c8dg and c8ds, 483 days post-implant; Figure 3B and 4B,C site c3bs-c3a, 496 days; Figure 5 site c33, 390 days.
Session dates differ by 13 days; the paper's days post-implant differ by 13
T-0913: subject fields {'age': 'P11Y', 'description': 'Two adult female rhesus monkeys were used for in vivo valida', 'sex': 'F', 'species': 'Macaca mulatta', 'subject_id': 'Monkey T'}
devices {'device_ecephys': 'Neuralynx recording'}; electrode columns ['group', 'group_name', 'id', 'location']
T-0926: subject fields {'age': 'P11Y', 'description': 'Two adult female rhesus monkeys were used for in vivo valida', 'sex': 'F', 'species': 'Macaca mulatta', 'subject_id': 'Monkey T'}
devices {'device_ecephys': 'Neuralynx recording'}; electrode columns ['ADMaxValue', 'DSPHighCutFilterEnabled', 'DSPLowCutFilterEnabled', 'DspDelayCompensation', 'DspFilterDelay_µs', 'DspHighCutFilterType', 'DspHighCutFrequency', 'DspHighCutNumTaps', 'DspLowCutFilterType', 'DspLowCutFrequency', 'DspLowCutNumTaps', 'FileType', 'InputRange', 'NumADChannels', 'RecordSize', 'channel_name', 'group', 'group_name', 'id', 'input_inverted', 'location', 'recording_closed', 'recording_opened']
P-0813: subject fields {'age': 'P8Y', 'description': 'Two adult female rhesus monkeys were used for in vivo valida', 'sex': 'F', 'species': 'Macaca mulatta', 'subject_id': 'Monkey P'}
devices {'device_ecephys': 'Neuralynx recording'}; electrode columns ['channel_name', 'group', 'group_name', 'id', 'location']
Every mention of 'implant' in the dandiset metadata and in the general section of the files, with context:
dandiset metadata: '-n92ldne8nv5b/v1", "https://www.protocols.io/view/baseplate-implantation-for-two-part-chamber-system-kqdg32b91v25/v1", "https:/'; date in context: False
dandiset metadata: 'www.protocols.io/view/craniotomy-and-second-and-third-phase-implantation-x54v92wd4l3e/v1", "https://www.protocols.io/view/chron'; date in context: False
dandiset metadata: 'ps://www.protocols.io/view/chronic-multi-electrode-assembly-implantation-with-bp2l62m95gqe/v1"], "schemaKey": "Dandiset", "iden'; date in context: False
T-0913: experiment_description: '-assisted micro-invasive probes (s-\\xc2\\xb5IPs) chronically implanted in the striatum of task-performing rhesus monkeys. The s-'; date in context: False
T-0913: experiment_description: 'ield potentials with micrometer spatial resolution. Chronic implants demonstrated stable, multi-modal neural activity measureme'; date in context: False
T-0913: experiment_description: "ents from micro-invasive devices over extended periods post-implantation.'"; date in context: False
T-0913: keywords: "[b'fast scan cyclic voltammetry' b'micro-invasive neural implant'\n b'dopamine' b'striatum']"; date in context: False
T-0913: subject/description: 'in vivo validation of the s-\\xc2\\xb5IP device, with chronic implants targeting the striatum for simultaneous FSCV and EPhys rec'; date in context: False
T-0926: experiment_description: '-assisted micro-invasive probes (s-\\xc2\\xb5IPs) chronically implanted in the striatum of task-performing rhesus monkeys. The s-'; date in context: False
T-0926: experiment_description: 'ield potentials with micrometer spatial resolution. Chronic implants demonstrated stable, multi-modal neural activity measureme'; date in context: False
T-0926: experiment_description: "ents from micro-invasive devices over extended periods post-implantation.'"; date in context: False
T-0926: extracellular_ephys/ElectrodeGroup (description): 'ectrodes from the EPhys system (Neuralynx, Digital Lynx SX) implanted in the striatum.'; date in context: False
T-0926: keywords: "[b'fast scan cyclic voltammetry' b'micro-invasive neural implant'\n b'dopamine' b'striatum']"; date in context: False
T-0926: subject/description: 'in vivo validation of the s-\\xc2\\xb5IP device, with chronic implants targeting the striatum for simultaneous FSCV and EPhys rec'; date in context: False
P-0813: experiment_description: '-assisted micro-invasive probes (s-\\xc2\\xb5IPs) chronically implanted in the striatum of task-performing rhesus monkeys. The s-'; date in context: False
P-0813: experiment_description: 'ield potentials with micrometer spatial resolution. Chronic implants demonstrated stable, multi-modal neural activity measureme'; date in context: False
P-0813: experiment_description: "ents from micro-invasive devices over extended periods post-implantation.'"; date in context: False
P-0813: extracellular_ephys/ElectrodeGroup (description): 'ectrodes from the EPhys system (Neuralynx, Digital Lynx SX) implanted in the striatum.'; date in context: False
P-0813: keywords: "[b'fast scan cyclic voltammetry' b'micro-invasive neural implant'\n b'dopamine' b'striatum']"; date in context: False
P-0813: subject/description: 'in vivo validation of the s-\\xc2\\xb5IP device, with chronic implants targeting the striatum for simultaneous FSCV and EPhys rec'; date in context: False
Files of the paper's Zenodo deposit: ['09262024_c3bs-c3a_data_trial35bigrewardexample.mat', '09262024_tr_nlx_c3bs-c3a.mat', '09262024_trialids_value object start_big rew and success_tr1-1501.mat', '09262024_trialids_value object start_left and success_tr1-1501.mat', '09262024_trialids_value object start_right and success_tr1-1501.mat', '09262024_trialids_value object start_small rew and success_tr1-1501.mat', '09262024_trlist.mat', 'FSCV_trlists_c8dg_09132024_firsthalf.mat', 'cl3_session126_csc12_100_spikes.mat', 'session126_events.mat']
Two adult female rhesus monkeys and three sessions, as the dandiset describes. Monkey P's site is named
cl3 in the file, c33 in the Figure 5 legend. The files' related_publications point to the bioRxiv
preprint (doi:10.1101/2025.01.30.635139), not the article. The subject, device and electrode fields hold no implantation date, and no mention of implantation in the
general section of the files or in the dandiset metadata comes with a date, so the days post-implant cannot be checked; the two Monkey T sessions are 13 days apart, which
matches the 483 and 496 days of the legends. Monkey P's file has no trials and no events. The paper's Zenodo deposit holds ten files (listed above): trial
lists, LFP, spikes and events of the example sessions, and processed FSCV of c8dg only (the one FSCV file is named for c8dg), and none of its file names points
to a list of implanted sensors. The conversion
code also reads eye tracking channels, which are not in the files.
ev13, evtab13 = mp.load_events(files["T-0913"])
ev26, evtab26 = mp.load_events(files["T-0926"])
for key, tab in (("T-0913", evtab13), ("T-0926", evtab26)):
f = files[key]
tr = f["intervals/trials"]
st, mid, sp = tr["start_time"][:], tr["midpoint_time"][:], tr["stop_time"][:]
ev = ev13 if key == "T-0913" else ev26
tags = {x.decode() for x in tr["tags"][:]}
print(f"{key}: {len(mid)} trials; distinct trial-start codes (128) {len(ev[128])}; trial midpoints that are a "
f"trial-start code {np.isin(np.round(mid, 6), ev[128]).sum()}; midpoint - start {np.unique(np.round(mid - st, 6))} s; "
f"stop - start {np.unique(np.round(sp - st, 6))} s; tags {tags}")
print(f" trial-to-trial interval {np.percentile(np.diff(mid), [0, 50, 100]).round(2)} s; stored event times "
f"{tab.n_stored.sum()}, distinct {tab.n_distinct.sum()} (ratio {tab.n_stored.sum() / tab.n_distinct.sum():.2f})")
print(tab.to_string(index=False))
T-0913: 1468 trials; distinct trial-start codes (128) 1468; trial midpoints that are a trial-start code 1468; midpoint - start [30.] s; stop - start [60.] s; tags {'initial cue display'}
trial-to-trial interval [ 6.53 13.56 165.7 ] s; stored event times 85665, distinct 17686 (ratio 4.84)
code label n_stored n_distinct
18 end trial 7099 1468
21 feedback 1318 303
23 value object start 1508 323
24 central cue end 4704 977
28 onedr 1745 358
30 left cue reward condition 1 3207 639
31 left cue reward condition 2 2340 471
35 right cue reward condition 1 3207 652
36 right cue reward condition 2 4085 816
50 central cue fixation started 6938 1404
51 value cue fixation started 4686 977
53 forced trial valuecue fix started 1495 323
60 error fixation break target 564 104
61 error fixation break initial 9 2
62 error fixation never started 324 62
65 error choice value initial fix break 5 1
66 error choice fixbreak right 946 195
67 error choice fixbreak left 708 147
69 error forced trial value cue fix break 97 20
81 image 1 6411 1285
88 image 8 6428 1293
101 big reward 2802 643
105 reward delivery end 4119 937
115 transient value condition 1745 358
117 left condition 807 166
118 right condition 938 192
119 forced trial condition 1745 358
120 choice trial condition 5547 1110
121 right choice chosen code 1245 277
122 left choice chosen code 1601 357
128 intended trial start 7292 1468
T-0926: 1501 trials; distinct trial-start codes (128) 1501; trial midpoints that are a trial-start code 1501; midpoint - start [30.] s; stop - start [60.] s; tags {'initial cue display'}
trial-to-trial interval [ 6.75 13.5 16.31] s; stored event times 59702, distinct 12434 (ratio 4.80)
code label n_stored n_distinct
18 end trial 7378 1501
21 feedback 5085 1172
23 value object start 6720 1372
50 central cue fixation started 7109 1411
51 value cue fixation started 6253 1303
60 error central cue fixation break 220 39
61 error central cue initial fixation break 40 7
62 error fixation never started 464 83
65 error value cue initial fixation break 336 69
69 error value cue fixation break 651 131
100 small reward 2639 602
101 big reward 2476 570
105 reward delivery end 5117 1172
117 left condition 3846 756
118 right condition 3761 745
128 intended trial start 7607 1501
The trials tables hold only times: each row is a 60 s window centred on a trial-start code (128), with no condition columns. The events tables store each event once for every trial window that contains it, so each time appears about five times (trials start every 13 s or so). The distinct times are used from here on: one trial-start code per trial in both sessions. The conditions below are decoded from the event codes.
d26 = mp.direction_task_trials(files["T-0926"])
s13 = mp.shape_task_trials(files["T-0913"])
print("Direction task, 09/26/2024:", len(d26), "trials;", int(d26.success.sum()), "successful (a reward code)")
print(pd.crosstab(d26.success, [d26.big, d26.left], rownames=["success"], colnames=["big", "left"]).to_string())
print("trials with both or neither side code:", int((d26.left == d26.right).sum()),
"; with both reward codes:", int((d26.big & d26.small).sum()))
ok = d26[d26.success]
print("value cue - trial start (s):", np.percentile(ok.value_cue - ok.start, [0, 50, 100]).round(3),
"; feedback - value cue (s):", np.percentile(ok.feedback - ok.value_cue, [0, 50, 100]).round(3))
side_big = np.where(ok.big == ok.left, "large on left", "large on right")
runs = np.diff(np.r_[0, np.where(side_big[1:] != side_big[:-1])[0] + 1, len(side_big)])
print("direction task: blocks of successful trials with the same side-reward mapping:", len(runs),
"; their lengths:", runs.tolist())
print("\nShape task, 09/13/2024:", len(s13), "trials; forced", int(s13.forced.sum()), "; choice", int(s13.choice.sum()),
"; both", int((s13.forced & s13.choice).sum()))
print("forced trials with feedback (successful):", int(s13.success.sum()), "; choice trials with feedback:",
int((s13.choice & ~np.isnan(s13.feedback)).sum()))
print("Successful forced trials by big-reward code (101) and cue condition codes 35 and 36:")
okf = s13[s13.success]
print(pd.crosstab(okf.big, [okf.cue_reward_1, okf.cue_reward_2], rownames=["code 101"], colnames=["code 35", "code 36"]).to_string())
print("value cue - trial start (s):", np.percentile(okf.value_cue - okf.start, [0, 50, 100]).round(3),
"; feedback - value cue (s):", np.percentile(okf.feedback - okf.value_cue, [0, 50, 100]).round(3))
Direction task, 09/26/2024: 1501 trials; 1172 successful (a reward code) big False True left False True False True success False 167 162 0 0 True 290 312 288 282 trials with both or neither side code: 0 ; with both reward codes: 0 value cue - trial start (s): [1.549 1.799 3.566] ; feedback - value cue (s): [3.715 3.799 4.582] direction task: blocks of successful trials with the same side-reward mapping: 16 ; their lengths: [77, 84, 53, 70, 88, 86, 88, 93, 65, 56, 65, 50, 61, 87, 75, 74] Shape task, 09/13/2024: 1468 trials; forced 358 ; choice 1110 ; both 0 forced trials with feedback (successful): 303 ; choice trials with feedback: 0 Successful forced trials by big-reward code (101) and cue condition codes 35 and 36: code 35 False True code 36 True False code 101 False 139 0 True 0 164 value cue - trial start (s): [1.549 1.85 3.566] ; feedback - value cue (s): [3.748 3.782 4.032]
Direction task: every successful trial has exactly one reward code and one side code, and large and small rewards occur on both sides, with the mapping switching in blocks, as the Methods describe. Shape task: forced and choice trials are exclusive, and only forced trials carry the feedback code (21), so success can be read only for forced trials, the only ones the paper shows. The session has no small-reward code (100). In successful forced trials the big-reward code (101) occurs exactly when the condition code 35 does, and never with 36. Two independent codes agree on the reward size, so a successful forced trial without 101 is taken as small. The codes' labels ("right cue reward condition 1/2") do not say this; the meaning is inferred from the data. The value cue fixation lasts 3.4 to 3.8 s (value cue to feedback) in both sessions.
fs = mp.fscv_table(files["T-0913"])
tab = files["T-0913"]["processing/fscv/trial_aligned_fscv_table"]
print("trial-aligned FSCV rows:", len(fs["trial"]), "; trial numbers", fs["trial"].min(), "to", fs["trial"].max(),
"; contiguous:", bool(np.all(np.diff(fs["trial"]) == 1)), "; good:", int(fs["good"].sum()),
"; rate:", np.unique(fs["rate"]), "; samples per row:", fs["da"].shape[1])
for name in ("da", "ph", "m", "iox"):
x = fs[name]
print(f" {name}: range {np.nanmin(x):.3g} to {np.nanmax(x):.3g}, NaN fraction {np.isnan(x).mean():.3f}; "
f"description: {tab[name].attrs['description']!r}; unit attribute: {tab[name].attrs.get('unit', 'none')}")
rows13 = fs["trial"] - 1
cue13 = (s13.value_cue.values - s13.start.values)[rows13]
sel = s13.success.values[rows13]
at_cue = fs["da"][rows13[sel] - rows13[0], np.round(300 + cue13[sel] * 10).astype(int)]
print("successful forced trials with a FSCV row:", int(sel.sum()), "; large", int(s13.big.values[rows13][sel].sum()),
"; small", int(s13.small.values[rows13][sel].sum()))
print(f"stored da at the value cue in these trials: median |da| {np.nanmedian(np.abs(at_cue)):.2f}, exactly 0 in "
f"{int(np.sum(at_cue == 0))}; median |da| over all their samples {np.nanmedian(np.abs(fs['da'][sel])):.2f}")
trial-aligned FSCV rows: 734 ; trial numbers 1 to 734 ; contiguous: True ; good: 734 ; rate: [10.] ; samples per row: 601 da: range -57.5 to 65.9, NaN fraction 0.027; description: 'PCA extracted dopamine concentration time series.'; unit attribute: none ph: range -3.29 to 2.31, NaN fraction 0.000; description: 'PH change time series.'; unit attribute: none m: range -0.785 to 0.926, NaN fraction 0.000; description: 'Motion artifact time series.'; unit attribute: none iox: range -1.47 to 1.93, NaN fraction 0.000; description: 'Measured oxidation current at 0.6 V (approximate voltage at which dopamine oxidizes).'; unit attribute: none successful forced trials with a FSCV row: 168 ; large 87 ; small 81 stored da at the value cue in these trials: median |da| 6.70, exactly 0 in 0; median |da| over all their samples 5.92
The trial-aligned FSCV table has rows only for the first 734 of the 1468 trials (the source file is named
FSCV_trlists_c8dg_09132024_firsthalf.mat in the conversion notes), all marked good. Of the paper's
forced trials, 168 successful ones (87 large, 81 small) have dopamine here. The columns carry no units.
The stored da is not zero at the value cue, so it is not referenced there as the Methods describe; the
analysis re-references it at the value cue (claims file).
raw = mp.raw_fscv(files["T-0913"])
v = raw["volts"]
rise = slice(0, int(np.argmax(v)) + 1)
dt = np.median(np.diff(raw["t"][:214]))
print("samples per scan:", raw["current"].shape[1], "; scans:", raw["current"].shape[0], "; sample interval",
round(dt * 1e6, 2), "µs")
print(f"applied voltage {v.min():.3f} to {v.max():.3f} V (first 50 scans {raw['volts_first'].min():.3f} to "
f"{raw['volts_first'].max():.3f}); stored attributes {raw['excitation_attrs']}")
print("applied voltage data attributes:", raw["excitation_data_attrs"])
print("current data attributes:", raw["current_data_attrs"])
print("scan start intervals (s):", np.unique(np.round(np.diff(raw["scan_start"]), 6)))
mid13 = s13.start.values
print(f"raw FSCV spans {raw['scan_start'][0]:.3f} to {raw['scan_start'][-1] + 0.1:.3f} s; trials (trial starts) span "
f"{mid13[0]:.3f} to {mid13[-1]:.3f} s; first scan - first trial start = {raw['scan_start'][0] - mid13[0]:.6f} s")
print(f"current: unit attribute {raw['unit']!r}, values {raw['current'].min():.1f} to {raw['current'].max():.1f}; "
f"other attributes {raw['attrs']}")
samples per scan: 214 ; scans: 241200 ; sample interval 40.0 µs
applied voltage -0.347 to 1.339 V (first 50 scans -0.347 to 1.340); stored attributes {'comments': 'no comments', 'description': 'FSCV excitation series for applied voltage waveform.', 'namespace': 'ndx-fscv', 'neurodata_type': 'FSCVExcitationSeries', 'scan_frequency': np.float64(10.0), 'sweep_rate': np.float64(400.0), 'waveform_shape': 'Triangle'}
applied voltage data attributes: {'conversion': np.float64(1.0), 'offset': np.float64(0.0), 'resolution': np.float64(-1.0), 'unit': 'volts'}
current data attributes: {'conversion': np.float64(1.0), 'offset': np.float64(0.0), 'resolution': np.float64(-1.0), 'unit': 'amperes'}
scan start intervals (s): [0.1]
raw FSCV spans 14120.168 to 38240.168 s; trials (trial starts) span 14120.168 to 32910.208 s; first scan - first trial start = 0.000000 s
current: unit attribute 'amperes', values -1526.5 to 1557.0; other attributes {'comments': 'no comments', 'current_to_voltage_factor': np.float64(0.00499), 'description': 'FSCV response current before background subtraction.', 'namespace': 'ndx-fscv', 'neurodata_type': 'FSCVResponseSeries'}
The applied waveform is a triangle from -0.347 to 1.340 V, as the dandiset's description says (-0.35 to
1.34 V); the Methods say -0.4 to 1.3 V (claim fscv-waveform). The data of both series have conversion 1 and offset 0, so the voltage is read as stored. The raw current is stored as amperes, but
its values are in the hundreds to thousands, so the unit attribute is wrong (they are likely nA). The scan
times are an exact 0.1 s grid that begins at the first trial start and runs 24120 s, past the last trial
(18790 s later). The conversion notes say they were interpolated between trial-start codes of the two
systems; a perfect grid starting exactly at the first trial start suggests they were not. The next cell
tests whether the raw current shares the trial clock.
k06 = int(np.argmin(np.abs(v[rise] - 0.6)))
c06 = raw["current"][:, k06]
hp = signal.butter(2, 0.5, btype="high", fs=10, output="sos")
c06_hp = signal.sosfiltfilt(hp, c06)
print(f"current at {v[k06]:.3f} V (sample {k06}) compared with the stored iox, both high-passed at 0.5 Hz; "
"search over the whole raw record:")
for r in [0, 50, 200, 400, 600, 733]:
y = signal.sosfiltfilt(hp, fs["iox"][r])
i0 = int(np.searchsorted(raw["scan_start"], fs["start"][r]))
k, rr = mp.best_lag(y, c06_hp)
r0 = np.corrcoef(c06_hp[i0:i0 + 601], y)[0, 1]
print(f" row {r} (trial {fs['trial'][r]}): r at the stored time {r0:.3f}; best r {rr:.3f}, raw {fs['start'][r] - raw['scan_start'][k]:+.1f} s from the stored time")
def raw_contrast(shift):
t = s13.value_cue.values[s13.success.values] + shift
i = np.searchsorted(raw["scan_start"], t)
ok = (i > 0) & (i + 31 < len(c06))
resp = np.array([c06[j + 5:j + 30].mean() - c06[j] for j in i[ok]])
big = s13.big.values[s13.success.values][ok]
return resp[big].mean() - resp[~big].mean(), stats.mannwhitneyu(resp[big], resp[~big]).pvalue, int(ok.sum())
print("\nRaw current at 0.6 V, mean over 0.5 to 3 s after the value cue minus at the cue, large minus small "
"(all successful forced trials), at the stored times and shifted:")
for shift in [0, -60, -30, -10, -5, 5, 10, 30, 60]:
d, p, n = raw_contrast(shift)
print(f" shift {shift:+4d} s: difference {d:+.3f}, Mann-Whitney p = {p:.3f}, n = {n}")
top = int(np.argmax(v))
print("\nThe same search over the whole raw record at neighbouring potentials (rising sweep, and -0.2 V on the falling sweep):")
for target in (0.5, 0.7, -0.2):
kv = int(np.argmin(np.abs(v[rise] - target))) if target > 0 else top + int(np.argmin(np.abs(v[top:] - target)))
cv = signal.sosfiltfilt(hp, raw["current"][:, kv])
found = []
for r in [0, 200, 400, 733]:
k, rr = mp.best_lag(signal.sosfiltfilt(hp, fs["iox"][r]), cv)
found.append(f"row {r}: best r {rr:.3f} at {fs['start'][r] - raw['scan_start'][k]:+.0f} s")
print(f" {v[kv]:+.3f} V: " + "; ".join(found))
print("\nCriterion of fig3a-redox-peaks computed at the stored times (not deciding, the clock is not established):")
for name, mask in (("large", s13.success & s13.big), ("small", s13.success & s13.small)):
i = np.searchsorted(raw["scan_start"], s13.value_cue[mask].values)
i = i[(i > 0) & (i + 20 < len(raw["current"]))]
vg = np.mean([raw["current"][k + 5:k + 20].mean(0) - raw["current"][k] for k in i], 0)
a_, b_ = int(np.argmax(vg[:top + 1])), top + int(np.argmin(vg[top:]))
print(f" {name}-reward forced trials (n {len(i)}), mean of 0.5 to 2 s minus the scan at the cue: largest on the rising "
f"sweep {vg[a_]:+.3f} at {v[a_]:+.3f} V; most negative on the falling sweep {vg[b_]:+.3f} at {v[b_]:+.3f} V; "
f"at 0.6 V {vg[k06]:+.3f}; SD over the scan {vg.std():.3f}")
current at 0.608 V (sample 61) compared with the stored iox, both high-passed at 0.5 Hz; search over the whole raw record: row 0 (trial 1): r at the stored time -0.058; best r 0.172, raw -14678.2 s from the stored time row 50 (trial 51): r at the stored time -0.001; best r 0.202, raw -11607.4 s from the stored time row 200 (trial 201): r at the stored time -0.017; best r 0.186, raw -1421.0 s from the stored time row 400 (trial 401): r at the stored time -0.014; best r 0.177, raw -14810.5 s from the stored time row 600 (trial 601): r at the stored time -0.014; best r 0.164, raw +3754.1 s from the stored time row 733 (trial 734): r at the stored time -0.021; best r 0.172, raw -1793.7 s from the stored time Raw current at 0.6 V, mean over 0.5 to 3 s after the value cue minus at the cue, large minus small (all successful forced trials), at the stored times and shifted: shift +0 s: difference -0.072, Mann-Whitney p = 0.456, n = 303 shift -60 s: difference +0.027, Mann-Whitney p = 0.318, n = 302 shift -30 s: difference +0.059, Mann-Whitney p = 0.033, n = 302 shift -10 s: difference +0.059, Mann-Whitney p = 0.020, n = 302 shift -5 s: difference +0.024, Mann-Whitney p = 0.426, n = 302 shift +5 s: difference -0.020, Mann-Whitney p = 0.703, n = 303 shift +10 s: difference +0.032, Mann-Whitney p = 0.159, n = 303 shift +30 s: difference -0.006, Mann-Whitney p = 0.978, n = 303 shift +60 s: difference +0.059, Mann-Whitney p = 0.054, n = 303 The same search over the whole raw record at neighbouring potentials (rising sweep, and -0.2 V on the falling sweep):
+0.496 V: row 0: best r 0.194 at -15369 s; row 200: best r 0.223 at -10644 s; row 400: best r 0.184 at -4815 s; row 733: best r 0.181 at -1388 s
+0.704 V: row 0: best r 0.178 at -13293 s; row 200: best r 0.200 at -9926 s; row 400: best r 0.195 at -4815 s; row 733: best r 0.185 at -7942 s -0.192 V: row 0: best r 0.182 at -13681 s; row 200: best r 0.185 at -10644 s; row 400: best r 0.190 at -6771 s; row 733: best r 0.174 at +5424 s Criterion of fig3a-redox-peaks computed at the stored times (not deciding, the clock is not established): large-reward forced trials (n 164), mean of 0.5 to 2 s minus the scan at the cue: largest on the rising sweep +0.039 at -0.240 V; most negative on the falling sweep -0.129 at +1.339 V; at 0.6 V -0.033; SD over the scan 0.042 small-reward forced trials (n 139), mean of 0.5 to 2 s minus the scan at the cue: largest on the rising sweep +0.033 at +0.528 V; most negative on the falling sweep -0.063 at +0.112 V; at 0.6 V +0.028; SD over the scan 0.022
The raw current at 0.6 V does not match the stored oxidation current iox at the stored times (|r| below
0.06), and searched over the whole raw record the best matches are weak (r 0.16 to 0.20) at offsets
that differ from row to row, so no single offset links the two. The raw current shows no large-against-small
contrast at the stored value cue times (p = 0.46), while two of the eight shifted times give p below 0.05 (-10 s, p = 0.020; -30 s, p = 0.033), so a
contrast with the task at some offset is not ruled out. The same
search at 0.50 and 0.70 V and at -0.19 V on the falling sweep gives the same picture (best r 0.17 to 0.22, at
offsets that differ from row to row). Computed at the stored times, the background-subtracted voltammogram of
large-reward trials has its largest rising-sweep value at -0.24 V, not near 0.6 V (where it is -0.033), and
its most negative value at the switching potential (1.34 V), not near -0.2 V; its values (SD 0.042 over the
scan) are small. This is not a test of the claim, since the clock is not established. The raw
FSCV therefore cannot be shown to share the trial clock, and may not be the recording the trial-aligned
series came from. It is not used for any claim (fig3a-redox-peaks is blocked).
lfp_t, lfp_x, lfp_conv = mp.load_lfp(files["T-0926"])
d = np.diff(lfp_t)
print(f"bipolar LFP: {len(lfp_t)} samples, {lfp_t[0]:.3f} to {lfp_t[-1]:.3f} s; median interval {np.median(d) * 1e3:.3f} ms; "
f"intervals off 1 ms by more than 1 µs: {int(np.sum(np.abs(d - 1e-3) > 1e-6))}, smallest {d.min():.2e} s")
print(f"stored values {np.percentile(lfp_x, [0.5, 50, 99.5])} ; conversion attribute {lfp_conv}; unit "
f"{files['T-0926']['processing/ecephys/FilteredEphys/bipolar_lfp_series/data'].attrs['unit']}")
print(f"stored values x conversion: SD {np.std(lfp_x) * lfp_conv:.3g} V")
jn = np.where(np.abs(d - 1e-3) > 1e-6)[0]
dx = np.abs(np.diff(lfp_x))
print(f"joins (intervals off 1 ms): {len(jn)}, intervals {d[jn].min() * 1e3:.2g} to {d[jn].max() * 1e3:.3f} ms, longer than 1 ms: "
f"{int(np.sum(d[jn] > 1e-3 + 1e-6))}")
print(f"step between neighbouring samples: median {np.median(dx) * 1e6:.2f} µV over all samples, {np.median(dx[jn]) * 1e6:.2f} µV at the "
f"joins; largest at a join {dx[jn].max() * 1e6:.2f} µV, 99.9th percentile over all samples {np.percentile(dx, 99.9) * 1e6:.2f} µV")
okj = d26[d26.success]
rel = lfp_t[jn][None, :] - okj.value_cue.values[:, None]
has_join = np.any((rel >= -0.5) & (rel < 1.0), axis=1)
print("joins after the preceding trial start (s), percentiles 0, 5, 50, 95, 100:",
np.percentile(lfp_t[jn] - d26.start.values[np.searchsorted(d26.start.values, lfp_t[jn]) - 1], [0, 5, 50, 95, 100]).round(2))
print("successful trials with a join within -0.5 to 1 s of the value cue:", int(has_join.sum()), "of", len(okj))
print(pd.crosstab(has_join, [okj.big.values, okj.left.values], rownames=["join near cue"], colnames=["big", "left"]).to_string())
st26 = d26.start.values
print("trial windows covered by the LFP:", int(np.sum((st26 - 30 >= lfp_t[0] - 1e-6) & (st26 + 30 <= lfp_t[-1] + 1e-3))), "of", len(st26))
e = files["T-0926"]["acquisition/electrical_series"]
print(f"raw ephys: {e['data'].shape} at {e['starting_time'].attrs['rate']} Hz from starting_time {e['starting_time'][()]}; "
f"duration {e['data'].shape[0] / 32000:.1f} s; recording opened and closed "
f"{[x.decode() for x in files['T-0926']['general/extracellular_ephys/electrodes/recording_opened'][:1]]} "
f"{[x.decode() for x in files['T-0926']['general/extracellular_ephys/electrodes/recording_closed'][:1]]}; "
f"no timestamps stored: {'timestamps' not in e}")
lp = signal.butter(4, 100, fs=1000, output="sos")
tb, xb = mp.raw_bipolar_1khz(files["T-0926"], 0.0, e["data"].shape[0] / 32000)
xb = signal.sosfiltfilt(lp, xb)
print(f"raw c3bs - c3a averaged to 1 kHz: {len(xb)} samples")
for trial in (2, 84, 300, 700, 1100, 1500):
vc = d26.start[trial] + 2.0
j = np.searchsorted(lfp_t, vc - 20)
ref = signal.sosfiltfilt(lp, lfp_x[j:j + 40000])
i_same = np.searchsorted(tb, lfp_t[j])
r_same = np.corrcoef(xb[i_same:i_same + 40000], ref)[0, 1] if i_same + 40000 <= len(xb) else np.nan
k, rr = mp.best_lag(ref, xb)
scale = np.std(xb[k:k + 40000]) / np.std(ref)
print(f"trial {trial}: LFP against raw (both low-passed 100 Hz), whole raw searched: r at the same stored time {r_same:.3f}; "
f"best r {rr:.3f} with the raw {lfp_t[j] - tb[k]:.3f} s earlier than its stored time; raw SD / LFP SD {scale:.3f}")
bipolar LFP: 19247100 samples, 12919.054 to 32165.400 s; median interval 1.000 ms; intervals off 1 ms by more than 1 µs: 1497, smallest 1.82e-12 s stored values [-3.20552381e-04 -7.94835903e-06 4.55821742e-04] ; conversion attribute 3.0517578125e-08; unit volts
stored values x conversion: SD 4.06e-12 V joins (intervals off 1 ms): 1497, intervals 1.8e-09 to 0.999 ms, longer than 1 ms: 0
step between neighbouring samples: median 5.17 µV over all samples, 2.73 µV at the joins; largest at a join 20.28 µV, 99.9th percentile over all samples 25.99 µV joins after the preceding trial start (s), percentiles 0, 5, 50, 95, 100: [1.000e-02 9.100e-01 1.000e+01 1.125e+01 1.572e+01] successful trials with a join within -0.5 to 1 s of the value cue: 118 of 1172 big False True left False True False True join near cue False 262 276 257 259 True 28 36 31 23 trial windows covered by the LFP: 1500 of 1501 raw ephys: (722978367, 2) at 32000.0 Hz from starting_time 0.0; duration 22593.1 s; recording opened and closed ['2024-09-26 09:01:38'] ['2024-09-26 16:02:14']; no timestamps stored: True
raw c3bs - c3a averaged to 1 kHz: 22593073 samples
trial 2: LFP against raw (both low-passed 100 Hz), whole raw searched: r at the same stored time -0.472; best r 0.982 with the raw 9831.349 s earlier than its stored time; raw SD / LFP SD 1.019
trial 84: LFP against raw (both low-passed 100 Hz), whole raw searched: r at the same stored time -0.023; best r 0.924 with the raw 9831.350 s earlier than its stored time; raw SD / LFP SD 1.021
trial 300: LFP against raw (both low-passed 100 Hz), whole raw searched: r at the same stored time -0.170; best r 0.981 with the raw 9831.354 s earlier than its stored time; raw SD / LFP SD 1.018
trial 700: LFP against raw (both low-passed 100 Hz), whole raw searched: r at the same stored time 0.340; best r 0.979 with the raw 9831.360 s earlier than its stored time; raw SD / LFP SD 1.020
trial 1100: LFP against raw (both low-passed 100 Hz), whole raw searched: r at the same stored time nan; best r 0.983 with the raw 9831.366 s earlier than its stored time; raw SD / LFP SD 1.015
trial 1500: LFP against raw (both low-passed 100 Hz), whole raw searched: r at the same stored time nan; best r 0.954 with the raw 9831.371 s earlier than its stored time; raw SD / LFP SD 1.027
The LFP covers every trial window but the last (whose window runs past the end of the series), at 1 kHz with small irregularities where the conversion joined trial windows. Its stored values are of the order of 1e-4, a plausible LFP in volts; the stored conversion factor (3.05e-8) would make them picovolts, so the conversion attribute is wrong and the values are used as stored (volts). The 1497 joins are all intervals shorter than 1 ms (overlaps of the trial windows, none longer than 1 ms), and the signal does not step at them: the change between the two samples at a join (median 2.73 µV, largest 20.28 µV) is no larger than between other neighbouring samples (median 5.17 µV, 99.9th percentile 25.99 µV). 118 of the 1172 successful trials have a join within -0.5 to 1 s of the value cue, spread over the four reward and side conditions (23 to 36 each). The raw ephys has a start time of 0 and a rate, and no timestamps. Searched over the whole record, the raw bipolar signal matches the stored LFP closely (r above 0.9, similar SD) when it is read 9831.35 s earlier than its stored time, the same offset at every trial tested, to within 0.03 s from the first trial to the last. The raw ephys of this session is therefore shifted from the trial clock by 9831 s; its stored starting time is wrong. The claims use only the LFP, which the conversion built from the trial-aligned segments, so it is on the trial clock by construction; the next cell checks that its responses are locked to the cue.
beta = mp.beta_power(lfp_x)
beta_env = mp.beta_power(lfp_x, smooth_s=0) # unsmoothed envelope, for variants
ok26 = d26[d26.success]
vc26 = ok26.value_cue.values
def cue_drop(shift):
pre = mp.window_mean(lfp_t, beta, vc26 + shift, -0.5, 0.0)
post = mp.window_mean(lfp_t, beta, vc26 + shift, 0.0, 0.5)
return np.median(post / pre), stats.wilcoxon(post, pre).pvalue
print("Beta power, successful trials, median of (0 to 0.5 s) / (-0.5 to 0 s) around the value cue:")
for shift in [0, -1.5, -1.0, 1.0, 1.5, 3.0]:
ratio, p = cue_drop(shift)
print(f" cue shifted {shift:+.1f} s: ratio {ratio:.3f}, Wilcoxon p = {p:.2g}")
okf_rows = np.where(sel)[0]
x = fs["da"][okf_rows]
c = cue13[sel]
def da_rise(shift):
post = mp.fscv_window_mean(x, c + shift, 0.0, 3.0)
pre = mp.fscv_window_mean(x, c + shift, -1.0, 0.0)
return np.nanmedian(post - pre), stats.wilcoxon(post - pre, nan_policy="omit").pvalue
print("\nStored da, successful forced trials with FSCV, median of mean(0 to 3 s) - mean(-1 to 0 s) around the value cue:")
for shift in [0, -10, -5, 5, 10]:
dd, p = da_rise(shift)
print(f" cue shifted {shift:+d} s: {dd:+.2f}, Wilcoxon p = {p:.2g}")
Beta power, successful trials, median of (0 to 0.5 s) / (-0.5 to 0 s) around the value cue:
cue shifted +0.0 s: ratio 0.495, Wilcoxon p = 1.2e-148
cue shifted -1.5 s: ratio 1.891, Wilcoxon p = 2.1e-78
cue shifted -1.0 s: ratio 2.205, Wilcoxon p = 6.5e-164
cue shifted +1.0 s: ratio 1.264, Wilcoxon p = 2.2e-27
cue shifted +1.5 s: ratio 0.950, Wilcoxon p = 0.005
cue shifted +3.0 s: ratio 1.014, Wilcoxon p = 0.98 Stored da, successful forced trials with FSCV, median of mean(0 to 3 s) - mean(-1 to 0 s) around the value cue: cue shifted +0 s: +0.71, Wilcoxon p = 0.00019 cue shifted -10 s: -1.55, Wilcoxon p = 1.9e-09 cue shifted -5 s: -0.20, Wilcoxon p = 0.72 cue shifted +5 s: -1.58, Wilcoxon p = 3.8e-10 cue shifted +10 s: +0.39, Wilcoxon p = 0.0098
At the true value cue times beta power over 0 to 0.5 s is half of that over -0.5 to 0 s (ratio 0.495); with the cue moved by 1.5 or 3 s the ratio is 0.95 and 1.01, and moved back by 1 or 1.5 s it is above 1. The stored dopamine rises after the true value cue (+0.71, p = 0.00019); moved by 5 or 10 s the change is negative or smaller. The trial-aligned series are consistent with sharing the trial clock; this is shown by responses locked to the cue, not proven.
fP = files["P-0813"]
tsum = mp.timestamp_summary(fP)
print(f"Monkey P raw: {tsum['n']} samples, timestamps {tsum['first']:.6f} to {tsum['last']:.3f} s, median interval "
f"{tsum['median_dt'] * 1e6:.3f} µs ({1 / tsum['median_dt']:.1f} Hz); gaps {len(tsum['gaps'])}, totalling "
f"{sum(g[1] for g in tsum['gaps']):.3f} s")
units = mp.ragged(fP, "units", "spike_times")
names = [x.decode() for x in fP["units/unit_name"][:]]
thr = mp.ragged(fP, "processing/ecephys/thresholded_units", "spike_times")
thr_names = [x.decode() for x in fP["processing/ecephys/thresholded_units/unit_name"][:]]
for n_, s in zip(names, units):
print(f"units table, unit {n_}: {len(s)} spikes, {s.min():.4f} to {s.max():.3f} s, sorted {bool(np.all(np.diff(s) >= 0))}, "
f"repeated times {len(s) - len(np.unique(s))}, below 1 s: {int(np.sum(s < 1))}")
for n_, s in zip(thr_names, thr):
print(f"thresholded, unit {n_}: {len(s)} spikes, {s.min():.4f} to {s.max():.3f} s, sorted {bool(np.all(np.diff(s) >= 0))}, "
f"below 6800 s: {int(np.sum(s < 6800))}; shared with units-table unit 0: {int(np.isin(s, units[0]).sum())}, "
f"unit 1: {int(np.isin(s, units[1]).sum())}")
print("units table description:", fP["units"].attrs["description"][:160], "...")
z_t, z_ttl, z_str = mp.zenodo_events("session126_events.mat")
offset_P = z_t[0]
print(f"\nZenodo session126_events.mat: {len(z_t)} events from {z_t[0]:.5f} to {z_t[-1]:.3f} s; first event {str(z_str[0])!r}")
print(f"its {len(np.unique(z_ttl))} distinct TTL values and their counts (no names are given):", {int(a): int(b) for a, b in zip(*np.unique(z_ttl, return_counts=True))})
wf, wf_ends = mp.thresholded_waveforms(fP)
w1 = wf[wf_ends[1]:wf_ends[2]]
n10 = 30000 * 600
for start_s in (1000, 7000):
t_raw, x_raw = mp.raw_slice(fP, int(start_s * 30000), n10)
xf = mp.spike_band(x_raw, 30000.0)
print(f"\nraw window {t_raw[0]:.1f} to {t_raw[-1]:.1f} s: noise (MAD/0.6745, 300-6000 Hz) {np.median(np.abs(xf)) / 0.6745:.2f} µV")
for off in (0.0, 6800.0, offset_P):
m, n_ = mp.spike_triggered_mean(t_raw, xf, units[1], off)
print(f" unit 1 spikes moved back {off:.5f} s: {n_} in window, trough of the spike-triggered mean {np.nanmin(m):.2f} µV")
print("trough of the stored mean waveform of unit 1:", round(float(w1.mean(0).min()), 2), "µV")
Monkey P raw: 359560287 samples, timestamps 0.000000 to 11985.419 s, median interval 33.333 µs (30000.0 Hz); gaps 9, totalling 0.076 s units table, unit 0: 30765 spikes, 0.0003 to 18785.219 s, sorted True, repeated times 3, below 1 s: 7 units table, unit 1: 31247 spikes, 6800.7506 to 18784.397 s, sorted True, repeated times 0, below 1 s: 0 thresholded, unit -1: 10885 spikes, 0.0003 to 18784.620 s, sorted False, below 6800 s: 7; shared with units-table unit 0: 10885, unit 1: 0 thresholded, unit 0: 19880 spikes, 6800.6664 to 18785.219 s, sorted True, below 6800 s: 0; shared with units-table unit 0: 19880, unit 1: 0 thresholded, unit 1: 31247 spikes, 6800.7506 to 18784.397 s, sorted True, below 6800 s: 0; shared with units-table unit 0: 0, unit 1: 31247 units table description: Spike-sorted units from Plexon Offline Sorter. Units from Plexon channel(s) csc12_100 of csc12.plx. Spike times were regenerated from csc12.plx with neo includi ...
Zenodo session126_events.mat: 20812 events from 6800.66578 to 18786.083 s; first event 'Starting Recording'
its 36 distinct TTL values and their counts (no names are given): {0: 826, 2: 1, 3: 824, 4: 932, 5: 824, 6: 211, 7: 614, 8: 824, 10: 613, 11: 593, 12: 113, 13: 480, 14: 344, 18: 53, 19: 24, 20: 825, 22: 33, 23: 20, 24: 390, 25: 435, 26: 824, 27: 825, 29: 434, 30: 390, 31: 414, 32: 410, 33: 1853, 34: 1854, 35: 825, 36: 824, 37: 109, 38: 109, 45: 270, 46: 210, 100: 1660, 110: 822}
raw window 1000.0 to 1600.0 s: noise (MAD/0.6745, 300-6000 Hz) 5.25 µV unit 1 spikes moved back 0.00000 s: 0 in window, trough of the spike-triggered mean nan µV unit 1 spikes moved back 6800.00000 s: 2251 in window, trough of the spike-triggered mean -0.46 µV unit 1 spikes moved back 6800.66578 s: 2251 in window, trough of the spike-triggered mean -0.29 µV
raw window 7000.0 to 7600.0 s: noise (MAD/0.6745, 300-6000 Hz) 5.41 µV unit 1 spikes moved back 0.00000 s: 1745 in window, trough of the spike-triggered mean -0.32 µV unit 1 spikes moved back 6800.00000 s: 1721 in window, trough of the spike-triggered mean -0.51 µV unit 1 spikes moved back 6800.66578 s: 1719 in window, trough of the spike-triggered mean -26.43 µV trough of the stored mean waveform of unit 1: -30.64 µV
The raw signal is sampled at 30 kHz with 9 short gaps; its timestamps run from 0 s. The spike times run from about 6800 to 18785 s. Moved back by 6800.666 s (the time of the first event, "Starting Recording", in the session's events file on Zenodo, read above), the raw signal around the spikes of unit 1 averages to a spike with a trough close to that of the stored waveforms; without the shift, or with a shift of 6800 s, there is nothing. The spike tables and the raw signal of this file are on different clocks, 6800.666 s apart, and the file does not say so. In the units table, unit 0 merges the thresholded spikes labelled 0 (unsorted) and -1 (poor quality), and it holds 7 spike times below 1 s and 3 repeated times; the thresholded spikes labelled -1 are not in time order. Only unit 1 is the sorted unit. The Zenodo events file holds 20812 events with 36 distinct TTL values and no names. They do not follow the Monkey T codes: the trial-start code 128 does not occur, and the value cue code 23 occurs 20 times, against about 825 for several others. The trials of this session cannot be decoded from them.
fscv-waveform¶
FSCV applied a triangular waveform from -0.4 to 1.3 V at 400 V/s, every 100 ms (10 Hz).
step = np.median(np.abs(np.diff(v[rise]))) / dt
scan_hz = 1 / np.median(np.diff(raw["scan_start"]))
measured = {"lowest (V)": v.min(), "highest (V)": v.max(), "sweep rate (V/s)": step, "scan frequency (Hz)": scan_hz}
expected = {"lowest (V)": (-0.4, 0.04), "highest (V)": (1.3, 0.13), "sweep rate (V/s)": (400, 40), "scan frequency (Hz)": (10, 1)}
for k_, val in measured.items():
e_, tol = expected[k_]
print(f"{k_}: {val:.4f}; paper {e_} +- {tol}; within tolerance: {abs(val - e_) <= tol}")
i_ = meta["description"].index("triangular")
print("dandiset description:", meta["description"][i_ - 20:i_ + 70])
lowest (V): -0.3470; paper -0.4 +- 0.04; within tolerance: False highest (V): 1.3395; paper 1.3 +- 0.13; within tolerance: True sweep rate (V/s): 400.4639; paper 400 +- 40; within tolerance: True scan frequency (Hz): 10.0000; paper 10 +- 1; within tolerance: True dandiset description: as performed with a triangular waveform (−0.35 to 1.34 V, 400 V/s) applied at 10 Hz, with
Result: partially reproduced. Three of the four values are within tolerance: the highest potential (1.340 V), the sweep rate and the scan frequency. The lowest potential is -0.347 V, outside -0.4 +- 0.04 V; the dandiset's own description (-0.35 V) matches the data.
ephys-acquisition¶
EPhys was sampled at 30 kHz in monkey P and 32 kHz in monkey T, and band-pass filtered from 0.1 to 7500 Hz.
el = files["T-0926"]["general/extracellular_ephys/electrodes"]
vals = {"rate, monkey P (Hz)": (1 / tsum["median_dt"], 30000, 3000),
"rate, monkey T (Hz)": (float(e["starting_time"].attrs["rate"]), 32000, 3200),
"low cut, monkey T (Hz)": (float(el["DspLowCutFrequency"][0]), 0.1, 0.01),
"high cut, monkey T (Hz)": (float(el["DspHighCutFrequency"][0]), 7500, 750)}
for k_, (val, e_, tol) in vals.items():
print(f"{k_}: {val:.4g}; paper {e_} +- {tol}; within tolerance: {abs(val - e_) <= tol}")
print("both monkey T channels:", el["DspLowCutFrequency"][:], el["DspHighCutFrequency"][:], "; input range (µV):", el["InputRange"][:])
notes_P = files["P-0813"]["general/notes"][()].decode()
print(f"monkey P notes ({len(notes_P)} characters):", notes_P)
print("filter settings in monkey P's notes or electrode columns:",
re.findall(r"(?i)[^,{]*(?:filter|cut|dsp)[^,}]*", notes_P),
[c for c in files["P-0813"]["general/extracellular_ephys/electrodes"] if re.search("(?i)filter|cut|dsp", c)])
i_ = meta["description"].index("7500")
print("dandiset description:", meta["description"][i_ - 120:i_ + 10])
rate, monkey P (Hz): 3e+04; paper 30000 +- 3000; within tolerance: True
rate, monkey T (Hz): 3.2e+04; paper 32000 +- 3200; within tolerance: True
low cut, monkey T (Hz): 0.1; paper 0.1 +- 0.01; within tolerance: True
high cut, monkey T (Hz): 8000; paper 7500 +- 750; within tolerance: True
both monkey T channels: [0.1 0.1] [8000 8000] ; input range (µV): [1000 1000]
monkey P notes (654 characters): {"FileVersion": "3.4", "TimeClosed": "2018/08/13 15:32:25", "input_inverted": "False", "ProbeName": "-OriginalFileName \"C:\\CheetahData\\2018-08-13_10-19-56\\CSC12.ncs\"", "channel_names": "['CSC12']", "TimeCreated": "2018/08/13 10:20:08", "RecordSize": "1044", "NumADChannels": "1", "FileUUID": "c2b561b7-088c-4ac2-b2f1-8d900cd15b49", "FileType": "NCS", "recording_closed": "2018-08-13 15:32:25", "ReferenceChannel": "\"Source 01 Reference 1\"", "channel_ids": "[11]", "OriginalFileName": "C:\\CheetahData\\2018-08-13_10-19-56\\CSC12.ncs", "ADMaxValue": "32767", "sampling_rate": "30000.0", "bit_to_microVolt": "[0.030517578125]", "InputRange": "1000"}
filter settings in monkey P's notes or electrode columns: [] []
dandiset description: ough a unity-gain headstage (HS-36, ±1 mV input range) at 30 kHz (Monkey P) or 32 kHz (Monkey T), bandpass filtered 0.1–7500 Hz. F
Result: reproduced. The sampling rates are 30000 and 32000 Hz and the acquisition filter's low cut is 0.1 Hz. The high cut is 8000 Hz, inside the 10% tolerance of the 7500 Hz stated by the paper and the dandiset description, but not equal to it. Monkey P's file stores no filter settings: its notes (printed in full) and its electrode columns name no filter or cut-off. The dandiset description states 0.1 to 7500 Hz.
shape-task-structure¶
In the shape task, forced trials made up 25% of all trials, and the cue-reward association was reversed after blocks of 50 to 300 trials.
share = 100 * s13.forced.mean()
blocks = mp.reward_blocks(s13)
print(f"forced trials: {int(s13.forced.sum())} of {len(s13)} = {share:.1f}% (paper 25 +- 2.5)")
print(blocks.to_string(index=False))
complete = blocks[blocks.complete]
print("complete blocks:", len(complete), "; lengths", complete.n_trials.tolist(), "; all within 45 to 330:",
bool(complete.n_trials.between(45, 330).all()))
fh = blocks[blocks.first_trial <= 734]
print("blocks starting in the first 734 trials (those with FSCV rows):", len(fh))
forced trials: 358 of 1468 = 24.4% (paper 25 +- 2.5)
first_trial last_trial n_trials large_reward_image complete
1 61 61 1 False
62 233 172 8 True
234 386 153 1 True
387 507 121 8 True
508 578 71 1 True
579 772 194 8 True
773 865 93 1 True
866 1001 136 8 True
1002 1274 273 1 True
1275 1468 194 8 False
complete blocks: 8 ; lengths [172, 153, 121, 71, 194, 93, 136, 273] ; all within 45 to 330: True
blocks starting in the first 734 trials (those with FSCV rows): 6
Result: reproduced. Forced trials are 24.4% of the 1468 trials (25 +- 2.5%). The eight complete blocks are 71 to 273 trials long, inside the widened range of 45 to 330 trials.
zenodo-link¶
The dandiset reproduces the processed data the paper deposited on Zenodo; the trial selections and the beta power of the paper's example trial computed from the dandiset agree with the paper's own files.
lists = {"big rew": d26.success & d26.big, "small rew": d26.success & d26.small,
"left": d26.success & d26.left, "right": d26.success & d26.right}
zen = {}
for name in lists:
zen[name] = mp.zenodo_mat(f"09262024_trialids_value object start_{name} and success_tr1-1501.mat")["trids_out"].ravel().astype(int)
for name, mask in lists.items():
mine, theirs = set(d26.index[mask]), set(zen[name])
both = len(mine & theirs)
print(f"{name}: Zenodo {len(theirs)}, here {len(mine)}, both {both} ({100 * both / max(len(mine), len(theirs)):.2f}% of the "
f"larger); only Zenodo {sorted(theirs - mine)}, only here {sorted(mine - theirs)}")
ex = mp.zenodo_mat("09262024_c3bs-c3a_data_trial35bigrewardexample.mat")
bp_paper = ex["trial35bigrewardataBetaPower"].ravel()
t_paper = ex["trial35bigrewardataTS"].ravel()
print("Zenodo example: its 'data' variable equals its time variable:", bool(np.allclose(ex["trial35bigrewardata"].ravel(), t_paper)),
"; times", t_paper[0], "to", t_paper[-1], "s")
trial35 = int(zen["big rew"][34])
vc = d26.value_cue[trial35]
i0 = np.searchsorted(lfp_t, vc - 3)
k, rr = mp.best_lag(bp_paper, beta[i0 - 100:i0 + 100 + len(bp_paper)])
lag_ms = k - 100
seg = beta[i0:i0 + len(bp_paper)]
print(f"35th large-reward trial = trial {trial35}; r at lag 0 {np.corrcoef(seg, bp_paper)[0, 1]:.3f}; best r {rr:.3f} at {lag_ms} ms; "
f"median ratio here / Zenodo {np.median(seg / bp_paper):.3f}")
fig, ax = plt.subplots(figsize=(8, 2.5))
ax.plot(t_paper, bp_paper * 1e12, color="0.5", lw=2, label="Zenodo file")
ax.plot(t_paper, seg * 1e12, color=SMALL, lw=1, label="computed here")
ax.set_xlabel("time from value cue (s)"); ax.set_ylabel("beta power (µV²)"); ax.legend(frameon=False)
ax.set_title(f"Trial {trial35}, beta power"); plt.show()
print(f"plotted range: Zenodo {bp_paper.min() * 1e12:.1f} to {bp_paper.max() * 1e12:.1f} µV², here {seg.min() * 1e12:.1f} to {seg.max() * 1e12:.1f}")
big rew: Zenodo 569, here 570, both 569 (99.82% of the larger); only Zenodo [], only here [1501] small rew: Zenodo 602, here 602, both 602 (100.00% of the larger); only Zenodo [], only here [] left: Zenodo 594, here 594, both 594 (100.00% of the larger); only Zenodo [], only here [] right: Zenodo 577, here 578, both 577 (99.83% of the larger); only Zenodo [], only here [1501]
Zenodo example: its 'data' variable equals its time variable: True ; times -3.0 to 12.0 s 35th large-reward trial = trial 84; r at lag 0 0.942; best r 0.943 at 6 ms; median ratio here / Zenodo 0.943
plotted range: Zenodo 36.0 to 1559.6 µV², here 14.1 to 2184.5
Result: reproduced. The four trial lists agree for 99.82% to 100% of their trials; the only difference is trial 1501, the session's last trial, which the Zenodo lists leave out. The paper's example "trial35" is the 35th trial of its large-reward list, trial 84, and its beta power correlates with the one computed here at r = 0.943, at a lag of 6 ms, with a median amplitude ratio of 0.943. The Zenodo example file's raw-data variable is a copy of its time variable, so its LFP trace itself cannot be compared.
fig3b-beta-bursts¶
Single-trial bipolar LFP from site c3bs-c3a shows beta bursts (one to a few cycles of high beta power) whose frequency changes around the cue and other task events.
in_trial = np.zeros(len(lfp_t), bool)
for s_ in st26:
in_trial[np.searchsorted(lfp_t, s_ - 30):np.searchsorted(lfp_t, s_ + 30)] = True
big26 = d26[d26.success & d26.big]
vcb = big26.value_cue.values
min_len = int(mp.BURST_MIN_S * mp.LFP_RATE)
for pct in (75, 50, 90):
thr_ = np.percentile(beta[in_trial], pct)
on_t = lfp_t[mp.burst_onsets(beta, thr_, min_len)]
pre = np.array([np.sum((on_t >= c_ - 1) & (on_t < c_)) for c_ in vcb])
post = np.array([np.sum((on_t >= c_) & (on_t < c_ + 1)) for c_ in vcb])
w = stats.wilcoxon(post, pre)
print(f"threshold at the {pct}th percentile ({thr_ * 1e12:.1f} µV²): bursts per s, -1 to 0 s {pre.mean():.3f}, 0 to 1 s "
f"{post.mean():.3f}; trials with more after / fewer after / equal {np.sum(post > pre)}/{np.sum(post < pre)}/"
f"{np.sum(post == pre)}; Wilcoxon p = {w.pvalue:.2g}; n = {len(vcb)} trials" + (" <- deciding" if pct == 75 else ""))
thr_env = np.percentile(beta_env[in_trial], 75)
on_env = lfp_t[mp.burst_onsets(beta_env, thr_env, min_len)]
pre = np.array([np.sum((on_env >= c_ - 1) & (on_env < c_)) for c_ in vcb])
post = np.array([np.sum((on_env >= c_) & (on_env < c_ + 1)) for c_ in vcb])
print(f"variant, unsmoothed envelope, 75th percentile ({thr_env * 1e12:.1f} µV²), at least 50 ms: bursts per s, -1 to 0 s "
f"{pre.mean():.3f}, 0 to 1 s {post.mean():.3f}; Wilcoxon p = {stats.wilcoxon(post, pre).pvalue:.2g}")
thr75 = np.percentile(beta[in_trial], 75)
for name_, x_, th_ in (("smoothed power", beta, thr75), ("unsmoothed envelope", beta_env, thr_env)):
e_ = np.diff(np.r_[False, x_ > th_, False].astype(int))
runs_s = np.where(e_ == -1)[0] - np.where(e_ == 1)[0]
print(f"{name_}, runs above its 75th percentile: {len(runs_s)}, {np.mean(runs_s >= min_len) * 100:.1f}% last at least 50 ms; "
f"median length {np.median(runs_s):.0f} ms")
on_t = lfp_t[mp.burst_onsets(beta, thr75, min_len)]
edges = np.arange(-3, 8.01, 0.25)
rate = np.array([np.histogram(on_t - c_, edges)[0] for c_ in vcb]).mean(0) / 0.25
fb = np.median(big26.feedback.values - vcb)
print("burst rate per 0.25 s bin from the value cue (bursts per s):")
print(pd.Series(rate.round(2), index=edges[:-1]).to_string())
print(f"median feedback time {fb:.2f} s after the value cue; burst rate range {rate.min():.2f} to {rate.max():.2f}")
threshold at the 75th percentile (356.2 µV²): bursts per s, -1 to 0 s 1.512, 0 to 1 s 1.207; trials with more after / fewer after / equal 123/240/207; Wilcoxon p = 9.5e-10; n = 570 trials <- deciding
threshold at the 50th percentile (159.2 µV²): bursts per s, -1 to 0 s 0.695, 0 to 1 s 1.421; trials with more after / fewer after / equal 348/50/172; Wilcoxon p = 2.5e-47; n = 570 trials
threshold at the 90th percentile (678.6 µV²): bursts per s, -1 to 0 s 1.475, 0 to 1 s 0.661; trials with more after / fewer after / equal 59/349/162; Wilcoxon p = 2.4e-45; n = 570 trials
variant, unsmoothed envelope, 75th percentile (350.1 µV²), at least 50 ms: bursts per s, -1 to 0 s 3.074, 0 to 1 s 2.105; Wilcoxon p = 6.8e-34 smoothed power, runs above its 75th percentile: 17265, 95.1% last at least 50 ms; median length 206 ms
unsmoothed envelope, runs above its 75th percentile: 46636, 67.1% last at least 50 ms; median length 67 ms burst rate per 0.25 s bin from the value cue (bursts per s): -3.00 0.46 -2.75 0.58 -2.50 0.46 -2.25 0.49 -2.00 0.55 -1.75 0.29 -1.50 1.16 -1.25 2.29 -1.00 1.79 -0.75 1.61 -0.50 1.30 -0.25 1.35 0.00 0.51 0.25 0.65 0.50 1.86 0.75 1.81 1.00 1.32 1.25 1.38 1.50 1.33 1.75 1.40 2.00 1.50 2.25 1.33 2.50 1.41 2.75 1.36 3.00 1.49 3.25 1.37 3.50 1.34 3.75 1.49 4.00 1.38 4.25 1.21 4.50 1.59 4.75 1.50 5.00 1.39 5.25 1.46 5.50 1.41 5.75 1.11 6.00 0.34 6.25 0.31 6.50 0.39 6.75 0.32 7.00 0.29 7.25 0.18 7.50 0.21 7.75 0.18 median feedback time 3.78 s after the value cue; burst rate range 0.18 to 2.29
Result: reproduced. With the threshold at the 75th percentile, large-reward trials have 1.51 bursts per s in the second before the value cue and 1.21 in the second after it (240 trials fewer, 123 more; Wilcoxon p = 9.5e-10, 570 trials). The rate changes around the cue under the other thresholds too, though the direction depends on the threshold: at the 50th percentile it rises (0.70 to 1.42 per s) and at the 90th it falls (1.48 to 0.66), both with p below 1e-40. In 0.25 s bins the rate drops to 0.51 per s just after the cue and falls again from 6 s after it (0.39 per s and below), after the feedback, which comes at a median of 3.78 s. On the smoothed power 95.1% of runs above the threshold last at least 50 ms (median 206 ms), so the minimum length does little; on the unsmoothed envelope (67.1% of runs last 50 ms) the rate also falls after the cue (3.07 to 2.11 per s, p = 6.8e-34). At the deciding threshold the fall in bursts after the cue coincides with the drop in beta power after the cue (power ratio 0.495 in the data check), so it is not separate evidence from the suppression tested in fig4b-beta-reward.
fig4a-dopamine-reward¶
Trial-averaged dopamine at site c8dg is higher after large-reward than after small-reward value cues, in successful forced trials of the shape task.
big_f = s13.big.values[rows13][sel]
x_da = fs["da"][okf_rows]
x_re = mp.rereference(fs["da"], okf_rows, c)
x_iox = mp.rereference(fs["iox"], okf_rows, c)
def contrast(x, a, b, label, deciding=False):
m = mp.fscv_window_mean(x, c, a, b)
L, S = m[big_f], m[~big_f]
L, S = L[~np.isnan(L)], S[~np.isnan(S)]
p = stats.mannwhitneyu(L, S).pvalue
print(f"{label}, {a} to {b} s: large median {np.median(L):.2f} (n {len(L)}), small median {np.median(S):.2f} (n {len(S)}), "
f"difference {np.median(L) - np.median(S):+.2f}, Mann-Whitney p = {p:.2g}" + (" <- deciding" if deciding else ""))
contrast(x_re, 0, 3, "da re-referenced at the cue", True)
contrast(x_re, 0, 1.5, "da re-referenced at the cue")
contrast(x_re, 1.5, 3, "da re-referenced at the cue")
contrast(x_re, -2, 0, "da re-referenced at the cue (pre-cue check)")
contrast(x_da, 0, 3, "da as stored")
contrast(x_da, -2, 0, "da as stored (pre-cue check)")
contrast(x_iox, 0, 3, "iox re-referenced at the cue")
x_base = x_da - mp.fscv_window_mean(x_da, c, -1, 0)[:, None]
contrast(x_base, 0, 3, "da minus its mean over -1 to 0 s")
contrast(x_base, 1.5, 3, "da minus its mean over -1 to 0 s")
for label, x_ in (("re-referenced at the cue", x_re), ("minus its mean over -1 to 0 s", x_base)):
m = mp.fscv_window_mean(x_, c, 0, 3)[big_f]
m = m[~np.isnan(m)]
print(f"c8dg half of fig3a-dual-channel, large-reward trials, mean 0 to 3 s of da {label}: above zero in {np.sum(m > 0)} of "
f"{len(m)} ({100 * np.mean(m > 0):.1f}%)")
first = fs["trial"][okf_rows] <= 367
for half, mask in (("trials 1-367", first), ("trials 368-734", ~first)):
m = mp.fscv_window_mean(x_re[mask], c[mask], 0, 3)
b_ = big_f[mask]
print(f"{half}: large {np.nanmedian(m[b_]):.2f} (n {b_.sum()}), small {np.nanmedian(m[~b_]):.2f} (n {(~b_).sum()}), "
f"p = {stats.mannwhitneyu(m[b_][~np.isnan(m[b_])], m[~b_][~np.isnan(m[~b_])]).pvalue:.2g}")
print("NaN samples in the deciding window:", int(np.isnan(mp.fscv_aligned(x_re, c, 0, 3)[1]).sum()), "of",
mp.fscv_aligned(x_re, c, 0, 3)[1].size)
tt, A = mp.fscv_aligned(x_re, c, 3, 8)
fig, ax = plt.subplots(figsize=(7, 3))
for mask, col, lab in ((big_f, LARGE, "large"), (~big_f, SMALL, "small")):
mu = np.nanmean(A[mask], 0); se = np.nanstd(A[mask], 0) / np.sqrt(mask.sum())
ax.plot(tt, mu, color=col, lw=2, label=f"{lab} (n {mask.sum()})"); ax.fill_between(tt, mu - se, mu + se, color=col, alpha=0.2, lw=0)
ax.axvline(0, color="0.6", lw=1); ax.set_xlabel("time from value cue (s)"); ax.set_ylabel("re-referenced da")
ax.legend(frameon=False); ax.set_title("c8dg, successful forced trials, mean +- SE"); plt.show()
print(f"plotted means range {np.nanmin([np.nanmean(A[big_f], 0), np.nanmean(A[~big_f], 0)]):.2f} to "
f"{np.nanmax([np.nanmean(A[big_f], 0), np.nanmean(A[~big_f], 0)]):.2f}")
da re-referenced at the cue, 0 to 3 s: large median 1.28 (n 85), small median 0.19 (n 80), difference +1.09, Mann-Whitney p = 0.023 <- deciding da re-referenced at the cue, 0 to 1.5 s: large median 0.35 (n 85), small median -0.01 (n 80), difference +0.35, Mann-Whitney p = 0.41 da re-referenced at the cue, 1.5 to 3 s: large median 1.99 (n 84), small median 0.71 (n 80), difference +1.28, Mann-Whitney p = 0.0072 da re-referenced at the cue (pre-cue check), -2 to 0 s: large median -0.36 (n 85), small median 0.37 (n 80), difference -0.73, Mann-Whitney p = 0.064 da as stored, 0 to 3 s: large median -3.48 (n 87), small median -3.79 (n 81), difference +0.31, Mann-Whitney p = 0.3
da as stored (pre-cue check), -2 to 0 s: large median -4.98 (n 87), small median -3.57 (n 81), difference -1.41, Mann-Whitney p = 0.85 iox re-referenced at the cue, 0 to 3 s: large median -0.01 (n 87), small median 0.00 (n 81), difference -0.01, Mann-Whitney p = 0.41 da minus its mean over -1 to 0 s, 0 to 3 s: large median 1.38 (n 85), small median 0.10 (n 81), difference +1.28, Mann-Whitney p = 0.022 da minus its mean over -1 to 0 s, 1.5 to 3 s: large median 2.38 (n 84), small median 0.44 (n 81), difference +1.94, Mann-Whitney p = 0.0062 c8dg half of fig3a-dual-channel, large-reward trials, mean 0 to 3 s of da re-referenced at the cue: above zero in 56 of 85 (65.9%) c8dg half of fig3a-dual-channel, large-reward trials, mean 0 to 3 s of da minus its mean over -1 to 0 s: above zero in 58 of 85 (68.2%) trials 1-367: large 1.62 (n 40), small 0.08 (n 49), p = 0.013 trials 368-734: large 1.02 (n 47), small 0.70 (n 32), p = 0.51 NaN samples in the deciding window: 133 of 5208
plotted means range -4.47 to 2.14
Result: reproduced. Over 0 to 3 s after the value cue, the dopamine re-referenced at the cue has a median of 1.28 in 85 large-reward trials and 0.19 in 80 small-reward trials (Mann-Whitney p = 0.023). The difference is in 1.5 to 3 s (1.99 against 0.71, p = 0.0072), not in 0 to 1.5 s (p = 0.41). It is weak: before the cue the two differ in the other direction with p = 0.064; the stored da without re-referencing gives no difference (p = 0.3), nor does the re-referenced oxidation current (p = 0.41); and only the first half of the FSCV trials shows it (trials 1 to 367 p = 0.013, trials 368 to 734 p = 0.51). With a baseline of the mean over -1 to 0 s in place of the single sample at the cue, the result is the same (1.38 against 0.10, p = 0.022). Only the first 734 of the session's 1468 trials have FSCV in the dandiset. For the c8dg half of fig3a-dual-channel (blocked, as c8ds is missing): the mean over 0 to 3 s is above zero in 56 of 85 large-reward trials (65.9%).
fig4b-beta-reward¶
Beta-band LFP power at c3bs-c3a in the direction task is more suppressed after large-reward than after small-reward value cues, in a brief window after the cue.
def beta_contrast(maskA, maskB, nameA, nameB):
va, vb = d26.value_cue[maskA].values, d26.value_cue[maskB].values
for a, b, norm, deciding in ((0, 0.5, False, True), (0, 1, False, False), (0.25, 0.75, False, False),
(0, 0.5, True, False), (-0.5, 0, False, False)):
A_ = mp.window_mean(lfp_t, beta, va, a, b); B_ = mp.window_mean(lfp_t, beta, vb, a, b)
if norm:
A_ = A_ / mp.window_mean(lfp_t, beta, va, -0.5, 0); B_ = B_ / mp.window_mean(lfp_t, beta, vb, -0.5, 0)
scale = 1 if norm else 1e12
p = stats.mannwhitneyu(A_, B_).pvalue
print(f"{a} to {b} s{' / (-0.5 to 0 s)' if norm else ' (µV²)'}: {nameA} median {np.median(A_) * scale:.3f} (n {len(A_)}), "
f"{nameB} median {np.median(B_) * scale:.3f} (n {len(B_)}), ratio {np.median(A_) / np.median(B_):.3f}, "
f"Mann-Whitney p = {p:.2g}" + (" <- deciding" if deciding else ""))
for a, b in ((0, 0.5), (0.125, 0.5)):
A_ = mp.window_mean(lfp_t, beta_env, va, a, b); B_ = mp.window_mean(lfp_t, beta_env, vb, a, b)
print(f"unsmoothed envelope, {a} to {b} s (µV²): {nameA} median {np.median(A_) * 1e12:.3f}, {nameB} median "
f"{np.median(B_) * 1e12:.3f}, ratio {np.median(A_) / np.median(B_):.3f}, Mann-Whitney p = {stats.mannwhitneyu(A_, B_).pvalue:.2g}")
bins = []
for a in np.arange(-5, 10) / 10:
A_ = mp.window_mean(lfp_t, beta_env, va, a, a + 0.1); B_ = mp.window_mean(lfp_t, beta_env, vb, a, a + 0.1)
bins.append(f"{a:+.1f} s {np.median(A_) / np.median(B_):.2f} (p {stats.mannwhitneyu(A_, B_).pvalue:.1g})")
print(f"unsmoothed envelope, {nameA} / {nameB} median ratio in 0.1 s bins from the value cue:\n " + "; ".join(bins))
tt, PA = mp.aligned(lfp_t, beta, va, 3, 8); _, PB = mp.aligned(lfp_t, beta, vb, 3, 8)
fig, ax = plt.subplots(figsize=(7, 3))
for P, col, lab in ((PA, LARGE, nameA), (PB, SMALL, nameB)):
mu = P.mean(0) * 1e12; se = P.std(0) / np.sqrt(len(P)) * 1e12
ax.plot(tt, mu, color=col, lw=2, label=f"{lab} (n {len(P)})"); ax.fill_between(tt, mu - se, mu + se, color=col, alpha=0.2, lw=0)
ax.axvline(0, color="0.6", lw=1); ax.set_xlabel("time from value cue (s)"); ax.set_ylabel("beta power (µV²)")
ax.legend(frameon=False); ax.set_title("c3bs-c3a, mean +- SE"); plt.show()
print(f"plotted means range {min(PA.mean(0).min(), PB.mean(0).min()) * 1e12:.2f} to {max(PA.mean(0).max(), PB.mean(0).max()) * 1e12:.2f} µV²")
beta_contrast(d26.success & d26.big, d26.success & d26.small, "large", "small")
for side in ("left", "right"):
m_ = d26.success & d26[side]
print(f"{side} trials only:", end=" ")
A_ = mp.window_mean(lfp_t, beta, d26.value_cue[m_ & d26.big].values, 0, 0.5)
B_ = mp.window_mean(lfp_t, beta, d26.value_cue[m_ & d26.small].values, 0, 0.5)
print(f"large {np.median(A_) * 1e12:.3f}, small {np.median(B_) * 1e12:.3f} µV², p = {stats.mannwhitneyu(A_, B_).pvalue:.2g}")
0 to 0.5 s (µV²): large median 275.125 (n 570), small median 327.206 (n 602), ratio 0.841, Mann-Whitney p = 1e-08 <- deciding
0 to 1 s (µV²): large median 323.846 (n 570), small median 311.146 (n 602), ratio 1.041, Mann-Whitney p = 0.09
0.25 to 0.75 s (µV²): large median 213.249 (n 570), small median 213.613 (n 602), ratio 0.998, Mann-Whitney p = 0.66
0 to 0.5 s / (-0.5 to 0 s): large median 0.455 (n 570), small median 0.553 (n 602), ratio 0.822, Mann-Whitney p = 1.1e-06
-0.5 to 0 s (µV²): large median 613.142 (n 570), small median 609.639 (n 602), ratio 1.006, Mann-Whitney p = 0.82
unsmoothed envelope, 0 to 0.5 s (µV²): large median 268.885, small median 322.255, ratio 0.834, Mann-Whitney p = 1.9e-08
unsmoothed envelope, 0.125 to 0.5 s (µV²): large median 163.923, small median 244.299, ratio 0.671, Mann-Whitney p = 4.8e-19
unsmoothed envelope, large / small median ratio in 0.1 s bins from the value cue: -0.5 s 1.05 (p 0.2); -0.4 s 0.95 (p 0.5); -0.3 s 0.97 (p 0.6); -0.2 s 0.95 (p 0.9); -0.1 s 1.02 (p 1); +0.0 s 1.01 (p 1); +0.1 s 0.85 (p 0.001); +0.2 s 0.59 (p 1e-18); +0.3 s 0.60 (p 1e-24); +0.4 s 0.98 (p 0.6); +0.5 s 1.32 (p 5e-06); +0.6 s 1.39 (p 3e-07); +0.7 s 1.35 (p 5e-06); +0.8 s 1.49 (p 9e-11); +0.9 s 1.39 (p 7e-12)
plotted means range 101.93 to 682.38 µV² left trials only:
large 244.540, small 290.286 µV², p = 0.0016 right trials only:
large 301.329, small 394.225 µV², p = 1e-07
Result: reproduced. Over 0 to 0.5 s after the value cue, the median beta power is 275 µV² in 570 large-reward trials and 327 µV² in 602 small-reward trials (ratio 0.84, Mann-Whitney p = 1e-08), with no difference before the cue (-0.5 to 0 s, ratio 1.006, p = 0.82). Relative to each trial's pre-cue power the result holds (ratio 0.82, p = 1.1e-06), and it holds on each side (left p = 0.0016, right p = 1e-07). It is confined to a brief window: over 0 to 1 s large is slightly higher (ratio 1.04, p = 0.09), and over 0.25 to 0.75 s there is no difference (p = 0.66). The 0.25 s smoothing does not carry the difference from before the cue: on the unsmoothed envelope it holds over 0 to 0.5 s (ratio 0.834, p = 1.9e-08) and over 0.125 to 0.5 s (ratio 0.671, p = 4.8e-19), and in 0.1 s bins large and small do not differ before the cue (ratios 0.95 to 1.05, p 0.2 or above), large is lower from 0.1 to 0.4 s (ratios 0.85, 0.59, 0.60) and higher from 0.5 s on (1.32 to 1.49). Beta power is lower after large-reward cues at 0.1 to 0.4 s and higher from 0.5 s on. These bins fit a suppression that is deeper after large-reward cues, and equally an earlier suppression with an earlier rebound; they do not tell the two apart. The criterion, lower power over 0 to 0.5 s, is met.
fig4c-beta-side¶
With the reward fixed at large, beta-band LFP power is more suppressed when the cue is on the left (contralateral) than on the right (ipsilateral).
beta_contrast(d26.success & d26.big & d26.left, d26.success & d26.big & d26.right, "left", "right")
A_ = mp.window_mean(lfp_t, beta, d26.value_cue[d26.success & d26.small & d26.left].values, 0, 0.5)
B_ = mp.window_mean(lfp_t, beta, d26.value_cue[d26.success & d26.small & d26.right].values, 0, 0.5)
print(f"small reward, 0 to 0.5 s: left {np.median(A_) * 1e12:.3f}, right {np.median(B_) * 1e12:.3f} µV², p = {stats.mannwhitneyu(A_, B_).pvalue:.2g}")
0 to 0.5 s (µV²): left median 244.540 (n 282), right median 301.329 (n 288), ratio 0.812, Mann-Whitney p = 1.6e-06 <- deciding
0 to 1 s (µV²): left median 307.014 (n 282), right median 341.442 (n 288), ratio 0.899, Mann-Whitney p = 0.019
0.25 to 0.75 s (µV²): left median 214.378 (n 282), right median 210.208 (n 288), ratio 1.020, Mann-Whitney p = 0.75
0 to 0.5 s / (-0.5 to 0 s): left median 0.389 (n 282), right median 0.547 (n 288), ratio 0.710, Mann-Whitney p = 2.1e-09
-0.5 to 0 s (µV²): left median 654.502 (n 282), right median 597.658 (n 288), ratio 1.095, Mann-Whitney p = 0.013
unsmoothed envelope, 0 to 0.5 s (µV²): left median 238.983, right median 296.141, ratio 0.807, Mann-Whitney p = 1.7e-06
unsmoothed envelope, 0.125 to 0.5 s (µV²): left median 139.268, right median 202.483, ratio 0.688, Mann-Whitney p = 3.2e-17
unsmoothed envelope, left / right median ratio in 0.1 s bins from the value cue: -0.5 s 1.17 (p 0.06); -0.4 s 1.16 (p 0.05); -0.3 s 1.24 (p 0.03); -0.2 s 1.32 (p 0.1); -0.1 s 0.95 (p 1); +0.0 s 0.94 (p 0.6); +0.1 s 0.54 (p 5e-12); +0.2 s 0.55 (p 6e-15); +0.3 s 0.91 (p 0.06); +0.4 s 0.98 (p 0.4); +0.5 s 1.14 (p 0.2); +0.6 s 1.13 (p 0.3); +0.7 s 1.19 (p 0.04); +0.8 s 1.10 (p 0.1); +0.9 s 0.90 (p 0.3)
plotted means range 97.74 to 699.87 µV²
small reward, 0 to 0.5 s: left 290.286, right 394.225 µV², p = 1.5e-10
Result: reproduced. Among large-reward trials, the median beta power over 0 to 0.5 s after the value cue is 245 µV² for left (contralateral) cues and 301 µV² for right cues (282 and 288 trials, ratio 0.81, Mann-Whitney p = 1.6e-06). Before the cue, left is slightly higher (ratio 1.095, p = 0.013), so relative to the pre-cue power the difference is larger (ratio 0.71, p = 2.1e-09). Over 0 to 1 s the difference holds (p = 0.019); over 0.25 to 0.75 s it does not (p = 0.75). Small-reward trials show the same side difference (p = 1.5e-10). On the unsmoothed envelope left is lower over 0 to 0.5 s (ratio 0.807, p = 1.7e-06) and over 0.125 to 0.5 s (ratio 0.688, p = 3.2e-17), most in the bins from 0.1 to 0.3 s (ratios 0.54 and 0.55). The files do not record the implanted hemisphere, so left as contralateral rests on the Methods.
fig5bc-unit-isolation¶
The s-µIP recorded a spike-sorted unit in monkey P, shown by its average waveform and its interspike interval histogram.
isi = np.diff(units[1]) * 1e3
mwf = w1.mean(0); sdwf = w1.std(0)
t_raw, x_raw = mp.raw_slice(fP, 7000 * 30000, n10)
noise = np.median(np.abs(mp.spike_band(x_raw, 30000.0))) / 0.6745
print(f"unit 1: {len(units[1])} spikes over {units[1][-1] - units[1][0]:.0f} s ({len(units[1]) / (units[1][-1] - units[1][0]):.2f} Hz)")
print(f"ISIs under 1 ms: {np.sum(isi < 1)} of {len(isi)} ({100 * np.mean(isi < 1):.3f}%); under 2 ms {100 * np.mean(isi < 2):.3f}%; "
f"median ISI {np.median(isi):.1f} ms")
print(f"mean waveform trough {mwf.min():.2f} µV (SD at the trough {sdwf[mwf.argmin()]:.2f}), peak {mwf.max():.2f} µV; "
f"raw noise {noise:.2f} µV; trough / noise {abs(mwf.min()) / noise:.2f}")
for n_, s in zip(names, units):
i_ = np.diff(s) * 1e3
print(f" units-table unit {n_}: ISIs under 1 ms {100 * np.mean(i_ < 1):.2f}%")
fig, axs = plt.subplots(1, 2, figsize=(9, 3))
tw = (np.arange(48) - 8) / 30.0
axs[0].plot(tw, mwf, color=SMALL, lw=2); axs[0].fill_between(tw, mwf - sdwf, mwf + sdwf, color=SMALL, alpha=0.2, lw=0)
axs[0].set_xlabel("time from threshold crossing (ms)"); axs[0].set_ylabel("µV"); axs[0].set_title("unit 1, mean +- SD")
axs[1].hist(isi, bins=np.arange(0, 101, 1), color="0.4"); axs[1].set_xlabel("ISI (ms)"); axs[1].set_ylabel("count")
axs[1].set_title("ISI, 1 ms bins"); plt.tight_layout(); plt.show()
print(f"ISIs above 100 ms, outside the plotted range: {np.sum(isi > 100)}; waveform mean +- SD range {(mwf - sdwf).min():.1f} to {(mwf + sdwf).max():.1f} µV")
unit 1: 31247 spikes over 11984 s (2.61 Hz) ISIs under 1 ms: 1 of 31246 (0.003%); under 2 ms 0.032%; median ISI 90.3 ms mean waveform trough -30.64 µV (SD at the trough 5.23), peak 17.93 µV; raw noise 5.41 µV; trough / noise 5.66 units-table unit 0: ISIs under 1 ms 11.65% units-table unit 1: ISIs under 1 ms 0.00%
ISIs above 100 ms, outside the plotted range: 14698; waveform mean +- SD range -35.9 to 23.3 µV
Result: reproduced. Plexon unit 1 has 31247 spikes (2.61 Hz). One of its 31246 interspike intervals is under 1 ms (0.003%), and the trough of its mean waveform (-30.6 µV) is 5.66 times the noise level of the raw signal (5.41 µV). The units table's unit 0, the merged unsorted and poor-quality spikes, has 11.65% of intervals under 1 ms.
Summary¶
| Claim | Outcome |
|---|---|
insertion-yield |
Blocked: data missing |
chronic-duration |
Blocked: metadata insufficient |
fscv-waveform |
Partially reproduced (lowest potential -0.347 V, not -0.4 V) |
ephys-acquisition |
Reproduced |
shape-task-structure |
Reproduced |
zenodo-link (with outside data) |
Reproduced |
fig3a-dual-channel |
Blocked: data missing |
fig3a-redox-peaks |
Blocked: metadata insufficient |
fig3b-beta-bursts |
Reproduced |
fig4a-dopamine-reward |
Reproduced, weakly |
fig4b-beta-reward |
Reproduced |
fig4c-beta-side |
Reproduced |
fig5a-spike-reward |
Blocked: data missing |
fig5bc-unit-isolation |
Reproduced |
A failure to reproduce here may be a failure of this attempt, not of the paper: see the reasoning and the differences from the paper's method.
How the analysis differs from the paper's. The paper's dopamine and beta code is not public; the beta
filter and envelope are common choices, which match the paper's example trace at r = 0.94. Trial
conditions are decoded from event codes; in the shape task small reward is inferred from the absence of the
big-reward code, which the cue condition codes confirm. Dopamine is re-referenced at the value cue from the
stored PCA output, without the original currents. Only the first 734 trials of the shape-task session have
dopamine. The paper gives no statistics or windows for Figures 3 to 5; the tests and windows here were fixed
before those analyses ran. Outside data come only from the paper's Zenodo deposit: the trial lists and example
trace of zenodo-link, and Monkey P's events file, used in the data check for its clock offset.
What was found about the dandiset.
- Each event is stored about five times in the events tables, once per overlapping 60 s trial window.
- The trials tables hold no conditions; the shape task has no small-reward code, and the labels of codes 35 and 36 do not say they code reward size.
- The trial-aligned FSCV table covers only the first 734 of 1468 trials; its columns have no units; the stored
dais not referenced at the value cue. - The raw FSCV current is labelled amperes with values in the hundreds; its scan times are a regular grid from the first trial start; it cannot be matched to the trial-aligned series at any offset tried, or to the task at the stored times.
- The bipolar LFP's conversion factor is wrong (values are already in volts).
- The raw ephys of 09/26/2024 starts at 0 s but belongs 9831.35 s later on the trial clock.
- Monkey P's spike times are 6800.666 s later than its raw signal's clock; its units table counts the unsorted and poor-quality spikes as unit 0; the poor-quality spikes are not in time order; the file has no trials or events, so Figure 5A cannot be made (the Zenodo events of this session are TTL values without names).
- The second FSCV site of Figure 3A (c8ds) and the implant dates are not in the dandiset; Monkey P's site is
cl3in the file andc33in the paper; the files cite the preprint; the acquisition high cut is 8000 Hz, not the 7500 Hz the description states.
print(f"total run time {(time.time() - started) / 60:.1f} min")
total run time 2.4 min