Reproducing Fujisawa et al. (2008) from dandiset 000067¶
Paper. Fujisawa S, Amarasingham A, Harrison MT, Buzsáki G (2008). Behavior-dependent short-term assembly dynamics in the medial prefrontal cortex. Nature Neuroscience 11, 823-833. doi:10.1038/nn.2134. This attempt used the full text of the author manuscript in PubMed Central (PMC2562676): Results, figure legends and Methods. The supplementary information could not be retrieved (PubMed Central and the publisher returned a reCAPTCHA page), so measures it defines are fixed in the claims file instead. This is the first revision, after the verifier's round-1 report: the variants it asked for are added in each section, and traj-diverge was compared on the x coordinate under a criterion added before that comparison was made; after the third round of verification it is blocked, as the x coordinate cannot show a lateral separation.
Dandiset. DANDI:000067, version 0.250624.0419. The paper (2008) cites no DANDI version. There are two published versions, 0.210812.1457 and 0.250624.0419; the data check below shows that both, and the draft, hold the same 28 assets (same asset and blob ids). The later one differs only in metadata (contributor identifiers and affiliations, a fuller link to the paper). The most recent published version is pinned.
What the dandiset holds. 28 sessions of three rats (EE, FF, GG): spike times of every sorted cluster, the LFP, in most files the wideband signal, and in ten files a trial table and position. The units table has no channel, shank, brain region or cell type. The recordings combine a 64-site prefrontal probe and a 32-site CA1 probe. The paper's four rats and its example session (EE.188) are not all here.
Outside inputs. The channel map (which channels are the prefrontal probe, the CA1 probe and the nose-poke sensors) comes from the CRCNS pfc-2 data description (version 0.7, doi:10.6080/K01V5BWK), which documents the same recordings; the link claim tests it against the data. Claims marked "with labels" use it. The analysis itself runs on the dandiset's data only.
Randomness. Each claim's section draws from its own generator, seeded by a fixed number and the session's index.
| Claim | Data | Status |
|---|---|---|
| channel-region-link | with labels | tested |
| rats-and-units | with labels | tested |
| task-performance | dandiset | blocked: data missing (only correct trials are in the files) |
| fig2c-active-fraction | with labels | tested |
| fig2-field-size | with labels | tested |
| fig2b-trajectory-selective | with labels | tested |
| traj-diverge | dandiset | blocked: data missing (no y coordinate; the x coordinate is compared) |
| odor-nosepoke | with labels | tested |
| fig3-connection-rates | with labels | tested |
| fig3-cell-classes | with labels | tested |
| fig3f-distance | with labels | tested |
| fig4-position-dependence | with labels | tested; blocked: methods underspecified (the test does not hold its false-positive rate) |
| fig6-facilitation-depression | with labels | tested |
| fig7-supralinear | with labels | tested |
| fig3-trajectory-by-class | with labels | tested |
| ee188-example-session | dandiset | blocked: data missing (session EE.188 is not in the dandiset) |
import time
import warnings
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from concurrent.futures import ThreadPoolExecutor
from dandi.dandiapi import DandiAPIClient
from scipy import stats
import pfc
pd.set_option("display.width", 220)
pd.set_option("display.max_columns", 40)
pd.set_option("display.max_rows", 300)
warnings.filterwarnings("ignore", category=RuntimeWarning)
started = time.time()
urls = pfc.asset_urls()
with ThreadPoolExecutor(8) as pool:
files = pd.DataFrame(list(pool.map(lambda p: pfc.describe_file(p, urls[p]), sorted(urls))))
wm_paths = [p for p in sorted(urls) if any(f"ses-{s.replace('.', '-')}_" in p for s in pfc.WM_SESSIONS)]
sessions = pfc.load_sessions(wm_paths, workers=9)
by_name = {s["session"]: s for s in sessions}
analysed = [s for s in sessions if s["has_raw"]]
print(f"{len(files)} files described and {len(sessions)} working memory sessions loaded from dandiset "
f"{pfc.DANDISET_ID} version {pfc.DANDISET_VERSION} in {time.time() - started:.0f} s")
28 files described and 9 working memory sessions loaded from dandiset 000067 version 0.250624.0419 in 418 s
Data check¶
Each check of the description against the data is stated under its output.
with DandiAPIClient() as client:
versions = [v.identifier for v in client.get_dandiset(pfc.DANDISET_ID).get_versions()]
blobs = {}
for v in versions:
ds = client.get_dandiset(pfc.DANDISET_ID, v)
blobs[v] = sorted((a.path, a.identifier, a.size) for a in ds.get_assets())
print("Versions:", versions)
for v in versions:
print(f" {v}: {len(blobs[v])} assets, {sum(b[2] for b in blobs[v]) / 1e9:.1f} GB, "
f"identical to {pfc.DANDISET_VERSION}: {blobs[v] == blobs[pfc.DANDISET_VERSION]}")
show = files.drop(columns=["lfp_head", "spike_head", "path"])
print("\nEvery file (raw_s, lfp_s: duration of the wideband and LFP series in seconds):")
print(show.to_string(index=False))
print("\nSubjects:", files.subject.value_counts().to_dict())
print("Channel locations:", sorted({l for ls in files.locations for l in ls}),
"; electrode x, y, z all NaN in every file:", bool(files.xyz_all_nan.all()),
"; electrode groups per file:", sorted(files.groups.unique()))
print("Unit table columns:", sorted({c for cs in files.unit_columns for c in cs}))
dup_lfp = files[files.lfp_head.duplicated(keep=False) & (files.lfp_s > 2)]
dup_spk = files[(files.spike_head != "") & files.spike_head.duplicated(keep=False)]
print("\nFiles sharing the first second of LFP (channel 1):", dup_lfp.session.tolist() or "none")
print("Files sharing their first 50 spike times:", dup_spk.session.tolist() or "none")
Versions: ['draft', '0.210812.1457', '0.250624.0419']
draft: 28 assets, 94.6 GB, identical to 0.250624.0419: True
0.210812.1457: 28 assets, 94.6 GB, identical to 0.250624.0419: True
0.250624.0419: 28 assets, 94.6 GB, identical to 0.250624.0419: True
Every file (raw_s, lfp_s: duration of the wideband and LFP series in seconds):
session subject start channels locations groups xyz_all_nan raw_s lfp_s units unit_columns spikes trials position
EE.042 EE 2008-06-22 104 [unknown] 13 True 0.0000 993.0000 0 [] 0 0 False
EE.043 EE 2008-06-22 104 [unknown] 13 True 1.6384 1.6384 0 [] 0 0 False
EE.044 EE 2008-06-22 104 [unknown] 13 True 1.3000 1.3000 0 [] 0 0 False
EE.045 EE 2008-06-22 104 [unknown] 13 True 2340.4544 2340.4544 0 [] 0 0 False
EE.046 EE 2008-06-22 104 [unknown] 13 True 140.0832 140.0832 0 [] 0 0 False
EE.047 EE 2008-06-22 104 [unknown] 13 True 1210.3680 1210.3680 0 [] 0 0 False
EE.048 EE 2008-06-22 104 [unknown] 13 True 1035.6000 1035.6000 0 [] 0 0 False
EE.049 EE 2008-06-22 104 [unknown] 13 True 1788.8000 1788.8000 118 [spike_times] 3246543 19 True
EE.050 EE 2008-06-22 104 [unknown] 13 True 1876.4000 1876.4000 0 [] 0 0 False
EE.080 EE 2008-06-27 104 [unknown] 13 True 0.0000 875.3152 362 [spike_times] 567663 0 False
EE.081 EE 2008-06-27 104 [unknown] 13 True 0.0000 174.4896 362 [spike_times] 123764 0 False
EE.082 EE 2008-06-27 104 [unknown] 13 True 0.0000 1334.0000 362 [spike_times] 938551 0 False
EE.083 EE 2008-06-27 104 [unknown] 13 True 0.0000 1359.8720 358 [spike_times] 2709136 0 False
EE.084 EE 2008-06-27 104 [unknown] 13 True 0.0000 1704.7552 362 [spike_times] 1101699 0 False
EE.085 EE 2008-06-27 104 [unknown] 13 True 0.0000 689.3568 362 [spike_times] 434497 0 False
EE.086 EE 2008-06-27 104 [unknown] 13 True 0.0000 151.9616 356 [spike_times] 114859 0 False
EE.087 EE 2008-06-27 104 [unknown] 13 True 0.0000 448.6000 362 [spike_times] 292328 0 False
EE.088 EE 2008-06-27 104 [unknown] 13 True 0.0000 2165.1456 350 [spike_times] 4507778 31 True
EE.089 EE 2008-06-27 104 [unknown] 13 True 0.0000 2057.8304 362 [spike_times] 1504438 0 False
EE.196 EE 2008-07-10 104 [unknown] 13 True 2844.2624 2844.2624 161 [spike_times] 2771822 48 True
EE.210 EE 2008-07-11 104 [unknown] 13 True 2411.8000 2411.8000 170 [spike_times] 1680867 0 False
FF.103 FF 2008-11-14 104 [unknown] 13 True 2866.8000 2866.8000 78 [spike_times] 5662607 39 True
FF.130 FF 2008-11-16 104 [unknown] 13 True 2030.2000 2030.2000 90 [spike_times] 4340159 46 True
FF.167 FF 2008-11-19 104 [unknown] 13 True 1977.8000 1977.8000 64 [spike_times] 821839 43 True
GG.063 GG 2008-03-30 101 [unknown] 13 True 2850.3000 2850.3000 74 [spike_times] 3662835 28 True
GG.069 GG 2008-03-30 101 [unknown] 13 True 2375.8000 2375.8000 76 [spike_times] 4244461 40 True
GG.097 GG 2008-04-01 101 [unknown] 13 True 1750.0000 1750.0000 69 [spike_times] 526073 38 True
GG.149 GG 2008-04-06 101 [unknown] 13 True 2846.2000 2846.2000 63 [spike_times] 721786 31 True
Subjects: {'EE': 21, 'GG': 4, 'FF': 3}
Channel locations: ['unknown'] ; electrode x, y, z all NaN in every file: True ; electrode groups per file: [np.int64(13)]
Unit table columns: ['spike_times']
Files sharing the first second of LFP (channel 1): none
Files sharing their first 50 spike times: none
All versions hold the same assets, so the pinned version's content is the content of every version. The dandiset has three rats; the paper recorded four. Ten files have a trial table and position, and twenty have units. EE.042 and the ten sessions of 2008-06-27 (EE.080 to EE.089, EE.088 among them) have no wideband signal, and the test sessions EE.043 and EE.044 last under two seconds. Channel locations are "unknown", electrode coordinates are missing, and the 104 (EE, FF) or 101 (GG) channels are cut into 13 groups of 8 called "shank1" to "shank13", although the probes have 8 and 4 shanks and some channels are sensors. The units table holds spike times only. No two files share their first second of LFP or their first spike times, so no recording is deposited twice; the sessions of one day of rat EE were sorted together (the CRCNS description says so), so their units share numbers but not spikes.
The CRCNS description lists 13 working memory sessions. Of these, EE.152, EE.165 and EE.188 are not in the dandiset and EE.210 has no trial table, which leaves the nine loaded here. GG.063, a spontaneous alternation session in the CRCNS table, also has a trial table; it is not used.
rows = []
for s in sessions:
tr = s["trials"]
t = pfc.frame_times(s)
pos, trial = pfc.positions_at(s, t)
in_run = np.zeros(len(t), bool)
in_trial = np.zeros(len(t), bool)
first_d, last_d = [], []
for a0, a, b in zip(tr.nosepoke_start, tr.run_start, tr.run_end):
m = (t >= a) & (t <= b)
in_run |= m
in_trial |= (t >= a0) & (t <= b)
d = s["lin"][m]
if (d > 0).any():
first_d.append(d[d > 0][0])
last_d.append(d[d > 0][-1])
rows.append({
"session": s["session"], "trials": len(tr), "left": int((tr.side == "Left").sum()),
"right": int((tr.side == "Right").sum()),
"poke_s_median": np.median(tr.run_start - tr.nosepoke_start),
"run_s_median": np.median(tr.run_end - tr.run_start),
"reward_after_stop": bool((tr.run_end > tr.run_start).all()),
"first_d_in_run": np.median(first_d), "last_d_in_run": np.median(last_d),
"run_frames_d>0": np.mean(s["lin"][in_run] > 0), "frames_d>0_outside_trials": np.mean(s["lin"][~in_trial] > 0),
"lin_t0": s["lin_t0"], "lin_min": s["lin"].min(), "lin_max": s["lin"].max(), "lin_unit": s["lin_unit"],
"x1_range": (s["x1"].min(), s["x1"].max()), "x2_all_zero": bool(np.all(s["x2"] == 0)),
"pos_unit": s["pos_unit"], "duration_s": round(s["duration"], 1)})
print(pd.DataFrame(rows).round(3).to_string(index=False))
print("\nTrial column descriptions:", sessions[0]["trial_descriptions"])
session trials left right poke_s_median run_s_median reward_after_stop first_d_in_run last_d_in_run run_frames_d>0 frames_d>0_outside_trials lin_t0 lin_min lin_max lin_unit x1_range x2_all_zero pos_unit duration_s
EE.049 19 9 10 3.282 9.894 True 1.018 1.971 0.850 0.000 0.0 0.0 2.0 meters (-1.0, 564.0) True meters 1788.8
EE.088 31 11 20 2.657 5.005 True 1.047 1.890 0.986 0.003 0.0 0.0 2.0 meters (-1.0, 284.0) True meters 2165.1
EE.196 48 21 27 4.551 4.593 True 1.023 1.967 0.971 0.003 0.0 0.0 2.0 meters (-1.0, 592.2105) True meters 2844.3
FF.103 39 19 20 0.842 6.579 True 1.026 1.954 0.990 0.003 0.0 0.0 2.0 meters (-1.0, 605.4951) True meters 2866.8
FF.130 46 23 23 1.366 6.037 True 1.024 1.941 0.986 0.004 0.0 0.0 2.0 meters (-1.0, 610.0) True meters 2030.2
FF.167 43 21 22 4.385 5.869 True 1.027 1.970 0.981 0.004 0.0 0.0 2.0 meters (-1.0, 605.0) True meters 1977.8
GG.069 40 20 20 10.714 9.375 True 1.027 1.977 0.994 0.001 0.0 0.0 2.0 meters (-1.0, 611.0) True meters 2375.8
GG.097 38 20 18 9.694 6.040 True 1.024 1.984 0.994 0.005 0.0 0.0 2.0 meters (-1.0, 619.0) True meters 1750.0
GG.149 31 16 15 5.815 5.721 True 1.024 1.986 0.996 0.002 0.0 0.0 2.0 meters (-1.0, 619.9849) True meters 2846.2
Trial column descriptions: {'start_time': 'Start time of epoch, in seconds', 'stop_time': 'Stop time of epoch, in seconds', 'reward_time': 'Time when subject began consuming reward.', 'left_or_right': 'Time when subject began consuming reward.'}
The trial table's start_time, stop_time and reward_time are the three columns of the CRCNS
SessionNP matrix: the start of the first nose poke, the end of the last nose poke (the paper's
trial start, position 0) and the start of reward consumption (the paper's trial end, position 1).
The data agree: the first tracked linearized position after stop_time is 1.02 to 1.05 and the
last before reward_time is 1.89 to 1.99 (medians over trials), the position is above 0 in 85%
(EE.049) to 99.6% of run frames, and in at most 0.5% of frames outside the trials (from the first
nose poke to the reward). stop_time is therefore not the end of the trial, and the description of
left_or_right ("Time when subject began consuming reward") is wrong. The linearized position starts at time
lin_t0 (0 s in every file, so frame k is at k / 39.0625 s) and is
stated in metres but runs from 0 to 2 (0 outside trials), and the two "(x,y) coordinates" series
are both x coordinates, of the two head LEDs (the conversion code takes columns 1 and 3 of the
position file); they are stated in metres but are in pixels, and the second LED's is zero
throughout. There is no y coordinate; x1 runs along the maze (see traj-diverge). The
session durations match the CRCNS table (for example EE.049, 29.81 min). The trial table holds
correct trials only (CRCNS description), so the task performance cannot be computed.
units = pd.concat([pfc.assign_units(s) for s in analysed], ignore_index=True)
print("Units by session and class (PFC and CA1 split into single units and multi-unit clusters):")
cls = np.where(units.region.isin(["PFC", "CA1"]), units.region + np.where(units.mua, " multi", " single"), units.region)
print(pd.crosstab(units.session, cls, margins=True).to_string())
print("\nUnits whose spike times were not in order (sorted on loading):",
{s["session"]: s["unsorted_units"] for s in sessions if s["unsorted_units"]})
for s in analysed:
for i in s["unsorted_units"]:
st = s["trains"][i]
r = units[(units.session == s["session"]) & (units.unit == i)].iloc[0]
print(f" {s['session']} unit {i}: {len(st)} spikes, {r.rate_hz:.1f} Hz, ISIs < 2 ms {r.isi_lt2ms:.3f}, region {r.region}")
print("\nFraction of ISIs < 2 ms, all units of the eight sessions (bins):")
print(pd.cut(units.isi_lt2ms, [-0.001, 0.001, 0.005, 0.01, 0.02, 0.05, 0.1, 0.2, 0.5, 1]).value_counts(sort=False).to_string())
print("\nMulti-unit rule (ISIs) against the conversion's ordering rule (first unit of a shank block), probe units:")
probe = units[units.region.isin(["PFC", "CA1"])]
print(pd.crosstab(probe.mua.map({True: "ISI multi-unit", False: "ISI single"}),
probe.first_of_block.map({True: "first of block", False: "not first"})).to_string())
print("\nMedian rate (Hz) and amplitude (uV) by class:")
print(probe.groupby(["region", "mua"]).agg(n=("unit", "size"), rate_hz=("rate_hz", "median"),
amp_uv=("amp_uv", "median"), isi_lt2ms=("isi_lt2ms", "median")).round(3).to_string())
sens = units[units.region == "sensor"]
print("\nSensor clusters (largest on a sensor channel, under 200 uV on the probes):", len(sens),
"; their sensor-channel amplitude (uV) ranges", round(sens.aux_amp_uv.min()), "to", round(sens.aux_amp_uv.max()))
Units by session and class (PFC and CA1 split into single units and multi-unit clusters):
col_0 CA1 multi CA1 single PFC multi PFC single sensor unassigned All
session
EE.049 4 31 9 70 1 3 118
EE.196 4 15 9 132 1 0 161
FF.103 4 11 11 50 2 0 78
FF.130 4 9 10 67 0 0 90
FF.167 2 8 10 43 1 0 64
GG.069 7 41 9 18 1 0 76
GG.097 4 25 8 19 9 4 69
GG.149 4 26 8 18 6 1 63
All 33 166 74 417 21 8 719
Units whose spike times were not in order (sorted on loading): {'EE.049': [1], 'EE.088': [1], 'EE.196': [1], 'FF.103': [1], 'FF.130': [1], 'GG.069': [1]}
EE.049 unit 1: 200169 spikes, 111.9 Hz, ISIs < 2 ms 0.400, region PFC
EE.196 unit 1: 133213 spikes, 46.8 Hz, ISIs < 2 ms 0.493, region PFC
FF.103 unit 1: 362233 spikes, 126.4 Hz, ISIs < 2 ms 0.311, region PFC
FF.130 unit 1: 135295 spikes, 66.6 Hz, ISIs < 2 ms 0.393, region PFC
GG.069 unit 1: 197824 spikes, 83.3 Hz, ISIs < 2 ms 0.327, region PFC
Fraction of ISIs < 2 ms, all units of the eight sessions (bins):
isi_lt2ms
(-0.001, 0.001] 418
(0.001, 0.005] 174
(0.005, 0.01] 12
(0.01, 0.02] 2
(0.02, 0.05] 5
(0.05, 0.1] 5
(0.1, 0.2] 20
(0.2, 0.5] 75
(0.5, 1.0] 8
Multi-unit rule (ISIs) against the conversion's ordering rule (first unit of a shank block), probe units:
first_of_block first of block not first
mua
ISI multi-unit 96 11
ISI single 35 548
Median rate (Hz) and amplitude (uV) by class:
n rate_hz amp_uv isi_lt2ms
region mua
CA1 False 166 1.599 449.587 0.001
True 33 28.629 224.278 0.151
PFC False 417 1.047 538.150 0.000
True 74 55.622 237.975 0.363
Sensor clusters (largest on a sensor channel, under 200 uV on the probes): 21 ; their sensor-channel amplitude (uV) ranges 24 to 9155
Units are placed by their mean waveform. The fraction of inter-spike intervals under 2 ms is bimodal: 604 of 719 units have at most 1%, 103 have more than 10%, and 12 lie between. The units above 2% are taken as multi-unit clusters. Of the 107 such units on the probes, 96 are the first unit of a shank block, which is where the conversion's reader puts cluster 1 (multi-unit activity), and 35 first units are single units by the ISI rule, so the two rules mostly agree. Multi-unit clusters fire at a median of 29 (CA1) and 56 Hz (PFC) with about half the amplitude of single units. 21 units are clusters of the sensor signals (up to 9 mV on a sensor channel, under 200 uV on the probes). In six sessions the spike times of unit 1 are not in order, which breaks the format; they are sorted on loading, and in the analysed sessions these units are multi-unit clusters. From here on, PFC single units are the units placed on the prefrontal probe that are not multi-unit clusters.
rows = []
for s in analysed:
th = s["theta_ratio"]
sub = s["subject"]
ca1 = np.array(pfc.CA1_CHANNELS[sub])
rows.append({"session": s["session"], "theta_pfc_1_64": np.median(th[:64]),
"theta_crcns_ca1_65_96": np.median(th[list(pfc.CRCNS_CA1)]),
"theta_ca1_used": np.median(th[ca1]), "ca1_used": f"{ca1[0] + 1}-{ca1[-1] + 1}",
"theta_channels_65_69": np.round(th[64:69], 2).tolist()})
print("Theta ratio (6-10 Hz over 1-40 Hz power, raw signal in windows during trials), medians (1-based channels):")
print(pd.DataFrame(rows).round(3).to_string(index=False))
rows = []
for s in sessions:
crcns = pfc.sensor_check(s, pfc.CRCNS_NOSEPOKE)
used = pfc.sensor_check(s)
act = np.where(used.side == "Right", used.right_poke, used.left_poke) / np.where(used.side == "Right", used.right_run, used.left_run)
rows.append({"session": s["session"], "trials": len(used),
"agree_crcns_97_98": (crcns.side == crcns.sensor_side).mean(),
"agree_used": (used.side == used.sensor_side).mean(),
"used_channels_r_l": tuple(c + 1 for c in pfc.NOSEPOKE_CHANNELS[s["subject"]]),
"poke_over_run_activity": np.median(act)})
print("\nNose-poke sensors: fraction of trials whose more active sensor during the nose poke is the trial's side:")
print(pd.DataFrame(rows).round(3).to_string(index=False))
Theta ratio (6-10 Hz over 1-40 Hz power, raw signal in windows during trials), medians (1-based channels): session theta_pfc_1_64 theta_crcns_ca1_65_96 theta_ca1_used ca1_used theta_channels_65_69 EE.049 0.274 0.513 0.513 65-96 [0.5099999904632568, 0.5299999713897705, 0.5299999713897705, 0.5, 0.5400000214576721] EE.196 0.211 0.534 0.534 65-96 [0.5400000214576721, 0.5400000214576721, 0.550000011920929, 0.5600000023841858, 0.5699999928474426] FF.103 0.194 0.440 0.440 65-96 [0.4399999976158142, 0.4399999976158142, 0.44999998807907104, 0.44999998807907104, 0.44999998807907104] FF.130 0.202 0.543 0.543 65-96 [0.5600000023841858, 0.5600000023841858, 0.5699999928474426, 0.5600000023841858, 0.5600000023841858] FF.167 0.215 0.504 0.504 65-96 [0.5199999809265137, 0.5299999713897705, 0.5199999809265137, 0.5099999904632568, 0.5199999809265137] GG.069 0.176 0.376 0.391 70-101 [0.10000000149011612, 0.12999999523162842, 0.20999999344348907, 0.25, 0.05999999865889549] GG.097 0.163 0.365 0.367 70-101 [0.07999999821186066, 0.07000000029802322, 0.23000000417232513, 0.1899999976158142, 0.05999999865889549] GG.149 0.238 0.439 0.446 70-101 [0.09000000357627869, 0.10000000149011612, 0.2199999988079071, 0.1899999976158142, 0.05999999865889549]
Nose-poke sensors: fraction of trials whose more active sensor during the nose poke is the trial's side: session trials agree_crcns_97_98 agree_used used_channels_r_l poke_over_run_activity EE.049 19 1.000 1.000 (97, 98) 784.776978 EE.088 31 0.968 0.968 (97, 98) 647.422974 EE.196 48 1.000 1.000 (97, 98) 687.692993 FF.103 39 1.000 1.000 (97, 98) 1151.776978 FF.130 46 1.000 1.000 (97, 98) 908.862976 FF.167 43 1.000 1.000 (97, 98) 241.934998 GG.069 40 0.350 1.000 (66, 65) 84.031998 GG.097 38 0.421 1.000 (66, 65) 127.758003 GG.149 31 0.419 1.000 (66, 65) 153.852997
The CA1 channels have a theta ratio 1.9 to 2.7 times that of the prefrontal channels in every session (0.37 to 0.54 against 0.16 to 0.27), as a hippocampal probe should. In rats EE and FF the CRCNS map holds: CA1 on 65 to 96, the right and left sensors on 97 and 98 (agreement 0.97 to 1.00). In rat GG it does not: channels 65 to 69 carry little theta (0.06 to 0.25), channels 97 and 98 are not sensors (agreement 0.35 to 0.42), the CA1 probe is on 70 to 101, and the nose-poke sensors are on 65 and 66, with 66 active on trials labelled Right (agreement 1.00). The order of the two sensors, or the side labels, of rat GG is therefore reversed relative to rats EE and FF; the data cannot tell which. The left-right tests below do not depend on which side is called left. The prefrontal probe is channels 1 to 64 in all three rats.
rows = []
for s in analysed:
u = units[units.session == s["session"]]
ok = (u.sta_spikes >= 5) & u.region.isin(["PFC", "CA1"])
true_amp = u.amp_uv[ok]
chans = pfc.neural_channels(s["subject"])
shifted = s["sta_amp_shifted"][:, chans].max(1)[u.unit[ok].values]
rows.append({"session": s["session"], "units": int(ok.sum()), "amp_uv_median": np.median(true_amp),
"amp_uv_spikes_moved_3ms": np.median(shifted), "ratio": np.median(shifted / true_amp.values)})
print("Mean waveform amplitude at the spike times and with every spike moved 3 ms later:")
print(pd.DataFrame(rows).round(3).to_string(index=False))
edges = pd.concat([pfc.sensor_edges(s).assign(session=s["session"]) for s in sessions])
def crossing(df, col, rising):
v = df[col].values
half = (v.max() + v.min()) / 2
k = np.flatnonzero(v > half)
return df.lag_s.values[k[0]] if rising else df.lag_s.values[k[-1]]
print("\nTime (s) at which the trial's sensor activity crosses half its range, relative to the table's times:")
def at(d, col, lag):
return d[col].values[np.argmin(np.abs(d.lag_s.values - lag))]
print(edges.groupby("session").apply(lambda d: pd.Series({
"onset_vs_poke_start": crossing(d, "at_poke_start", True),
"offset_vs_poke_end": crossing(d, "at_poke_end", False),
"level_0.5s_before_start": at(d, "at_poke_start", -0.5), "level_0.5s_after_start": at(d, "at_poke_start", 0.5),
"level_0.5s_before_end": at(d, "at_poke_end", -0.5), "level_0.5s_after_end": at(d, "at_poke_end", 0.5)})).round(3).to_string())
Mean waveform amplitude at the spike times and with every spike moved 3 ms later: session units amp_uv_median amp_uv_spikes_moved_3ms ratio EE.049 107 501.360 28.398 0.071 EE.196 159 494.404 37.749 0.065 FF.103 76 474.253 24.132 0.052 FF.130 89 526.238 21.742 0.038 FF.167 63 530.358 23.771 0.047 GG.069 73 371.981 22.179 0.065 GG.097 54 493.294 48.180 0.086 GG.149 56 414.796 32.271 0.078
Time (s) at which the trial's sensor activity crosses half its range, relative to the table's times:
onset_vs_poke_start offset_vs_poke_end level_0.5s_before_start level_0.5s_after_start level_0.5s_before_end level_0.5s_after_end
session
EE.049 -0.001 -0.002 0.001 0.489 0.689 0.001
EE.088 -0.001 -0.002 0.001 0.393 0.692 0.001
EE.196 -0.001 -0.002 0.001 0.689 0.842 0.001
FF.103 -0.001 -0.002 0.002 0.407 0.456 0.001
FF.130 -0.001 -0.002 0.001 0.365 0.285 0.001
FF.167 -0.001 -0.002 0.002 0.365 0.389 0.001
GG.069 -0.001 -0.002 0.002 0.439 0.233 0.002
GG.097 -0.001 -0.002 0.005 0.374 0.170 0.001
GG.149 -0.001 -0.002 0.003 0.418 0.451 0.002
The spikes share the clock of the wideband signal: the mean waveform at the spike times is 12 to 26 times larger than with every spike moved 3 ms later. The trial table shares the clock of the recorded signals: in every session the trial's sensor turns on within 2 ms of the start of the first nose poke and off within 2 ms of the end of the last, and it is off 0.5 s before the start and 0.5 s after the end. The table's times were evidently taken from these sensor signals.
offsets = np.arange(-20, 8, 2.0)
dec = []
for s in analysed:
single = units[(units.session == s["session"]) & (units.region == "PFC") & ~units.mua].unit.tolist()
d = pfc.decode_side(s, single, offsets)
d["session"] = s["session"]
dec.append(d)
dec = pd.concat(dec)
table = dec.pivot(index="session", columns="offset_s", values="accuracy")
print("Leave-one-out decoding of the trial side from PFC single units, 1 s windows starting at run start + offset (s):")
print(table.round(2).to_string())
far = table[[c for c in table.columns if c <= -10]].mean(axis=1)
near = table[[c for c in table.columns if 0 <= c <= 4]].mean(axis=1)
print("\nMean accuracy, windows from -20 to -10 s:", far.round(2).to_dict())
print("Mean accuracy, windows from 0 to 4 s:", near.round(2).to_dict())
print("Sessions where 0-4 s exceeds -20 to -10 s:", int((near > far).sum()), "of", len(near))
Leave-one-out decoding of the trial side from PFC single units, 1 s windows starting at run start + offset (s):
offset_s -20.0 -18.0 -16.0 -14.0 -12.0 -10.0 -8.0 -6.0 -4.0 -2.0 0.0 2.0 4.0 6.0
session
EE.049 0.53 0.37 0.37 0.58 0.63 0.79 0.79 0.53 0.47 0.89 0.68 0.79 0.79 0.79
EE.196 0.77 0.67 0.60 0.79 0.54 0.75 0.69 0.81 0.85 0.92 0.96 1.00 0.96 0.92
FF.103 0.56 0.59 0.49 0.56 0.64 0.72 0.59 0.46 0.59 0.36 0.69 0.72 0.90 0.90
FF.130 0.63 0.59 0.61 0.74 0.41 0.48 0.57 0.70 0.76 0.89 0.91 0.96 0.91 0.96
FF.167 0.47 0.53 0.37 0.35 0.42 0.53 0.56 0.63 0.81 0.88 0.77 0.95 0.93 0.88
GG.069 0.52 0.45 0.50 0.52 0.62 0.65 0.55 0.45 0.48 0.55 0.55 0.75 0.70 0.75
GG.097 0.61 0.53 0.68 0.63 0.50 0.66 0.68 0.63 0.68 0.66 0.42 0.37 0.66 0.84
GG.149 0.39 0.52 0.52 0.52 0.58 0.48 0.55 0.39 0.58 0.52 0.55 0.58 0.77 0.87
Mean accuracy, windows from -20 to -10 s: {'EE.049': 0.54, 'EE.196': 0.69, 'FF.103': 0.59, 'FF.130': 0.58, 'FF.167': 0.45, 'GG.069': 0.55, 'GG.097': 0.6, 'GG.149': 0.5}
Mean accuracy, windows from 0 to 4 s: {'EE.049': 0.75, 'EE.196': 0.97, 'FF.103': 0.77, 'FF.130': 0.93, 'FF.167': 0.88, 'GG.069': 0.67, 'GG.097': 0.48, 'GG.149': 0.63}
Sessions where 0-4 s exceeds -20 to -10 s: 7 of 8
The side of the trial is decoded from PFC spikes better from 0 to 4 s after the run starts than 10 to 20 s before it, between trials, in 7 of 8 sessions (for example EE.196, 0.97 against 0.69). The z-scoring uses the training trials only. In GG.097 (19 units) the accuracy rises only at 6 s (0.84). With the sensor check above, this is consistent with spikes, trials and positions sharing one clock in every session; for GG.097 it is shown by the sensors and the waveforms, and only weakly by the spikes.
channel-region-link¶
The dandiset's recording channels can be assigned to the prefrontal probe, the CA1 probe and the nose-poke sensors by the channel map of the CRCNS pfc-2 description of the same recordings, and the data agree with that map.
rows = []
for s in sessions:
th = s.get("theta_ratio")
crcns = pfc.sensor_check(s, pfc.CRCNS_NOSEPOKE)
used = pfc.sensor_check(s)
rows.append({"session": s["session"], "rat": s["subject"],
"theta_check_crcns": None if th is None else bool(np.median(th[list(pfc.CRCNS_CA1)]) > np.median(th[:64])),
"theta_check_corrected": None if th is None else bool(np.median(th[list(pfc.CA1_CHANNELS[s["subject"]])]) > np.median(th[:64])),
"sensor_agree_crcns": (crcns.side == crcns.sensor_side).mean(),
"sensor_agree_corrected": (used.side == used.sensor_side).mean()})
link = pd.DataFrame(rows)
print(link.round(3).to_string(index=False))
crcns_ok = link.theta_check_crcns.dropna().all() and (link.sensor_agree_crcns >= 0.9).all()
corr_ok = link.theta_check_corrected.dropna().all() and (link.sensor_agree_corrected >= 0.9).all()
print("\nBoth checks hold in every session with the CRCNS numbers:", bool(crcns_ok))
print("Both checks hold in every session with the per-rat correction:", bool(corr_ok))
session rat theta_check_crcns theta_check_corrected sensor_agree_crcns sensor_agree_corrected EE.049 EE True True 1.000 1.000 EE.088 EE None None 0.968 0.968 EE.196 EE True True 1.000 1.000 FF.103 FF True True 1.000 1.000 FF.130 FF True True 1.000 1.000 FF.167 FF True True 1.000 1.000 GG.069 GG True True 0.350 1.000 GG.097 GG True True 0.421 1.000 GG.149 GG True True 0.419 1.000 Both checks hold in every session with the CRCNS numbers: False Both checks hold in every session with the per-rat correction: True
Result: partially reproduced. Channels 1 to 64 carry the prefrontal probe in all three rats, and in rats EE and FF the CRCNS map holds in every session: the CA1 channels have the higher theta ratio and the sensors agree with the trial side on 0.97 to 1.00 of trials. In rat GG the CRCNS sensor channels agree on 0.35 to 0.42 of trials. With the sensors on channels 65 and 66 (right on 66) and the CA1 probe on 70 to 101, both checks hold in every session (agreement 1.00).
rats-and-units¶
633 well-isolated units were recorded from the mPFC of four rats.
maps = {s["session"]: None for s in analysed}
peak = {}
for s in analysed:
single = units[(units.session == s["session"]) & (units.region == "PFC") & ~units.mua].unit.tolist()
m = pfc.rate_maps(s, single)
maps[s["session"]] = m
for u, mm in m.items():
peak[(s["session"], u)] = max(mm["Left"].max(), mm["Right"].max())
units["peak_hz"] = [peak.get((a, b), np.nan) for a, b in zip(units.session, units.unit)]
pfc_units = units[units.region == "PFC"]
counts = pd.DataFrame({
"a_single": pfc_units[~pfc_units.mua].groupby("session").size(),
"b_single_peak_1hz": pfc_units[~pfc_units.mua & (pfc_units.peak_hz >= 1)].groupby("session").size(),
"c_all_on_probe": pfc_units.groupby("session").size()}).fillna(0).astype(int)
counts.loc["total"] = counts.sum()
print(counts.to_string())
print("\nRats in the dandiset:", files.subject.nunique(), "(", sorted(files.subject.unique()), "); paper: 4")
print("Session EE.088 (no wideband signal) has", len(by_name["EE.088"]["trains"]), "units that cannot be placed by waveform.")
probe1 = pfc_units[~pfc_units.mua & ~pfc_units.first_of_block]
print("Variant (d): PFC single units, also leaving out the first unit of each shank block:", len(probe1))
# Variant: place units by the conversion's block order. Units are stored shank by shank with
# cluster 1 (multi-unit) first, so a block starts at each multi-unit cluster (by the ISI rule)
# that does not follow another; blocks 1 to 8 are taken as the prefrontal shanks.
def block_rule(isi_lt2ms):
mua = np.asarray(isi_lt2ms) > pfc.MAX_ISI_VIOLATION
start = mua & ~np.r_[False, mua[:-1]]
start[0] = True
block = np.cumsum(start) - 1
return block, mua
rows = []
for s in analysed:
u = units[units.session == s["session"]].sort_values("unit")
block, mua = block_rule(u.isi_lt2ms.values)
known = u.region.isin(["PFC", "CA1"]).values & (u.shank.values >= 0)
pred_pfc = block < 8
rows.append({"session": s["session"], "blocks": int(block.max() + 1), "units_placed": int(known.sum()),
"same_shank_as_waveform": np.mean(block[known] == u.shank.values[known]),
"pfc_agreement": np.mean(pred_pfc[known] == (u.region.values[known] == "PFC")),
"pfc_single_by_waveform": int(((u.region == "PFC") & ~u.mua).sum()),
"pfc_single_by_blocks": int((pred_pfc & ~mua).sum())})
print("\nBlock rule against the waveform placement:")
print(pd.DataFrame(rows).round(3).to_string(index=False))
s88 = by_name["EE.088"]
isi88 = [float(np.mean(np.diff(st) < 0.002)) if len(st) > 1 else 0.0 for st in s88["trains"]]
block, mua = block_rule(isi88)
n88 = int(((block < 8) & ~mua).sum())
print(f"EE.088 by the block rule: {int(block.max() + 1)} blocks, {n88} PFC single units; reading (a) with them: "
f"{int(counts.loc['total', 'a_single']) + n88}")
a_single b_single_peak_1hz c_all_on_probe session EE.049 70 52 79 EE.196 132 105 141 FF.103 50 41 61 FF.130 67 65 77 FF.167 43 41 53 GG.069 18 13 27 GG.097 19 17 27 GG.149 18 12 26 total 417 346 491 Rats in the dandiset: 3 ( ['EE', 'FF', 'GG'] ); paper: 4 Session EE.088 (no wideband signal) has 350 units that cannot be placed by waveform. Variant (d): PFC single units, also leaving out the first unit of each shank block: 416 Block rule against the waveform placement: session blocks units_placed same_shank_as_waveform pfc_agreement pfc_single_by_waveform pfc_single_by_blocks EE.049 13 114 1.000 1.000 70 73 EE.196 12 160 1.000 1.000 132 132 FF.103 10 76 0.250 0.816 50 60 FF.130 11 90 0.789 0.956 67 69 FF.167 10 63 0.016 0.889 43 46 GG.069 9 75 0.027 0.360 18 59 GG.097 9 56 0.018 0.536 19 42 GG.149 11 56 0.250 0.696 18 33 EE.088 by the block rule: 12 blocks, 214 PFC single units; reading (a) with them: 631
Result: not reproduced. The dandiset holds three rats, against the paper's four. The eight sessions hold 417 PFC single units (346 with a peak of at least 1 Hz, 491 counting multi-unit clusters), all outside 633 plus or minus 63. Leaving out also the first unit of each shank block gives 416. The 350 units of EE.088 cannot be placed by waveform. Placed by the conversion's block order, a rule that matches the waveform placement in both other EE sessions (agreement 1.000) but not in FF and GG (0.360 to 0.956), EE.088 adds 214 PFC single units, and reading (a) becomes 631, within tolerance; the number of rats still fails.
task-performance¶
Rats performed the odor-place matching task at 92% correct (91.9% in the recording sessions).
Blocked: the trial table holds correct trials only (CRCNS description of SessionNP). The cell below counts what the files hold outside the listed trials: visits to the far end of the central arm found from x1 (laps), and nose-poke bouts in the sensor channels. Neither says whether an unlisted lap was an error, an unrewarded lap or a return, so they give a bound, not the performance.
# Not a test of the claim (it stays blocked): bounds on unlisted trials, asked for by the verifier.
rows = []
for s in sessions:
lp = pfc.laps(s)
pb = pfc.poke_bouts(s)
rows.append({"session": s["session"], "listed_trials": len(s["trials"]),
"laps": len(lp), "laps_in_listed_runs": int(lp.in_run.sum()), "laps_outside": int((~lp.in_run).sum()),
"poke_bouts": len(pb), "bouts_in_listed_pokes": int(pb.listed.sum()), "bouts_outside": int((~pb.listed).sum())})
perf = pd.DataFrame(rows)
print(perf.to_string(index=False))
tot = perf[["listed_trials", "laps_outside", "bouts_outside"]].sum()
print(f"\nAll sessions: {tot.listed_trials} listed trials, {tot.laps_outside} laps and {tot.bouts_outside} nose-poke bouts outside them;"
f" if every lap outside a listed run were an error trial, performance would be {100 * tot.listed_trials / (tot.listed_trials + tot.laps_outside):.1f}%;"
f" if every unlisted nose-poke bout were one, {100 * tot.listed_trials / (tot.listed_trials + tot.bouts_outside):.1f}%")
session listed_trials laps laps_in_listed_runs laps_outside poke_bouts bouts_in_listed_pokes bouts_outside EE.049 19 32 19 13 34 26 8 EE.088 31 42 31 11 62 48 14 EE.196 48 68 48 20 97 80 17 FF.103 39 52 39 13 58 48 10 FF.130 46 52 46 6 58 53 5 FF.167 43 47 43 4 72 68 4 GG.069 40 48 40 8 136 125 11 GG.097 38 49 38 11 161 129 32 GG.149 31 45 31 14 109 63 46 All sessions: 335 listed trials, 100 laps and 147 nose-poke bouts outside them; if every lap outside a listed run were an error trial, performance would be 77.0%; if every unlisted nose-poke bout were one, 69.5%
Still blocked. Every listed run holds one lap, and the files hold 100 laps and 147 nose-poke bouts outside the 335 listed trials. If every lap outside a listed run were an error trial, performance would be 77.0%; if every unlisted nose-poke bout were one, 69.5%. Neither is a bound in one direction: the files do not say which laps or bouts were error trials, nor whether every error trial produced one, so the 91.9% can be neither confirmed nor excluded.
fig2c-active-fraction¶
About 10% (layer 2/3) and 20% (layer 5) of the recorded neurons fire at least one spike in any 100 ms window, a fraction that is relatively constant along the maze.
rows = []
for s in analysed:
single = units[(units.session == s["session"]) & (units.region == "PFC") & ~units.mua].unit.tolist()
af = pfc.active_fraction(s, single)
bins = pd.cut(af.pos, np.linspace(0, 1, 11), include_lowest=True)
bm = af.groupby(bins, observed=True).fraction.mean()
rows.append({"session": s["session"], "rat": s["subject"], "units": len(single), "windows": len(af),
"mean_pct": 100 * af.fraction.mean(), "bin_min_pct": 100 * bm.min(), "bin_max_pct": 100 * bm.max(),
"max_over_min": bm.max() / bm.min()})
act = pd.DataFrame(rows)
print(act.round(2).to_string(index=False))
per_rat = act.groupby("rat").mean_pct.median()
print("\nPer rat, median over sessions (%):", per_rat.round(2).to_dict())
layer = {r: ("2/3" if 9 <= v <= 11 else "5" if 18 <= v <= 22 else "none") for r, v in per_rat.items()}
print("Layer reading matched per rat:", layer)
part1 = all(v != "none" for v in layer.values()) and sum(v == "5" for v in layer.values()) <= 1
part2 = bool((act.max_over_min <= 2).all())
print("Every rat matches a layer value, at most one at 20%:", part1, "; every session's max/min over position bins <= 2:", part2)
session rat units windows mean_pct bin_min_pct bin_max_pct max_over_min
EE.049 EE 70 2346 9.67 8.45 11.12 1.32
EE.196 EE 132 2315 10.73 7.37 13.43 1.82
FF.103 FF 50 2655 14.25 12.94 15.05 1.16
FF.130 FF 67 2955 20.73 18.90 22.74 1.20
FF.167 FF 43 2650 19.68 17.87 21.75 1.22
GG.069 GG 18 3781 7.61 5.88 9.81 1.67
GG.097 GG 19 2302 5.75 4.38 8.24 1.88
GG.149 GG 18 1851 6.16 3.62 9.28 2.57
Per rat, median over sessions (%): {'EE': 10.2, 'FF': 19.68, 'GG': 6.16}
Layer reading matched per rat: {'EE': '2/3', 'FF': '5', 'GG': 'none'}
Every rat matches a layer value, at most one at 20%: False ; every session's max/min over position bins <= 2: False
Result: not reproduced. The per-rat medians are 10.2% (EE), 19.7% (FF) and 6.2% (GG): EE matches the layer 2/3 value and FF the layer 5 value, but GG matches neither. Along the track the fraction varies by a factor of 1.16 to 1.88 in seven sessions and 2.57 in GG.149, above the threshold of 2. Neither part of the criterion holds for all rats, though both hold for rats EE and FF.
fig2-field-size¶
Neurons are active over similar standardized distances of the maze, 0.27 (62 cm) plus or minus 0.17 (39 cm) (mean and s.d.), from the 50% boundaries of the peak firing rate.
rows = []
for s in analysed:
for u, mm in maps[s["session"]].items():
pk = {sd: mm[sd].max() for sd in mm}
if max(pk.values()) < pfc.MIN_PEAK_HZ:
continue
pref = max(pk, key=pk.get)
rows.append({"session": s["session"], "rat": s["subject"], "unit": u,
"a_contiguous_pref": pfc.field_size(mm[pref]), "c_all_bins_pref": pfc.field_size(mm[pref], False),
"b_left": pfc.field_size(mm["Left"]), "b_right": pfc.field_size(mm["Right"])})
fields = pd.DataFrame(rows)
b_vals = np.r_[fields.b_left, fields.b_right]
summary = pd.DataFrame({
"a": [fields.a_contiguous_pref.mean(), fields.a_contiguous_pref.std()],
"b": [np.nanmean(b_vals), np.nanstd(b_vals, ddof=1)],
"c": [fields.c_all_bins_pref.mean(), fields.c_all_bins_pref.std()]}, index=["mean", "sd"])
print("Units with a peak of at least 1 Hz:", len(fields), "of", int(counts.loc["total", "a_single"]), "PFC single units")
print(summary.round(3).to_string())
for k in summary:
m, sd = summary.loc["mean", k], summary.loc["sd", k]
print(f"reading {k}: mean within 0.27 +- 0.17: {abs(m - 0.27) <= 0.17}; s.d. within 0.17 +- 0.017: {abs(sd - 0.17) <= 0.017}")
own = np.r_[[max(mm["Left"].max(), 0) for s in analysed for u, mm in maps[s["session"]].items() if max(mm["Left"].max(), mm["Right"].max()) >= pfc.MIN_PEAK_HZ],
[max(mm["Right"].max(), 0) for s in analysed for u, mm in maps[s["session"]].items() if max(mm["Left"].max(), mm["Right"].max()) >= pfc.MIN_PEAK_HZ]]
b_own = b_vals[own >= pfc.MIN_PEAK_HZ]
print(f"Reading b: {np.isnan(b_vals).sum()} sides without spikes left out (NaN); variant with only sides whose own peak is at least 1 Hz: "
f"{len(b_own)} values, mean {np.nanmean(b_own):.3f}, s.d. {np.nanstd(b_own, ddof=1):.3f}")
print("\nBy rat, reading a:")
print(fields.groupby("rat").a_contiguous_pref.agg(["mean", "std", "size"]).round(3).to_string())
print("By session, reading a:")
print(fields.groupby("session").a_contiguous_pref.agg(["mean", "std", "size"]).round(3).to_string())
sub = []
for thr in (0.0, 0.5, 1.0, 2.0, 5.0):
vals = []
for s in analysed:
for u, mm in maps[s["session"]].items():
pk = {sd: mm[sd].max() for sd in mm}
if max(pk.values()) >= thr and max(pk.values()) > 0:
vals.append(pfc.field_size(mm[max(pk, key=pk.get)]))
sub.append({"min_peak_hz": thr, "units": len(vals), "mean": np.mean(vals), "sd": np.std(vals, ddof=1)})
print("\nReading a with other peak-rate thresholds:")
print(pd.DataFrame(sub).round(3).to_string(index=False))
fig, ax = plt.subplots(figsize=(5, 3))
ax.hist(fields.a_contiguous_pref, bins=np.linspace(0, 1, 21), color="0.4")
ax.set_xlabel("field size (normalized length), reading a")
ax.set_ylabel("units")
plt.show()
print("Plotted range:", fields.a_contiguous_pref.min(), "to", fields.a_contiguous_pref.max(), "(all within the axis)")
Units with a peak of at least 1 Hz: 346 of 417 PFC single units
a b c
mean 0.270 0.281 0.331
sd 0.196 0.206 0.220
reading a: mean within 0.27 +- 0.17: True; s.d. within 0.17 +- 0.017: False
reading b: mean within 0.27 +- 0.17: True; s.d. within 0.17 +- 0.017: False
reading c: mean within 0.27 +- 0.17: True; s.d. within 0.17 +- 0.017: False
Reading b: 0 sides without spikes left out (NaN); variant with only sides whose own peak is at least 1 Hz: 610 values, mean 0.292, s.d. 0.214
By rat, reading a:
mean std size
rat
EE 0.221 0.163 157
FF 0.322 0.223 147
GG 0.272 0.164 42
By session, reading a:
mean std size
session
EE.049 0.261 0.128 52
EE.196 0.202 0.175 105
FF.103 0.345 0.233 41
FF.130 0.331 0.236 65
FF.167 0.285 0.191 41
GG.069 0.346 0.171 13
GG.097 0.265 0.177 17
GG.149 0.202 0.106 12
Reading a with other peak-rate thresholds:
min_peak_hz units mean sd
0.0 417 0.254 0.189
0.5 387 0.263 0.192
1.0 346 0.270 0.196
2.0 263 0.287 0.212
5.0 139 0.330 0.243
Plotted range: 0.04 to 1.0 (all within the axis)
Result: partially reproduced. The mean field size is 0.270 over 346 units, the paper's value, but the s.d. is 0.196, outside 0.17 plus or minus 0.017. Readings b and c give means of 0.281 and 0.331 with s.d. of 0.206 and 0.220; no side without spikes entered reading b, and with only sides whose own peak is at least 1 Hz it gives 0.292 and 0.214. For peak thresholds from 0 to 5 Hz the s.d. is 0.189 to 0.243. Per rat the means are 0.221 to 0.322.
fig2b-trajectory-selective¶
The fraction of neurons with significantly different firing on left and right trials is 16% (layer 2/3) and 18% (layer 5) in the central arm (positions 0 to 0.3), where the trajectories overlap, and 22% (layer 2/3) and 40% (layer 5) in the side arms (0.5 to 1).
rows = []
for si, s in enumerate(analysed):
rng = np.random.default_rng([2008, 2, si])
is_left = (s["trials"].side == "Left").values
single = units[(units.session == s["session"]) & (units.region == "PFC") & ~units.mua].unit.tolist()
for u in single:
pos, tr = pfc.positions_at(s, s["trains"][u])
dens = pfc.trial_densities(pos, tr, len(is_left))
sig, _ = pfc.lr_test(dens, is_left, rng)
sig_sd, _ = pfc.lr_test(dens, is_left, rng, kind="sd")
sig0, _ = pfc.lr_test(dens, rng.permutation(is_left), rng)
rows.append({"session": s["session"], "rat": s["subject"], "unit": u, "peak_hz": peak[(s["session"], u)],
"central": sig[pfc.GRID <= 0.3].any(), "side": sig[pfc.GRID >= 0.5].any(), "any": sig.any(),
"central_sd": sig_sd[pfc.GRID <= 0.3].any(), "side_sd": sig_sd[pfc.GRID >= 0.5].any(),
"central_chance": sig0[pfc.GRID <= 0.3].any(), "side_chance": sig0[pfc.GRID >= 0.5].any(),
"any_chance": sig0.any()})
lr = pd.DataFrame(rows)
sel = lr[lr.peak_hz >= pfc.MIN_PEAK_HZ]
cols = ["central", "side", "central_chance", "side_chance", "central_sd", "side_sd"]
print("Units with a peak of at least 1 Hz; percent selective, by rat:")
print((100 * sel.groupby("rat")[cols].mean()).round(1).join(sel.groupby("rat").size().rename("units")).to_string())
print("\nBy session:")
print((100 * sel.groupby("session")[cols].mean()).round(1).join(sel.groupby("session").size().rename("units")).to_string())
pooled = 100 * sel[cols].mean()
print("\nPooled (%):", pooled.round(1).to_dict())
rat_ok = {}
for rat, g in sel.groupby("rat"):
c, sd = 100 * g.central.mean(), 100 * g.side.mean()
rat_ok[rat] = {"central_in_range": 14.4 <= c <= 19.8, "side_layer": "2/3" if 19.8 <= sd <= 24.2 else "5" if 36 <= sd <= 44 else "none"}
print("Per rat against the criterion:", rat_ok)
table = [[int(sel.central.sum()), int((~sel.central).sum())], [int(sel.central_chance.sum()), int((~sel.central_chance).sum())]]
print("Central arm, real against chance: Fisher one-sided p =", f"{stats.fisher_exact(table, alternative='greater')[1]:.3g}")
n10, n01 = int((sel.central & ~sel.central_chance).sum()), int((~sel.central & sel.central_chance).sum())
print(f"Variant, McNemar exact test (paired on units): {n10} units significant only with real labels, {n01} only relabelled, "
f"p = {stats.binomtest(n10, n10 + n01, 0.5, alternative='greater').pvalue:.3g}")
print("\nAll PFC single units (no peak threshold), percent:", (100 * lr[["central", "side", "any"]].mean()).round(1).to_dict(), "of", len(lr))
Units with a peak of at least 1 Hz; percent selective, by rat:
central side central_chance side_chance central_sd side_sd units
rat
EE 22.9 51.6 0.6 1.3 24.8 56.7 157
FF 24.5 59.9 1.4 0.7 23.8 63.3 147
GG 26.2 52.4 0.0 0.0 23.8 59.5 42
By session:
central side central_chance side_chance central_sd side_sd units
session
EE.049 1.9 34.6 0.0 3.8 13.5 34.6 52
EE.196 33.3 60.0 1.0 0.0 30.5 67.6 105
FF.103 7.3 29.3 0.0 0.0 7.3 34.1 41
FF.130 24.6 70.8 1.5 0.0 27.7 73.8 65
FF.167 41.5 73.2 2.4 2.4 34.1 75.6 41
GG.069 69.2 46.2 0.0 0.0 61.5 46.2 13
GG.097 5.9 41.2 0.0 0.0 5.9 58.8 17
GG.149 8.3 75.0 0.0 0.0 8.3 75.0 12
Pooled (%): {'central': 24.0, 'side': 55.2, 'central_chance': 0.9, 'side_chance': 0.9, 'central_sd': 24.3, 'side_sd': 59.8}
Per rat against the criterion: {'EE': {'central_in_range': np.False_, 'side_layer': 'none'}, 'FF': {'central_in_range': np.False_, 'side_layer': 'none'}, 'GG': {'central_in_range': np.False_, 'side_layer': 'none'}}
Central arm, real against chance: Fisher one-sided p = 7.34e-24
Variant, McNemar exact test (paired on units): 82 units significant only with real labels, 2 only relabelled, p = 1.85e-22
All PFC single units (no peak threshold), percent: {'central': 21.3, 'side': 48.2, 'any': 56.6} of 417
Result: partially reproduced. Units that fire differently on left and right trials are common in the central arm, where the paths overlap: 24.0% of 346 units pooled (22.9 to 26.2% per rat), against 0.9% with relabelled trials (Fisher p = 7.3e-24; McNemar, 82 against 2 discordant units, p = 1.9e-22). Every rat is above the paper's 16 to 18%, and in the side arms 51.6 to 59.9% of units per rat are selective, above both the 22% and the 40% of the paper. Per session the central-arm fraction ranges from 1.9% to 69.2%. The s.d. band gives 24.3% and 59.8%.
traj-diverge¶
Left and right trajectories overlap in the central arm and begin to differ significantly at position 0.3.
The two position series in the files are the x coordinates of the two LEDs (the second is zero throughout); there is no y coordinate. The verifier asked whether x1 separates the sides. The criterion (added in the first revision, before this comparison): per trial the mean x1 in 50 bins of 0.02, Welch t-test of left against right per bin, divergence point the first bin from which every bin up to 0.6 differs at p < 0.01.
bins = np.linspace(0, 1, 51)
rows, curves = [], {}
for s in sessions:
X = pfc.x_by_position(s)
L = (s["trials"].side == "Left").values
p = np.array([stats.ttest_ind(X[L, j], X[~L, j], equal_var=False, nan_policy="omit").pvalue for j in range(50)])
sig = p < 0.01
upto = np.flatnonzero(bins[:-1] < 0.6)
div = next((bins[j] for j in upto if sig[j:upto[-1] + 1].all()), np.nan)
ml, mr = np.nanmean(X[L], 0), np.nanmean(X[~L], 0)
curves[s["session"]] = (ml, mr)
rows.append({"session": s["session"], "divergence": div,
"sig_bins_central_0_0.3": f"{int(sig[:15].sum())}/15", "sig_bins_side_0.5_1": f"{int(sig[25:].sum())}/25",
"max_abs_diff_px_0_0.3": np.nanmax(np.abs(ml - mr)[:15]), "max_abs_diff_px_0.5_1": np.nanmax(np.abs(ml - mr)[25:]),
"x1_range_px": np.nanmax(np.r_[ml, mr]) - np.nanmin(np.r_[ml, mr])})
div = pd.DataFrame(rows)
print(div.round(3).to_string(index=False))
print("\nMean x1 (px) by side at positions 0, 0.1, ..., 0.9:")
for k, (ml, mr) in curves.items():
print(f" {k} Left ", np.round(ml[::5]).tolist(), "\n " + " " * len(k) + " Right", np.round(mr[::5]).tolist())
side_ok = (div["sig_bins_side_0.5_1"].str.split("/").str[0].astype(int) > 0).mean()
med = div.divergence.median()
print(f"\nSessions where x1 differs somewhere in 0.5 to 1: {side_ok:.2f}; sessions with a divergence point: {div.divergence.notna().sum()} of {len(div)};"
f" median divergence point {med:.2f} (criterion: 0.27 to 0.33)")
fig, ax = plt.subplots(figsize=(5, 3))
ml, mr = curves["FF.130"]
ax.plot(bins[:-1] + 0.01, ml, label="Left")
ax.plot(bins[:-1] + 0.01, mr, label="Right")
ax.set_xlabel("normalized position")
ax.set_ylabel("x1 (px), FF.130")
ax.legend()
plt.show()
print("Plotted range (px):", round(np.nanmin(np.r_[ml, mr])), "to", round(np.nanmax(np.r_[ml, mr])), "(all within the axes)")
session divergence sig_bins_central_0_0.3 sig_bins_side_0.5_1 max_abs_diff_px_0_0.3 max_abs_diff_px_0.5_1 x1_range_px
EE.049 NaN 11/15 18/25 12.356 9.566 402.051
EE.088 0.48 0/15 21/25 0.586 10.094 203.197
EE.196 0.50 1/15 24/25 1.199 40.443 413.954
FF.103 0.50 0/15 23/25 1.704 15.592 385.024
FF.130 0.48 11/15 22/25 5.177 48.694 385.640
FF.167 0.46 0/15 21/25 0.959 57.290 380.379
GG.069 0.58 14/15 21/25 12.008 6.794 390.922
GG.097 0.58 12/15 21/25 17.000 7.180 388.714
GG.149 0.06 12/15 24/25 7.200 14.338 387.734
Mean x1 (px) by side at positions 0, 0.1, ..., 0.9:
EE.049 Left [559.0, 472.0, 387.0, 301.0, 207.0, 163.0, 214.0, 300.0, 387.0, 471.0]
Right [558.0, 483.0, 397.0, 311.0, 226.0, 157.0, 204.0, 292.0, 378.0, 464.0]
EE.088 Left [nan, 248.0, 203.0, 157.0, 110.0, 79.0, 104.0, 150.0, 196.0, 240.0]
Right [nan, 248.0, 202.0, 157.0, 111.0, 82.0, 113.0, 159.0, 205.0, 249.0]
EE.196 Left [582.0, 503.0, 417.0, 330.0, 245.0, 170.0, 188.0, 273.0, 359.0, 445.0]
Right [581.0, 502.0, 417.0, 332.0, 246.0, 174.0, 226.0, 312.0, 399.0, 485.0]
FF.103 Left [571.0, 493.0, 408.0, 323.0, 236.0, 195.0, 255.0, 340.0, 424.0, 510.0]
Right [570.0, 493.0, 409.0, 323.0, 240.0, 186.0, 241.0, 325.0, 409.0, 494.0]
FF.130 Left [571.0, 498.0, 415.0, 333.0, 251.0, 213.0, 275.0, 358.0, 441.0, 523.0]
Right [576.0, 495.0, 413.0, 330.0, 248.0, 190.0, 226.0, 310.0, 393.0, 475.0]
FF.167 Left [575.0, 498.0, 416.0, 334.0, 249.0, 220.0, 291.0, 373.0, 455.0, 532.0]
Right [575.0, 498.0, 417.0, 335.0, 254.0, 195.0, 235.0, 317.0, 399.0, 481.0]
GG.069 Left [597.0, 514.0, 430.0, 347.0, 261.0, 211.0, 259.0, 342.0, 426.0, 510.0]
Right [599.0, 525.0, 441.0, 358.0, 275.0, 213.0, 255.0, 338.0, 421.0, 505.0]
GG.097 Left [597.0, 512.0, 430.0, 347.0, 266.0, 209.0, 253.0, 335.0, 417.0, 500.0]
Right [597.0, 529.0, 447.0, 364.0, 280.0, 215.0, 247.0, 328.0, 411.0, 493.0]
GG.149 Left [592.0, 514.0, 433.0, 350.0, 266.0, 204.0, 246.0, 326.0, 409.0, 492.0]
Right [591.0, 520.0, 439.0, 357.0, 275.0, 211.0, 259.0, 339.0, 422.0, 503.0]
Sessions where x1 differs somewhere in 0.5 to 1: 1.00; sessions with a divergence point: 8 of 9; median divergence point 0.49 (criterion: 0.27 to 0.33)
Plotted range (px): 190 to 576 (all within the axes)
# Analysis behind the proposed change of outcome (verifier round 2): where the x1 difference
# falls relative to the turn of the path in x. x1 falls along the central arm to a minimum and
# rises again; if the sides differ in x1 only from about that minimum on, x1 marks where the
# path turns back in x, not where the two paths separate laterally.
rows = []
for k, (ml, mr) in curves.items():
r = div[div.session == k].iloc[0]
j = int(round(0.3 / 0.02))
rows.append({"session": k, "x1_min_pos_left": bins[np.nanargmin(ml)] + 0.01, "x1_min_pos_right": bins[np.nanargmin(mr)] + 0.01,
"divergence": r.divergence, "diff_px_at_0.3": abs(ml[j] - mr[j]), "x1_drop_px_0_to_0.3": np.nanmean(np.r_[ml[0], mr[0]]) - np.nanmean(np.r_[ml[j], mr[j]])})
turn = pd.DataFrame(rows)
turn["divergence_minus_turn"] = turn.divergence - turn[["x1_min_pos_left", "x1_min_pos_right"]].mean(axis=1)
print(turn.round(3).to_string(index=False))
print(f"\nMedian position of the x1 minimum: {turn[['x1_min_pos_left', 'x1_min_pos_right']].stack().median():.2f};"
f" median divergence minus turn: {turn.divergence_minus_turn.median():+.2f} ({turn.divergence_minus_turn.notna().sum()} sessions);"
f" median left-right difference at 0.3, {turn['diff_px_at_0.3'].median():.1f} px, against a median fall of x1 from 0 to 0.3 of {turn['x1_drop_px_0_to_0.3'].median():.0f} px")
session x1_min_pos_left x1_min_pos_right divergence diff_px_at_0.3 x1_drop_px_0_to_0.3 divergence_minus_turn EE.049 0.49 0.51 NaN 10.581 252.127 NaN EE.088 0.49 0.49 0.48 0.189 NaN -0.01 EE.196 0.53 0.51 0.50 1.563 250.086 -0.02 FF.103 0.49 0.51 0.50 0.153 247.468 0.00 FF.130 0.47 0.51 0.48 2.808 241.996 -0.01 FF.167 0.47 0.51 0.46 0.325 240.818 -0.03 GG.069 0.49 0.51 0.58 11.627 245.539 0.08 GG.097 0.49 0.51 0.58 17.073 241.701 0.08 GG.149 0.51 0.51 0.06 7.103 238.042 -0.45 Median position of the x1 minimum: 0.51; median divergence minus turn: -0.01 (8 sessions); median left-right difference at 0.3, 2.8 px, against a median fall of x1 from 0 to 0.3 of 244 px
Result: blocked (data missing). x1 runs along the maze: it falls from about 570 px at the start to about 190 px near position 0.5 and rises again on the way to the reward, on both sides. It differs between left and right trials somewhere in 0.5 to 1 in all 9 sessions, by at most 7 to 57 px, against a range of about 390 px (200 px in EE.088). The median divergence point is 0.49 (8 of 9 sessions have one; 0.46 to 0.58 in seven, 0.06 in GG.149), above 0.33, which the written criterion would rate as partially reproduced. The central arm (0 to 0.3) also differs in 11 to 14 of 15 bins in EE.049, FF.130 and the three GG sessions, by at most 5 to 17 px, so by x1 alone the paths are not identical in the central arm in those sessions either. The cell above shows what the x1 difference measures: in every session with a divergence point but GG.149 it falls within 0.08 of the position where x1 reaches its minimum (median -0.01; the minimum is at 0.47 to 0.53), so x1 begins to differ where the path turns back in x, and at position 0.3 the sides differ by a median 2.8 px while x1 falls by a median 244 px from 0 to 0.3. A lateral separation of the paths would lie along the y coordinate, which is not in the dandiset, so x1 cannot test where the trajectories diverge, and the criterion's partial rating would follow from any x1 divergence past the turn. Whether the paths overlap in the central arm cannot be judged from x1 either. The verifiers of rounds 2 and 3 found the same, and the outcome is blocked (data missing).
odor-nosepoke¶
About one quarter of the neurons fire at significantly different rates (P < 0.05 per cell) during nose poking for the two odors.
rows = []
for si, s in enumerate(analysed):
rng = np.random.default_rng([2008, 8, si])
tr = s["trials"]
dur = (tr.run_start - tr.nosepoke_start).values
is_left = (tr.side == "Left").values
single = units[(units.session == s["session"]) & (units.region == "PFC") & ~units.mua].unit.tolist()
for u in single:
st = s["trains"][u]
c = np.searchsorted(st, tr.run_start.values) - np.searchsorted(st, tr.nosepoke_start.values)
if c.sum() < 10:
continue
r = c / dur
sh = rng.permutation(is_left)
rows.append({"session": s["session"], "rat": s["subject"], "unit": u,
"p": stats.mannwhitneyu(r[is_left], r[~is_left]).pvalue,
"p_t": stats.ttest_ind(r[is_left], r[~is_left]).pvalue,
"p_chance": stats.mannwhitneyu(r[sh], r[~sh]).pvalue})
odor = pd.DataFrame(rows)
n_single = int(counts.loc["total", "a_single"])
print("Units tested (at least 10 spikes in nose pokes):", len(odor), "of", n_single)
print("Nose-poke duration (s): median", round(float(np.median(np.concatenate([(s['trials'].run_start - s['trials'].nosepoke_start).values for s in analysed]))), 2))
frac = 100 * (odor.p < 0.05).mean()
chance = 100 * (odor.p_chance < 0.05).mean()
print(f"Mann-Whitney P < 0.05: {frac:.1f}% (chance, relabelled trials: {chance:.1f}%)")
print(f"Variant, t-test P < 0.05: {100 * (odor.p_t < 0.05).mean():.1f}%")
print("By rat (%):", (100 * odor.groupby("rat").p.apply(lambda p: (p < 0.05).mean())).round(1).to_dict())
print("By session (%):", (100 * odor.groupby("session").p.apply(lambda p: (p < 0.05).mean())).round(1).to_dict())
table = [[int((odor.p < 0.05).sum()), int((odor.p >= 0.05).sum())], [int((odor.p_chance < 0.05).sum()), int((odor.p_chance >= 0.05).sum())]]
print("Real against chance: Fisher one-sided p =", f"{stats.fisher_exact(table, alternative='greater')[1]:.3g}")
print("Within 20 to 30%:", 20 <= frac <= 30)
n10, n01 = int(((odor.p < 0.05) & (odor.p_chance >= 0.05)).sum()), int(((odor.p >= 0.05) & (odor.p_chance < 0.05)).sum())
print(f"Variant, McNemar exact test (paired on units): {n10} against {n01}, p = {stats.binomtest(n10, n10 + n01, 0.5, alternative='greater').pvalue:.3g}")
Units tested (at least 10 spikes in nose pokes): 397 of 417
Nose-poke duration (s): median 4.66
Mann-Whitney P < 0.05: 37.8% (chance, relabelled trials: 4.0%)
Variant, t-test P < 0.05: 28.5%
By rat (%): {'EE': 34.9, 'FF': 50.3, 'GG': 11.1}
By session (%): {'EE.049': 18.3, 'EE.196': 42.9, 'FF.103': 54.0, 'FF.130': 59.4, 'FF.167': 32.6, 'GG.069': 22.2, 'GG.097': 0.0, 'GG.149': 11.8}
Real against chance: Fisher one-sided p = 5.51e-35
Within 20 to 30%: False
Variant, McNemar exact test (paired on units): 142 against 8, p = 3.9e-33
Result: partially reproduced. 37.8% of 397 units fire at different rates during the nose poke for the two odors (Mann-Whitney P < 0.05), against 4.0% with relabelled trials (Fisher p = 5.5e-35; McNemar, 142 against 8 discordant units, p = 3.9e-33), so odor-dependent firing is common, but the fraction is above the 20 to 30% of the criterion. A t-test gives 28.5%. Per rat the fraction is 34.9% (EE), 50.3% (FF) and 11.1% (GG).
fig3-connection-rates¶
Of all ordered pairs of simultaneously recorded neurons, 0.79% had short-latency, narrow peaks or troughs in their cross-correlograms (0.55% excitatory, 0.24% inhibitory), found with a jitter test.
t0 = time.time()
conn, directed = {}, []
for si, s in enumerate(analysed):
rng = np.random.default_rng([2008, 9, si])
single = units[(units.session == s["session"]) & (units.region == "PFC") & ~units.mua].unit.tolist()
c = pfc.connections(s, single, rng, check_sample=100)
conn[s["session"]] = c
d = pfc.directed(c)
sh = units[units.session == s["session"]].set_index("unit").shank
d["session"], d["rat"] = s["session"], s["subject"]
d["dist_um"] = 200 * np.abs(sh[d.pre].values - sh[d.post].values)
directed.append(d)
directed = pd.concat(directed, ignore_index=True)
directed["conn"] = directed.exc | directed.inh
rows = []
for name, c in conn.items():
d = directed[directed.session == name]
chk = c[c.check]
rows.append({"session": name, "units": len(set(c.a) | set(c.b)), "ordered_pairs": len(d),
"screened_unordered": int(c[~c.check].screened.sum()), "exc": int(d.exc.sum()),
"exc_wider_than_2ms": int((d.exc_any & ~d.exc).sum()), "inh": int(d.inh.sum()),
"check_sample": len(chk),
"check_sample_significant": int(chk[["exc_ab", "exc_ba", "inh_ab", "inh_ba"]].fillna(False).astype(bool).any(axis=1).sum())})
per = pd.DataFrame(rows)
print(per.to_string(index=False))
n = len(directed)
rates = {"total": 100 * directed.conn.mean(), "excitatory": 100 * directed.exc.mean(), "inhibitory": 100 * directed.inh.mean()}
print(f"\n{n} ordered pairs: connections {rates['total']:.3f}%, excitatory {rates['excitatory']:.3f}%, inhibitory {rates['inhibitory']:.3f}%")
for k, (v, tol) in {"total": (0.79, 0.079), "excitatory": (0.55, 0.055), "inhibitory": (0.24, 0.024)}.items():
print(f" {k}: {rates[k]:.3f} against {v} +- {tol}: within {abs(rates[k] - v) <= tol}")
cs = set(map(tuple, directed.loc[directed.conn, ["session", "pre", "post"]].values))
recip = sum(1 for (se, a, b) in cs if (se, b, a) in cs) // 2
print(f"Reciprocally connected unordered pairs: {recip}, {100 * recip / (n / 2):.3f}% of unordered pairs (paper 0.17%)")
print("\nBy rat (%):")
print((100 * directed.groupby("rat")[["conn", "exc", "inh"]].mean()).round(3).join(directed.groupby("rat").size().rename("pairs")).to_string())
same = directed.dist_um == 0
print("\nSame shank (%):", (100 * directed[same][["exc", "inh"]].mean()).round(3).to_dict(),
"; different shanks (%):", (100 * directed[~same][["exc", "inh"]].mean()).round(3).to_dict())
print("Variant, peaks of any width counted as excitatory:", f"{100 * directed.exc_any.mean():.3f}%")
print(f"time {time.time() - t0:.0f} s")
session units ordered_pairs screened_unordered exc exc_wider_than_2ms inh check_sample check_sample_significant
EE.049 70 4830 86 29 0 2 100 0
EE.196 132 17292 430 156 7 47 100 0
FF.103 50 2450 81 26 0 4 100 0
FF.130 67 4422 160 51 6 17 100 0
FF.167 43 1806 81 14 0 4 100 0
GG.069 18 306 12 6 0 0 100 0
GG.097 19 342 4 2 0 0 100 0
GG.149 18 306 3 1 0 0 100 0
31754 ordered pairs: connections 1.131%, excitatory 0.898%, inhibitory 0.233%
total: 1.131 against 0.79 +- 0.079: within False
excitatory: 0.898 against 0.55 +- 0.055: within False
inhibitory: 0.233 against 0.24 +- 0.024: within True
Reciprocally connected unordered pairs: 40, 0.252% of unordered pairs (paper 0.17%)
By rat (%):
conn exc inh pairs
rat
EE 1.058 0.836 0.221 22122
FF 1.337 1.049 0.288 8678
GG 0.943 0.943 0.000 954
Same shank (%): {'exc': 4.014, 'inh': 0.588} ; different shanks (%): {'exc': 0.348, 'inh': 0.17}
Variant, peaks of any width counted as excitatory: 0.938%
time 57 s
# Troughs that mirror a peak in the other direction (verifier's point B1).
opp = directed.inh & directed.opposite_peak
print(f"Inhibitory connections: {int(directed.inh.sum())}; of them, with a peak in the other direction of the same pair: {int(opp.sum())}")
clean = {"total": 100 * (directed.exc | directed.inh_clean).mean(), "excitatory": 100 * directed.exc.mean(),
"inhibitory": 100 * directed.inh_clean.mean()}
print("Variant without those troughs (%):", {k: round(v, 3) for k, v in clean.items()},
"; within tolerance:", {k: abs(clean[k] - v) <= 0.1 * v for k, v in {"total": 0.79, "excitatory": 0.55, "inhibitory": 0.24}.items()})
cs2 = set(map(tuple, directed.loc[directed.exc | directed.inh_clean, ["session", "pre", "post"]].values))
recip2 = sum(1 for (se, a, b) in cs2 if (se, b, a) in cs2) // 2
print(f"Reciprocally connected unordered pairs without them: {recip2} ({100 * recip2 / (n / 2):.3f}%)")
# Synthetic pairs: a drives b with a 1 to 4 ms delay; no inhibition anywhere.
rng = np.random.default_rng([2008, 9, 99])
T, rows = 1800.0, []
for k in range(30):
ra, rb = rng.uniform(5, 20, 2)
a = np.sort(rng.uniform(0, T, rng.poisson(ra * T)))
m = rng.random(len(a)) < 0.05
b = np.sort(np.r_[rng.uniform(0, T, rng.poisson(rb * T)), a[m] + rng.uniform(0.001, 0.004, m.sum())])
syn = pfc.directed(pfc.connections({"trains": [a, b]}, [0, 1], rng))
rows.append({"exc_a_to_b": bool(syn.exc_any[(syn.pre == 0)].any()), "trough_b_to_a": bool(syn.inh[(syn.pre == 1)].any()),
"trough_b_to_a_clean": bool(syn.inh_clean[(syn.pre == 1)].any())})
syn = pd.DataFrame(rows)
print(f"\nSynthetic excitatory pairs (5 to 20 Hz, transmission 0.05): {int(syn.exc_a_to_b.sum())} of 30 detected; "
f"a trough in the other direction in {int(syn.trough_b_to_a.sum())} (criterion's rule), {int(syn.trough_b_to_a_clean.sum())} without mirrored troughs")
Inhibitory connections: 74; of them, with a peak in the other direction of the same pair: 34
Variant without those troughs (%): {'total': np.float64(1.023), 'excitatory': np.float64(0.898), 'inhibitory': np.float64(0.126)} ; within tolerance: {'total': np.False_, 'excitatory': np.False_, 'inhibitory': np.False_}
Reciprocally connected unordered pairs without them: 10 (0.063%)
Synthetic excitatory pairs (5 to 20 Hz, transmission 0.05): 30 of 30 detected; a trough in the other direction in 21 (criterion's rule), 0 without mirrored troughs
Result: partially reproduced. Of 31754 ordered pairs of PFC single units, 1.131% are connected, 0.898% excitatory and 0.233% inhibitory. The inhibitory rate matches the paper's 0.24%; the excitatory rate is 1.6 times the paper's 0.55%, so the total is above 0.79%. Excitatory connections are found in 4.01% of same-shank pairs and 0.35% of pairs on different shanks. None of the 800 unscreened pairs given the Monte Carlo test was significant. 40 unordered pairs (0.252%) are connected both ways, against the paper's 0.17%. The inhibitory match is fragile: 34 of the 74 inhibitory connections are troughs in the reverse direction of a pair with a peak, and on 30 synthetic excitatory pairs without any inhibition the criterion's rule finds such a trough in 21. Without those troughs the rates are 1.023%, 0.898% and 0.126%, none within tolerance, and 10 pairs (0.063%) are connected both ways.
fig3-cell-classes¶
By their connections, 32.5% of units are putative pyramidal cells and 12.5% putative interneurons, a ratio of 2.82 plus or minus 0.51 (s.d.), and 50.7% of the postsynaptic targets of putative excitatory cells themselves inhibit other neurons.
rows = []
for s in analysed:
d = directed[directed.session == s["session"]]
single = units[(units.session == s["session"]) & (units.region == "PFC") & ~units.mua].unit.tolist()
for u in single:
e, i = d[d.pre == u].exc.any(), d[d.pre == u].inh.any()
rows.append({"session": s["session"], "rat": s["subject"], "unit": u, "pyr": e and not i, "int": i and not e, "both": e and i})
cells = pd.DataFrame(rows)
print(f"PFC single units {len(cells)}: putative pyramidal {100 * cells.pyr.mean():.1f}%, putative interneurons "
f"{100 * cells.int.mean():.1f}%, both kinds of outgoing connection {int(cells.both.sum())} units")
g = cells.groupby("session")[["pyr", "int", "both"]].sum()
g["ratio"] = g.pyr / g.int.replace(0, np.nan)
print(g.to_string())
gr = cells.groupby("rat")[["pyr", "int"]].sum()
ratios = {"a_mean_of_sessions": g.ratio.mean(), "a_sd": g.ratio.std(), "b_pooled": cells.pyr.sum() / cells.int.sum(),
"c_mean_of_rats": (gr.pyr / gr.int.replace(0, np.nan)).mean()}
print("Ratios:", {k: round(v, 2) for k, v in ratios.items()}, "; per rat:", (gr.pyr / gr.int.replace(0, np.nan)).round(2).to_dict())
inter = set(map(tuple, cells[cells.int][["session", "unit"]].values))
targets = directed[directed.exc][["session", "post"]].drop_duplicates()
tfrac = 100 * np.mean([(a, b) in inter for a, b in targets.values])
inhib_any = set(map(tuple, directed[directed.inh][["session", "pre"]].values))
tfrac_any = 100 * np.mean([(a, b) in inhib_any for a, b in targets.values])
print(f"Targets of excitatory connections: {len(targets)} units, {tfrac:.1f}% putative interneurons "
f"({tfrac_any:.1f}% with any inhibitory outgoing connection)")
checks = {"pyramidal": abs(100 * cells.pyr.mean() - 32.5) <= 3.25, "interneurons": abs(100 * cells.int.mean() - 12.5) <= 1.25,
"ratio_a": abs(ratios["a_mean_of_sessions"] - 2.82) <= 0.51, "ratio_b": abs(ratios["b_pooled"] - 2.82) <= 0.51,
"ratio_c": abs(ratios["c_mean_of_rats"] - 2.82) <= 0.51, "targets": abs(tfrac - 50.7) <= 5.07}
print("Within tolerance:", checks)
print("\nInhibitory connections by the presynaptic unit's class:",
directed[directed.inh].merge(cells, left_on=["session", "pre"], right_on=["session", "unit"]).pipe(
lambda d: {"interneuron": int(d.int.sum()), "both": int(d.both.sum())}))
cl = []
for s in analysed:
d = directed[directed.session == s["session"]]
for u in units[(units.session == s["session"]) & (units.region == "PFC") & ~units.mua].unit:
e, i = d[d.pre == u].exc.any(), d[d.pre == u].inh_clean.any()
cl.append({"session": s["session"], "rat": s["subject"], "unit": u, "pyr": e and not i, "int": i and not e, "both": e and i})
cells_clean = pd.DataFrame(cl)
gc = cells_clean.groupby("session")[["pyr", "int"]].sum()
inter_c = set(map(tuple, cells_clean[cells_clean.int][["session", "unit"]].values))
print(f"\nVariant without troughs that mirror a peak: pyramidal {100 * cells_clean.pyr.mean():.1f}%, interneurons "
f"{100 * cells_clean.int.mean():.1f}% ({int(cells_clean.int.sum())} units), both {int(cells_clean.both.sum())}; "
f"pooled ratio {cells_clean.pyr.sum() / max(cells_clean.int.sum(), 1):.2f}; "
f"targets that are interneurons {100 * np.mean([(a, b) in inter_c for a, b in targets.values]):.1f}%")
print("Inhibitory connections on the same shank:", int((directed.inh & (directed.dist_um == 0)).sum()), "of", int(directed.inh.sum()))
PFC single units 417: putative pyramidal 37.9%, putative interneurons 3.1%, both kinds of outgoing connection 21 units
pyr int both ratio
session
EE.049 21 1 1 21.000000
EE.196 73 4 14 18.250000
FF.103 18 3 0 6.000000
FF.130 31 3 4 10.333333
FF.167 8 2 2 4.000000
GG.069 4 0 0 NaN
GG.097 2 0 0 NaN
GG.149 1 0 0 NaN
Ratios: {'a_mean_of_sessions': np.float64(11.92), 'a_sd': np.float64(7.46), 'b_pooled': np.float64(12.15), 'c_mean_of_rats': np.float64(12.96)} ; per rat: {'EE': 18.8, 'FF': 7.12, 'GG': nan}
Targets of excitatory connections: 124 units, 7.3% putative interneurons (20.2% with any inhibitory outgoing connection)
Within tolerance: {'pyramidal': np.False_, 'interneurons': np.False_, 'ratio_a': np.False_, 'ratio_b': np.False_, 'ratio_c': np.False_, 'targets': np.False_}
Inhibitory connections by the presynaptic unit's class: {'interneuron': 20, 'both': 54}
Variant without troughs that mirror a peak: pyramidal 38.6%, interneurons 2.4% (10 units), both 18; pooled ratio 16.10; targets that are interneurons 4.8% Inhibitory connections on the same shank: 28 of 74
Result: not reproduced. 37.9% of 417 units are putative pyramidal cells and 3.1% putative interneurons, against 32.5% and 12.5%, and 21 units have outgoing connections of both kinds. The pyramidal to interneuron ratio is 11.92 plus or minus 7.46 over sessions (12.15 pooled), against 2.82 plus or minus 0.51. Of the 124 targets of excitatory connections, 7.3% are putative interneurons (20.2% have some inhibitory connection), against 50.7%. 54 of the 74 inhibitory connections come from units that also excite other units. Without the troughs that mirror a peak, 38.6% are putative pyramidal cells and 2.4% (10 units) interneurons, a pooled ratio of 16.10, and 4.8% of targets are interneurons, further from the paper.
fig3f-distance¶
Most connected pairs were recorded on the same shank, and the probability of a connection falls with the distance between the neurons, though connections were found up to 1200 micrometres apart.
same = directed.dist_um == 0
tab = [[int(directed.conn[same].sum()), int((~directed.conn[same]).sum())],
[int(directed.conn[~same].sum()), int((~directed.conn[~same]).sum())]]
odds, p_f = stats.fisher_exact(tab)
prob = directed.groupby("dist_um").conn.agg(["mean", "sum", "size"])
prob["mean"] *= 100
# Check that the order of channel groups is the order of the shanks: a unit's second-largest
# shank amplitude should most often be on a neighbouring group if it is.
second = []
for s in analysed:
u = units[(units.session == s["session"]) & (units.region == "PFC") & ~units.mua & (units.sta_spikes >= 5)]
amp = s["sta_amp"][u.unit.values][:, :64].reshape(len(u), 8, 8).max(2)
order = np.argsort(-amp, axis=1)
second.append(np.abs(order[:, 1] - order[:, 0]))
second = np.concatenate(second)
expected = np.mean([np.mean(np.abs(np.delete(np.arange(8), k) - k) == 1) for k in range(8)])
print(f"Shank order check: second-largest shank is a neighbour of the largest for {100 * np.mean(second == 1):.1f}% of "
f"{len(second)} units (random order: {100 * expected:.1f}%); distribution of |difference|:",
pd.Series(second).value_counts().sort_index().to_dict())
print("Connection probability (%) by distance:")
print(prob.round(3).to_string())
rho, p_s = stats.spearmanr(prob.index, prob["mean"])
print(f"\nSame shank {100 * directed.conn[same].mean():.2f}% against different shanks {100 * directed.conn[~same].mean():.2f}%: "
f"Fisher p = {p_f:.3g}")
print(f"Spearman correlation over distances 0 to 1400 um: rho = {rho:.3f}, p = {p_s:.3g}")
print("Largest distance of a connection (um):", int(directed[directed.conn].dist_um.max()),
"; connections at 1000 um or more:", int((directed.conn & (directed.dist_um >= 1000)).sum()))
print(f"Share of connections on the same shank: {100 * (directed.conn & same).sum() / directed.conn.sum():.1f}%")
for sub in ("exc", "inh"):
pr = directed.groupby("dist_um")[sub].mean()
r, p = stats.spearmanr(pr.index, pr.values)
print(f"Variant, {sub} only: rho = {r:.3f}, p = {p:.3g}")
conn_d = directed.dist_um[directed.conn]
unconn_d = directed.dist_um[~directed.conn]
print(f"Variant, pair-level Mann-Whitney test that connected pairs are closer: p = {stats.mannwhitneyu(conn_d, unconn_d, alternative='less').pvalue:.3g}")
pr = prob[prob.index <= 1200]
r, p = stats.spearmanr(pr.index, pr["mean"])
print(f"Variant, distances 0 to 1200 um only: rho = {r:.3f}, p = {p:.3g}")
Shank order check: second-largest shank is a neighbour of the largest for 47.7% of 407 units (random order: 25.0%); distribution of |difference|: {1: 194, 2: 102, 3: 39, 4: 35, 5: 17, 6: 14, 7: 6}
Connection probability (%) by distance:
mean sum size
dist_um
0 4.603 219 4758
200 1.159 93 8022
400 0.316 19 6010
600 0.239 12 5024
800 0.193 8 4152
1000 0.168 4 2374
1200 0.181 2 1104
1400 0.645 2 310
Same shank 4.60% against different shanks 0.52%: Fisher p = 1.13e-89
Spearman correlation over distances 0 to 1400 um: rho = -0.619, p = 0.102
Largest distance of a connection (um): 1400 ; connections at 1000 um or more: 8
Share of connections on the same shank: 61.0%
Variant, exc only: rho = -0.571, p = 0.139
Variant, inh only: rho = -0.405, p = 0.32
Variant, pair-level Mann-Whitney test that connected pairs are closer: p = 2.27e-87
Variant, distances 0 to 1200 um only: rho = -0.964, p = 0.000454
Result: not reproduced. Connections are far more likely on the same shank (4.60% of pairs against 0.52%, Fisher p = 1.1e-89; 61.0% of all connections), and 8 connections are 1000 micrometres or more apart (the largest 1400). The Spearman correlation over all distances is -0.619 with p = 0.102, which fails the criterion: the 1400 micrometre bin (2 of 310 pairs, 0.645%) is above the bins from 600 to 1200. Over 0 to 1200 micrometres the correlation is -0.964 (p = 0.00045), and a test over pairs finds connected pairs closer than unconnected ones (Mann-Whitney p = 2.3e-87). The fall with distance is present; the outcome comes from the eight-point test fixed beforehand. The distance assumes that the order of the channel groups is the order of the shanks; the second-largest shank of a unit is a neighbour of its largest for 47.7% of 407 units, against 25.0% for a random order, which is consistent with that.
fig4-position-dependence¶
Out of 343 excitatory monosynaptic pairs, 67 showed position dependence of their short-latency coincidences in the maze (P < 0.01), beyond what the firing rates explain.
The Methods' rank-based global band can rarely reach a one-sided level of 0.01 here: with 1000 surrogates and a smooth statistic on 101 points, more than 1% of the data sets reach the most extreme rank somewhere on the track for most pairs, and then no rank band reaches the level and nothing can be significant (the cell below counts these pairs). The main test uses bands of the surrogate mean plus c times its s.d., with c searched as the Methods search the band; the rank version is shown as a variant.
t0 = time.time()
exc = directed[directed.exc].reset_index(drop=True)
rows = []
for k, r in exc.iterrows():
s = by_name[r.session]
rng = np.random.default_rng([2008, 12, k])
is_left = (s["trials"].side == "Left").values
both = pfc.coincidence_test(s, r.pre, r.post, rng, kinds=("sd", "rank", "efficacy"))
left = pfc.coincidence_test(s, r.pre, r.post, rng, is_left)
right = pfc.coincidence_test(s, r.pre, r.post, rng, ~is_left)
rows.append({"session": r.session, "pre": r.pre, "post": r.post, "coincidences": both["n_coincidences"],
"a_sig": both["sd_sig"], "a_fraction": both["sd_fraction"],
"b_sig": left["sd_sig"] or right["sd_sig"],
"a_rank": both["rank_sig"], "b_rank": left["rank_sig"] or right["rank_sig"],
"rank_level": both["rank_level"], "a_efficacy": both["efficacy_sig"]})
pos_dep = pd.DataFrame(rows)
print("Excitatory pairs:", len(pos_dep), "; coincidences in runs per pair: median", pos_dep.coincidences.median(),
", range", pos_dep.coincidences.min(), "to", pos_dep.coincidences.max())
for kcol, name in (("a_sig", "a, left and right together"), ("b_sig", "b, left or right"),
("a_rank", "variant, rank band, together"), ("b_rank", "variant, rank band, left or right")):
print(f"{name}: {int(pos_dep[kcol].sum())} pairs, {100 * pos_dep[kcol].mean():.1f}%")
print(f"Rank band: level 0 (no band reaches 0.01) in {int((pos_dep.rank_level == 0).sum())} of {len(pos_dep)} pairs")
print(f"Variant, null of constant efficacy along the track: {int(pos_dep.a_efficacy.sum())} pairs, {100 * pos_dep.a_efficacy.mean():.1f}%")
print("Within 19.5 +- 1.95%:", {k: abs(100 * pos_dep[k].mean() - 19.5) <= 1.95 for k in ("a_sig", "b_sig")})
print("Fraction of the track significant, pairs significant under a:", pos_dep.a_fraction[pos_dep.a_sig].round(3).tolist())
print("By session (% under a):", (100 * pos_dep.groupby("session").a_sig.mean()).round(1).to_dict())
for lo in (50, 100):
sub = pos_dep[pos_dep.coincidences >= lo]
print(f"Pairs with at least {lo} coincidences: {len(sub)}, significant under a: {100 * sub.a_sig.mean() if len(sub) else float('nan'):.1f}%")
print(f"time {time.time() - t0:.0f} s")
Excitatory pairs: 285 ; coincidences in runs per pair: median 7.0 , range 0 to 442
a, left and right together: 10 pairs, 3.5%
b, left or right: 5 pairs, 1.8%
variant, rank band, together: 4 pairs, 1.4%
variant, rank band, left or right: 6 pairs, 2.1%
Rank band: level 0 (no band reaches 0.01) in 183 of 285 pairs
Variant, null of constant efficacy along the track: 0 pairs, 0.0%
Within 19.5 +- 1.95%: {'a_sig': np.False_, 'b_sig': np.False_}
Fraction of the track significant, pairs significant under a: [0.03, 0.089, 0.158, 0.535, 0.05, 0.109, 0.079, 0.04, 0.416, 0.05]
By session (% under a): {'EE.049': 3.4, 'EE.196': 3.2, 'FF.103': 0.0, 'FF.130': 7.8, 'FF.167': 0.0, 'GG.069': 0.0, 'GG.097': 0.0, 'GG.149': 0.0}
Pairs with at least 50 coincidences: 28, significant under a: 21.4%
Pairs with at least 100 coincidences: 13, significant under a: 7.7%
time 20 s
# Calibration (verifier's point B3): synthetic posts driven by real presynaptic trains with a
# constant efficacy everywhere, so no position dependence of efficacy; and, as a positive control,
# posts driven only on one part of the track.
rng = np.random.default_rng([2008, 12, 999])
rows = []
for k, r in exc.sample(60, random_state=1).reset_index(drop=True).iterrows():
s = by_name[r.session]
a = s["trains"][r.pre]
rate_b = len(s["trains"][r.post]) / s["duration"]
pa, _ = pfc.positions_at(s, a)
for eff in (0.05, 0.2, "0.2 at 0.5 to 0.8 only"):
if isinstance(eff, str): # positive control: efficacy limited to one part of the track
m = (rng.random(len(a)) < 0.2) & (pa >= 0.5) & (pa <= 0.8)
else:
m = rng.random(len(a)) < eff
b = np.sort(np.r_[rng.uniform(0, s["duration"], rng.poisson(rate_b * s["duration"])), a[m] + rng.uniform(0.001, 0.004, m.sum())])
fake = {**s, "trains": [a, b]}
t = pfc.coincidence_test(fake, 0, 1, rng, kinds=("sd", "efficacy"))
rows.append({"efficacy": eff, "coincidences": t["n_coincidences"], "sd_sig": t["sd_sig"], "efficacy_sig": t["efficacy_sig"]})
cal = pd.DataFrame(rows)
print("Synthetic pairs, fraction flagged as position dependent:")
for eff, g in cal.groupby(cal.efficacy.astype(str), sort=False):
big = g[g.coincidences >= 50]
print(f" efficacy {eff}: {len(g)} pairs, median {g.coincidences.median():.0f} coincidences; criterion's test {100 * g.sd_sig.mean():.1f}%"
f" ({100 * big.sd_sig.mean():.1f}% of {len(big)} with at least 50); constant-efficacy null {100 * g.efficacy_sig.mean():.1f}%"
f" ({100 * big.efficacy_sig.mean():.1f}%)")
Synthetic pairs, fraction flagged as position dependent: efficacy 0.05: 60 pairs, median 17 coincidences; criterion's test 30.0% (75.0% of 16 with at least 50); constant-efficacy null 0.0% (0.0%) efficacy 0.2: 60 pairs, median 54 coincidences; criterion's test 75.0% (100.0% of 31 with at least 50); constant-efficacy null 0.0% (0.0%) efficacy 0.2 at 0.5 to 0.8 only: 60 pairs, median 18 coincidences; criterion's test 51.7% (100.0% of 16 with at least 50); constant-efficacy null 6.7% (25.0%)
Result: blocked (methods underspecified). Of 285 excitatory pairs, 10 (3.5%) show significant excess coincidences with left and right trials together and 5 (1.8%) with each side tested alone, against the paper's 19.5% (67 of 343). The criterion's test, however, compares the coincidences with jittered trains, so it detects any excess of coincidences somewhere on the track, not a change of efficacy along it: of synthetic pairs with a constant efficacy, it flags 30.0% (efficacy 0.05) and 75.0% (efficacy 0.2), against the paper's P < 0.01. The procedure here is therefore not the paper's, whose thinning step and null are in the supplement. With a null of constant efficacy, as a variant, no real pair is significant (0 of 285); that test flags none of the constant-efficacy synthetic pairs, but only 6.7% of the pairs driven only on part of the track (25.0% of the 16 with at least 50 coincidences), so its power at these spike counts is low. The rank band reaches the 0.01 level in only 102 of 285 pairs. The pairs have a median of 7 coincidences during runs (0 to 442).
fig6-facilitation-depression¶
About equal fractions of connected pairs show firing pattern-dependent depression (12.7%) and facilitation (10.7%), comparing the transmission of the first spike of a train with that of later spikes.
t0 = time.time()
rows = []
for k, r in exc.iterrows():
s = by_name[r.session]
rng = np.random.default_rng([2008, 13, k])
rows.append({"session": r.session, "pre": r.pre, "post": r.post,
**pfc.first_vs_later(s["trains"][r.pre], s["trains"][r.post], rng)})
fl = pd.DataFrame(rows)
ok = fl[(fl.n_first >= 20) & (fl.n_later >= 20)].copy()
ok["depressing"] = (ok.p < 0.10) & (ok["diff"] > 0)
ok["facilitating"] = (ok.p < 0.10) & (ok["diff"] < 0)
dep, fac = 100 * ok.depressing.mean(), 100 * ok.facilitating.mean()
print(f"Excitatory pairs {len(fl)}, with at least 20 first and 20 later spikes {len(ok)}")
print(f"Depressing {int(ok.depressing.sum())} ({dep:.1f}%), facilitating {int(ok.facilitating.sum())} ({fac:.1f}%)")
print("Within tolerance:", {"depressing": abs(dep - 12.7) <= 1.27, "facilitating": abs(fac - 10.7) <= 1.07})
for name, col in (("depressing", "depressing"), ("facilitating", "facilitating")):
p = stats.binomtest(int(ok[col].sum()), len(ok), 0.05, alternative="greater").pvalue
print(f"{name} above 5%: binomial p = {p:.3g}")
print("Median transmission per presynaptic spike: first", round(ok.eff_first.median(), 4), ", later", round(ok.eff_later.median(), 4),
"; Wilcoxon of first minus later: p =", f"{stats.wilcoxon(ok['diff']).pvalue:.3g}")
print("By session:", ok.groupby("session")[["depressing", "facilitating"]].sum().to_dict())
for lo in (50, 200):
sub = ok[(ok.n_first >= lo) & (ok.n_later >= lo)]
print(f"Pairs with at least {lo} of each: {len(sub)}, depressing {100 * sub.depressing.mean():.1f}%, facilitating {100 * sub.facilitating.mean():.1f}%")
print(f"time {time.time() - t0:.0f} s")
Excitatory pairs 285, with at least 20 first and 20 later spikes 284
Depressing 14 (4.9%), facilitating 44 (15.5%)
Within tolerance: {'depressing': np.False_, 'facilitating': np.False_}
depressing above 5%: binomial p = 0.56
facilitating above 5%: binomial p = 3.5e-11
Median transmission per presynaptic spike: first 0.0068 , later 0.0099 ; Wilcoxon of first minus later: p = 6.08e-08
By session: {'depressing': {'EE.049': 0, 'EE.196': 7, 'FF.103': 3, 'FF.130': 2, 'FF.167': 2, 'GG.069': 0, 'GG.097': 0, 'GG.149': 0}, 'facilitating': {'EE.049': 7, 'EE.196': 26, 'FF.103': 2, 'FF.130': 7, 'FF.167': 1, 'GG.069': 1, 'GG.097': 0, 'GG.149': 0}}
Pairs with at least 50 of each: 279, depressing 5.0%, facilitating 15.4%
Pairs with at least 200 of each: 246, depressing 5.3%, facilitating 15.9%
time 36 s
Result: not reproduced. Of 284 excitatory pairs with enough spikes, 4.9% are depressing and 15.5% facilitating, against 12.7% and 10.7%. Facilitation is above the 5% the test gives by chance (binomial p = 3.5e-11) and depression is not (p = 0.56). Later spikes transmit more than first spikes (median excess 0.0099 against 0.0068 per presynaptic spike, Wilcoxon p = 6.1e-08). Pairs with at least 50 or 200 spikes of each kind give the same fractions within 0.5 points.
fig7-supralinear¶
Coincident firing of two or more presynaptic neurons within 5 ms on a common target is more effective than the sum of their separate effects (supralinear), while spikes spread over windows longer than 10 ms add linearly.
rows = []
for name in exc.session.unique():
s = by_name[name]
e = exc[exc.session == name]
for post, g in e.groupby("post"):
if len(g) < 2:
continue
pres = [s["trains"][p] for p in g.pre]
for w in (0.005, 0.02):
cg = pfc.coincidence_gain(pres, s["trains"][post], s["duration"], w).set_index("n_partners")
row = {"session": name, "post": post, "partners": len(g), "window_ms": int(w * 1000)}
for n_p in (1, 2, 3):
row[f"excess_{n_p}"] = cg.excess.get(n_p, np.nan)
row[f"events_{n_p}"] = int(cg.events.get(n_p, 0))
rows.append(row)
sat = pd.DataFrame(rows)
print("Satellites (targets with at least two excitatory partners):", int((sat.window_ms == 5).sum()))
for w in (5, 20):
g = sat[(sat.window_ms == w) & (sat.events_2 >= 20)].copy()
g["supra"] = g.excess_2 - 2 * g.excess_1
p = stats.wilcoxon(g.supra, alternative="greater").pvalue if len(g) else np.nan
print(f"\n{w} ms events: {len(g)} satellites with at least 20 two-partner events; supralinear in {int((g.supra > 0).sum())} "
f"({100 * (g.supra > 0).mean():.0f}%); median excess one partner {g.excess_1.median():.4f}, two partners "
f"{g.excess_2.median():.4f}; Wilcoxon (two > twice one) p = {p:.3g}")
three = g[g.events_3 >= 20]
if len(three):
print(f" with at least 20 three-partner events: {len(three)}; median excess three partners {three.excess_3.median():.4f}")
print(sat[sat.window_ms == 5].round(4).to_string(index=False))
Satellites (targets with at least two excitatory partners): 41 5 ms events: 40 satellites with at least 20 two-partner events; supralinear in 15 (38%); median excess one partner 0.0145, two partners 0.0197; Wilcoxon (two > twice one) p = 0.993 with at least 20 three-partner events: 16; median excess three partners 0.0276 20 ms events: 41 satellites with at least 20 two-partner events; supralinear in 12 (29%); median excess one partner 0.0129, two partners 0.0189; Wilcoxon (two > twice one) p = 0.998 with at least 20 three-partner events: 22; median excess three partners 0.0334 session post partners window_ms excess_1 events_1 excess_2 events_2 excess_3 events_3 EE.049 11 8 5 0.0207 19939 0.0361 1076 0.0201 23 EE.049 20 6 5 0.0180 28501 0.0361 2281 0.0724 83 EE.049 60 3 5 0.0237 6789 -0.0056 119 NaN 0 EE.049 68 3 5 0.0252 4754 -0.0131 40 NaN 0 EE.196 5 6 5 0.0036 97509 0.0013 31796 0.0062 2544 EE.196 6 2 5 0.0036 7366 -0.0020 41 NaN 0 EE.196 7 11 5 0.0280 46157 0.0509 3291 0.0127 97 EE.196 8 5 5 0.0104 73642 0.0269 2941 -0.0530 29 EE.196 14 2 5 0.0023 4478 -0.0052 4 NaN 0 EE.196 15 11 5 0.0158 39174 0.0258 2853 0.0581 66 EE.196 26 3 5 0.0034 25528 0.0045 749 -0.0169 5 EE.196 31 11 5 0.0165 35631 0.0223 2261 0.0065 51 EE.196 37 3 5 0.0192 5406 0.0586 98 NaN 0 EE.196 50 25 5 0.0440 86954 0.0919 16481 0.1140 1501 EE.196 53 2 5 0.0013 3534 -0.0049 31 NaN 0 EE.196 63 6 5 0.0050 45526 0.0124 3290 0.0216 128 EE.196 79 2 5 0.0021 8426 0.0185 50 NaN 0 EE.196 84 10 5 0.0301 26451 0.0279 1296 0.0336 30 EE.196 105 2 5 0.0325 6133 -0.0001 21 NaN 0 EE.196 106 5 5 0.0564 16446 0.0540 439 0.3219 5 EE.196 117 2 5 0.0023 18606 0.0085 235 NaN 0 EE.196 118 5 5 0.0333 21895 0.0610 572 -0.0562 3 EE.196 125 2 5 0.0022 19030 0.0039 415 NaN 0 EE.196 127 3 5 0.0068 53167 0.0150 1047 -0.0137 6 EE.196 137 6 5 0.0199 34934 0.0297 1740 0.1331 28 FF.103 21 2 5 0.0101 11984 0.0120 150 NaN 0 FF.103 30 2 5 0.0046 8438 -0.0033 71 NaN 0 FF.103 34 5 5 0.0015 40443 0.0036 1913 0.0155 30 FF.103 42 5 5 0.0101 23920 0.0244 819 -0.0391 9 FF.103 50 2 5 0.0094 16334 0.0557 206 NaN 0 FF.130 5 4 5 0.0282 10943 0.0475 242 -0.0434 1 FF.130 6 5 5 0.0132 41033 0.0112 2820 0.1258 53 FF.130 21 6 5 0.0335 23015 0.0756 1501 0.1219 31 FF.130 22 2 5 0.0041 51832 0.0135 1107 NaN 0 FF.130 34 3 5 0.0161 10605 0.0131 495 -0.0111 4 FF.130 53 10 5 0.0282 33355 0.0319 2603 0.0899 66 FF.130 66 3 5 0.0118 15470 0.0196 2757 -0.0047 37 FF.130 67 2 5 0.0057 13241 0.0127 135 NaN 0 FF.167 15 2 5 0.0095 7786 0.0547 45 NaN 0 FF.167 48 2 5 0.0204 17317 0.0290 523 NaN 0 GG.069 17 3 5 0.0275 7129 0.0199 146 NaN 0
# Variant (verifier's point B2): two-spike events against a linear prediction from each
# partner's one-spike events, response window from 1 ms after the first spike to 4 ms after the last.
rng = np.random.default_rng([2008, 14, 0])
T, check = 1800.0, []
for kind in ("linear", "supralinear"):
for k in range(10):
a, b = (np.sort(rng.uniform(0, T, rng.poisson(15 * T))) for _ in range(2))
common = np.sort(rng.uniform(0, T, rng.poisson(3 * T)))
a, b = np.sort(np.r_[a, common]), np.sort(np.r_[b, common + rng.uniform(0, 0.004, len(common))])
post = [rng.uniform(0, T, rng.poisson(5 * T))]
for pre in (a, b):
m = rng.random(len(pre)) < 0.05
post.append(pre[m] + rng.uniform(0.001, 0.004, m.sum()))
if kind == "supralinear":
lags, i, j = pfc.near_lags(a, b, 0.005)
m = rng.random(len(i)) < 0.2
post.append(np.maximum(a[i[m]], b[j[m]]) + rng.uniform(0.001, 0.004, m.sum()))
post = np.sort(np.concatenate(post))
lg = pfc.linear_gain([a, b], post)
cg = pfc.coincidence_gain([a, b], post, T).set_index("n_partners")
check.append({"kind": kind, "linear_variant": lg["excess_2"] - lg["linear_2"], "criterion_measure": cg.excess[2] - 2 * cg.excess[1]})
check = pd.DataFrame(check)
print("Synthetic satellites (two partners at 15 Hz, transmission 0.05 each): two-partner effect minus linear prediction, mean (s.d.)")
print(check.groupby("kind").agg(["mean", "std"]).round(4).to_string())
rows = []
for name in exc.session.unique():
s = by_name[name]
e = exc[exc.session == name]
for post, g in e.groupby("post"):
if len(g) < 2:
continue
for w in (0.005, 0.02):
lg = pfc.linear_gain([s["trains"][p] for p in g.pre], s["trains"][post], w)
rows.append({"session": name, "post": post, "window_ms": int(w * 1000), **lg})
lin = pd.DataFrame(rows)
for w in (5, 20):
g = lin[(lin.window_ms == w) & (lin.events_2 >= 20)].copy()
g["supra"] = g.excess_2 - g.linear_2
print(f"\nVariant, {w} ms events: {len(g)} satellites with at least 20 two-spike events; above the linear prediction in "
f"{int((g.supra > 0).sum())} ({100 * (g.supra > 0).mean():.0f}%); median response {g.excess_2.median():.4f} against prediction "
f"{g.linear_2.median():.4f}; Wilcoxon (above) p = {stats.wilcoxon(g.supra, alternative='greater').pvalue:.3g}")
Synthetic satellites (two partners at 15 Hz, transmission 0.05 each): two-partner effect minus linear prediction, mean (s.d.)
linear_variant criterion_measure
mean std mean std
kind
linear 0.0016 0.0046 -0.0278 0.0030
supralinear 0.1983 0.0041 0.1729 0.0047
Variant, 5 ms events: 40 satellites with at least 20 two-spike events; above the linear prediction in 4 (10%); median response 0.0015 against prediction 0.0239; Wilcoxon (above) p = 1 Variant, 20 ms events: 40 satellites with at least 20 two-spike events; above the linear prediction in 12 (30%); median response 0.0078 against prediction 0.0213; Wilcoxon (above) p = 0.998
# Checks of the variant's local baseline (verifier round 2): the post rate in the flanks of
# one-partner and two-partner events, and the variant with the post unit's mean rate over the
# session as baseline in place of the flanks.
rows = []
for name in exc.session.unique():
s = by_name[name]
e = exc[exc.session == name]
for post, g in e.groupby("post"):
if len(g) < 2:
continue
pres = [s["trains"][p] for p in g.pre]
lg = pfc.linear_gain(pres, s["trains"][post], 0.005, duration=s["duration"])
rows.append({"session": name, "post": post, "mean_rate_hz": len(s["trains"][post]) / s["duration"], **lg})
glob = pd.DataFrame(rows)
glob = glob[glob.events_2 >= 20].copy()
glob["supra"] = glob.excess_2 - glob.linear_2
print(f"{len(glob)} satellites with at least 20 two-spike events (5 ms)")
print(f"Flank rate, median over satellites: one-partner events {glob.flank_hz_1.median():.2f} Hz, two-partner events {glob.flank_hz_2.median():.2f} Hz,"
f" mean rate {glob.mean_rate_hz.median():.2f} Hz; two above one in {int((glob.flank_hz_2 > glob.flank_hz_1).sum())} of {len(glob)},"
f" ratio two/one median {(glob.flank_hz_2 / glob.flank_hz_1).median():.2f}; Wilcoxon p = {stats.wilcoxon(glob.flank_hz_2 - glob.flank_hz_1).pvalue:.3g}")
print(f"Global baseline: one-partner response median {glob.excess_1_mean.median():.4f}; two-partner response median {glob.excess_2.median():.4f}"
f" against linear prediction {glob.linear_2.median():.4f}; above the prediction in {int((glob.supra > 0).sum())} of {len(glob)}"
f" ({100 * (glob.supra > 0).mean():.0f}%); Wilcoxon (above) p = {stats.wilcoxon(glob.supra, alternative='greater').pvalue:.3g},"
f" (below) p = {stats.wilcoxon(glob.supra, alternative='less').pvalue:.3g}")
infl_2 = (glob.flank_hz_2 - glob.mean_rate_hz) * (glob.sep_ms_median + 3) / 1000
infl_pred = 2 * (glob.flank_hz_1 - glob.mean_rate_hz) * 0.003
print(f"Background the global baseline leaves in (flank rate minus mean rate, times the window), median: two-partner response {infl_2.median():.4f},"
f" linear prediction {infl_pred.median():.4f}; larger for the two-partner response in {int((infl_2 > infl_pred).sum())} of {len(glob)}")
loc = lin[(lin.window_ms == 5) & (lin.events_2 >= 20)]
print(f"Same satellites, local baseline: two-partner response median {loc.excess_2.median():.4f} against {loc.linear_2.median():.4f}")
synth = []
rng = np.random.default_rng([2008, 14, 1])
for kind in ("linear", "supralinear"):
for k in range(10):
a, b = (np.sort(rng.uniform(0, T, rng.poisson(15 * T))) for _ in range(2))
common = np.sort(rng.uniform(0, T, rng.poisson(3 * T)))
a, b = np.sort(np.r_[a, common]), np.sort(np.r_[b, common + rng.uniform(0, 0.004, len(common))])
post = [rng.uniform(0, T, rng.poisson(5 * T))]
for pre in (a, b):
m = rng.random(len(pre)) < 0.05
post.append(pre[m] + rng.uniform(0.001, 0.004, m.sum()))
if kind == "supralinear":
lags, i, j = pfc.near_lags(a, b, 0.005)
m = rng.random(len(i)) < 0.2
post.append(np.maximum(a[i[m]], b[j[m]]) + rng.uniform(0.001, 0.004, m.sum()))
post = np.sort(np.concatenate(post))
lg = pfc.linear_gain([a, b], post, duration=T)
synth.append({"kind": kind, "global_variant": lg["excess_2"] - lg["linear_2"]})
print("Synthetic satellites (as above), global baseline: two-partner effect minus linear prediction, mean (s.d.)")
print(pd.DataFrame(synth).groupby("kind").agg(["mean", "std"]).round(4).to_string())
print(glob[["session", "post", "events_2", "flank_hz_1", "flank_hz_2", "mean_rate_hz", "excess_1_mean", "excess_2", "linear_2"]].round(4).to_string(index=False))
40 satellites with at least 20 two-spike events (5 ms) Flank rate, median over satellites: one-partner events 7.26 Hz, two-partner events 9.44 Hz, mean rate 5.97 Hz; two above one in 28 of 40, ratio two/one median 1.15; Wilcoxon p = 0.000219 Global baseline: one-partner response median 0.0180; two-partner response median 0.0135 against linear prediction 0.0337; above the prediction in 9 of 40 (22%); Wilcoxon (above) p = 0.999, (below) p = 0.000584 Background the global baseline leaves in (flank rate minus mean rate, times the window), median: two-partner response 0.0125, linear prediction 0.0074; larger for the two-partner response in 28 of 40 Same satellites, local baseline: two-partner response median 0.0015 against 0.0239
Synthetic satellites (as above), global baseline: two-partner effect minus linear prediction, mean (s.d.)
global_variant
mean std
kind
linear -0.0021 0.0068
supralinear 0.1991 0.0077
session post events_2 flank_hz_1 flank_hz_2 mean_rate_hz excess_1_mean excess_2 linear_2
EE.049 11 1055 9.9653 12.1090 7.8885 0.0262 0.0556 0.0411
EE.049 20 2254 10.3255 11.9122 8.0775 0.0194 0.0343 0.0347
EE.049 60 117 11.5982 11.5385 10.4372 0.0247 -0.0195 0.0488
EE.049 68 40 5.8424 5.0000 4.3974 0.0256 -0.0243 0.0491
EE.196 5 31012 7.6448 7.6478 7.8639 0.0100 -0.0010 0.0072
EE.196 6 38 1.4515 0.6579 0.6698 0.0035 -0.0035 0.0070
EE.196 7 3256 25.3982 29.7681 24.1054 0.0539 0.0498 0.0691
EE.196 8 2899 19.6971 22.6457 18.1558 0.0281 0.0095 0.0416
EE.196 15 2813 7.0194 8.3985 5.9295 0.0219 0.0285 0.0344
EE.196 26 731 5.9010 5.5404 5.6830 0.0036 -0.0030 0.0074
EE.196 31 2181 5.4707 6.6598 4.3885 0.0163 0.0336 0.0322
EE.196 37 98 3.7061 12.2449 0.8653 0.0206 0.0870 0.0438
EE.196 50 16053 32.8046 39.4584 26.3073 0.0565 0.1263 0.0868
EE.196 53 31 1.2263 1.6129 1.6461 0.0013 -0.0096 0.0026
EE.196 63 3245 5.9455 6.2018 5.8634 0.0115 0.0143 0.0157
EE.196 79 50 0.7067 0.0000 0.4852 0.0018 0.0370 0.0037
EE.196 84 1287 12.8802 15.3263 11.1997 0.0314 0.0363 0.0587
EE.196 105 21 21.2228 19.0476 16.2963 0.0301 0.0041 0.0601
EE.196 106 433 36.3734 38.2217 27.1069 0.0521 0.0574 0.1043
EE.196 117 234 1.8158 2.2436 1.4236 0.0029 0.0046 0.0059
EE.196 118 570 24.8819 29.7368 19.2651 0.0351 0.0555 0.0634
EE.196 125 401 0.5959 0.4988 0.2960 0.0021 0.0034 0.0041
EE.196 127 1010 6.0281 8.1188 4.5966 0.0110 0.0132 0.0200
EE.196 137 1685 30.8388 32.3887 28.2210 0.0356 0.0100 0.0365
FF.103 21 150 6.2396 7.1667 4.9135 0.0080 0.0127 0.0160
FF.103 30 71 1.4552 0.7042 1.0890 0.0050 -0.0063 0.0100
FF.103 34 1903 4.7641 5.3468 6.0053 0.0077 -0.0047 0.0082
FF.103 42 818 13.2164 14.7311 13.3019 0.0130 -0.0001 0.0196
FF.103 50 206 2.6509 10.0728 0.8490 0.0165 0.0632 0.0330
FF.130 5 242 16.7491 16.2190 14.7897 0.0392 0.0345 0.0512
FF.130 6 2805 21.5249 20.2496 21.6594 0.0102 -0.0134 0.0220
FF.130 21 1499 26.9057 30.7205 24.7813 0.0312 0.0655 0.0590
FF.130 22 1099 2.5721 3.9581 1.8116 0.0063 0.0132 0.0126
FF.130 34 480 4.4667 3.9063 3.7248 0.0158 0.0099 0.0318
FF.130 53 2594 16.9582 19.1981 15.8615 0.0311 0.0205 0.0552
FF.130 66 2736 4.0384 5.7931 1.5836 0.0117 0.0309 0.0276
FF.130 67 135 11.2015 9.6296 13.3046 0.0063 -0.0166 0.0125
FF.167 15 44 4.2947 6.2500 4.0004 0.0329 0.0470 0.0657
FF.167 48 520 11.7462 16.9231 8.2703 0.0254 0.0499 0.0508
GG.069 17 146 7.4986 9.2466 4.8177 0.0361 0.0138 0.0628
Result: not reproduced. Of 40 satellites with at least 20 two-partner events, 15 (38%) are supralinear by the criterion's measure. The median excess transmission is 0.0145 after one partner and 0.0197 after two, less than twice the single effect, and the Wilcoxon test gives p = 0.993. With 20 ms events 29% are supralinear (p = 0.998). The criterion's measure is biased towards sublinear: on synthetic satellites with linear summation it gives -0.0278 (s.d. 0.0030). The variant against a linear prediction is unbiased there (0.0016, s.d. 0.0046) and detects supralinear input (0.198); on the data, 4 of 40 satellites (10%) are above the linear prediction (median response 0.0015 against 0.0239, p = 1), and with 20 ms events 12 of 40 (30%). The variant's local baseline is not neutral on the data: the post unit fires faster in the flanks of two-partner events than of one-partner events (median 9.44 against 7.26 Hz, higher in 28 of 40 satellites, Wilcoxon p = 0.00022), which lowers the two-partner response. With the post unit's mean rate over the session as baseline (unbiased on the synthetic satellites, -0.0021, s.d. 0.0068; detects supralinear input, 0.199), the two-partner response is 0.0135 against a linear prediction of 0.0337, above it in 9 of 40 satellites (22%), and below it by the Wilcoxon test (p = 0.00058). That baseline leaves in more background for two-partner events than for the prediction (median 0.0125 against 0.0074), which favours a supralinear result, so it does not explain the shortfall. Under none of the three measures is summation supralinear; with both baselines the two-partner response is below the response to one partner (0.0135 against 0.0180 with the global baseline), which these measures do not explain. The paper's satellites converge on interneurons; here any target with two excitatory partners is used.
fig3-trajectory-by-class¶
Large fractions of putative pyramidal cells (72.2% in layer 2/3, 84.7% in layer 5) and of interneurons (58.3% and 83%) are trajectory specific.
m = cells.merge(lr[["session", "unit", "any", "any_chance"]], on=["session", "unit"])
rows = []
for rat, g in m.groupby("rat"):
rows.append({"rat": rat, "pyramidal": int(g.pyr.sum()), "pyr_specific_pct": 100 * g[g.pyr]["any"].mean() if g.pyr.any() else np.nan,
"interneurons": int(g.int.sum()), "int_specific_pct": 100 * g[g.int]["any"].mean() if g.int.any() else np.nan,
"chance_pct": 100 * g.any_chance.mean()})
byclass = pd.DataFrame(rows)
print(byclass.round(1).to_string(index=False))
def layer_match(p, i):
out = []
for name, (vp, vi) in {"2/3": (72.2, 58.3), "5": (84.7, 83.0)}.items():
if abs(p - vp) <= 0.1 * vp and abs(i - vi) <= 0.1 * vi:
out.append(name)
return out
print("Layer readings matched per rat:", {r.rat: layer_match(r.pyr_specific_pct, r.int_specific_pct) for r in byclass.itertuples()})
mc = cells_clean.merge(lr[["session", "unit", "any"]], on=["session", "unit"])
print("Variant with the classes without mirrored troughs (% specific, n):",
{rat: {"pyr": (round(100 * g[g.pyr]["any"].mean(), 1) if g.pyr.any() else None, int(g.pyr.sum())),
"int": (round(100 * g[g.int]["any"].mean(), 1) if g.int.any() else None, int(g.int.sum()))} for rat, g in mc.groupby("rat")})
print("All units, trajectory specific (%):", round(100 * m["any"].mean(), 1), "; chance:", round(100 * m.any_chance.mean(), 1))
rat pyramidal pyr_specific_pct interneurons int_specific_pct chance_pct
EE 94 58.5 5 60.0 2.5
FF 57 66.7 8 62.5 1.9
GG 7 14.3 0 NaN 0.0
Layer readings matched per rat: {'EE': [], 'FF': ['2/3'], 'GG': []}
Variant with the classes without mirrored troughs (% specific, n): {'EE': {'pyr': (np.float64(57.3), 96), 'int': (np.float64(50.0), 4)}, 'FF': {'pyr': (np.float64(67.2), 58), 'int': (np.float64(66.7), 6)}, 'GG': {'pyr': (np.float64(14.3), 7), 'int': (None, 0)}}
All units, trajectory specific (%): 56.6 ; chance: 1.9
Result: not reproduced. In rats EE and FF, 58.5% and 66.7% of putative pyramidal cells and 60.0% and 62.5% of putative interneurons are trajectory specific, far above the chance of 2.5% and 1.9%; FF matches the layer 2/3 values within 10% and EE matches neither layer. Rat GG has 7 putative pyramidal cells (14.3% specific) and no putative interneuron, so the criterion, which asks this of every rat, fails. With the classes recomputed without the troughs that mirror a peak, EE has 4 interneurons (50.0% specific) and FF 6 (66.7%); GG still has none.
ee188-example-session¶
In the example layer 2/3 session, 117 neurons gave 1172 pairs with 48 excitatory and 30 inhibitory connections, 39 of the 117 cells formed one connected hub, and the example pairs show position-dependent, depressing, facilitating and supralinear transmission.
Blocked: the CRCNS description names EE.188 as the session of Figures 2a and 3d, and the dandiset holds no file of that session (nor of EE.152 or EE.165).
print(f"Notebook ran in {time.time() - started:.0f} s")
Notebook ran in 622 s
Summary¶
A failure to reproduce here may be a failure of this attempt, not of the paper: see the reasoning and the differences from the paper's method.
| Claim | Data | Outcome |
|---|---|---|
| channel-region-link | with labels | partially reproduced |
| rats-and-units | with labels | not reproduced |
| task-performance | dandiset | blocked: data missing |
| fig2c-active-fraction | with labels | not reproduced |
| fig2-field-size | with labels | partially reproduced |
| fig2b-trajectory-selective | with labels | partially reproduced |
| traj-diverge | dandiset | blocked: data missing |
| odor-nosepoke | with labels | partially reproduced |
| fig3-connection-rates | with labels | partially reproduced |
| fig3-cell-classes | with labels | not reproduced |
| fig3f-distance | with labels | not reproduced |
| fig4-position-dependence | with labels | blocked: methods underspecified |
| fig6-facilitation-depression | with labels | not reproduced |
| fig7-supralinear | with labels | not reproduced |
| fig3-trajectory-by-class | with labels | not reproduced |
| ee188-example-session | dandiset | blocked: data missing |
How the analysis differs from the paper's. Units are placed on the prefrontal probe by their mean waveform in the wideband signal, with the CRCNS channel map corrected for rat GG, and multi-unit and sensor clusters are left out by rules chosen here; the paper's own selection of well-isolated units is not known. Eight sessions of the dandiset are used (the paper does not list its sessions); EE.088 is left out because its units cannot be placed by waveform. With EE.088 placed by the conversion's block order the unit count is 631, within the paper's 633 plus or minus 63; the count of rats fails because the fourth rat is not in the dandiset. Units are not matched across sessions. The layer of each rat is not known, so values by layer are matched rat by rat. The cross-correlogram bin (1 ms), the number of permutations, the Poisson screen before the jitter test, the measure of transmission for first and later spikes, the definition of satellites and the 1 Hz rule for position analyses are choices made here, as the supplementary information was not available. The test of position dependence uses a band of mean plus c times s.d. because the rank band reaches 0.01 in only 102 of 285 pairs, and without the paper's thinning step it flags pairs of constant efficacy, so that claim is blocked. Inhibitory connections include troughs that mirror a peak in the other direction, as the criterion was written; the rates and classes without them are shown. The trajectories are compared on the x coordinate only, which marks the turn of the path in x and not a lateral separation, so traj-diverge is blocked.
What was found about the dandiset. It holds three of the paper's four rats, and not the
example session EE.188 (nor EE.152 and EE.165). The units table has no channel, shank, region,
cell type or multi-unit flag, and mixes prefrontal units, CA1 units, multi-unit clusters and
clusters of the sensor signals; channel locations are "unknown" and the electrode groups are cut
every eight channels whatever the probe. In rat GG the channel order differs from the CRCNS map,
and the nose-poke sensors are in reverse order relative to the side labels. In six files one
unit's spike times are out of order. The trial table's stop_time is the end of the nose poke
(the paper's trial start), left_or_right carries the description of reward_time, and only
correct trials are listed. The two position series are the x coordinates of the two LEDs, in
pixels though stated in metres, the second all zero; there is no y coordinate. 100 laps
outside the listed runs cannot be classed. The linearized
position is stated in metres but runs from 0 to 2. EE.088 and the other sessions of its day have
no wideband signal.