Reproducing Knierim, McNaughton & Poe (2000) from dandiset 001754¶
Paper. Knierim JJ, McNaughton BL, Poe GR (2000). Three-dimensional spatial selectivity of hippocampal neurons during space flight. Nature Neuroscience 3:209-210. doi:10.1038/72910. The paper is not open access. This attempt used only its abstract and the titles of its three figures, from the publisher's page. Four claims come from the dandiset's own description, which summarises the paper's findings; they could not be checked against the paper's text.
Dandiset. DANDI:001754, version draft. This
version was asked for. It is a repeat attempt: the first attempt, 001754-knierim-2000-01, read
version 0.260728.1352, and the data owners have since replaced all eight files in the draft with a
new conversion that addresses its findings. The fixes are not yet published, so the draft is read.
A draft can change: the draft read here was last modified on 2026-10-14 at 21:55:09 UTC (the data
check prints the time it has when the notebook runs). The three published versions (0.260227.0029,
0.260227.0031, 0.260728.1352) hold the same blobs; every file of the draft is a new blob, and the
dandiset's description is unchanged. The data check compares the draft with 0.260728.1352 where the
analysis uses the difference.
What the dandiset holds. Three rats on the Neurolab mission (STS-90, April 1998): a preflight recording of each, flight day 4 (FD4) of each, and flight day 9 (FD9) of Rats 1 and 2. Each file has sorted units with the original spike files merged into each and their windows, task epochs (BL baseline, ES Escher staircase, MC Magic Carpet) with their original names, and the 64 x 64 rate and occupancy maps of the original analysis for ES and MC, each linked to its unit. The flight files also hold the video position, but only at spike times and only within the rat's own ES and MC epochs. There is no EEG and no head direction.
Conversion code read. catalystneuro/mcnaughton2000-to-nwb at commit 734688253c17a32e26d84c57f41fbe48ffbcd009 (2026-10-14), which made the draft's files; the first attempt read b538d9e. The new commit recovers the original long file names (ES.cell1.rmap, BL2.cell5) from each tetrode folder's WS_FTP.LOG, links each rate map to its unit, names the epochs as the originals did, records each unit's spike files and windows, restricts each flight rat's position to its own ES and MC epochs (merging the shared tracker's samples from both rats' spike files), stores the lost position and unvisited bins as NaN, leaves out rate maps that are identical copies of another map, and drops a stray spike from before the window in 48 spike files. The paper's own analysis code is not linked and was not found.
| Claim | Source | Tested |
|---|---|---|
flight-spatial-code: a robust spatial code on the Escher staircase in flight |
Abstract | Yes |
three-orthogonal-surfaces: the code covers each of the three surfaces |
Abstract | Blocked: metadata insufficient |
si-comparable-to-preflight: flight spatial information comparable to preflight |
Dandiset description | Yes |
fd4-abnormal-selectivity: abnormal selectivity on flight day 4 |
Dandiset description | Yes |
fig3-stability-across-sessions: place fields stable across sessions |
Figure 3 title | Yes |
theta-normal: normal theta in flight |
Dandiset description | Blocked: data missing |
ripples-normal: normal sharp waves and ripples in flight |
Dandiset description | Blocked: data missing |
In [1]:
import json
import re
import time
import warnings
from collections import Counter
from concurrent.futures import ThreadPoolExecutor
import h5py
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import remfile
import requests
from scipy import stats
import neurolab as nl
pd.set_option("display.width", 200)
pd.set_option("display.max_columns", 30)
pd.set_option("display.max_rows", 300)
warnings.filterwarnings("ignore", category=RuntimeWarning)
OLD_VERSION = "0.260728.1352" # the version of the earlier attempt, 001754-knierim-2000-01
started = time.time()
sessions = nl.load_all()
by_label = {nl.label(s): s for s in sessions}
print(f"loaded {len(sessions)} files from dandiset {nl.DANDISET_ID} version {nl.DANDISET_VERSION} "
f"in {time.time() - started:.0f} s")
loaded 8 files from dandiset 001754 version draft in 7 s
Data check¶
Each check of the metadata against the values is stated under its output. The status of each finding of the first attempt on this version is given at the end of the data check.
In [2]:
api = "https://api.dandiarchive.org/api/dandisets/001754/versions"
listing = requests.get(f"{api}/", timeout=60).json()["results"]
print("Versions and when each was last modified:")
for v in listing:
print(f" {v['version']:15s} {v['modified']}")
print("The draft read here was last modified", next(v["modified"] for v in listing if v["version"] == "draft"),
"; when this attempt was made it was 2026-10-14T21:55:09Z, and a later time means the draft has changed since.")
blobs = {}
for v in [x["version"] for x in listing]:
assets = requests.get(f"{api}/{v}/assets/", params={"page_size": 100}, timeout=60).json()["results"]
blobs[v] = sorted((a["path"], a["size"], a.get("blob")) for a in assets)
published = [v for v in blobs if v != "draft"]
print("\nThe published versions hold the same files with the same blobs:", all(blobs[v] == blobs[published[0]] for v in published))
old_blobs = {p: (size, blob) for p, size, blob in blobs[OLD_VERSION]}
print(f"Files of the draft, with their size in {OLD_VERSION} and whether the draft's file is a new blob:")
for path, size, blob in blobs["draft"]:
print(f" {path} {size / 1e6:.2f} MB (was {old_blobs[path][0] / 1e6:.2f} MB) new blob: {blob != old_blobs[path][1]}")
meta = {v: requests.get(f"{api}/{v}/", timeout=60).json() for v in ["draft", OLD_VERSION]}
print(f"The dandiset description is the same in the draft and {OLD_VERSION}:", meta["draft"]["description"] == meta[OLD_VERSION]["description"])
print("\nEvery group (with its neurodata type) and dataset (shape, type, description) of every file, unfiltered:")
types = {}
for path_url in nl.asset_urls():
tree = nl.h5_tree(path_url)
print(f"\n{path_url[0]}: {len(tree)} groups and datasets")
for line in tree:
print(" ", line)
if " group " in line and line.split(" group ")[1]:
types.setdefault(line.split(" group ")[1], set()).add(path_url[0])
print("\nNeurodata types of groups, and the number of files holding each:", {t: len(v) for t, v in sorted(types.items())})
Versions and when each was last modified: draft 2026-10-15T14:09:09.218930Z 0.260227.0029 2026-02-27T00:29:28.150330Z 0.260227.0031 2026-02-27T00:31:34.105096Z 0.260728.1352 2026-07-28T13:52:39.253571Z The draft read here was last modified 2026-10-15T14:09:09.218930Z ; when this attempt was made it was 2026-10-14T21:55:09Z, and a later time means the draft has changed since.
The published versions hold the same files with the same blobs: True Files of the draft, with their size in 0.260728.1352 and whether the draft's file is a new blob: sub-Rat1/sub-Rat1_ses-19980414T125300_ecephys.nwb 2.54 MB (was 2.45 MB) new blob: True sub-Rat1/sub-Rat1_ses-19980420T095700_behavior+ecephys.nwb 3.28 MB (was 3.59 MB) new blob: True sub-Rat1/sub-Rat1_ses-19980425T124500_behavior+ecephys.nwb 7.52 MB (was 8.88 MB) new blob: True sub-Rat2/sub-Rat2_ses-19980413T163700_ecephys.nwb 2.55 MB (was 2.47 MB) new blob: True sub-Rat2/sub-Rat2_ses-19980420T095700_behavior+ecephys.nwb 3.10 MB (was 3.24 MB) new blob: True sub-Rat2/sub-Rat2_ses-19980425T124500_behavior+ecephys.nwb 5.41 MB (was 4.32 MB) new blob: True sub-Rat3/sub-Rat3_ses-19980414T134900_ecephys.nwb 3.53 MB (was 3.46 MB) new blob: True sub-Rat3/sub-Rat3_ses-19980420T152800_behavior+ecephys.nwb 4.14 MB (was 4.04 MB) new blob: True The dandiset description is the same in the draft and 0.260728.1352: False Every group (with its neurodata type) and dataset (shape, type, description) of every file, unfiltered:
sub-Rat1/sub-Rat1_ses-19980414T125300_ecephys.nwb: 72 groups and datasets acquisition/ group analysis/ group file_create_date (1,) object general/ group general/devices/ group general/devices/tetrode_array/ group Device general/experiment_description () object general/experimenter (3,) object general/extracellular_ephys/ group general/extracellular_ephys/TT11/ group ElectrodeGroup general/extracellular_ephys/TT3/ group ElectrodeGroup general/extracellular_ephys/TT4/ group ElectrodeGroup general/extracellular_ephys/TT5/ group ElectrodeGroup general/extracellular_ephys/TT6/ group ElectrodeGroup general/extracellular_ephys/TT7/ group ElectrodeGroup general/extracellular_ephys/TT8/ group ElectrodeGroup general/extracellular_ephys/TT9/ group ElectrodeGroup general/extracellular_ephys/electrodes/ group ElectrodesTable general/extracellular_ephys/electrodes/group (32,) object Reference to the ElectrodeGroup. general/extracellular_ephys/electrodes/group_name (32,) object Name of the ElectrodeGroup. general/extracellular_ephys/electrodes/id (32,) int64 general/extracellular_ephys/electrodes/location (32,) object Location of the electrode (channel). general/institution () object general/keywords (8,) object general/lab () object general/related_publications (1,) object general/subject/ group Subject general/subject/age () object general/subject/description () object general/subject/sex () object general/subject/species () object general/subject/strain () object general/subject/subject_id () object identifier () object intervals/ group intervals/epochs/ group TimeIntervals intervals/epochs/epoch_name (2,) object Epoch name from the original file names (e.g. BL2, ES, ES1) intervals/epochs/id (2,) int64 intervals/epochs/session_type (2,) object Task type: BL (Baseline), ES (Escher Staircase), or MC (Magic Carpet) intervals/epochs/session_type_description (2,) object Description of the task type intervals/epochs/start_time (2,) float64 Start time of epoch, in seconds intervals/epochs/stop_time (2,) float64 Stop time of epoch, in seconds processing/ group processing/ecephys/ group ProcessingModule processing/ecephys/rate_maps/ group DynamicTable processing/ecephys/rate_maps/archived_file_name (62,) object Windows 8.3 short name under which the file is archived in NASA OSDR O processing/ecephys/rate_maps/cluster_id (62,) int64 Cluster ID from Xclust sorting processing/ecephys/rate_maps/epoch_name (62,) object Epoch name from the original file name (e.g. ES, ES2, MC) processing/ecephys/rate_maps/id (62,) int64 processing/ecephys/rate_maps/occupancy_map (62, 64, 64) int32 64x64 occupancy map: the number of 60 Hz video samples in each bin. Sa processing/ecephys/rate_maps/rate_map (62, 64, 64) float32 64x64 spatial firing rate map from the original analysis. The original processing/ecephys/rate_maps/source_file (62,) object Original file name (e.g. ES.cell5.rmap) processing/ecephys/rate_maps/tetrode (62,) object Tetrode name (e.g. TT0) processing/ecephys/rate_maps/unit (62,) int64 Unit in the units table (same tetrode and cluster) session_description () object session_start_time () object stimulus/ group stimulus/presentation/ group stimulus/templates/ group timestamps_reference_time () object units/ group Units units/cluster_id (35,) int64 Cluster ID from Xclust sorting units/electrodes (140,) int64 the electrodes that each spike unit came from units/electrodes_index (35,) uint8 Index for VectorData 'electrodes' units/id (35,) int64 units/obs_intervals (62, 2) float64 the observation intervals for each unit units/obs_intervals_index (35,) uint8 Index for VectorData 'obs_intervals' units/source_files (62,) object Original spike files merged into this unit, as tetrode/long name (e.g. units/source_files_index (35,) uint8 Index for VectorData 'source_files' units/spike_times (21992,) float64 the spike times for each unit in seconds units/spike_times_index (35,) uint16 Index for VectorData 'spike_times' units/tetrode (35,) object Tetrode name (e.g. TT0)
sub-Rat1/sub-Rat1_ses-19980420T095700_behavior+ecephys.nwb: 78 groups and datasets acquisition/ group analysis/ group file_create_date (1,) object general/ group general/devices/ group general/devices/tetrode_array/ group Device general/experiment_description () object general/experimenter (3,) object general/extracellular_ephys/ group general/extracellular_ephys/TT0/ group ElectrodeGroup general/extracellular_ephys/TT1/ group ElectrodeGroup general/extracellular_ephys/TT10/ group ElectrodeGroup general/extracellular_ephys/TT11/ group ElectrodeGroup general/extracellular_ephys/TT2/ group ElectrodeGroup general/extracellular_ephys/TT7/ group ElectrodeGroup general/extracellular_ephys/TT8/ group ElectrodeGroup general/extracellular_ephys/TT9/ group ElectrodeGroup general/extracellular_ephys/electrodes/ group ElectrodesTable general/extracellular_ephys/electrodes/group (32,) object Reference to the ElectrodeGroup. general/extracellular_ephys/electrodes/group_name (32,) object Name of the ElectrodeGroup. general/extracellular_ephys/electrodes/id (32,) int64 general/extracellular_ephys/electrodes/location (32,) object Location of the electrode (channel). general/institution () object general/keywords (8,) object general/lab () object general/related_publications (1,) object general/subject/ group Subject general/subject/age () object general/subject/description () object general/subject/sex () object general/subject/species () object general/subject/strain () object general/subject/subject_id () object identifier () object intervals/ group intervals/epochs/ group TimeIntervals intervals/epochs/epoch_name (5,) object Epoch name from the original file names (e.g. BL2, ES, ES1) intervals/epochs/id (5,) int64 intervals/epochs/session_type (5,) object Task type: BL (Baseline), ES (Escher Staircase), or MC (Magic Carpet) intervals/epochs/session_type_description (5,) object Description of the task type intervals/epochs/start_time (5,) float64 Start time of epoch, in seconds intervals/epochs/stop_time (5,) float64 Stop time of epoch, in seconds processing/ group processing/behavior/ group ProcessingModule processing/behavior/position/ group Position processing/behavior/position/spatial_series/ group SpatialSeries processing/behavior/position/spatial_series/data (20444, 2) float64 pixels processing/behavior/position/spatial_series/reference_frame () object processing/behavior/position/spatial_series/timestamps (20444,) float64 seconds processing/ecephys/ group ProcessingModule processing/ecephys/rate_maps/ group DynamicTable processing/ecephys/rate_maps/archived_file_name (65,) object Windows 8.3 short name under which the file is archived in NASA OSDR O processing/ecephys/rate_maps/cluster_id (65,) int64 Cluster ID from Xclust sorting processing/ecephys/rate_maps/epoch_name (65,) object Epoch name from the original file name (e.g. ES, ES2, MC) processing/ecephys/rate_maps/id (65,) int64 processing/ecephys/rate_maps/occupancy_map (65, 64, 64) int32 64x64 occupancy map: the number of 60 Hz video samples in each bin. Sa processing/ecephys/rate_maps/rate_map (65, 64, 64) float32 64x64 spatial firing rate map from the original analysis. The original processing/ecephys/rate_maps/source_file (65,) object Original file name (e.g. ES.cell5.rmap) processing/ecephys/rate_maps/tetrode (65,) object Tetrode name (e.g. TT0) processing/ecephys/rate_maps/unit (65,) int64 Unit in the units table (same tetrode and cluster) session_description () object session_start_time () object stimulus/ group stimulus/presentation/ group stimulus/templates/ group timestamps_reference_time () object units/ group Units units/cluster_id (47,) int64 Cluster ID from Xclust sorting units/electrodes (188,) int64 the electrodes that each spike unit came from units/electrodes_index (47,) uint8 Index for VectorData 'electrodes' units/id (47,) int64 units/obs_intervals (196, 2) float64 the observation intervals for each unit units/obs_intervals_index (47,) uint8 Index for VectorData 'obs_intervals' units/source_files (196,) object Original spike files merged into this unit, as tetrode/long name (e.g. units/source_files_index (47,) uint8 Index for VectorData 'source_files' units/spike_times (39971,) float64 the spike times for each unit in seconds units/spike_times_index (47,) uint16 Index for VectorData 'spike_times' units/tetrode (47,) object Tetrode name (e.g. TT0)
sub-Rat1/sub-Rat1_ses-19980425T124500_behavior+ecephys.nwb: 78 groups and datasets acquisition/ group analysis/ group file_create_date (1,) object general/ group general/devices/ group general/devices/tetrode_array/ group Device general/experiment_description () object general/experimenter (3,) object general/extracellular_ephys/ group general/extracellular_ephys/TT0/ group ElectrodeGroup general/extracellular_ephys/TT1/ group ElectrodeGroup general/extracellular_ephys/TT10/ group ElectrodeGroup general/extracellular_ephys/TT11/ group ElectrodeGroup general/extracellular_ephys/TT2/ group ElectrodeGroup general/extracellular_ephys/TT6/ group ElectrodeGroup general/extracellular_ephys/TT8/ group ElectrodeGroup general/extracellular_ephys/TT9/ group ElectrodeGroup general/extracellular_ephys/electrodes/ group ElectrodesTable general/extracellular_ephys/electrodes/group (32,) object Reference to the ElectrodeGroup. general/extracellular_ephys/electrodes/group_name (32,) object Name of the ElectrodeGroup. general/extracellular_ephys/electrodes/id (32,) int64 general/extracellular_ephys/electrodes/location (32,) object Location of the electrode (channel). general/institution () object general/keywords (8,) object general/lab () object general/related_publications (1,) object general/subject/ group Subject general/subject/age () object general/subject/description () object general/subject/sex () object general/subject/species () object general/subject/strain () object general/subject/subject_id () object identifier () object intervals/ group intervals/epochs/ group TimeIntervals intervals/epochs/epoch_name (4,) object Epoch name from the original file names (e.g. BL2, ES, ES1) intervals/epochs/id (4,) int64 intervals/epochs/session_type (4,) object Task type: BL (Baseline), ES (Escher Staircase), or MC (Magic Carpet) intervals/epochs/session_type_description (4,) object Description of the task type intervals/epochs/start_time (4,) float64 Start time of epoch, in seconds intervals/epochs/stop_time (4,) float64 Stop time of epoch, in seconds processing/ group processing/behavior/ group ProcessingModule processing/behavior/position/ group Position processing/behavior/position/spatial_series/ group SpatialSeries processing/behavior/position/spatial_series/data (131057, 2) float64 pixels processing/behavior/position/spatial_series/reference_frame () object processing/behavior/position/spatial_series/timestamps (131057,) float64 seconds processing/ecephys/ group ProcessingModule processing/ecephys/rate_maps/ group DynamicTable processing/ecephys/rate_maps/archived_file_name (64,) object Windows 8.3 short name under which the file is archived in NASA OSDR O processing/ecephys/rate_maps/cluster_id (64,) int64 Cluster ID from Xclust sorting processing/ecephys/rate_maps/epoch_name (64,) object Epoch name from the original file name (e.g. ES, ES2, MC) processing/ecephys/rate_maps/id (64,) int64 processing/ecephys/rate_maps/occupancy_map (64, 64, 64) int32 64x64 occupancy map: the number of 60 Hz video samples in each bin. Sa processing/ecephys/rate_maps/rate_map (64, 64, 64) float32 64x64 spatial firing rate map from the original analysis. The original processing/ecephys/rate_maps/source_file (64,) object Original file name (e.g. ES.cell5.rmap) processing/ecephys/rate_maps/tetrode (64,) object Tetrode name (e.g. TT0) processing/ecephys/rate_maps/unit (64,) int64 Unit in the units table (same tetrode and cluster) session_description () object session_start_time () object stimulus/ group stimulus/presentation/ group stimulus/templates/ group timestamps_reference_time () object units/ group Units units/cluster_id (39,) int64 Cluster ID from Xclust sorting units/electrodes (156,) int64 the electrodes that each spike unit came from units/electrodes_index (39,) uint8 Index for VectorData 'electrodes' units/id (39,) int64 units/obs_intervals (139, 2) float64 the observation intervals for each unit units/obs_intervals_index (39,) uint8 Index for VectorData 'obs_intervals' units/source_files (139,) object Original spike files merged into this unit, as tetrode/long name (e.g. units/source_files_index (39,) uint8 Index for VectorData 'source_files' units/spike_times (243517,) float64 the spike times for each unit in seconds units/spike_times_index (39,) uint32 Index for VectorData 'spike_times' units/tetrode (39,) object Tetrode name (e.g. TT0)
sub-Rat2/sub-Rat2_ses-19980413T163700_ecephys.nwb: 70 groups and datasets acquisition/ group analysis/ group file_create_date (1,) object general/ group general/devices/ group general/devices/tetrode_array/ group Device general/experiment_description () object general/experimenter (3,) object general/extracellular_ephys/ group general/extracellular_ephys/TT0/ group ElectrodeGroup general/extracellular_ephys/TT1/ group ElectrodeGroup general/extracellular_ephys/TT10/ group ElectrodeGroup general/extracellular_ephys/TT2/ group ElectrodeGroup general/extracellular_ephys/TT7/ group ElectrodeGroup general/extracellular_ephys/TT8/ group ElectrodeGroup general/extracellular_ephys/electrodes/ group ElectrodesTable general/extracellular_ephys/electrodes/group (24,) object Reference to the ElectrodeGroup. general/extracellular_ephys/electrodes/group_name (24,) object Name of the ElectrodeGroup. general/extracellular_ephys/electrodes/id (24,) int64 general/extracellular_ephys/electrodes/location (24,) object Location of the electrode (channel). general/institution () object general/keywords (8,) object general/lab () object general/related_publications (1,) object general/subject/ group Subject general/subject/age () object general/subject/description () object general/subject/sex () object general/subject/species () object general/subject/strain () object general/subject/subject_id () object identifier () object intervals/ group intervals/epochs/ group TimeIntervals intervals/epochs/epoch_name (2,) object Epoch name from the original file names (e.g. BL2, ES, ES1) intervals/epochs/id (2,) int64 intervals/epochs/session_type (2,) object Task type: BL (Baseline), ES (Escher Staircase), or MC (Magic Carpet) intervals/epochs/session_type_description (2,) object Description of the task type intervals/epochs/start_time (2,) float64 Start time of epoch, in seconds intervals/epochs/stop_time (2,) float64 Stop time of epoch, in seconds processing/ group processing/ecephys/ group ProcessingModule processing/ecephys/rate_maps/ group DynamicTable processing/ecephys/rate_maps/archived_file_name (59,) object Windows 8.3 short name under which the file is archived in NASA OSDR O processing/ecephys/rate_maps/cluster_id (59,) int64 Cluster ID from Xclust sorting processing/ecephys/rate_maps/epoch_name (59,) object Epoch name from the original file name (e.g. ES, ES2, MC) processing/ecephys/rate_maps/id (59,) int64 processing/ecephys/rate_maps/occupancy_map (59, 64, 64) int32 64x64 occupancy map: the number of 60 Hz video samples in each bin. Sa processing/ecephys/rate_maps/rate_map (59, 64, 64) float32 64x64 spatial firing rate map from the original analysis. The original processing/ecephys/rate_maps/source_file (59,) object Original file name (e.g. ES.cell5.rmap) processing/ecephys/rate_maps/tetrode (59,) object Tetrode name (e.g. TT0) processing/ecephys/rate_maps/unit (59,) int64 Unit in the units table (same tetrode and cluster) session_description () object session_start_time () object stimulus/ group stimulus/presentation/ group stimulus/templates/ group timestamps_reference_time () object units/ group Units units/cluster_id (33,) int64 Cluster ID from Xclust sorting units/electrodes (132,) int64 the electrodes that each spike unit came from units/electrodes_index (33,) uint8 Index for VectorData 'electrodes' units/id (33,) int64 units/obs_intervals (59, 2) float64 the observation intervals for each unit units/obs_intervals_index (33,) uint8 Index for VectorData 'obs_intervals' units/source_files (59,) object Original spike files merged into this unit, as tetrode/long name (e.g. units/source_files_index (33,) uint8 Index for VectorData 'source_files' units/spike_times (36664,) float64 the spike times for each unit in seconds units/spike_times_index (33,) uint16 Index for VectorData 'spike_times' units/tetrode (33,) object Tetrode name (e.g. TT0)
sub-Rat2/sub-Rat2_ses-19980420T095700_behavior+ecephys.nwb: 76 groups and datasets acquisition/ group analysis/ group file_create_date (1,) object general/ group general/devices/ group general/devices/tetrode_array/ group Device general/experiment_description () object general/experimenter (3,) object general/extracellular_ephys/ group general/extracellular_ephys/TT12/ group ElectrodeGroup general/extracellular_ephys/TT13/ group ElectrodeGroup general/extracellular_ephys/TT14/ group ElectrodeGroup general/extracellular_ephys/TT17/ group ElectrodeGroup general/extracellular_ephys/TT19/ group ElectrodeGroup general/extracellular_ephys/TT22/ group ElectrodeGroup general/extracellular_ephys/electrodes/ group ElectrodesTable general/extracellular_ephys/electrodes/group (24,) object Reference to the ElectrodeGroup. general/extracellular_ephys/electrodes/group_name (24,) object Name of the ElectrodeGroup. general/extracellular_ephys/electrodes/id (24,) int64 general/extracellular_ephys/electrodes/location (24,) object Location of the electrode (channel). general/institution () object general/keywords (8,) object general/lab () object general/related_publications (1,) object general/subject/ group Subject general/subject/age () object general/subject/description () object general/subject/sex () object general/subject/species () object general/subject/strain () object general/subject/subject_id () object identifier () object intervals/ group intervals/epochs/ group TimeIntervals intervals/epochs/epoch_name (5,) object Epoch name from the original file names (e.g. BL2, ES, ES1) intervals/epochs/id (5,) int64 intervals/epochs/session_type (5,) object Task type: BL (Baseline), ES (Escher Staircase), or MC (Magic Carpet) intervals/epochs/session_type_description (5,) object Description of the task type intervals/epochs/start_time (5,) float64 Start time of epoch, in seconds intervals/epochs/stop_time (5,) float64 Stop time of epoch, in seconds processing/ group processing/behavior/ group ProcessingModule processing/behavior/position/ group Position processing/behavior/position/spatial_series/ group SpatialSeries processing/behavior/position/spatial_series/data (27120, 2) float64 pixels processing/behavior/position/spatial_series/reference_frame () object processing/behavior/position/spatial_series/timestamps (27120,) float64 seconds processing/ecephys/ group ProcessingModule processing/ecephys/rate_maps/ group DynamicTable processing/ecephys/rate_maps/archived_file_name (57,) object Windows 8.3 short name under which the file is archived in NASA OSDR O processing/ecephys/rate_maps/cluster_id (57,) int64 Cluster ID from Xclust sorting processing/ecephys/rate_maps/epoch_name (57,) object Epoch name from the original file name (e.g. ES, ES2, MC) processing/ecephys/rate_maps/id (57,) int64 processing/ecephys/rate_maps/occupancy_map (57, 64, 64) int32 64x64 occupancy map: the number of 60 Hz video samples in each bin. Sa processing/ecephys/rate_maps/rate_map (57, 64, 64) float32 64x64 spatial firing rate map from the original analysis. The original processing/ecephys/rate_maps/source_file (57,) object Original file name (e.g. ES.cell5.rmap) processing/ecephys/rate_maps/tetrode (57,) object Tetrode name (e.g. TT0) processing/ecephys/rate_maps/unit (57,) int64 Unit in the units table (same tetrode and cluster) session_description () object session_start_time () object stimulus/ group stimulus/presentation/ group stimulus/templates/ group timestamps_reference_time () object units/ group Units units/cluster_id (72,) int64 Cluster ID from Xclust sorting units/electrodes (288,) int64 the electrodes that each spike unit came from units/electrodes_index (72,) uint16 Index for VectorData 'electrodes' units/id (72,) int64 units/obs_intervals (194, 2) float64 the observation intervals for each unit units/obs_intervals_index (72,) uint8 Index for VectorData 'obs_intervals' units/source_files (194,) object Original spike files merged into this unit, as tetrode/long name (e.g. units/source_files_index (72,) uint8 Index for VectorData 'source_files' units/spike_times (31766,) float64 the spike times for each unit in seconds units/spike_times_index (72,) uint16 Index for VectorData 'spike_times' units/tetrode (72,) object Tetrode name (e.g. TT0)
sub-Rat2/sub-Rat2_ses-19980425T124500_behavior+ecephys.nwb: 76 groups and datasets acquisition/ group analysis/ group file_create_date (1,) object general/ group general/devices/ group general/devices/tetrode_array/ group Device general/experiment_description () object general/experimenter (3,) object general/extracellular_ephys/ group general/extracellular_ephys/TT12/ group ElectrodeGroup general/extracellular_ephys/TT13/ group ElectrodeGroup general/extracellular_ephys/TT14/ group ElectrodeGroup general/extracellular_ephys/TT15/ group ElectrodeGroup general/extracellular_ephys/TT19/ group ElectrodeGroup general/extracellular_ephys/TT22/ group ElectrodeGroup general/extracellular_ephys/electrodes/ group ElectrodesTable general/extracellular_ephys/electrodes/group (24,) object Reference to the ElectrodeGroup. general/extracellular_ephys/electrodes/group_name (24,) object Name of the ElectrodeGroup. general/extracellular_ephys/electrodes/id (24,) int64 general/extracellular_ephys/electrodes/location (24,) object Location of the electrode (channel). general/institution () object general/keywords (8,) object general/lab () object general/related_publications (1,) object general/subject/ group Subject general/subject/age () object general/subject/description () object general/subject/sex () object general/subject/species () object general/subject/strain () object general/subject/subject_id () object identifier () object intervals/ group intervals/epochs/ group TimeIntervals intervals/epochs/epoch_name (4,) object Epoch name from the original file names (e.g. BL2, ES, ES1) intervals/epochs/id (4,) int64 intervals/epochs/session_type (4,) object Task type: BL (Baseline), ES (Escher Staircase), or MC (Magic Carpet) intervals/epochs/session_type_description (4,) object Description of the task type intervals/epochs/start_time (4,) float64 Start time of epoch, in seconds intervals/epochs/stop_time (4,) float64 Stop time of epoch, in seconds processing/ group processing/behavior/ group ProcessingModule processing/behavior/position/ group Position processing/behavior/position/spatial_series/ group SpatialSeries processing/behavior/position/spatial_series/data (90084, 2) float64 pixels processing/behavior/position/spatial_series/reference_frame () object processing/behavior/position/spatial_series/timestamps (90084,) float64 seconds processing/ecephys/ group ProcessingModule processing/ecephys/rate_maps/ group DynamicTable processing/ecephys/rate_maps/archived_file_name (76,) object Windows 8.3 short name under which the file is archived in NASA OSDR O processing/ecephys/rate_maps/cluster_id (76,) int64 Cluster ID from Xclust sorting processing/ecephys/rate_maps/epoch_name (76,) object Epoch name from the original file name (e.g. ES, ES2, MC) processing/ecephys/rate_maps/id (76,) int64 processing/ecephys/rate_maps/occupancy_map (76, 64, 64) int32 64x64 occupancy map: the number of 60 Hz video samples in each bin. Sa processing/ecephys/rate_maps/rate_map (76, 64, 64) float32 64x64 spatial firing rate map from the original analysis. The original processing/ecephys/rate_maps/source_file (76,) object Original file name (e.g. ES.cell5.rmap) processing/ecephys/rate_maps/tetrode (76,) object Tetrode name (e.g. TT0) processing/ecephys/rate_maps/unit (76,) int64 Unit in the units table (same tetrode and cluster) session_description () object session_start_time () object stimulus/ group stimulus/presentation/ group stimulus/templates/ group timestamps_reference_time () object units/ group Units units/cluster_id (54,) int64 Cluster ID from Xclust sorting units/electrodes (216,) int64 the electrodes that each spike unit came from units/electrodes_index (54,) uint8 Index for VectorData 'electrodes' units/id (54,) int64 units/obs_intervals (179, 2) float64 the observation intervals for each unit units/obs_intervals_index (54,) uint8 Index for VectorData 'obs_intervals' units/source_files (179,) object Original spike files merged into this unit, as tetrode/long name (e.g. units/source_files_index (54,) uint8 Index for VectorData 'source_files' units/spike_times (52170,) float64 the spike times for each unit in seconds units/spike_times_index (54,) uint16 Index for VectorData 'spike_times' units/tetrode (54,) object Tetrode name (e.g. TT0)
sub-Rat3/sub-Rat3_ses-19980414T134900_ecephys.nwb: 74 groups and datasets acquisition/ group analysis/ group file_create_date (1,) object general/ group general/devices/ group general/devices/tetrode_array/ group Device general/experiment_description () object general/experimenter (3,) object general/extracellular_ephys/ group general/extracellular_ephys/TT0/ group ElectrodeGroup general/extracellular_ephys/TT1/ group ElectrodeGroup general/extracellular_ephys/TT10/ group ElectrodeGroup general/extracellular_ephys/TT11/ group ElectrodeGroup general/extracellular_ephys/TT2/ group ElectrodeGroup general/extracellular_ephys/TT4/ group ElectrodeGroup general/extracellular_ephys/TT5/ group ElectrodeGroup general/extracellular_ephys/TT6/ group ElectrodeGroup general/extracellular_ephys/TT8/ group ElectrodeGroup general/extracellular_ephys/TT9/ group ElectrodeGroup general/extracellular_ephys/electrodes/ group ElectrodesTable general/extracellular_ephys/electrodes/group (40,) object Reference to the ElectrodeGroup. general/extracellular_ephys/electrodes/group_name (40,) object Name of the ElectrodeGroup. general/extracellular_ephys/electrodes/id (40,) int64 general/extracellular_ephys/electrodes/location (40,) object Location of the electrode (channel). general/institution () object general/keywords (8,) object general/lab () object general/related_publications (1,) object general/subject/ group Subject general/subject/age () object general/subject/description () object general/subject/sex () object general/subject/species () object general/subject/strain () object general/subject/subject_id () object identifier () object intervals/ group intervals/epochs/ group TimeIntervals intervals/epochs/epoch_name (2,) object Epoch name from the original file names (e.g. BL2, ES, ES1) intervals/epochs/id (2,) int64 intervals/epochs/session_type (2,) object Task type: BL (Baseline), ES (Escher Staircase), or MC (Magic Carpet) intervals/epochs/session_type_description (2,) object Description of the task type intervals/epochs/start_time (2,) float64 Start time of epoch, in seconds intervals/epochs/stop_time (2,) float64 Stop time of epoch, in seconds processing/ group processing/ecephys/ group ProcessingModule processing/ecephys/rate_maps/ group DynamicTable processing/ecephys/rate_maps/archived_file_name (86,) object Windows 8.3 short name under which the file is archived in NASA OSDR O processing/ecephys/rate_maps/cluster_id (86,) int64 Cluster ID from Xclust sorting processing/ecephys/rate_maps/epoch_name (86,) object Epoch name from the original file name (e.g. ES, ES2, MC) processing/ecephys/rate_maps/id (86,) int64 processing/ecephys/rate_maps/occupancy_map (86, 64, 64) int32 64x64 occupancy map: the number of 60 Hz video samples in each bin. Sa processing/ecephys/rate_maps/rate_map (86, 64, 64) float32 64x64 spatial firing rate map from the original analysis. The original processing/ecephys/rate_maps/source_file (86,) object Original file name (e.g. ES.cell5.rmap) processing/ecephys/rate_maps/tetrode (86,) object Tetrode name (e.g. TT0) processing/ecephys/rate_maps/unit (86,) int64 Unit in the units table (same tetrode and cluster) session_description () object session_start_time () object stimulus/ group stimulus/presentation/ group stimulus/templates/ group timestamps_reference_time () object units/ group Units units/cluster_id (48,) int64 Cluster ID from Xclust sorting units/electrodes (192,) int64 the electrodes that each spike unit came from units/electrodes_index (48,) uint8 Index for VectorData 'electrodes' units/id (48,) int64 units/obs_intervals (87, 2) float64 the observation intervals for each unit units/obs_intervals_index (48,) uint8 Index for VectorData 'obs_intervals' units/source_files (87,) object Original spike files merged into this unit, as tetrode/long name (e.g. units/source_files_index (48,) uint8 Index for VectorData 'source_files' units/spike_times (48000,) float64 the spike times for each unit in seconds units/spike_times_index (48,) uint16 Index for VectorData 'spike_times' units/tetrode (48,) object Tetrode name (e.g. TT0)
sub-Rat3/sub-Rat3_ses-19980420T152800_behavior+ecephys.nwb: 85 groups and datasets
acquisition/ group
analysis/ group
file_create_date (1,) object
general/ group
general/devices/ group
general/devices/tetrode_array/ group Device
general/experiment_description () object
general/experimenter (3,) object
general/extracellular_ephys/ group
general/extracellular_ephys/TT0/ group ElectrodeGroup
general/extracellular_ephys/TT10/ group ElectrodeGroup
general/extracellular_ephys/TT11/ group ElectrodeGroup
general/extracellular_ephys/TT13/ group ElectrodeGroup
general/extracellular_ephys/TT16/ group ElectrodeGroup
general/extracellular_ephys/TT17/ group ElectrodeGroup
general/extracellular_ephys/TT18/ group ElectrodeGroup
general/extracellular_ephys/TT19/ group ElectrodeGroup
general/extracellular_ephys/TT20/ group ElectrodeGroup
general/extracellular_ephys/TT22/ group ElectrodeGroup
general/extracellular_ephys/TT23/ group ElectrodeGroup
general/extracellular_ephys/TT3/ group ElectrodeGroup
general/extracellular_ephys/TT4/ group ElectrodeGroup
general/extracellular_ephys/TT5/ group ElectrodeGroup
general/extracellular_ephys/TT8/ group ElectrodeGroup
general/extracellular_ephys/electrodes/ group ElectrodesTable
general/extracellular_ephys/electrodes/group (60,) object Reference to the ElectrodeGroup.
general/extracellular_ephys/electrodes/group_name (60,) object Name of the ElectrodeGroup.
general/extracellular_ephys/electrodes/id (60,) int64
general/extracellular_ephys/electrodes/location (60,) object Location of the electrode (channel).
general/institution () object
general/keywords (8,) object
general/lab () object
general/related_publications (1,) object
general/subject/ group Subject
general/subject/age () object
general/subject/description () object
general/subject/sex () object
general/subject/species () object
general/subject/strain () object
general/subject/subject_id () object
identifier () object
intervals/ group
intervals/epochs/ group TimeIntervals
intervals/epochs/epoch_name (3,) object Epoch name from the original file names (e.g. BL2, ES, ES1)
intervals/epochs/id (3,) int64
intervals/epochs/session_type (3,) object Task type: BL (Baseline), ES (Escher Staircase), or MC (Magic Carpet)
intervals/epochs/session_type_description (3,) object Description of the task type
intervals/epochs/start_time (3,) float64 Start time of epoch, in seconds
intervals/epochs/stop_time (3,) float64 Stop time of epoch, in seconds
processing/ group
processing/behavior/ group ProcessingModule
processing/behavior/position/ group Position
processing/behavior/position/spatial_series/ group SpatialSeries
processing/behavior/position/spatial_series/data (50168, 2) float64 pixels
processing/behavior/position/spatial_series/reference_frame () object
processing/behavior/position/spatial_series/timestamps (50168,) float64 seconds
processing/ecephys/ group ProcessingModule
processing/ecephys/rate_maps/ group DynamicTable
processing/ecephys/rate_maps/archived_file_name (65,) object Windows 8.3 short name under which the file is archived in NASA OSDR O
processing/ecephys/rate_maps/cluster_id (65,) int64 Cluster ID from Xclust sorting
processing/ecephys/rate_maps/epoch_name (65,) object Epoch name from the original file name (e.g. ES, ES2, MC)
processing/ecephys/rate_maps/id (65,) int64
processing/ecephys/rate_maps/occupancy_map (65, 64, 64) int32 64x64 occupancy map: the number of 60 Hz video samples in each bin. Sa
processing/ecephys/rate_maps/rate_map (65, 64, 64) float32 64x64 spatial firing rate map from the original analysis. The original
processing/ecephys/rate_maps/source_file (65,) object Original file name (e.g. ES.cell5.rmap)
processing/ecephys/rate_maps/tetrode (65,) object Tetrode name (e.g. TT0)
processing/ecephys/rate_maps/unit (65,) int64 Unit in the units table (same tetrode and cluster)
session_description () object
session_start_time () object
stimulus/ group
stimulus/presentation/ group
stimulus/templates/ group
timestamps_reference_time () object
units/ group Units
units/cluster_id (47,) int64 Cluster ID from Xclust sorting
units/electrodes (188,) int64 the electrodes that each spike unit came from
units/electrodes_index (47,) uint8 Index for VectorData 'electrodes'
units/id (47,) int64
units/obs_intervals (66, 2) float64 the observation intervals for each unit
units/obs_intervals_index (47,) uint8 Index for VectorData 'obs_intervals'
units/source_files (66,) object Original spike files merged into this unit, as tetrode/long name (e.g.
units/source_files_index (47,) uint8 Index for VectorData 'source_files'
units/spike_times (58716,) float64 the spike times for each unit in seconds
units/spike_times_index (47,) uint16 Index for VectorData 'spike_times'
units/tetrode (47,) object Tetrode name (e.g. TT0)
Neurodata types of groups, and the number of files holding each: {'Device': 8, 'DynamicTable': 8, 'ElectrodeGroup': 8, 'ElectrodesTable': 8, 'Position': 5, 'ProcessingModule': 8, 'SpatialSeries': 5, 'Subject': 8, 'TimeIntervals': 8, 'Units': 8}
The published versions hold the same files; every file of the draft is a new blob, and the
dandiset description has not changed. The listing shows every group and dataset of every file. In
each file acquisition/, analysis/ and the stimulus/ groups are empty, and the neurodata types
are only those printed under the listing (devices, electrode groups and table, subject, epochs,
units, position and the rate-map table): there is no ElectricalSeries or LFP, so no file holds EEG
or any continuous voltage. The units table now has tetrode, cluster_id, spike_times,
obs_intervals, source_files and electrodes; the epochs table epoch_name, session_type and
its description; the rate-map table tetrode, source_file, archived_file_name, epoch_name,
cluster_id, unit, rate_map and occupancy_map; the position one series of x and y in pixels.
In [3]:
# The full descriptions of the tables and series the analysis uses, from the Rat1 FD9 file (the
# other files carry the same text, except the position description, which names each file's
# partner rat), with the subject's age and description from every file.
path_url = [p for p in nl.asset_urls() if p[0].startswith("sub-Rat1/sub-Rat1_ses-19980425")][0]
with h5py.File(remfile.File(path_url[1]), "r") as f:
print("experiment_description:", nl._text(f["general/experiment_description"][()]), "\n")
for name in ["units", "units/spike_times", "intervals/epochs", "processing/ecephys/rate_maps",
"processing/ecephys/rate_maps/rate_map", "processing/ecephys/rate_maps/occupancy_map",
"processing/behavior/position/spatial_series"]:
attrs = {k: nl._text(v) for k, v in f[name].attrs.items() if k in ("description", "resolution")}
print(f"{name}: {attrs}\n")
print("position reference_frame:", nl._text(f["processing/behavior/position/spatial_series/reference_frame"][()]), "\n")
for p, url in nl.asset_urls():
with h5py.File(remfile.File(url), "r") as f:
print(f"{p.split('/')[0]} {p.split('_ses-')[1][:8]}: age {nl._text(f['general/subject/age'][()])!r}; "
f"{nl._text(f['general/subject/description'][()])}")
experiment_description: Hippocampal place cell recordings from rats aboard the Space Shuttle during the Neurolab STS-90 mission (April–May 1998). Three rats with chronically implanted tetrode arrays in hippocampal area CA1 were trained to traverse a three-dimensional track ('Escher staircase') and a flat two-dimensional track ('Magic Carpet') for medial forebrain bundle stimulation reward. Recordings were made preflight on the ground and in-flight on Flight Day 4 and Flight Day 9. In baseline (BL) epochs the rat rested quietly in the Rodent Sleeping Pouch. Rat 1 and Rat 2 were recorded in the same acquisition session on each flight day and took turns: while one rat ran the Escher staircase or Magic Carpet, the other was in a baseline epoch. The source files hold sorted spikes and the original rate maps; no EEG or other continuous voltage was archived.
units: {'description': "Sorted units. Each tetrode was cluster-cut once in Xclust over the whole recording (one TTn.pos.parms source file per tetrode), and the original spike files hold the spikes of one cluster within one epoch. A unit merges the files of one cluster on one tetrode, so it is the same sorted cluster across epochs. Spike times were written with six significant digits, so they are resolved to 1 ms below 1000 s, 10 ms from 1000 to 9999 s, and 0.1 s from 10000 s; spikes in one unit can share a time. In 48 original files the first row is a spike from before the file's epoch window; spikes outside a file's window are left out."}
units/spike_times: {'description': 'the spike times for each unit in seconds', 'resolution': '0.1'}
intervals/epochs: {'description': "Task epochs within the recording session, one per epoch name in the original file names. Where the original files cut an epoch with different windows on different tetrodes, the epoch spans all of them; each unit's obs_intervals give the windows used for that unit."}
processing/ecephys/rate_maps: {'description': "Spatial firing rate maps and occupancy maps from the original analysis, one row per original .rmap file, linked to the unit whose cluster the file name gives. Maps are 64 x 64 bins over the tracker's 256 x 256 pixel image: bin [i, j] covers pos_y from 4i to 4i + 3 and pos_x from 4j to 4j + 3. Maps that were byte-identical copies of another map in the same session were left out."}
processing/ecephys/rate_maps/rate_map: {'description': "64x64 spatial firing rate map from the original analysis. The original documentation describes it as spikes in each bin divided by the bin's occupancy; the maps are smoothed and the kernel is not recorded in the source files. Unvisited bins are NaN (-999 in the source)."}
processing/ecephys/rate_maps/occupancy_map: {'description': '64x64 occupancy map: the number of 60 Hz video samples in each bin. Samples where the tracker lost the animal (pixel 0, 255) are counted in bin [63, 0].'}
processing/behavior/position/spatial_series: {'description': "Position from the video tracker, available in the source only at spike times (the pos_x and pos_y columns of the original spike files), restricted to this rat's ES and MC epochs. Samples where the tracker lost the animal, written as (0, 255) in the source, are NaN. Occupancy cannot be computed from these samples; the rate_maps table holds the original analysis's occupancy, counted from the 60 Hz tracker. This recording was shared with another rat (folder FD9RAT2), and one tracker served both; the samples here merge both rats' spike files. In each window one rat ran the task while the other rested in the sleeping pouch, so the trace in Rat1's ES and MC epochs is taken to be Rat1's position, as in the original rate maps. The source files do not name the tracked animal."}
position reference_frame: Pixel coordinates (0 to 255) of the video tracker image, as written in the original spike files. The camera view and the pixel size on the apparatus are not recorded in the source files.
sub-Rat1 19980414: age '/'; Adult male Fischer 344 rat (Rat1) from Taconic Biosciences, implanted 28 to 30 days before launch. The age is not recorded in the source files.
sub-Rat1 19980420: age '/'; Adult male Fischer 344 rat (Rat1) from Taconic Biosciences, implanted 28 to 30 days before launch. The age is not recorded in the source files.
sub-Rat1 19980425: age '/'; Adult male Fischer 344 rat (Rat1) from Taconic Biosciences, implanted 28 to 30 days before launch. The age is not recorded in the source files.
sub-Rat2 19980413: age '/'; Adult male Fischer 344 rat (Rat2) from Taconic Biosciences, implanted 28 to 30 days before launch. The age is not recorded in the source files.
sub-Rat2 19980420: age '/'; Adult male Fischer 344 rat (Rat2) from Taconic Biosciences, implanted 28 to 30 days before launch. The age is not recorded in the source files.
sub-Rat2 19980425: age '/'; Adult male Fischer 344 rat (Rat2) from Taconic Biosciences, implanted 28 to 30 days before launch. The age is not recorded in the source files.
sub-Rat3 19980414: age '/'; Adult male Fischer 344 rat (Rat3) from Taconic Biosciences, implanted 28 to 30 days before launch. The age is not recorded in the source files.
sub-Rat3 19980420: age '/'; Adult male Fischer 344 rat (Rat3) from Taconic Biosciences, implanted 28 to 30 days before launch. The age is not recorded in the source files.
The draft's descriptions now state what the first attempt had to find from the values: that unvisited rate-map bins are NaN (-999 in the source), that the maps are smoothed with an unrecorded kernel, that the occupancy counts 60 Hz video samples and counts the lost position in bin [63, 0], how rate-map bins map to pixels, that each tetrode was cluster-cut once so a unit is the same sorted cluster in every epoch, how the spike times are rounded, that Rat1 and Rat2 shared one tracker on each flight day, and that the camera view and pixel size are not recorded. Each of these is checked against the values below. The position description says the trace in a rat's ES and MC epochs "is taken to be" that rat's position and that the source files do not name the tracked animal. The age is still "/", now with a subject description saying the age is not recorded. The experiment description says no EEG was archived; the dandiset's own description still reports EEG.
In [4]:
# What changed from the version of the earlier attempt. Both versions are read: epochs, units,
# spikes, position and rate maps of each file. In the old files unvisited rate-map bins are -999,
# in the draft NaN; maps are compared with NaN read as -999.
def read_old(path_url):
with h5py.File(remfile.File(path_url[1]), "r") as f:
ep, u, rm = f["intervals/epochs"], f["units"], f["processing/ecephys/rate_maps"]
out = {"epochs": list(zip(ep["start_time"][:], ep["stop_time"][:], nl._column(ep, "session_type"))),
"units": list(zip(nl._column(u, "tetrode"), u["cluster_id"][:])), "spikes": nl._ragged(u, "spike_times"),
"maps": rm["rate_map"][:], "occ": rm["occupancy_map"][:], "map_types": nl._column(rm, "session_type"),
"map_tetrodes": nl._column(rm, "tetrode"), "map_files": nl._column(rm, "source_file")}
if "processing/behavior" in f:
ss = f["processing/behavior/position/spatial_series"]
out["t"], out["xy"] = ss["timestamps"][:], ss["data"][:]
return out
with ThreadPoolExecutor(8) as pool:
old_urls = nl.asset_urls(OLD_VERSION)
old_files = dict(zip([p for p, _ in old_urls], pool.map(read_old, old_urls)))
def overlapping(intervals):
iv = sorted(intervals)
return sum(1 for (a1, b1), (a2, b2) in zip(iv, iv[1:]) if a2 < b1)
rows, lost_maps = [], []
for s in sessions:
o = old_files[s["path"]]
old_spikes = {(t, int(c)): sp for (t, c), sp in zip(o["units"], o["spikes"])}
new_spikes = {(t, int(c)): sp for t, c, sp in zip(s["units"].tetrode, s["units"].cluster, s["units"].spikes)}
same_units = set(old_spikes) == set(new_spikes)
removed = [len(old_spikes[k]) - len(new_spikes[k]) for k in new_spikes]
subset = all(np.isin(new_spikes[k], old_spikes[k]).all() for k in new_spikes)
new_maps = {np.where(np.isnan(m), -999, m).astype(np.float32).tobytes() for m in s["rate_maps"]}
old_keys = [m.tobytes() for m in o["maps"]]
groups = Counter(old_keys)
for k, n in groups.items():
if k not in new_maps:
lost_maps.append((nl.label(s), n, [f"{o['map_tetrodes'][i]}/{o['map_files'][i]}" for i, x in enumerate(old_keys) if x == k]))
row = {"session": nl.label(s), "epochs_old": len(o["epochs"]), "epochs_new": len(s["epochs"]),
"overlapping_old": overlapping([(a, b) for a, b, _ in o["epochs"]]), "overlapping_new": overlapping(list(zip(s["epochs"].start, s["epochs"].stop))),
"same_units": same_units, "spikes_removed": int(sum(removed)), "units_losing_spikes": int(sum(r > 0 for r in removed)),
"new_spikes_subset_of_old": subset, "maps_old": len(o["maps"]), "distinct_maps_old": len(groups), "maps_new": len(s["rate_maps"]),
"old_distinct_maps_in_new": int(sum(k in new_maps for k in groups))}
if "t" in s:
row["position_old"], row["position_new"] = len(o["t"]), len(s["t"])
rows.append(row)
print(pd.DataFrame(rows).to_string(index=False))
print("\nOld distinct rate maps that are not in the draft (file, copies in the old file, tetrode/short name of each copy):")
for x in lost_maps:
print(" ", x)
print("\nPosition. For each draft flight file, its samples whose timestamp is in the old file of the same rat or of the")
print("other rat, and whether the positions agree there (within 2 px, both valid; the old files wrote (0, 255) for a lost position):")
rows = []
for s in sessions:
if "t" not in s:
continue
partner = {"Rat1": "Rat2", "Rat2": "Rat1"}.get(s["subject"])
row = {"session": nl.label(s), "samples_new": len(s["t"])}
for who, key in [("same_rat", s["path"]), ("other_rat", by_label[f"{partner} {s['day']}"]["path"] if partner else None)]:
if key is None:
continue
o = old_files[key]
common, i_new, i_old = np.intersect1d(s["t"], o["t"], return_indices=True)
old_valid = ~((o["xy"][i_old, 0] == 0) & (o["xy"][i_old, 1] == 255))
both = s["valid"][i_new] & old_valid
dist = np.hypot(*(s["xy"][i_new][both] - o["xy"][i_old][both]).T)
row[f"in_old_{who}_file"] = round(len(common) / len(s["t"]), 3)
row[f"{who}_within_2px"] = round(float((dist < 2).mean()), 3)
row[f"{who}_lost_agrees"] = round(float((~s["valid"][i_new] == ~old_valid).mean()), 3)
in_either = np.isin(s["t"], np.concatenate([old_files[s["path"]]["t"]] + ([old_files[by_label[f"{partner} {s['day']}"]["path"]]["t"]] if partner else [])))
row["in_neither_old_file"] = int((~in_either).sum())
row["their_time_range_s"] = f"{s['t'][~in_either].min():.0f}-{s['t'][~in_either].max():.0f}" if (~in_either).any() else ""
row["old_same_rat_samples_not_in_new"] = int((~np.isin(old_files[s["path"]]["t"], s["t"])).sum())
rows.append(row)
print(pd.DataFrame(rows).to_string(index=False))
session epochs_old epochs_new overlapping_old overlapping_new same_units spikes_removed units_losing_spikes new_spikes_subset_of_old maps_old distinct_maps_old maps_new old_distinct_maps_in_new position_old position_new
Rat1 preflight 2 2 0 0 True 3 3 True 62 62 62 62 NaN NaN
Rat1 FD4 6 5 1 0 True 16 14 True 67 66 65 65 34537.0 20444.0
Rat1 FD9 4 4 0 0 True 10 6 True 64 64 64 64 191065.0 131057.0
Rat2 preflight 2 2 0 0 True 5 5 True 59 59 59 59 NaN NaN
Rat2 FD4 5 5 0 0 True 3 3 True 65 60 57 57 25320.0 27120.0
Rat2 FD9 4 4 0 0 True 5 5 True 80 78 76 76 42858.0 90084.0
Rat3 preflight 2 2 0 0 True 5 5 True 86 86 86 86 NaN NaN
Rat3 FD4 3 3 0 0 True 1 1 True 65 65 65 65 48910.0 50168.0
Old distinct rate maps that are not in the draft (file, copies in the old file, tetrode/short name of each copy):
('Rat1 FD4', 2, ['TT1/ESCELL~1.RMA', 'TT11/ESCELL~1.RMA'])
('Rat2 FD4', 4, ['TT13/ES2BC4~1.RMA', 'TT13/ES7D18~1.RMA', 'TT19/ES3BCC~1.RMA', 'TT22/ESCELL~4.RMA'])
('Rat2 FD4', 2, ['TT13/ES3BCA~1.RMA', 'TT13/ESCELL~4.RMA'])
('Rat2 FD4', 2, ['TT19/ES3BC8~1.RMA', 'TT22/ESCELL~1.RMA'])
('Rat2 FD9', 2, ['TT13/ES3BCC~1.RMA', 'TT14/ESCELL~4.RMA'])
('Rat2 FD9', 2, ['TT13/MC2DC6~1.RMA', 'TT14/MC2DC2~1.RMA'])
Position. For each draft flight file, its samples whose timestamp is in the old file of the same rat or of the
other rat, and whether the positions agree there (within 2 px, both valid; the old files wrote (0, 255) for a lost position):
session samples_new in_old_same_rat_file same_rat_within_2px same_rat_lost_agrees in_old_other_rat_file other_rat_within_2px other_rat_lost_agrees in_neither_old_file their_time_range_s old_same_rat_samples_not_in_new
Rat1 FD4 20444 0.726 1.0 1.0 0.320 0.997 1.000 0 19691
Rat1 FD9 131057 0.917 1.0 1.0 0.143 0.996 0.999 0 70829
Rat2 FD4 27120 0.487 1.0 1.0 0.559 0.999 1.000 0 12107
Rat2 FD9 90084 0.267 1.0 1.0 0.786 0.999 1.000 44 10144-10151 18808
Rat3 FD4 50168 0.975 1.0 1.0 NaN NaN NaN 1259 6309-6348 1
The draft holds the same units as 0.260728.1352, and every spike of the draft is a spike of the earlier version: 1 to 16 spikes per file were removed, the stray first spikes the units' description mentions. Rat1 FD4 has five epochs in place of six, without the overlapping baseline. The rate maps of the draft are all among the earlier version's maps; what is missing are the six groups of identical maps, every copy of which is left out (14 maps), so these six maps are no longer in the dandiset. The draft's position samples are all at timestamps of the earlier files of the same rat or of the other rat, where they agree (99.6 to 100% within 2 pixels), except 44 in the last 7 s of Rat2 FD9 and 1259 in the last 40 s of Rat3 FD4, both within MC epochs. Each file now holds the samples of both rats' spike files within its own rat's task epochs and none elsewhere.
In [5]:
rows = []
for s in sessions:
ep = s["epochs"]
rows.append({"file": s["path"].split("/")[1][:31], "rat": s["subject"], "day": s["day"], "start": s["start"][:16],
"units": len(s["units"]), "spikes": int(sum(len(x) for x in s["units"]["spikes"])),
"tetrodes": s["units"]["tetrode"].nunique(),
"epochs": " ".join(f"{n}[{a:.0f}-{b:.0f}]" for a, b, n in ep[["start", "stop", "name"]].itertuples(index=False)),
"rate_maps": len(s["maps"]), "map_epochs": "/".join(sorted(s["maps"]["epoch_name"].unique())),
"position_samples": len(s["t"]) if "t" in s else 0})
info = pd.DataFrame(rows)
print(info.to_string(index=False))
print("\nFiles per rat and day:"); print(pd.crosstab(info.rat, info.day).to_string())
for s in sessions:
print(f"{nl.label(s):14s} {s['description']}")
print("\nTetrodes with units, per file:")
for s in sessions:
print(f" {nl.label(s):14s}", sorted(s["units"].tetrode.unique(), key=lambda t: int(t[2:])))
file rat day start units spikes tetrodes epochs rate_maps map_epochs position_samples sub-Rat1_ses-19980414T125300_ec Rat1 preflight 1998-04-14T12:53 35 21992 8 ES[276-741] MC[936-1349] 62 ES/MC 0 sub-Rat1_ses-19980420T095700_be Rat1 FD4 1998-04-20T09:57 47 39971 8 BL1[2525-3390] BL2[3850-4800] ES[5935-6420] MC[7020-7550] BL3[7850-8473] 65 ES/MC 20444 sub-Rat1_ses-19980425T124500_be Rat1 FD9 1998-04-25T12:45 39 243517 8 ES[3178-4789] BL1[4871-7535] MC[7888-9200] BL2[9220-10151] 64 ES/MC 131057 sub-Rat2_ses-19980413T163700_ec Rat2 preflight 1998-04-13T16:37 33 36664 6 ES[1013-1573] MC[1667-2234] 59 ES/MC 0 sub-Rat2_ses-19980420T095700_be Rat2 FD4 1998-04-20T09:57 72 31766 6 BL1[2525-3390] ES[3850-4800] BL2[5935-6420] BL3[7020-7550] MC[7850-8473] 57 ES/MC 27120 sub-Rat2_ses-19980425T124500_be Rat2 FD9 1998-04-25T12:45 54 52170 6 BL1[3178-4789] ES[4871-7535] BL2[7888-9200] MC[9220-10151] 76 ES/MC 90084 sub-Rat3_ses-19980414T134900_ec Rat3 preflight 1998-04-14T13:49 48 48000 10 ES[148-762] MC[1095-1481] 86 ES/MC 0 sub-Rat3_ses-19980420T152800_be Rat3 FD4 1998-04-20T15:28 47 58716 15 ES1[32-1345] ES2[2397-3729] MC[5713-6349] 65 ES1/ES2/MC 50168 Files per rat and day: day FD4 FD9 preflight rat Rat1 1 1 1 Rat2 1 1 1 Rat3 1 0 1 Rat1 preflight Preflight ground recording — Rat 1. Recorded 3 days before launch at Kennedy Space Center. Rat1 FD4 Flight Day 4 recording — Rat 1. Recorded in the same acquisition session as Rat 2; the two rats took turns on the task apparatus. Rat1 FD9 Flight Day 9 recording — Rat 1. Recorded in the same acquisition session as Rat 2; the two rats took turns on the task apparatus. Rat2 preflight Preflight ground recording — Rat 2. Recorded 4 days before launch at Kennedy Space Center. Rat2 FD4 Flight Day 4 recording — Rat 2. Recorded in the same acquisition session as Rat 1; the two rats took turns on the task apparatus. Rat2 FD9 Flight Day 9 recording — Rat 2. Recorded in the same acquisition session as Rat 1; the two rats took turns on the task apparatus. Rat3 preflight Preflight ground recording — Rat 3. Recorded 3 days before launch at Kennedy Space Center. Rat3 FD4 Flight Day 4 recording — Rat 3. Separate recording session; partial data recovered due to technical issues. Tetrodes with units, per file: Rat1 preflight ['TT3', 'TT4', 'TT5', 'TT6', 'TT7', 'TT8', 'TT9', 'TT11'] Rat1 FD4 ['TT0', 'TT1', 'TT2', 'TT7', 'TT8', 'TT9', 'TT10', 'TT11'] Rat1 FD9 ['TT0', 'TT1', 'TT2', 'TT6', 'TT8', 'TT9', 'TT10', 'TT11'] Rat2 preflight ['TT0', 'TT1', 'TT2', 'TT7', 'TT8', 'TT10'] Rat2 FD4 ['TT12', 'TT13', 'TT14', 'TT17', 'TT19', 'TT22'] Rat2 FD9 ['TT12', 'TT13', 'TT14', 'TT15', 'TT19', 'TT22'] Rat3 preflight ['TT0', 'TT1', 'TT2', 'TT4', 'TT5', 'TT6', 'TT8', 'TT9', 'TT10', 'TT11'] Rat3 FD4 ['TT0', 'TT3', 'TT4', 'TT5', 'TT8', 'TT10', 'TT11', 'TT13', 'TT16', 'TT17', 'TT18', 'TT19', 'TT20', 'TT22', 'TT23']
Three rats, eight files: preflight and FD4 for each rat, FD9 for Rats 1 and 2, as the description says. Rate maps exist only for ES and MC epochs, not for baseline; Rat3 FD4 has two ES epochs, ES1 and ES2, each with its own maps. Position exists only in the five flight files. No epochs overlap now. Rat2's preflight units are on tetrodes TT0 to TT10 and its flight units on TT12 to TT22, with no tetrode in common, while Rat1's and Rat3's preflight and flight tetrodes overlap. The conversion notes list whether the preflight folders are assigned to the right rats as an open question; nothing in the files settles it, so the per-rat comparisons with preflight below assume the assignment.
In [6]:
print("Position, flight files. Stated unit:", {s['xy_unit'] for s in sessions if 'xy' in s})
rows = []
for s in sessions:
if "xy" not in s:
continue
xy, v, t = s["xy"], s["valid"], s["t"]
own = np.unique(np.concatenate(list(s["units"]["spikes"])))
partner = {"Rat1": "Rat2", "Rat2": "Rat1"}.get(s["subject"])
other = np.unique(np.concatenate(list(by_label[f"{partner} {s['day']}"]["units"]["spikes"]))) if partner else np.array([])
task = np.zeros(len(t), bool)
own_task = np.zeros(len(own), bool)
for a, b, typ in s["epochs"][["start", "stop", "type"]].itertuples(index=False):
if typ != "BL":
task |= (t >= a) & (t <= b)
own_task |= (own >= a) & (own <= b)
rows.append({"session": nl.label(s), "samples": len(t), "lost_(NaN)": round(1 - v.mean(), 3),
"x_min": xy[v, 0].min(), "x_max": xy[v, 0].max(), "y_min": xy[v, 1].min(), "y_max": xy[v, 1].max(),
"t_first": t[0], "t_last": t[-1], "in_own_ES_MC_epochs": round(task.mean(), 3),
"at_own_spike_times": round(float(np.isin(t, own).mean()), 3),
"at_other_rat_spike_times": round(float(np.isin(t, other).mean()), 3) if partner else np.nan,
"at_either": round(float(np.isin(t, np.union1d(own, other)).mean()), 3),
"own_task_spike_times_with_sample": round(float(np.isin(own[own_task], t).mean()), 4)})
print(pd.DataFrame(rows).round(2).to_string(index=False))
Position, flight files. Stated unit: {'pixels'}
session samples lost_(NaN) x_min x_max y_min y_max t_first t_last in_own_ES_MC_epochs at_own_spike_times at_other_rat_spike_times at_either own_task_spike_times_with_sample Rat1 FD4 20444 0.12 0.0 255.0 64.00 255.00 5935.53 7549.60 1.0 0.73 0.32 1.0 1.0 Rat1 FD9 131057 0.09 0.0 255.0 64.00 255.00 3178.02 9199.50 1.0 0.92 0.14 1.0 1.0 Rat2 FD4 27120 0.13 0.0 251.8 65.00 255.00 3850.58 8472.99 1.0 0.49 0.56 1.0 1.0 Rat2 FD9 90084 0.06 0.0 255.0 64.00 254.90 4871.06 10151.00 1.0 0.27 0.79 1.0 1.0 Rat3 FD4 50168 0.11 0.0 255.0 54.33 254.94 32.17 6348.26 1.0 1.00 NaN 1.0 1.0
Position is in video pixels from 0 to 255 in x and y. Lost positions are NaN (6 to 13% of samples); these are dropped. Every sample lies in one of the file's own ES or MC epochs and at a spike time of one of the two rats sharing the tracker (in Rat3's file, of its own units); in Rat2's files most samples are at Rat1's spike times (0.56 on FD4, 0.79 on FD9). Every spike of the file's own units in its task epochs has a sample. The trace is sampled only when some unit fires, at irregular intervals, consistent with the description. No occupancy time can be built from it; the occupancy maps of the original analysis are used.
In [7]:
# The position lookup on synthetic input: samples at 1.0 s, 1.1 s (the lost position, NaN) and 1.5 s.
fake = {"t": np.array([1.0, 1.1, 1.5]), "xy": np.array([[40.0, 40.0], [np.nan, np.nan], [200.0, 200.0]])}
fake["valid"] = ~np.isnan(fake["xy"]).any(axis=1)
found, r, c = nl.PositionLookup(fake)(np.array([1.0, 1.1, 1.12, 1.4, 2.0]))
print("times 1.0, 1.1, 1.12, 1.4, 2.0 s; found:", found.tolist(), "; rows:", r.tolist(), "; columns:", c.tolist())
assert found.tolist() == [True, False, False, True, False] and r.tolist() == [10, 50] and c.tolist() == [10, 50]
times 1.0, 1.1, 1.12, 1.4, 2.0 s; found: [True, False, False, True, False] ; rows: [10, 50] ; columns: [10, 50]
The lookup gives a spike the position of the nearest sample within 0.25 s and drops the spike when that sample is the lost position (NaN in the draft). On the synthetic input, the spikes at 1.1 s (the lost sample itself) and at 1.12 s (nearest to it) are dropped, as is the spike at 2.0 s, 0.5 s from any sample.
In [8]:
print("Do the Rat1 and Rat2 files of a flight day share position samples, and which epochs do their samples fall in?")
for day in ["FD4", "FD9"]:
a, b = by_label[f"Rat1 {day}"], by_label[f"Rat2 {day}"]
print(f" {day}: timestamps in both files: {len(np.intersect1d(a['t'], b['t']))}")
for s in (a, b):
n = {f"{name} ({typ})": int(((s["t"] >= x) & (s["t"] <= y)).sum()) for x, y, typ, name in s["epochs"][["start", "stop", "type", "name"]].itertuples(index=False)}
print(f" {nl.label(s)} samples per own epoch: {n}")
print(" Rat1 epochs:", list(a["epochs"][["start", "stop", "name"]].itertuples(index=False, name=None)))
print(" Rat2 epochs:", list(b["epochs"][["start", "stop", "name"]].itertuples(index=False, name=None)))
Do the Rat1 and Rat2 files of a flight day share position samples, and which epochs do their samples fall in?
FD4: timestamps in both files: 0
Rat1 FD4 samples per own epoch: {'BL1 (BL)': 0, 'BL2 (BL)': 0, 'ES (ES)': 8657, 'MC (MC)': 11787, 'BL3 (BL)': 0}
Rat2 FD4 samples per own epoch: {'BL1 (BL)': 0, 'ES (ES)': 17325, 'BL2 (BL)': 0, 'BL3 (BL)': 0, 'MC (MC)': 9795}
Rat1 epochs: [(2525.0, 3390.0, 'BL1'), (3850.0, 4800.0, 'BL2'), (5935.0, 6420.0, 'ES'), (7020.0, 7550.0, 'MC'), (7850.0, 8473.0, 'BL3')]
Rat2 epochs: [(2525.0, 3390.0, 'BL1'), (3850.0, 4800.0, 'ES'), (5935.0, 6420.0, 'BL2'), (7020.0, 7550.0, 'BL3'), (7850.0, 8473.0, 'MC')]
FD9: timestamps in both files: 0
Rat1 FD9 samples per own epoch: {'ES (ES)': 65526, 'BL1 (BL)': 0, 'MC (MC)': 65531, 'BL2 (BL)': 0}
Rat2 FD9 samples per own epoch: {'BL1 (BL)': 0, 'ES (ES)': 67380, 'BL2 (BL)': 0, 'MC (MC)': 22704}
Rat1 epochs: [(3178.0, 4789.0, 'ES'), (4871.0, 7535.0, 'BL1'), (7888.0, 9200.0, 'MC'), (9220.0, 10151.0, 'BL2')]
Rat2 epochs: [(3178.0, 4789.0, 'BL1'), (4871.0, 7535.0, 'ES'), (7888.0, 9200.0, 'BL2'), (9220.0, 10151.0, 'MC')]
The Rat1 and Rat2 files of a flight day no longer share any timestamp: each holds samples only in its own rat's ES and MC epochs, and the epochs show that the rats took turns (when one rat is labelled ES or MC, the other is labelled BL). In 0.260728.1352 both files held the same trace over the whole session, presented as each rat's position; the draft assigns each window to one rat. The next checks ask whether the trace in each window is the position of the rat it is assigned to.
In [9]:
print("Distinct occupancy maps per epoch name, per file:", {nl.label(s): dict(Counter(n for n, _, _ in nl.map_occupancy_groups(s))) for s in sessions})
rows = []
for s in sessions:
if "xy" not in s:
continue
for a, b, typ, name in s["epochs"][["start", "stop", "type", "name"]].itertuples(index=False):
if typ == "BL":
continue
scores = nl.occupancy_match(s, a, b)
occ = nl.epoch_occupancy(s, name)
rows.append({"session": nl.label(s), "epoch": f"{name} {a:.0f}-{b:.0f}", "duration_s": b - a,
"match_own_name": round(scores[name], 2), "best_match": max(scores, key=scores.get),
"match_other_names": {k: round(v, 2) for k, v in scores.items() if k != name},
"occupancy_counts": int(occ.sum()), "counts_per_s": round(occ.sum() / (b - a), 1)})
print("Correlation of each binary occupancy map with the binary histogram of the epoch's position samples:")
print(pd.DataFrame(rows).to_string(index=False))
Distinct occupancy maps per epoch name, per file: {'Rat1 preflight': {'ES': 1, 'MC': 1}, 'Rat1 FD4': {'ES': 1, 'MC': 1}, 'Rat1 FD9': {'ES': 1, 'MC': 1}, 'Rat2 preflight': {'ES': 1, 'MC': 1}, 'Rat2 FD4': {'ES': 1, 'MC': 1}, 'Rat2 FD9': {'ES': 1, 'MC': 1}, 'Rat3 preflight': {'ES': 1, 'MC': 1}, 'Rat3 FD4': {'ES1': 1, 'ES2': 1, 'MC': 1}}
Correlation of each binary occupancy map with the binary histogram of the epoch's position samples:
session epoch duration_s match_own_name best_match match_other_names occupancy_counts counts_per_s
Rat1 FD4 ES 5935-6420 485.0 0.79 ES {'MC': 0.31} 26967 55.6
Rat1 FD4 MC 7020-7550 530.0 0.71 MC {'ES': 0.3} 25739 48.6
Rat1 FD9 ES 3178-4789 1611.0 0.88 ES {'MC': 0.55} 84975 52.7
Rat1 FD9 MC 7888-9200 1312.0 0.88 MC {'ES': 0.54} 59340 45.2
Rat2 FD4 ES 3850-4800 950.0 0.79 ES {'MC': 0.47} 53333 56.1
Rat2 FD4 MC 7850-8473 623.0 0.70 MC {'ES': 0.43} 34325 55.1
Rat2 FD9 ES 4871-7535 2664.0 0.82 ES {'MC': 0.51} 123177 46.2
Rat2 FD9 MC 9220-10151 931.0 0.75 MC {'ES': 0.5} 46561 50.0
Rat3 FD4 ES1 32-1345 1313.0 0.81 ES1 {'ES2': 0.67, 'MC': 0.32} 74454 56.7
Rat3 FD4 ES2 2397-3729 1332.0 0.83 ES2 {'ES1': 0.63, 'MC': 0.32} 65669 49.3
Rat3 FD4 MC 5713-6349 636.0 0.79 MC {'ES1': 0.34, 'ES2': 0.32} 24357 38.3
Each epoch name has one occupancy map in its file. For each ES and MC epoch, the position samples in the epoch cover the bins of the occupancy map of the same name best (0.70 to 0.88, against 0.30 to 0.67 for the other maps), including Rat3's ES1 and ES2. The occupancy totals are 38 to 57 counts per second of the epoch, the same order in every epoch and below the 60 Hz video rate the occupancy description now states; why not every frame is counted is not stated. The epoch names, times and position are therefore consistent with the original analysis. This is a check of position against occupancy; whether spikes and position share a clock is shown with spike-based maps below.
In [10]:
print("Rate-map values:")
r_all = np.concatenate([s["rate_maps"].ravel() for s in sessions])
o_all = np.concatenate([s["occupancy_maps"].ravel() for s in sessions])
print(" NaN bins:", round(float(np.isnan(r_all).mean()), 3), "; NaN exactly where occupancy is 0:", bool((np.isnan(r_all) == (o_all == 0)).all()),
"; negative values:", int((r_all < 0).sum()))
print(" max rate per file (Hz):", {nl.label(s): round(float(np.nanmax(s["rate_maps"])), 1) for s in sessions})
print("\nOccupancy in bin [63, 0], the bin of the lost position (0, 255), as a fraction of each occupancy map:")
for s in sessions:
print(" ", nl.label(s), [f"{name} {occ[nl.SENTINEL_BIN] / occ.sum():.3f}" for name, occ, _ in nl.map_occupancy_groups(s)])
maps = nl.rma_table(sessions)
print("\nRate maps per file and type: all, identical copies of another map in the same file, archived short names with CEL:")
print(maps.groupby(["rat", "day", "type"]).agg(maps=("info", "size"), duplicates=("duplicate", "sum"), cel_names=("cel_name", "sum")).to_string())
print("\nExamples of source_file and archived_file_name:", [f"{a} = {b}" for a, b in zip(sessions[2]["maps"].source_file[:8], sessions[2]["maps"].archived_file_name[:8])])
print("source_file names that do not read <epoch name>.cell<cluster id>.rmap:",
sum(f != f"{n}.cell{c}.rmap" for s in sessions for f, n, c in zip(s["maps"].source_file, s["maps"].epoch_name, s["maps"].cluster)))
print("Rate maps whose unit has a different tetrode or cluster id:",
sum((s["units"].tetrode[u] != t) or (s["units"].cluster[u] != c) for s in sessions for u, t, c in zip(s["maps"].unit, s["maps"].tetrode, s["maps"].cluster)))
print("Rate maps whose unit has no spike file of the map's epoch:",
sum(f"{t}/{n}.cell{c}" not in s["units"].source_files[u] for s in sessions for u, t, n, c in zip(s["maps"].unit, s["maps"].tetrode, s["maps"].epoch_name, s["maps"].cluster)))
pk = []
for s in sessions:
for i in range(len(s["maps"])):
occ = nl.clean_occupancy(s["occupancy_maps"][i])
rate = np.where((occ > 0) & np.isfinite(s["rate_maps"][i]), s["rate_maps"][i], -1.0)
j = np.unravel_index(rate.argmax(), rate.shape)
pk.append({"session": nl.label(s), "peak_rate": rate[j], "occupancy_at_peak": occ[j]})
pk = pd.DataFrame(pk)
print("\nOccupancy counts in the bin of each map's peak rate:")
print(pk.groupby("session").agg(maps=("peak_rate", "size"), max_peak_hz=("peak_rate", "max"), median_occupancy_at_peak=("occupancy_at_peak", "median"),
peaks_in_bins_below_5=("occupancy_at_peak", lambda o: int((o < nl.MIN_OCC).sum()))).round(1).to_string())
print("Five highest peaks:"); print(pk.nlargest(5, "peak_rate").round(1).to_string(index=False))
Rate-map values:
NaN bins: 0.607 ; NaN exactly where occupancy is 0: True ; negative values: 0
max rate per file (Hz): {'Rat1 preflight': 16.0, 'Rat1 FD4': 10.5, 'Rat1 FD9': 277.8, 'Rat2 preflight': 16.8, 'Rat2 FD4': 30.5, 'Rat2 FD9': 5.2, 'Rat3 preflight': 15.4, 'Rat3 FD4': 26.7}
Occupancy in bin [63, 0], the bin of the lost position (0, 255), as a fraction of each occupancy map:
Rat1 preflight ['ES 0.006', 'MC 0.005']
Rat1 FD4 ['ES 0.079', 'MC 0.179']
Rat1 FD9 ['ES 0.045', 'MC 0.168']
Rat2 preflight ['ES 0.007', 'MC 0.009']
Rat2 FD4 ['ES 0.175', 'MC 0.113']
Rat2 FD9 ['ES 0.043', 'MC 0.172']
Rat3 preflight ['ES 0.005', 'MC 0.014']
Rat3 FD4 ['ES1 0.040', 'ES2 0.051', 'MC 0.190']
Rate maps per file and type: all, identical copies of another map in the same file, archived short names with CEL:
maps duplicates cel_names
rat day type
Rat1 FD4 ES 33 0 22
MC 32 0 24
FD9 ES 26 0 23
MC 38 0 27
preflight ES 27 0 22
MC 35 0 26
Rat2 FD4 ES 31 0 15
MC 26 0 19
FD9 ES 39 0 23
MC 37 0 23
preflight ES 26 0 20
MC 33 0 21
Rat3 FD4 ES 42 0 37
MC 23 0 19
preflight ES 45 0 34
MC 41 0 29
Examples of source_file and archived_file_name: ['ES.cell1.rmap = ESCELL~1.RMA', 'ES.cell2.rmap = ESCELL~2.RMA', 'ES.cell3.rmap = ESCELL~3.RMA', 'MC.cell7.rmap = MC2DC0~1.RMA', 'MC.cell6.rmap = MC2DC2~1.RMA', 'MC.cell1.rmap = MCCELL~1.RMA', 'MC.cell2.rmap = MCCELL~2.RMA', 'MC.cell3.rmap = MCCELL~3.RMA']
source_file names that do not read <epoch name>.cell<cluster id>.rmap: 0
Rate maps whose unit has a different tetrode or cluster id: 0
Rate maps whose unit has no spike file of the map's epoch: 0
Occupancy counts in the bin of each map's peak rate:
maps max_peak_hz median_occupancy_at_peak peaks_in_bins_below_5
session
Rat1 FD4 65 10.5 4.0 36
Rat1 FD9 64 277.8 24.0 10
Rat1 preflight 62 16.0 15.5 11
Rat2 FD4 57 30.5 6.0 20
Rat2 FD9 76 5.2 15.0 25
Rat2 preflight 59 16.8 17.0 9
Rat3 FD4 65 25.9 14.0 20
Rat3 preflight 86 15.4 12.0 15
Five highest peaks:
session peak_rate occupancy_at_peak
Rat1 FD9 277.8 55.0
Rat2 FD4 30.5 1.0
Rat3 FD4 25.9 22.0
Rat3 FD4 21.6 22.0
Rat3 FD4 19.1 32.0
Unvisited rate-map bins are NaN, exactly where the occupancy is zero, as the description now
says, and no value is negative. The lost position is still counted in the original occupancy (bin
[63, 0], up to 19% of an occupancy map in flight, under 2% preflight), as the occupancy description
now says; that bin is left out here. No file holds two identical rate maps. Every source_file reads
<epoch name>.cell<cluster id>.rmap, every map's unit has the map's tetrode and cluster id, and
every map's unit has a spike file of the map's epoch; the old 8.3 short names are kept in
archived_file_name. The "CEL names only" variant below uses those short names, as the first
attempt did.
The highest peak rate, 277.8 Hz in a Rat1 FD9 map, lies in a bin with 55 occupancy counts, not a sparsely visited bin. Per file, 9 to 36 maps have their peak in a bin with fewer than 5 counts; spatial information over bins with at least 5 counts is shown as a variant of the rate-map claims.
In [11]:
# What apparatus the preflight epochs labelled ES and MC were recorded on: the dandiset description,
# and the extent of every distinct occupancy map (visited bins, rows and columns of the 64 x 64 grid).
desc = meta["draft"]["description"]
print("Dandiset description, sentences on the track:")
for sentence in desc.split(". "):
if "track" in sentence:
print(" ", sentence.strip())
print("\nEpoch names per file:", {nl.label(s): sorted(set(s["epochs"].name)) for s in sessions})
rows, footprint = [], {}
for s in sessions:
for name, occ, idx in nl.map_occupancy_groups(s):
v = nl.clean_occupancy(occ) > 0
footprint.setdefault((s["subject"], s["day"], nl.epoch_type(name)), v)
rr, cc = np.where(v)
rows.append({"session": nl.label(s), "epoch": name, "maps": len(idx), "visited_bins": int(v.sum()),
"rows": f"{rr.min()}-{rr.max()}", "columns": f"{cc.min()}-{cc.max()}", "occupancy_counts": int(occ.sum())})
print(pd.DataFrame(rows).to_string(index=False))
jaccard = lambda a, b: round(float((a & b).sum() / (a | b).sum()), 2)
rows = []
for rat in ["Rat1", "Rat2", "Rat3"]:
fp = lambda day, typ: footprint[(rat, day, typ)]
for day in ["FD4", "FD9"]:
if (rat, day, "ES") in footprint:
rows.append({"rat": rat, "flight_day": day, "preflight ES / flight ES": jaccard(fp("preflight", "ES"), fp(day, "ES")),
"preflight MC / flight MC": jaccard(fp("preflight", "MC"), fp(day, "MC")),
"preflight ES / preflight MC": jaccard(fp("preflight", "ES"), fp("preflight", "MC")),
"flight ES / flight MC": jaccard(fp(day, "ES"), fp(day, "MC"))})
print("\nOverlap of visited bins (Jaccard index; for Rat3 FD4 the ES1 map):")
print(pd.DataFrame(rows).to_string(index=False))
Dandiset description, sentences on the track:
Rats were implanted with multi-electrode recording arrays and trained to traverse a three-dimensional track called the “Escher staircase,” in which three 90-degree yaw turns were interleaved with three 90-degree pitch turns, creating a closed loop after only 270 degrees of total yaw rotation
Rats received medial forebrain bundle stimulation as a reward for traversing the track.
Epoch names per file: {'Rat1 preflight': ['ES', 'MC'], 'Rat1 FD4': ['BL1', 'BL2', 'BL3', 'ES', 'MC'], 'Rat1 FD9': ['BL1', 'BL2', 'ES', 'MC'], 'Rat2 preflight': ['ES', 'MC'], 'Rat2 FD4': ['BL1', 'BL2', 'BL3', 'ES', 'MC'], 'Rat2 FD9': ['BL1', 'BL2', 'ES', 'MC'], 'Rat3 preflight': ['ES', 'MC'], 'Rat3 FD4': ['ES1', 'ES2', 'MC']}
session epoch maps visited_bins rows columns occupancy_counts
Rat1 preflight ES 27 1040 24-59 6-54 27895
Rat1 preflight MC 35 827 21-61 7-47 24774
Rat1 FD4 ES 33 1337 16-62 13-60 26967
Rat1 FD4 MC 32 1803 15-63 0-63 25739
Rat1 FD9 ES 26 2136 15-63 0-63 84975
Rat1 FD9 MC 38 2662 15-63 0-63 59340
Rat2 preflight ES 26 1342 23-58 6-56 33564
Rat2 preflight MC 33 993 20-60 8-50 34019
Rat2 FD4 ES 31 1854 16-63 4-63 53333
Rat2 FD4 MC 26 2184 15-63 0-63 34325
Rat2 FD9 ES 39 1686 15-63 9-59 123177
Rat2 FD9 MC 37 2421 15-63 0-63 46561
Rat3 preflight ES 45 1446 23-59 4-54 36835
Rat3 preflight MC 41 1057 21-63 7-49 23127
Rat3 FD4 ES1 19 1251 15-60 10-58 74454
Rat3 FD4 ES2 23 1121 15-57 15-58 65669
Rat3 FD4 MC 23 2062 15-63 0-63 24357
Overlap of visited bins (Jaccard index; for Rat3 FD4 the ES1 map):
rat flight_day preflight ES / flight ES preflight MC / flight MC preflight ES / preflight MC flight ES / flight MC
Rat1 FD4 0.27 0.27 0.31 0.42
Rat1 FD9 0.37 0.30 0.31 0.67
Rat2 FD4 0.47 0.37 0.46 0.61
Rat2 FD9 0.49 0.39 0.46 0.61
Rat3 FD4 0.45 0.35 0.46 0.43
The dandiset's description says the preflight recordings were made on a flat rectangular track, but the preflight files label their epochs ES and MC and hold ES and MC rate maps, as the flight files do; the draft's files describe MC as a flat two-dimensional track and say nothing of a different apparatus preflight. The preflight occupancy maps cover 827 to 1446 bins within rows 20 to 63 and columns 4 to 56; the flight maps reach row 15 or 16 and cover 1121 to 2662 bins. The visited bins of preflight ES and flight ES overlap with a Jaccard index of 0.27 to 0.49, below the overlap of flight ES with flight MC (0.42 to 0.67). The camera view on the ground and in the shuttle may differ, so the footprints cannot show which apparatus was used. What the preflight epochs labelled ES and MC were recorded on is not established; if the description is right, the comparisons with preflight below compare the Escher staircase in flight with another apparatus on the ground.
In [12]:
# For each flight unit and task epoch, the correlation of the unit's spike-count map with the
# spike-count prediction (rate x occupancy) of the rate map the dandiset links to it, and of every
# other map of the same tetrode and epoch.
rows = []
for s in sessions:
if "xy" not in s:
continue
lookup = nl.PositionLookup(s)
for a, b, typ, name in s["epochs"][["start", "stop", "type", "name"]].itertuples(index=False):
if typ == "BL":
continue
for ui, u in s["units"].iterrows():
sp = u.spikes[(u.spikes >= a) & (u.spikes <= b)]
found, r, c = lookup(sp)
if found.sum() < nl.MIN_SPIKES:
continue
counts = nl.count_map(r, c)
cand = s["maps"][(s["maps"].tetrode == u.tetrode) & (s["maps"].epoch_name == name)]
pred = {i: np.nan_to_num(s["rate_maps"][i]) * nl.clean_occupancy(s["occupancy_maps"][i]) for i in cand.index}
rs = {i: np.corrcoef(p.ravel(), counts.ravel())[0, 1] for i, p in pred.items()}
linked = list(cand.index[cand.unit == ui])
row = {"session": nl.label(s), "type": typ, "epoch": name, "tetrode": u.tetrode, "cluster": u.cluster, "spikes": int(found.sum()),
"maps_on_tetrode": len(cand), "has_linked_map": bool(linked)}
if linked:
i = linked[0]
rate_i, occ_i = s["rate_maps"][i], nl.clean_occupancy(s["occupancy_maps"][i])
positive = (np.nan_to_num(rate_i) > 0) & (occ_i > 0)
row |= {"r_linked": rs[i], "linked_is_best": max(rs, key=rs.get) == i,
"r_linked_x_y_swapped": np.corrcoef(pred[i].ravel(), counts.T.ravel())[0, 1],
"best_other_r": max([v for k, v in rs.items() if k != i], default=np.nan),
"positive_bins_without_spikes": float((positive & (counts == 0)).sum() / positive.sum()),
"rate_x_occupancy_per_spike": float(pred[i].sum() / found.sum())}
rows.append(row)
link = pd.DataFrame(rows)
print(f"Unit-epochs with at least {nl.MIN_SPIKES} spikes at a known position: {len(link)}; with a linked rate map: {int(link.has_linked_map.sum())}")
lk = link[link.has_linked_map]
print(lk.groupby(["session", "epoch"]).agg(n=("r_linked", "size"), median_r_linked=("r_linked", "median"), median_r_swapped=("r_linked_x_y_swapped", "median"),
median_best_other_r=("best_other_r", "median"), linked_is_best=("linked_is_best", "mean")).round(2).to_string())
print(f"\nAll: median r with the linked map {lk.r_linked.median():.2f}; the linked map is the best of its tetrode and epoch in "
f"{int(lk.linked_is_best.sum())} of {len(lk)}; r > 0.5 in {(lk.r_linked > 0.5).sum()}")
print("Unit-epochs without a linked map:"); print(link[~link.has_linked_map][["session", "epoch", "tetrode", "cluster", "spikes", "maps_on_tetrode"]].to_string(index=False))
good = lk[lk.r_linked > 0.5]
print(f"\nLinked maps with r > 0.5 ({len(good)}): fraction of bins with a positive rate where the unit has no spike, median "
f"{good.positive_bins_without_spikes.median():.2f}; rate x occupancy per spike, median "
f"{good.rate_x_occupancy_per_spike.median():.1f} (range {good.rate_x_occupancy_per_spike.min():.1f} to "
f"{good.rate_x_occupancy_per_spike.max():.1f})")
print(good.groupby("session").rate_x_occupancy_per_spike.median().round(1).to_string())
Unit-epochs with at least 50 spikes at a known position: 211; with a linked rate map: 211
n median_r_linked median_r_swapped median_best_other_r linked_is_best
session epoch
Rat1 FD4 ES 20 0.71 0.05 0.32 0.95
MC 18 0.38 0.01 0.22 0.94
Rat1 FD9 ES 25 0.76 0.07 0.21 1.00
MC 23 0.56 -0.01 0.29 1.00
Rat2 FD4 ES 17 0.47 0.05 0.24 1.00
MC 14 0.54 0.01 0.21 1.00
Rat2 FD9 ES 28 0.66 0.06 0.41 0.96
MC 11 0.48 0.01 0.18 1.00
Rat3 FD4 ES1 16 0.82 0.03 0.30 0.88
ES2 19 0.73 0.05 0.40 1.00
MC 20 0.18 -0.00 0.13 0.75
All: median r with the linked map 0.57; the linked map is the best of its tetrode and epoch in 201 of 211; r > 0.5 in 135
Unit-epochs without a linked map:
Empty DataFrame
Columns: [session, epoch, tetrode, cluster, spikes, maps_on_tetrode]
Index: []
Linked maps with r > 0.5 (135): fraction of bins with a positive rate where the unit has no spike, median 0.91; rate x occupancy per spike, median 22.0 (range 16.9 to 36.4)
session
Rat1 FD4 22.4
Rat1 FD9 21.6
Rat2 FD4 23.8
Rat2 FD9 22.0
Rat3 FD4 21.1
With position binned as x to columns and y to rows at 4 pixels per bin, as the rate-map table now
states, each of the 211 unit-epochs with 50 or more spikes at a known position has a linked rate
map, and that map is the best match among the maps of its tetrode and epoch for 201 of them (median
r 0.57 against a median of 0.13 to 0.41 for the best other map; -0.01 to 0.07 with x and y swapped).
So the rate maps are now linked to their units, as the table says, and were made from the same
spikes and positions. In 0.260728.1352 the cell_number column matched the cluster id for 84 of
211 unit-epochs.
For the 135 linked maps that correlate above 0.5, in the median map 91% of the bins with a positive rate hold no spike of the unit, so the maps are smoothed, as the description now says. Their summed rate x occupancy is about 22 per spike of the unit (median per file 21.1 to 23.8), so below it is used as a quantity proportional to a map's spike count. The preflight maps cannot be checked this way, as the preflight files have no position.
In [13]:
# Spike times: sorting, resolution and zero intervals by time range (the units' description says
# times have six significant digits), spikes outside each unit's obs_intervals, and the unit
# windows that differ from their epoch's.
ranges = [(0, 1000), (1000, 10000), (10000, np.inf)]
rows = []
for s in sessions:
allsp = np.concatenate(list(s["units"].spikes))
for lo, hi in ranges:
z = n = 0
for sp in s["units"].spikes:
w = sp[(sp >= lo) & (sp < hi)]
n += max(len(w) - 1, 0)
z += int((np.diff(w) == 0).sum())
distinct = np.unique(allsp[(allsp >= lo) & (allsp < hi)])
if len(distinct) > 1:
rows.append({"session": nl.label(s), "time_range_s": f"{lo}-{hi}", "spikes": int(((allsp >= lo) & (allsp < hi)).sum()),
"smallest_step_s": round(float(np.diff(distinct).min()), 4), "intervals": n, "zero_intervals": z,
"zero_fraction": round(z / max(n, 1), 3)})
print("Spike times are sorted in every unit:", all((np.diff(sp) >= 0).all() for s in sessions for sp in s["units"].spikes))
print(pd.DataFrame(rows).to_string(index=False))
rows = []
for s in sessions:
ep = {n: (a, b) for a, b, n in s["epochs"][["start", "stop", "name"]].itertuples(index=False)}
outside = differ = 0
for _, u in s["units"].iterrows():
inside = np.zeros(len(u.spikes), bool)
for (a, b), f in zip(u.obs_intervals, u.source_files):
inside |= (u.spikes >= a) & (u.spikes <= b)
differ += ep[f.split("/")[1].split(".")[0]] != (a, b)
outside += int((~inside).sum())
iv = sorted(zip(s["epochs"].start, s["epochs"].stop))
rows.append({"session": nl.label(s), "units": len(s["units"]), "spike_files": int(s["units"].source_files.map(len).sum()),
"spikes_outside_obs_intervals": outside, "unit_windows_differing_from_epoch": differ,
"overlapping_epochs": sum(1 for (a1, b1), (a2, b2) in zip(iv, iv[1:]) if a2 < b1),
"median_epochs_per_unit": float(s["units"].source_files.map(len).median())})
print("\nUnits, their spike files and windows:")
print(pd.DataFrame(rows).to_string(index=False))
s = by_label["Rat1 FD4"]
print("Rat1 FD4 unit windows in BL1 that differ from the epoch (2525-3390):",
Counter(tuple(w) for _, u in s["units"].iterrows() for w, f in zip(u.obs_intervals, u.source_files) if "/BL1." in f and tuple(w) != (2525.0, 3390.0)))
Spike times are sorted in every unit: True
session time_range_s spikes smallest_step_s intervals zero_intervals zero_fraction
Rat1 preflight 0-1000 13359 0.0010 13324 3 0.000
Rat1 preflight 1000-10000 8633 0.0100 8598 1034 0.120
Rat1 FD4 1000-10000 39971 0.0100 39924 3004 0.075
Rat1 FD9 1000-10000 238974 0.0100 238935 12204 0.051
Rat1 FD9 10000-inf 4543 0.1000 4510 2377 0.527
Rat2 preflight 1000-10000 36664 0.0100 36631 2884 0.079
Rat2 FD4 1000-10000 31766 0.0100 31694 4766 0.150
Rat2 FD9 1000-10000 51586 0.0100 51532 5080 0.099
Rat2 FD9 10000-inf 584 0.1000 556 197 0.354
Rat3 preflight 0-1000 30092 0.0010 30047 0 0.000
Rat3 preflight 1000-10000 17908 0.0100 17867 974 0.055
Rat3 FD4 0-1000 13451 0.0004 13432 0 0.000
Rat3 FD4 1000-10000 45265 0.0100 45219 2980 0.066
Units, their spike files and windows:
session units spike_files spikes_outside_obs_intervals unit_windows_differing_from_epoch overlapping_epochs median_epochs_per_unit
Rat1 preflight 35 62 0 0 0 2.0
Rat1 FD4 47 196 0 14 0 5.0
Rat1 FD9 39 139 0 0 0 4.0
Rat2 preflight 33 59 0 0 0 2.0
Rat2 FD4 72 194 0 0 0 1.0
Rat2 FD9 54 179 0 0 0 4.0
Rat3 preflight 48 87 0 0 0 2.0
Rat3 FD4 47 66 0 0 0 1.0
Rat1 FD4 unit windows in BL1 that differ from the epoch (2525-3390): Counter({(np.float64(2720.0), np.float64(3390.0)): 14})
Spike times are sorted. Their resolution follows the time, as the units' description says: 1 ms or
finer below 1000 s, 10 ms from 1000 s and 0.1 s from 10000 s. Zero intervals within a unit follow
the resolution: none or nearly none below 1000 s (0 to 3), 5 to 15% at 10 ms, and 35 to 53% at 0.1 s.
This explains the zero intervals the first attempt could not account for (Rat1's preflight MC epoch
crosses 1000 s). Within FD9 sessions, spikes after 10000 s (the last 151 s of each) are resolved
to 0.1 s only; of the task epochs only Rat2's FD9 MC epoch reaches it, and no claim is judged on it. No spike lies outside its unit's obs_intervals. Units
merge the per-epoch spike files of one cluster; 14 Rat1 FD4 units have their BL1 window at 2720 to
3390 s, inside the BL1 epoch (2525 to 3390 s), as the epochs' description explains. Every result here
except the stability claim is computed per epoch; the stability claim uses that a unit is the same
sorted cluster in both epochs, which the units' description now states.
In [14]:
# Clock and attribution contrast: spatial information of each rat's units in its own ES epoch and
# in the other rat's ES epoch of the same day. In each window the position comes from the file of
# the rat labelled as on the task, whose samples merge both rats' spike times, with that file's
# occupancy map for the epoch.
rng_check = np.random.default_rng(1754)
rows = []
for day in ["FD4", "FD9"]:
files = {rat: by_label[f"{rat} {day}"] for rat in ["Rat1", "Rat2"]}
for rat, other in [("Rat1", "Rat2"), ("Rat2", "Rat1")]:
for owner in (rat, other):
src = files[owner]
for a, b, typ, name in src["epochs"][["start", "stop", "type", "name"]].itertuples(index=False):
if typ != "ES":
continue
df = nl.epoch_units(files[rat], a, b, nl.epoch_occupancy(src, name), rng_check, n_shuffles=200, position=src)
d = df.dropna(subset=["info"])
rows.append({"day": day, "units_of": rat, "ES_epoch_of": owner, "window": f"{a:.0f}-{b:.0f}",
"units": len(d), "significant": int((d.p < 0.05).sum()),
"fraction": round(float((d.p < 0.05).mean()), 2) if len(d) else np.nan,
"binom_p_vs_0.05": stats.binomtest(int((d.p < 0.05).sum()), len(d), 0.05, alternative="greater").pvalue if len(d) else np.nan,
"median_info": round(float(d["info"].median()), 3) if len(d) else np.nan})
contrast = pd.DataFrame(rows)
print(contrast.to_string(index=False))
day units_of ES_epoch_of window units significant fraction binom_p_vs_0.05 median_info FD4 Rat1 Rat1 5935-6420 20 11 0.55 5.379601e-10 2.854 FD4 Rat1 Rat2 3850-4800 15 10 0.67 2.324454e-10 2.524 FD4 Rat2 Rat2 3850-4800 16 1 0.06 5.598733e-01 1.954 FD4 Rat2 Rat1 5935-6420 15 0 0.00 1.000000e+00 2.591 FD9 Rat1 Rat1 3178-4789 25 10 0.40 1.591912e-07 2.555 FD9 Rat1 Rat2 4871-7535 36 1 0.03 8.422208e-01 1.346 FD9 Rat2 Rat2 4871-7535 28 2 0.07 4.116914e-01 2.250 FD9 Rat2 Rat1 3178-4789 41 0 0.00 1.000000e+00 2.258
This contrast tests whether spikes and position share a clock and whose position the trace is. Each rat's units are tested (200 shuffles) in its own ES epoch and in the other rat's ES epoch of the same day, with the trace and occupancy of the file of the rat labelled as on the task. Rat1's units are selective in their own ES epochs (11 of 20 on FD4, 10 of 25 on FD9, against 5% expected by chance), so for Rat1 spikes and position share a clock. Rat2's units are no more selective than chance in their own ES epochs (1 of 16 on FD4, 2 of 28 on FD9) or in Rat1's (0 of 15, 0 of 41). Rat1's units are selective in Rat2's FD4 ES epoch (10 of 15), when Rat1 is labelled as on baseline, but not in Rat2's FD9 ES epoch (1 of 36). So, as in 0.260728.1352, a shared clock could not be established for Rat2's spikes and the position, and in Rat2's FD4 ES epoch the trace is consistent with Rat1's position, not Rat2's.
In [15]:
# Whose position is the trace in each task window? For every ES and MC window, the samples of the
# file of the rat on the task, the fraction of them at each rat's spike times, and each rat's
# spikes per second in the window.
rows = []
for day in ["FD4", "FD9"]:
f1, f2 = by_label[f"Rat1 {day}"], by_label[f"Rat2 {day}"]
spikes = {"Rat1": np.unique(np.concatenate(list(f1["units"].spikes))), "Rat2": np.unique(np.concatenate(list(f2["units"].spikes)))}
for owner, src in [("Rat1", f1), ("Rat2", f2)]:
for a, b, typ, name in src["epochs"][["start", "stop", "type", "name"]].itertuples(index=False):
if typ == "BL":
continue
w = (src["t"] >= a) & (src["t"] <= b)
tt = src["t"][w]
other_file = f2 if owner == "Rat1" else f1
row = {"day": day, "window": f"{owner} {name} {a:.0f}-{b:.0f}", "samples_per_s": round(w.sum() / (b - a), 1),
"other_file_samples": int(((other_file["t"] >= a) & (other_file["t"] <= b)).sum()),
"lost_fraction": round(float(1 - src["valid"][w].mean()), 3)}
for rat in ["Rat1", "Rat2"]:
row[f"at_{rat}_spikes"] = round(float(np.isin(tt, spikes[rat]).mean()), 3)
row[f"{rat}_spikes_per_s"] = round(((spikes[rat] >= a) & (spikes[rat] <= b)).sum() / (b - a), 1)
rows.append(row)
print(pd.DataFrame(rows).to_string(index=False))
day window samples_per_s other_file_samples lost_fraction at_Rat1_spikes Rat1_spikes_per_s at_Rat2_spikes Rat2_spikes_per_s FD4 Rat1 ES 5935-6420 17.8 0 0.050 0.613 10.9 0.434 7.7 FD4 Rat1 MC 7020-7550 22.2 0 0.166 0.809 18.0 0.236 5.2 FD4 Rat2 ES 3850-4800 18.2 0 0.149 0.500 9.1 0.548 10.0 FD4 Rat2 MC 7850-8473 15.7 0 0.103 0.662 10.4 0.380 6.0 FD9 Rat1 ES 3178-4789 40.7 0 0.038 0.877 35.7 0.197 8.0 FD9 Rat1 MC 7888-9200 49.9 0 0.138 0.958 47.9 0.090 4.5 FD9 Rat2 ES 4871-7535 25.3 0 0.027 0.751 19.0 0.305 7.7 FD9 Rat2 MC 9220-10151 24.4 0 0.151 0.892 21.8 0.154 3.8
In [16]:
# Movement of the trace in each task window (Rat3's for reference, as it had its own tracker):
# steps between consecutive valid samples, path length per second of the window, and the share of
# steps over 20 pixels. Then each of Rat1's and Rat2's units tested against the trace in its own and
# the other rat's MC windows, as for the ES windows above (200 shuffles).
rows = []
for s in sessions:
if "xy" not in s:
continue
for a, b, typ, name in s["epochs"][["start", "stop", "type", "name"]].itertuples(index=False):
if typ == "BL":
continue
w = (s["t"] >= a) & (s["t"] <= b) & s["valid"]
step = np.hypot(*np.diff(s["xy"][w], axis=0).T)
dt = np.diff(s["t"][w])
rows.append({"window": f"{nl.label(s)} {name} {a:.0f}-{b:.0f}", "valid_samples": int(w.sum()),
"path_px_per_s": round(float(step.sum() / (b - a)), 1), "median_step_px": round(float(np.median(step)), 2),
"steps_over_20px": round(float((step > 20).mean()), 4),
"px_per_s_over_gaps_0.5s_or_more": round(float(step[dt >= 0.5].sum() / dt[dt >= 0.5].sum()), 1) if (dt >= 0.5).any() else np.nan})
print(pd.DataFrame(rows).to_string(index=False))
rng_mc = np.random.default_rng(1755)
rows = []
for day in ["FD4", "FD9"]:
files = {rat: by_label[f"{rat} {day}"] for rat in ["Rat1", "Rat2"]}
for rat, other in [("Rat1", "Rat2"), ("Rat2", "Rat1")]:
for owner in (rat, other):
src = files[owner]
for a, b, typ, name in src["epochs"][["start", "stop", "type", "name"]].itertuples(index=False):
if typ != "MC":
continue
d = nl.epoch_units(files[rat], a, b, nl.epoch_occupancy(src, name), rng_mc, n_shuffles=200, position=src).dropna(subset=["info"])
k = int((d.p < 0.05).sum())
rows.append({"day": day, "units_of": rat, "MC_epoch_of": owner, "window": f"{a:.0f}-{b:.0f}", "units": len(d), "significant": k,
"fraction": round(k / len(d), 2) if len(d) else np.nan,
"binom_p_vs_0.05": stats.binomtest(k, len(d), 0.05, alternative="greater").pvalue if len(d) else np.nan})
print("\nEach rat's units against the trace in its own and the other rat's MC windows:")
print(pd.DataFrame(rows).to_string(index=False))
window valid_samples path_px_per_s median_step_px steps_over_20px px_per_s_over_gaps_0.5s_or_more Rat1 FD4 ES 5935-6420 8224 40.0 0.86 0.0100 27.1 Rat1 FD4 MC 7020-7550 9829 83.1 0.80 0.0339 40.5 Rat1 FD9 ES 3178-4789 63016 63.7 0.72 0.0065 25.2 Rat1 FD9 MC 7888-9200 56506 91.6 0.68 0.0133 20.3 Rat2 FD4 ES 3850-4800 14745 46.5 0.79 0.0193 22.8 Rat2 FD4 MC 7850-8473 8791 52.6 1.18 0.0263 32.9 Rat2 FD9 ES 4871-7535 65530 33.0 0.50 0.0058 23.7 Rat2 FD9 MC 9220-10151 19275 75.4 0.76 0.0282 42.4 Rat3 FD4 ES1 32-1345 15651 28.1 0.79 0.0114 13.0 Rat3 FD4 ES2 2397-3729 9670 20.5 0.98 0.0128 11.3 Rat3 FD4 MC 5713-6349 19458 72.6 0.64 0.0131 33.2
Each rat's units against the trace in its own and the other rat's MC windows: day units_of MC_epoch_of window units significant fraction binom_p_vs_0.05 FD4 Rat1 Rat1 7020-7550 18 2 0.11 0.226477 FD4 Rat1 Rat2 7850-8473 22 2 0.09 0.301849 FD4 Rat2 Rat2 7850-8473 14 0 0.00 1.000000 FD4 Rat2 Rat1 7020-7550 16 2 0.12 0.189240 FD9 Rat1 Rat1 7888-9200 23 0 0.00 1.000000 FD9 Rat1 Rat2 9220-10151 25 0 0.00 1.000000 FD9 Rat2 Rat2 9220-10151 11 1 0.09 0.431200 FD9 Rat2 Rat1 7888-9200 32 0 0.00 1.000000
Within each task window only the file of the rat on the task has samples, and they fall at both rats' spike times. Rat1's units fire more in Rat1's FD9 task windows (35.7 and 47.9 spikes per second in ES and MC, against 19.0 and 21.8 in Rat2's), while Rat2's units fire at about the same rate in Rat2's windows as in Rat1's (on FD9, 7.7 and 3.8 against 8.0 and 4.5; on FD4, 10.0 and 6.0 against 7.7 and 5.2). The files do not name the animal the tracker followed, and the draft says so; whose position the trace is in Rat2's task windows could not be established. The trace moves as a tracked animal would in every window: the median step between samples is 0.50 to 1.18 pixels, 0.6 to 3.4% of steps exceed 20 pixels, and the path length is 20.5 to 91.6 pixels per second (33.0 in Rat2's FD9 ES window, 63.7 in Rat1's, 28.1 and 20.5 in Rat3's two ES windows). It is a continuous track in Rat2's windows, but the movement does not show whose. In the MC windows neither rat's units are selective, in their own or the other rat's windows (0 to 2 of 11 to 32 units), so the MC windows cannot settle it either. The judged FD9 session of Rat2 is kept, as the criterion requires, with this caveat.
Findings of 001754-knierim-2000-01 on the draft, from the outputs above:
- The Rat1 and Rat2 files held the same trace, presented as each rat's position. Partly fixed. Each file now holds the trace only in its own rat's task epochs (no timestamp in common), and the position description states the shared tracker, the assumption that the trace in a rat's task epochs is that rat's, and that the files do not name the tracked animal. Rat2's units are still not selective in any window, and Rat1's are selective in Rat2's FD4 ES window (10 of 15), so whose position the trace is in Rat2's epochs is still not established.
- Position only at spike times, none preflight. Still present; the description now states it.
- The description places the preflight recordings on a flat rectangular track, while the files label them ES and MC. Still present; the dandiset description is unchanged.
cell_numberwas not a cell number, so rate maps could not be linked to units. Fixed: every map links to a unit with its tetrode and cluster, and the linked map is the best match for 201 of 211 unit-epochs.- The rate maps are smoothed with an unstated kernel. Partly fixed: the description states that they are smoothed and that the kernel is not recorded.
- -999 for unvisited bins and the lost position counted in bin [63, 0], both undocumented. Fixed: unvisited bins are NaN, lost positions are NaN, and the occupancy description states the bin [63, 0] count.
- Identical copies of rate maps, some under another tetrode. Fixed: no file holds two identical maps. The draft leaves out every copy in the six groups (14 maps), so those six maps, whose unit is ambiguous, are no longer in the dandiset.
- Units merge per-epoch spike files without saying a cluster id names one neuron across epochs, and
lack the spike-file columns. Fixed: the units' description says each tetrode was cut once, and
source_filesandobs_intervalsgive each unit's files and windows. - No EEG, though the description reports theta and ripples. Partly fixed: the files' experiment description says no EEG was archived; the dandiset description still reports EEG.
- Spike times at 10 ms with unexplained zero intervals, and two overlapping baseline epochs in Rat1 FD4. Fixed: the resolution is documented and the zero intervals follow it, and the epochs no longer overlap.
- The age is "/" and the pixel scale and camera view are not described. Partly fixed: the subject description and the position's reference frame now say the age, the camera view and the pixel size are not recorded, and the map geometry is given in pixels; the values themselves are still missing.
flight-spatial-code¶
Hippocampal neurons recorded on the space shuttle form a robust spatial code on the three-dimensional Escher staircase track in weightlessness. (Abstract)
Criterion (from the claims file): per in-flight ES epoch, a unit is significant when its spatial information exceeds the 95th percentile of 500 circular-shift shuffles. Judged on the two FD9 sessions: reproduced when the fraction of significant units is at least 0.5 in both.
In [17]:
rng_c1 = np.random.default_rng(20000301)
results = []
for s in sessions:
if "xy" not in s:
continue
for a, b, typ, name in s["epochs"][["start", "stop", "type", "name"]].itertuples(index=False):
if typ in ("ES", "MC"):
df = nl.epoch_units(s, a, b, nl.epoch_occupancy(s, name), rng_c1)
df["type"] = typ
df["epoch"] = f"{typ} {a:.0f}-{b:.0f}"
results.append(df)
units = pd.concat(results, ignore_index=True)
tested = units.dropna(subset=["info"]).copy()
tested["sig"] = tested.p < 0.05
def summary(d):
n, k = len(d), int(d.sig.sum())
return pd.Series({"units_in_file": np.nan, "tested": n, "significant": k, "fraction": round(k / n, 3) if n else np.nan,
"binom_p_vs_0.05": stats.binomtest(k, n, 0.05, alternative="greater").pvalue if n else np.nan,
"median_info": round(d["info"].median(), 3), "median_null95": round(d["null95"].median(), 3)})
per = tested.groupby(["type", "rat", "day", "epoch"]).apply(summary)
per["units_in_file"] = units.groupby(["type", "rat", "day", "epoch"]).size()
per["dropped_<50_spikes"] = per["units_in_file"] - per["tested"]
print(per.to_string())
units_in_file tested significant fraction binom_p_vs_0.05 median_info median_null95 dropped_<50_spikes
type rat day epoch
ES Rat1 FD4 ES 5935-6420 47 20.0 11.0 0.550 5.379601e-10 2.854 2.728 27.0
FD9 ES 3178-4789 39 25.0 10.0 0.400 1.591912e-07 2.555 2.603 14.0
Rat2 FD4 ES 3850-4800 72 16.0 1.0 0.062 5.598733e-01 1.954 2.420 56.0
FD9 ES 4871-7535 54 28.0 2.0 0.071 4.116914e-01 2.250 3.572 26.0
Rat3 FD4 ES 2397-3729 47 19.0 1.0 0.053 6.226464e-01 1.629 2.183 28.0
ES 32-1345 47 16.0 11.0 0.688 1.687114e-11 2.155 2.141 31.0
MC Rat1 FD4 MC 7020-7550 47 18.0 2.0 0.111 2.264774e-01 2.583 3.117 29.0
FD9 MC 7888-9200 39 23.0 0.0 0.000 1.000000e+00 2.478 4.058 16.0
Rat2 FD4 MC 7850-8473 72 14.0 1.0 0.071 5.123250e-01 2.597 3.228 58.0
FD9 MC 9220-10151 54 11.0 1.0 0.091 4.311999e-01 3.141 4.395 43.0
Rat3 FD4 MC 5713-6349 47 20.0 1.0 0.050 6.415141e-01 3.400 5.810 27.0
In [18]:
es = tested[tested.type == "ES"]
print("Judged sessions (FD9 ES):")
for rat in ["Rat1", "Rat2"]:
d = es[(es.rat == rat) & (es.day == "FD9")]
print(f" {rat} FD9: {d.sig.sum()} of {len(d)} significant, fraction {d.sig.mean():.3f}, "
f"binomial p {stats.binomtest(int(d.sig.sum()), len(d), 0.05, alternative='greater').pvalue:.2e}")
print("\nThreshold variants, fraction of significant ES units by session:")
var = []
for (rat, day, epoch), d in es.groupby(["rat", "day", "epoch"]):
row = {"session": f"{rat} {day} {epoch}", "n": len(d)}
for alpha in [0.05, 0.01]:
row[f"p<{alpha}"] = round(float((d.p < alpha).mean()), 3)
big = d[d.spikes_used >= 100]
row["n_>=100_spikes"] = len(big)
row["p<0.05,>=100_spikes"] = round(float((big.p < 0.05).mean()), 3) if len(big) else np.nan
row["info>0.5_bits"] = round(float((d["info"] > 0.5).mean()), 3)
var.append(row)
print(pd.DataFrame(var).to_string(index=False))
print("\nSpatial information of tested ES units, range:", round(es["info"].min(), 3), "to", round(es["info"].max(), 3), "bits/spike")
fig, ax = plt.subplots(figsize=(8, 3.5))
labels = []
for i, ((rat, day, epoch), d) in enumerate(es.groupby(["rat", "day", "epoch"])):
ax.scatter(np.full(len(d), i) + np.random.default_rng(i).uniform(-0.15, 0.15, len(d)), d["info"], s=10,
c=np.where(d.sig, "tab:blue", "tab:gray"))
ax.scatter(np.full(len(d), i) + 0.3, d["null95"], s=6, marker="_", c="k")
labels.append(f"{rat} {day}\n{epoch.split()[1]}")
ax.set_xticks(range(len(labels)), labels, fontsize=8)
ax.set_ylabel("spatial information (bits/spike)")
ax.set_title("ES units: blue significant, grey not; dashes: 95th percentile of shuffles")
plt.tight_layout()
plt.show()
Judged sessions (FD9 ES):
Rat1 FD9: 10 of 25 significant, fraction 0.400, binomial p 1.59e-07
Rat2 FD9: 2 of 28 significant, fraction 0.071, binomial p 4.12e-01
Threshold variants, fraction of significant ES units by session:
session n p<0.05 p<0.01 n_>=100_spikes p<0.05,>=100_spikes info>0.5_bits
Rat1 FD4 ES 5935-6420 20 0.550 0.350 16 0.438 1.00
Rat1 FD9 ES 3178-4789 25 0.400 0.320 24 0.417 0.92
Rat2 FD4 ES 3850-4800 16 0.062 0.000 13 0.000 1.00
Rat2 FD9 ES 4871-7535 28 0.071 0.071 26 0.077 1.00
Rat3 FD4 ES 2397-3729 19 0.053 0.000 19 0.053 1.00
Rat3 FD4 ES 32-1345 16 0.688 0.438 14 0.786 1.00
Spatial information of tested ES units, range: 0.276 to 6.07 bits/spike
In [19]:
# Per-unit results of the two judged sessions, with the Monte Carlo standard error of each p
# (500 shuffles) and the fraction of spikes dropped for lack of a position, real and shifted.
fd9 = es[es.day == "FD9"].copy()
fd9["mc_se"] = np.sqrt(fd9.p * (1 - fd9.p) / nl.N_SHUFFLES)
fd9["spikes_dropped"] = 1 - fd9.spikes_used / fd9.spikes_in_epoch
fd9["shuffle_spikes_dropped"] = 1 - fd9.null_spikes_used / fd9.spikes_in_epoch
cols = ["rat", "tetrode", "cluster", "spikes_in_epoch", "spikes_used", "spikes_dropped", "shuffle_spikes_dropped", "info", "null_mean", "null95", "p", "mc_se", "sig"]
print(fd9.sort_values(["rat", "p"])[cols].round(3).to_string(index=False))
for rat, d in fd9.groupby("rat"):
print(f"{rat} FD9: units with 0.01 < p < 0.1: {int(((d.p > 0.01) & (d.p < 0.1)).sum())}; "
f"within 2 Monte Carlo SE of 0.05: {int((np.abs(d.p - 0.05) < 2 * np.maximum(d.mc_se, np.sqrt(0.05 * 0.95 / nl.N_SHUFFLES))).sum())}")
print("\nRepeats of the judged sessions with 2000 shuffles and five other seeds:")
rows = []
for seed in range(5):
rng_rep = np.random.default_rng([20000301, seed])
for rat in ["Rat1", "Rat2"]:
s = by_label[f"{rat} FD9"]
a, b, name = s["epochs"][s["epochs"].type == "ES"][["start", "stop", "name"]].values[0]
d = nl.epoch_units(s, a, b, nl.epoch_occupancy(s, name), rng_rep, n_shuffles=2000).dropna(subset=["info"])
rows.append({"seed": seed, "rat": rat, "tested": len(d), "significant": int((d.p < 0.05).sum()),
"fraction": round(float((d.p < 0.05).mean()), 3)})
repeats = pd.DataFrame(rows)
print(repeats.pivot(index="seed", columns="rat", values=["significant", "fraction"]).to_string())
rat tetrode cluster spikes_in_epoch spikes_used spikes_dropped shuffle_spikes_dropped info null_mean null95 p mc_se sig Rat1 TT0 1 1714 1687 0.016 0.060 2.087 1.413 1.914 0.000 0.000 True Rat1 TT9 2 1180 1135 0.038 0.059 2.499 1.575 1.926 0.000 0.000 True Rat1 TT10 5 1735 1680 0.032 0.060 1.764 1.162 1.444 0.000 0.000 True Rat1 TT10 6 210 209 0.005 0.064 5.069 2.888 3.743 0.000 0.000 True Rat1 TT11 2 341 332 0.026 0.063 4.496 2.790 3.321 0.000 0.000 True Rat1 TT0 2 339 336 0.009 0.061 2.728 2.037 2.514 0.002 0.002 True Rat1 TT10 2 704 687 0.024 0.059 2.555 1.983 2.311 0.002 0.002 True Rat1 TT10 3 433 431 0.005 0.060 2.579 2.040 2.356 0.002 0.002 True Rat1 TT8 3 407 392 0.037 0.060 3.441 2.564 3.042 0.010 0.004 True Rat1 TT9 4 293 292 0.003 0.061 3.405 2.710 3.129 0.010 0.004 True Rat1 TT2 1 1100 1056 0.040 0.062 1.573 1.257 1.573 0.052 0.010 False Rat1 TT1 4 370 341 0.078 0.061 3.057 2.506 3.119 0.062 0.011 False Rat1 TT9 1 623 606 0.027 0.063 3.151 2.185 3.236 0.066 0.011 False Rat1 TT1 3 440 430 0.023 0.060 1.815 1.570 1.856 0.074 0.012 False Rat1 TT8 2 1237 1209 0.023 0.058 2.058 1.598 2.125 0.076 0.012 False Rat1 TT1 2 248 235 0.052 0.059 3.065 2.638 3.466 0.116 0.014 False Rat1 TT0 3 222 221 0.005 0.057 2.931 2.494 3.310 0.150 0.016 False Rat1 TT1 5 476 445 0.065 0.060 2.495 2.307 2.837 0.248 0.019 False Rat1 TT6 1 712 711 0.001 0.060 1.800 1.645 2.603 0.280 0.020 False Rat1 TT9 3 235 235 0.000 0.059 3.083 3.040 4.301 0.316 0.021 False Rat1 TT10 1 55 55 0.000 0.061 3.149 3.069 3.416 0.324 0.021 False Rat1 TT11 1 472 458 0.030 0.060 1.226 1.206 1.614 0.350 0.021 False Rat1 TT1 1 237 228 0.038 0.060 2.269 2.212 2.661 0.356 0.021 False Rat1 TT6 2 4837 4505 0.069 0.060 0.290 0.611 1.216 0.788 0.018 False Rat1 TT8 1 48565 46509 0.042 0.061 0.276 0.367 0.509 0.846 0.016 False Rat2 TT13 1 1040 1034 0.006 0.058 2.938 1.665 2.582 0.000 0.000 True Rat2 TT13 3 314 310 0.013 0.059 3.284 2.386 3.025 0.008 0.004 True Rat2 TT19 2 252 250 0.008 0.058 4.433 2.656 4.485 0.060 0.011 False Rat2 TT22 1 1353 1319 0.025 0.058 2.044 1.584 2.223 0.114 0.014 False Rat2 TT14 8 368 353 0.041 0.056 2.415 1.898 3.538 0.140 0.016 False Rat2 TT12 16 237 233 0.017 0.057 3.871 2.240 4.680 0.158 0.016 False Rat2 TT19 4 1965 1946 0.010 0.057 1.744 1.298 1.982 0.208 0.018 False Rat2 TT19 6 704 658 0.065 0.057 2.201 1.750 3.287 0.218 0.018 False Rat2 TT13 11 520 495 0.048 0.051 2.299 2.396 7.928 0.224 0.019 False Rat2 TT22 3 425 406 0.045 0.058 2.670 2.227 4.293 0.246 0.019 False Rat2 TT14 10 249 231 0.072 0.057 2.307 2.139 3.879 0.256 0.020 False Rat2 TT13 5 101 100 0.010 0.057 3.076 3.151 6.890 0.292 0.020 False Rat2 TT13 10 242 234 0.033 0.057 2.112 1.928 2.684 0.304 0.021 False Rat2 TT19 8 1467 1449 0.012 0.059 1.558 1.610 4.332 0.308 0.021 False Rat2 TT12 12 1224 1203 0.017 0.056 1.687 1.476 2.377 0.318 0.021 False Rat2 TT15 1 1003 991 0.012 0.058 1.770 1.622 2.765 0.344 0.021 False Rat2 TT13 8 186 162 0.129 0.061 3.577 3.605 5.833 0.368 0.022 False Rat2 TT12 13 2434 2407 0.011 0.058 1.518 1.374 2.654 0.378 0.022 False Rat2 TT13 9 193 190 0.016 0.061 2.468 2.522 4.156 0.392 0.022 False Rat2 TT19 5 51 51 0.000 0.058 3.222 3.368 4.562 0.532 0.022 False Rat2 TT19 1 170 153 0.100 0.056 2.498 2.672 3.647 0.546 0.022 False Rat2 TT12 15 2944 2854 0.031 0.058 1.212 1.241 1.774 0.588 0.022 False Rat2 TT12 14 553 530 0.042 0.059 1.786 2.018 3.061 0.608 0.022 False Rat2 TT15 2 309 295 0.045 0.058 1.486 2.059 3.606 0.672 0.021 False Rat2 TT14 7 5909 5778 0.022 0.058 0.787 1.017 1.550 0.696 0.021 False Rat2 TT15 3 143 141 0.014 0.056 1.877 2.528 3.912 0.826 0.017 False Rat2 TT19 7 205 198 0.034 0.058 1.355 1.778 2.366 0.906 0.013 False Rat2 TT14 9 53 52 0.019 0.057 2.588 3.877 6.382 0.956 0.009 False Rat1 FD9: units with 0.01 < p < 0.1: 5; within 2 Monte Carlo SE of 0.05: 3 Rat2 FD9: units with 0.01 < p < 0.1: 1; within 2 Monte Carlo SE of 0.05: 1 Repeats of the judged sessions with 2000 shuffles and five other seeds:
significant fraction rat Rat1 Rat2 Rat1 Rat2 seed 0 10.0 2.0 0.40 0.071 1 11.0 2.0 0.44 0.071 2 10.0 2.0 0.40 0.071 3 10.0 2.0 0.40 0.071 4 11.0 2.0 0.44 0.071
In [20]:
# Spikes dropped for lack of a position in every task epoch, and the count-matched variant, in which
# each shuffle is compared with the unit's information on a subsample of its own spikes of the same size.
tested["spikes_dropped"] = 1 - tested.spikes_used / tested.spikes_in_epoch
tested["shuffle_spikes_dropped"] = 1 - tested.null_spikes_used / tested.spikes_in_epoch
print(tested.groupby(["type", "rat", "day", "epoch"]).agg(units=("p", "size"), median_spikes_dropped=("spikes_dropped", "median"),
median_shuffle_spikes_dropped=("shuffle_spikes_dropped", "median"), max_shuffle_spikes_dropped=("shuffle_spikes_dropped", "max"),
empty_shuffles=("empty_shuffles", "sum"), median_info=("info", "median"), median_null_mean=("null_mean", "median")).round(3).to_string())
rng_cm = np.random.default_rng(20000302)
rows = []
for s in sessions:
if "xy" not in s:
continue
for a, b, typ, name in s["epochs"][["start", "stop", "type", "name"]].itertuples(index=False):
if typ != "ES":
continue
d = nl.epoch_units(s, a, b, nl.epoch_occupancy(s, name), rng_cm, match_counts=True).dropna(subset=["info"])
main_d = es[es.epoch == f"ES {a:.0f}-{b:.0f}"]
main_d = main_d[(main_d.rat == s["subject"]) & (main_d.day == s["day"])]
rows.append({"session": nl.label(s), "epoch": f"ES {a:.0f}-{b:.0f}", "tested": len(d), "significant_main": int(main_d.sig.sum()),
"significant_count_matched": int((d.p < 0.05).sum()), "fraction_count_matched": round(float((d.p < 0.05).mean()), 3),
"binom_p_vs_0.05": stats.binomtest(int((d.p < 0.05).sum()), len(d), 0.05, alternative="greater").pvalue})
print("\nCount-matched variant, ES epochs:")
print(pd.DataFrame(rows).to_string(index=False))
units median_spikes_dropped median_shuffle_spikes_dropped max_shuffle_spikes_dropped empty_shuffles median_info median_null_mean
type rat day epoch
ES Rat1 FD4 ES 5935-6420 20 0.046 0.084 0.088 0.0 2.854 2.307
FD9 ES 3178-4789 25 0.026 0.060 0.064 0.0 2.555 2.040
Rat2 FD4 ES 3850-4800 16 0.234 0.246 0.254 0.0 1.954 2.166
FD9 ES 4871-7535 28 0.021 0.058 0.061 0.0 2.250 2.039
Rat3 FD4 ES 2397-3729 19 0.058 0.215 0.221 0.0 1.629 1.646
ES 32-1345 16 0.040 0.205 0.208 0.0 2.155 1.667
MC Rat1 FD4 MC 7020-7550 18 0.206 0.281 0.289 0.0 2.583 2.556
FD9 MC 7888-9200 23 0.123 0.212 0.224 0.0 2.478 2.665
Rat2 FD4 MC 7850-8473 14 0.125 0.156 0.166 0.0 2.597 2.723
FD9 MC 9220-10151 11 0.133 0.155 0.162 0.0 3.141 3.288
Rat3 FD4 MC 5713-6349 20 0.084 0.242 0.256 0.0 3.400 3.954
Count-matched variant, ES epochs: session epoch tested significant_main significant_count_matched fraction_count_matched binom_p_vs_0.05 Rat1 FD4 ES 5935-6420 20 11 11 0.550 5.379601e-10 Rat1 FD9 ES 3178-4789 25 10 12 0.480 6.876528e-10 Rat2 FD4 ES 3850-4800 16 1 2 0.125 1.892403e-01 Rat2 FD9 ES 4871-7535 28 2 2 0.071 4.116914e-01 Rat3 FD4 ES 32-1345 16 11 11 0.688 1.687114e-11 Rat3 FD4 ES 2397-3729 19 1 1 0.053 6.226464e-01
In [21]:
# Calibration of the shuffle test at the real drop rates. In every flight ES epoch each unit's spikes
# are replaced by a train without spatial tuning and run through the same test (200 shuffles):
# (a) a Poisson train with the unit's mean count, uniform in the epoch; (b) the same, with each spike
# moved to the time of its nearest position sample, so that, as for real spikes, every spike has a
# sample at its own time and only shifted spikes are dropped for lack of one.
rng_cal = np.random.default_rng(20000303)
rows = []
for s in sessions:
if "xy" not in s:
continue
for a, b, typ, name in s["epochs"][["start", "stop", "type", "name"]].itertuples(index=False):
if typ != "ES":
continue
win = s["t"][(s["t"] >= a) & (s["t"] <= b)]
for kind in ["poisson", "poisson at sample times"]:
fake_spikes = []
for sp in s["units"].spikes:
t = np.sort(rng_cal.uniform(a, b, rng_cal.poisson(((sp >= a) & (sp <= b)).sum())))
if kind != "poisson":
i = np.clip(np.searchsorted(win, t), 1, len(win) - 1)
t = np.where(np.abs(win[i - 1] - t) <= np.abs(win[i] - t), win[i - 1], win[i])
fake_spikes.append(t)
fake = dict(s, units=s["units"].assign(spikes=fake_spikes))
d = nl.epoch_units(fake, a, b, nl.epoch_occupancy(s, name), rng_cal, n_shuffles=200).dropna(subset=["info"])
k = int((d.p < 0.05).sum())
ci = stats.binomtest(k, len(d)).proportion_ci() if len(d) else None
rows.append({"session": nl.label(s), "epoch": f"ES {a:.0f}-{b:.0f}", "null": kind, "tested": len(d), "false_positives": k,
"fraction": round(k / len(d), 3) if len(d) else np.nan,
"ci95": f"{ci.low:.3f}-{ci.high:.3f}" if ci else "",
"median_spikes_dropped": round(float((1 - d.spikes_used / d.spikes_in_epoch).median()), 3),
"median_shuffle_spikes_dropped": round(float((1 - d.null_spikes_used / d.spikes_in_epoch).median()), 3),
"median_null95": round(float(d.null95.median()), 3)})
cal = pd.DataFrame(rows)
print(cal.to_string(index=False))
pooled = {}
for kind, g in cal.groupby("null"):
k, n = int(g.false_positives.sum()), int(g.tested.sum())
pooled[kind] = k / n
print(f"{kind}: {k} of {n} units significant over all ES epochs ({k / n:.3f}); binomial p against 0.05: "
f"{stats.binomtest(k, n, 0.05, alternative='greater').pvalue:.3f}")
print("\nThe judged FD9 counts against the false-positive fraction of the same session and the pooled one (one-sided binomial p):")
for rat in ["Rat1", "Rat2"]:
d = es[(es.rat == rat) & (es.day == "FD9")]
k, n = int(d.sig.sum()), len(d)
for kind in pooled:
f_s = float(cal[(cal.session == f"{rat} FD9") & (cal.null == kind)].fraction.iloc[0])
print(f" {rat} FD9 {k} of {n}, null '{kind}': session {f_s:.3f}, p {stats.binomtest(k, n, f_s, alternative='greater').pvalue:.2g}; "
f"pooled {pooled[kind]:.3f}, p {stats.binomtest(k, n, pooled[kind], alternative='greater').pvalue:.2g}")
session epoch null tested false_positives fraction ci95 median_spikes_dropped median_shuffle_spikes_dropped median_null95 Rat1 FD4 ES 5935-6420 poisson 19 0 0.000 0.000-0.176 0.079 0.084 1.292 Rat1 FD4 ES 5935-6420 poisson at sample times 19 2 0.105 0.013-0.331 0.078 0.083 1.298 Rat1 FD9 ES 3178-4789 poisson 25 4 0.160 0.045-0.361 0.059 0.060 1.232 Rat1 FD9 ES 3178-4789 poisson at sample times 26 4 0.154 0.044-0.349 0.038 0.060 1.255 Rat2 FD4 ES 3850-4800 poisson 16 1 0.062 0.002-0.302 0.240 0.246 1.151 Rat2 FD4 ES 3850-4800 poisson at sample times 16 0 0.000 0.000-0.206 0.234 0.245 1.153 Rat2 FD9 ES 4871-7535 poisson 27 1 0.037 0.001-0.190 0.059 0.057 1.560 Rat2 FD9 ES 4871-7535 poisson at sample times 28 6 0.214 0.083-0.410 0.032 0.058 1.538 Rat3 FD4 ES 32-1345 poisson 16 1 0.062 0.002-0.302 0.211 0.206 1.005 Rat3 FD4 ES 32-1345 poisson at sample times 16 0 0.000 0.000-0.206 0.050 0.206 1.070 Rat3 FD4 ES 2397-3729 poisson 19 2 0.105 0.013-0.331 0.212 0.214 1.234 Rat3 FD4 ES 2397-3729 poisson at sample times 20 0 0.000 0.000-0.168 0.040 0.214 1.291 poisson: 9 of 122 units significant over all ES epochs (0.074); binomial p against 0.05: 0.158 poisson at sample times: 12 of 125 units significant over all ES epochs (0.096); binomial p against 0.05: 0.023 The judged FD9 counts against the false-positive fraction of the same session and the pooled one (one-sided binomial p): Rat1 FD9 10 of 25, null 'poisson': session 0.160, p 0.0035; pooled 0.074, p 5.5e-06 Rat1 FD9 10 of 25, null 'poisson at sample times': session 0.154, p 0.0026; pooled 0.096, p 5.6e-05 Rat2 FD9 2 of 28, null 'poisson': session 0.037, p 0.28; pooled 0.074, p 0.62 Rat2 FD9 2 of 28, null 'poisson at sample times': session 0.214, p 0.99; pooled 0.096, p 0.76
In [22]:
# The judged sessions on the earlier version's files (0.260728.1352), with which the first attempt
# ran: each rat's units and its own old file's position trace, in which a lost position is (0, 255).
# Same test (500 shuffles), the draft's occupancy maps, five seeds.
rows = []
for seed in range(5):
rng_old = np.random.default_rng([20000305, seed])
for rat in ["Rat1", "Rat2"]:
s = by_label[f"{rat} FD9"]
o = old_files[s["path"]]
xy = o["xy"].astype(float).copy()
lost = (xy[:, 0] == 0) & (xy[:, 1] == 255)
xy[lost] = np.nan
old_s = dict(s, t=o["t"], xy=xy, valid=~lost,
units=pd.DataFrame({"tetrode": [t for t, _ in o["units"]], "cluster": [int(c) for _, c in o["units"]], "spikes": o["spikes"]}))
a, b, name = s["epochs"][s["epochs"].type == "ES"][["start", "stop", "name"]].values[0]
d = nl.epoch_units(old_s, a, b, nl.epoch_occupancy(s, name), rng_old).dropna(subset=["info"])
rows.append({"seed": seed, "rat": rat, "tested": len(d), "significant": int((d.p < 0.05).sum()),
"median_spikes_dropped": round(float((1 - d.spikes_used / d.spikes_in_epoch).median()), 3),
"median_shuffle_spikes_dropped": round(float((1 - d.null_spikes_used / d.spikes_in_epoch).median()), 3)})
print(f"FD9 ES on {OLD_VERSION}:")
print(pd.DataFrame(rows).to_string(index=False))
FD9 ES on 0.260728.1352:
seed rat tested significant median_spikes_dropped median_shuffle_spikes_dropped
0 Rat1 25 12 0.026 0.074
0 Rat2 28 3 0.021 0.306
1 Rat1 25 11 0.026 0.074
1 Rat2 28 3 0.021 0.306
2 Rat1 25 11 0.026 0.074
2 Rat2 28 3 0.021 0.307
3 Rat1 25 11 0.026 0.075
3 Rat2 28 2 0.021 0.307
4 Rat1 25 12 0.026 0.074
4 Rat2 28 3 0.021 0.307
In [23]:
# Real spike trains decoupled from the position, to see what sets the upper tail of the shuffles.
# (a) In every flight ES epoch, each unit's spikes reversed in time within the epoch (t -> start +
# stop - t) and moved to the nearest position sample: the train keeps its bursts and slow changes in
# rate and loses its relation to the position. (b) In Rat3 FD4, each unit's ES1 spikes run against
# the ES2 trace and its ES2 spikes against the ES1 trace, at the same time from the epoch's start (the
# part beyond the shorter epoch is cut) and moved to the nearest sample. Same test, 200 shuffles.
# Temporal clumping of a train is the effective fraction of 10 s bins it fills, (sum n)^2 /
# sum(n^2) / bins: 1 for an even rate, small when the spikes fall in a few bins.
def to_samples(t, win):
i = np.clip(np.searchsorted(win, t), 1, len(win) - 1)
return np.sort(np.where(np.abs(win[i - 1] - t) <= np.abs(win[i] - t), win[i - 1], win[i]))
def spread(sp, a, b, w=10.0):
n = np.histogram(sp, np.arange(a, b + w, w))[0]
return n.sum() ** 2 / (n ** 2).sum() / len(n) if n.sum() else np.nan
rng_dec = np.random.default_rng(20000306)
cal_idx = {(r.session, r.epoch): r for r in cal[cal.null == "poisson at sample times"].itertuples()}
rows, dec = [], []
for s in sessions:
if "xy" not in s:
continue
for a, b, typ, name in s["epochs"][["start", "stop", "type", "name"]].itertuples(index=False):
if typ != "ES":
continue
win = s["t"][(s["t"] >= a) & (s["t"] <= b)]
rev = [to_samples(a + b - sp[(sp >= a) & (sp <= b)], win) for sp in s["units"].spikes]
d = nl.epoch_units(dict(s, units=s["units"].assign(spikes=rev)), a, b, nl.epoch_occupancy(s, name), rng_dec, n_shuffles=200)
epoch = f"ES {a:.0f}-{b:.0f}"
real = tested[(tested.rat == s["subject"]) & (tested.day == s["day"]) & (tested.epoch == epoch)].set_index(["tetrode", "cluster"])
real = real.assign(spread=[spread(sp[(sp >= a) & (sp <= b)], a, b) for sp in s["units"].set_index(["tetrode", "cluster"]).loc[real.index].spikes])
real["null95_reversed"] = d.set_index(["tetrode", "cluster"]).null95.reindex(real.index)
dec.append(real.reset_index().assign(session=nl.label(s), epoch=epoch))
k, n = int((d.p < 0.05).sum()), int(d.p.notna().sum())
rows.append({"session": nl.label(s), "epoch": epoch, "tested": len(real), "significant": int(real.sig.sum()),
"reversed_tested": n, "reversed_significant": k,
"median_null95_real": round(float(real.null95.median()), 3), "median_null95_reversed": round(float(d.null95.median()), 3),
"median_null95_poisson_at_samples": round(float(cal_idx[(nl.label(s), epoch)].median_null95), 3),
"median_spikes_used": float(real.spikes_used.median()), "median_spread": round(float(real.spread.median()), 3)})
print("(a) Time-reversed real trains, and the median 95th percentile of the shuffles for real, reversed and Poisson trains:")
dec_table = pd.DataFrame(rows)
print(dec_table.to_string(index=False))
print(f"Reversed trains: {dec_table.reversed_significant.sum()} of {dec_table.reversed_tested.sum()} units significant over all ES epochs")
r3 = by_label["Rat3 FD4"]
es3 = r3["epochs"][r3["epochs"].type == "ES"].sort_values("start")
e = es3[["start", "stop"]].values
length = min(e[0, 1] - e[0, 0], e[1, 1] - e[1, 0])
rows = []
for src, dst in [(0, 1), (1, 0)]:
win = r3["t"][(r3["t"] >= e[dst, 0]) & (r3["t"] <= e[dst, 1])]
moved = []
for sp in r3["units"].spikes:
x = sp[(sp >= e[src, 0]) & (sp <= e[src, 1])] - e[src, 0]
moved.append(to_samples(e[dst, 0] + x[x < length], win))
d = nl.epoch_units(dict(r3, units=r3["units"].assign(spikes=moved)), e[dst, 0], e[dst, 1],
nl.epoch_occupancy(r3, es3.name.values[dst]), rng_dec, n_shuffles=200).dropna(subset=["info"])
rows.append({"spikes_of": es3.name.values[src], "trace_of": es3.name.values[dst], "tested": len(d),
"significant": int((d.p < 0.05).sum()), "median_null95": round(float(d.null95.median()), 3)})
print("\n(b) Rat3 FD4, each epoch's spikes against the other epoch's trace:")
print(pd.DataFrame(rows).to_string(index=False))
D = pd.concat(dec, ignore_index=True).dropna(subset=["null95_reversed"])
X = np.column_stack([np.ones(len(D)), np.log(D.spikes_used), np.log(D.spread)])
coef, *_ = np.linalg.lstsq(X, np.log(D.null95), rcond=None)
r2 = 1 - np.var(np.log(D.null95) - X @ coef) / np.var(np.log(D.null95))
print(f"\nOver {len(D)} tested flight ES units: Spearman rho of the real null's 95th percentile with the reversed train's "
f"{stats.spearmanr(D.null95, D.null95_reversed).statistic:.3f}, with spikes used {stats.spearmanr(D.null95, D.spikes_used).statistic:.3f}, "
f"with spread {stats.spearmanr(D.null95, D.spread).statistic:.3f}")
print(f"log null95 = {coef[0]:.2f} + {coef[1]:.3f} log(spikes used) + {coef[2]:.3f} log(spread), R^2 {r2:.2f}")
print("\nUnits named in the round-2 verification (spread and 95th percentile in each epoch):")
print(D[D.set_index(["session", "tetrode", "cluster"]).index.isin([("Rat3 FD4", "TT4", 1), ("Rat3 FD4", "TT5", 1), ("Rat3 FD4", "TT3", 2), ("Rat2 FD9", "TT13", 11)])]
[["session", "epoch", "tetrode", "cluster", "spikes_used", "spread", "info", "null95", "null95_reversed", "p"]].round(3).to_string(index=False))
(a) Time-reversed real trains, and the median 95th percentile of the shuffles for real, reversed and Poisson trains: session epoch tested significant reversed_tested reversed_significant median_null95_real median_null95_reversed median_null95_poisson_at_samples median_spikes_used median_spread Rat1 FD4 ES 5935-6420 20 11 19 0 2.728 2.670 1.298 203.0 0.278 Rat1 FD9 ES 3178-4789 25 10 25 1 2.603 2.535 1.255 431.0 0.190 Rat2 FD4 ES 3850-4800 16 1 16 0 2.420 2.501 1.153 315.0 0.358 Rat2 FD9 ES 4871-7535 28 2 27 0 3.572 3.459 1.538 331.5 0.145 Rat3 FD4 ES 32-1345 16 11 16 0 2.141 2.239 1.070 301.0 0.222 Rat3 FD4 ES 2397-3729 19 1 19 1 2.183 2.378 1.291 455.0 0.307 Reversed trains: 2 of 122 units significant over all ES epochs
(b) Rat3 FD4, each epoch's spikes against the other epoch's trace:
spikes_of trace_of tested significant median_null95
ES1 ES2 16 0 3.081
ES2 ES1 19 0 1.600
Over 121 tested flight ES units: Spearman rho of the real null's 95th percentile with the reversed train's 0.995, with spikes used -0.601, with spread -0.874
log null95 = 0.97 + -0.139 log(spikes used) + -0.496 log(spread), R^2 0.80
Units named in the round-2 verification (spread and 95th percentile in each epoch):
session epoch tetrode cluster spikes_used spread info null95 null95_reversed p
Rat2 FD9 ES 4871-7535 TT13 11 495 0.064 2.299 7.928 8.237 0.224
Rat3 FD4 ES 32-1345 TT3 2 118 0.214 2.947 2.573 2.744 0.000
Rat3 FD4 ES 32-1345 TT4 1 1582 0.267 2.449 1.326 1.528 0.000
Rat3 FD4 ES 32-1345 TT5 1 8844 0.575 0.910 0.504 0.603 0.000
Rat3 FD4 ES 2397-3729 TT3 2 157 0.211 2.127 2.777 2.861 0.634
Rat3 FD4 ES 2397-3729 TT4 1 423 0.096 3.592 4.305 4.416 0.160
Rat3 FD4 ES 2397-3729 TT5 1 2221 0.307 1.472 2.393 2.458 0.338
In [24]:
print("Attribution contrast from the data check, ES epochs:")
print(contrast.to_string(index=False))
Attribution contrast from the data check, ES epochs: day units_of ES_epoch_of window units significant fraction binom_p_vs_0.05 median_info FD4 Rat1 Rat1 5935-6420 20 11 0.55 5.379601e-10 2.854 FD4 Rat1 Rat2 3850-4800 15 10 0.67 2.324454e-10 2.524 FD4 Rat2 Rat2 3850-4800 16 1 0.06 5.598733e-01 1.954 FD4 Rat2 Rat1 5935-6420 15 0 0.00 1.000000e+00 2.591 FD9 Rat1 Rat1 3178-4789 25 10 0.40 1.591912e-07 2.555 FD9 Rat1 Rat2 4871-7535 36 1 0.03 8.422208e-01 1.346 FD9 Rat2 Rat2 4871-7535 28 2 0.07 4.116914e-01 2.250 FD9 Rat2 Rat1 3178-4789 41 0 0.00 1.000000e+00 2.258
Result: not reproduced. On flight day 9, 10 of 25 Rat1 units (0.40) have significant spatial information on the Escher staircase (binomial p = 1.6e-7 against 0.05), and 2 of 28 Rat2 units (0.071, p = 0.41). By the criterion, Rat1 is above chance but below 0.5 and Rat2 is not above chance, so neither the rule for reproduced nor the rule for partially reproduced is met. Repeats with 2000 shuffles and five other seeds give 10 or 11 of 25 for Rat1 (0.40 to 0.44) and 2 of 28 for Rat2; the run of 200 shuffles in the data check gives 10 of 25 and 2 of 28; 5 Rat1 units have p between 0.01 and 0.1, 3 within two Monte Carlo standard errors of 0.05. On the files of 0.260728.1352 the same test gives 11 or 12 of 25 for Rat1 and 2 or 3 of 28 for Rat2 over five seeds; the first attempt recorded 12 of 25 (11 to 13 over seven runs). Every real spike has a sample at its own time in both versions; the draft's trace adds samples at the other rat's spike times in each window, so a shifted spike finds a sample within 0.25 s more often and the shuffle distributions change. Shifted spikes are now dropped about as often in Rat2's FD9 ES epoch (median 0.058) as in Rat1's (0.060), against 0.306 to 0.307 and 0.074 to 0.075 on the earlier version's files, so the drops no longer differ between the rats; the count-matched variant gives 12 of 25 (0.48) for Rat1 and 2 of 28 for Rat2, still below 0.5. The test was calibrated at the real drop rates by replacing each unit's spikes in every flight ES epoch with a Poisson train of the same mean count: 9 of 122 units are significant over all ES epochs (0.074, binomial p = 0.16 against 0.05), and 12 of 125 (0.096, p = 0.023) when each spike is moved to its nearest position sample, so that shifted spikes are dropped more often than the observed ones, as in the real data. Per session the false-positive fraction is 0 to 0.21; in the judged sessions it is 0.160 and 0.154 for Rat1 and 0.037 and 0.214 for Rat2. The criterion's chance level of 0.05 is thus below the false-positive fraction of some sessions. Against each session's own fractions, Rat1's 10 of 25 is still above chance (p = 0.0035 and 0.0026) and Rat2's 2 of 28 is not (p = 0.28 and 0.99), so the outcome is the same.
Poisson trains are not a good model of the real trains, whose spikes come in bursts and in periods of high and low rate. When each unit's real train is reversed in time within its epoch, so that it keeps this structure and loses its relation to the position, 2 of 122 units are significant over all flight ES epochs (1 of 25 in Rat1's FD9 session, 0 of 27 in Rat2's), and Rat3's ES1 trains against the ES2 trace and the reverse give 0 of 16 and 0 of 19. So with real temporal structure the test is conservative, not lenient. The 95th percentile of a unit's shuffles is set by its train: it is nearly the same for the reversed train (Spearman rho 0.995 over 121 units) and 1.7 to 2.3 times that of a Poisson train of the same count in each session (medians 2.141 to 3.572 against 1.070 to 1.538). It is higher the fewer 10 s periods hold a unit's spikes (rho -0.874 with the spread, the effective fraction of 10 s bins filled) and the fewer spikes it has (rho -0.601); the two explain 80% of its variance on a log scale. Rat2's FD9 units are the most clumped (median spread 0.145, against 0.190 for Rat1 FD9 and 0.222 to 0.358 elsewhere) and have the highest median 95th percentile (3.572); TT13 cluster 11 has a spread of 0.064 and a 95th percentile of 7.928. A unit whose spikes fall in a few periods needs much more information to pass, so the test has less power for such units, and this power was not measured. A low fraction of significant units, Rat2's especially, is therefore weak evidence against a spatial code. The criterion's rule still decides the outcome. Rat2's units are not selective with the trace in any window, and Rat1's are selective in Rat2's FD4 ES epoch (10 of 15), so Rat2's result may come from a trace that is not Rat2's position; the files do not show whose position it is. On flight day 4, 11 of 20 Rat1 units, 1 of 16 Rat2 units, and 11 of 16 and 1 of 19 Rat3 units in its ES1 and ES2 epochs are significant. In the MC epochs, which the claim does not cover, 0 to 2 units per session are significant.
three-orthogonal-surfaces¶
The spatial code covers each of the three orthogonal surfaces of the Escher staircase. (Abstract)
Blocked. The dandiset holds one two-dimensional pixel position from one camera; the NASA record of the study (OSD-968) says position was recorded with two synchronised cameras (printed below). The draft's position reference frame says the camera view and the pixel size on the apparatus are not recorded in the source files (data check), nothing in the files or the dandiset metadata describes the track layout in the image, and the conversion code describes the ES track only as a three-dimensional track with 90 degree yaw and pitch turns and the position only as video tracker pixel coordinates (lines printed below), so a pixel cannot be assigned to one of the three surfaces.
In [25]:
# The NASA record of the study (OSD-968): what it says of the cameras, its description of the files,
# and its file listing. Then every line of the conversion code that describes the track or the position.
osd_text = json.dumps(requests.get("https://osdr.nasa.gov/osdr/data/osd/meta/968", timeout=60).json(), ensure_ascii=False)
for word in ["camera", "Raw (Recording) Data Files", "EEG", "theta", "ripple", "LFP"]:
hits = [m.group(0) for m in re.finditer(r"[^.\"]{0,150}" + re.escape(word) + r"[^\"]{0,450}", osd_text)]
print(f"OSD-968 metadata, '{word}': {len(hits)} passage(s)")
for h in hits[:1]:
print(" ", h[:600])
files = requests.get("https://osdr.nasa.gov/osdr/data/osd/files/968", timeout=60).json()["studies"]["OSD-968"]["study_files"]
print(f"\nOSD-968 files ({len(files)}):")
for f in files:
print(f" {f['file_name']} {f['file_size'] / 1e6:.1f} MB")
repo, commit = "catalystneuro/mcnaughton2000-to-nwb", "734688253c17a32e26d84c57f41fbe48ffbcd009"
tree = requests.get(f"https://api.github.com/repos/{repo}/git/trees/{commit}?recursive=1", timeout=60).json()["tree"]
paths = [t["path"] for t in tree if t["type"] == "blob" and t["path"].endswith((".py", ".md"))]
print(f"\nConversion code at {commit[:7]}, every line on the track, camera, surfaces, coordinates or waveforms:")
for path in paths:
text = requests.get(f"https://raw.githubusercontent.com/{repo}/{commit}/{path}", timeout=60).text
for n, line in enumerate(text.splitlines(), 1):
if re.search(r"camera|surface|staircase|layout|geometry|reference_frame|coordinate|waveform", line, re.I):
print(f" {path}:{n}: {line.strip()[:160]}")
OSD-968 metadata, 'camera': 1 passage(s)
Position and orientation were recorded using two synchronized cameras integrated with the Data Acquisition System (DAS). Neural activity from chronically implanted tetrodes in the hippocampus was recorded simultaneously during behavioral navigation.
OSD-968 metadata, 'Raw (Recording) Data Files': 1 passage(s)
Raw (Recording) Data Files: These contain ASCII .CEL / .cell spike files with headers and spike event tables. These were also converted to .xls format (“Converted Raw Files”). The “Processed Data Files” compressed folder contains Analyzed Place Field Outputs as binary .RMA map files (fixed-length 32 KB). Additional legacy metadata files (FileMaker .fp5, .mcd, and Excel manifests) provided by the original PI team are also included for provenance. Please refer to the “READM
OSD-968 metadata, 'EEG': 0 passage(s)
OSD-968 metadata, 'theta': 0 passage(s)
OSD-968 metadata, 'ripple': 0 passage(s)
OSD-968 metadata, 'LFP': 0 passage(s)
OSD-968 files (27): OSD-968_metadata_OSD-968-ISA.zip 0.0 MB LSDS-173_Extracellular-Electrophysiology-Recording_Original_Raw_PREFLI3.zip 1.5 MB LSDS-173_Extracellular-Electrophysiology-Recording_Original_Raw_PREFLI2.zip 1.1 MB LSDS-173_Extracellular-Electrophysiology-Recording_Original_Raw_PREFLI1.zip 0.8 MB LSDS-173_Extracellular-Electrophysiology-Recording_Original_Raw_FD9RAT2.zip 2.1 MB LSDS-173_Extracellular-Electrophysiology-Recording_Original_Raw_FD9RAT1.zip 8.6 MB LSDS-173_Extracellular-Electrophysiology-Recording_Original_Raw_FD4RAT3.zip 2.1 MB LSDS-173_Extracellular-Electrophysiology-Recording_Original_Raw_FD4RAT2.zip 1.4 MB LSDS-173_Extracellular-Electrophysiology-Recording_Original_Raw_FD4RAT1.zip 1.7 MB LSDS-173_Extracellular-Electrophysiology-Recording_Original Processed Files.zip 272.3 MB LSDS-173_Extracellular-Electrophysiology-Recording_PREFLI3.nwb.zip 0.7 MB LSDS-173_Extracellular-Electrophysiology-Recording_PREFLI2.nwb.zip 0.5 MB LSDS-173_Extracellular-Electrophysiology-Recording_PREFLI1.nwb.zip 0.4 MB LSDS-173_Extracellular-Electrophysiology-Recording_NWB_Conversion_Notes.docx 0.0 MB LSDS-173_Extracellular-Electrophysiology-Recording_FD9RAT2.nwb.zip 0.8 MB LSDS-173_Extracellular-Electrophysiology-Recording_FD9RAT1.nwb.zip 1.9 MB LSDS-173_Extracellular-Electrophysiology-Recording_FD4RAT3.nwb.zip 0.7 MB LSDS-173_Extracellular-Electrophysiology-Recording_FD4RAT2.nwb.zip 0.6 MB LSDS-173_Extracellular-Electrophysiology-Recording_FD4RAT1.nwb.zip 0.6 MB LSDS-173_Extracellular-Electrophysiology-Recording_Converted_Raw_PREFLI3.zip 2.4 MB LSDS-173_Extracellular-Electrophysiology-Recording_Converted_Raw_PREFLI2.zip 3.0 MB LSDS-173_Extracellular-Electrophysiology-Recording_Converted_Raw_PREFLI1.zip 1.2 MB LSDS-173_Extracellular-Electrophysiology-Recording_Converted_Raw_FD9RAT2.zip 3.6 MB LSDS-173_Extracellular-Electrophysiology-Recording_Converted_Raw_FD9RAT1.zip 14.2 MB LSDS-173_Extracellular-Electrophysiology-Recording_Converted_Raw_FD4RAT3.zip 3.7 MB LSDS-173_Extracellular-Electrophysiology-Recording_Converted_Raw_FD4RAT2.zip 2.6 MB LSDS-173_Extracellular-Electrophysiology-Recording_Converted_Raw_FD4RAT1.zip 3.3 MB Conversion code at 7346882, every line on the track, camera, surfaces, coordinates or waveforms:
README.md:7: This converts data from the Neurolab Space Shuttle mission (April 1998), in which hippocampal CA1 place cells were recorded from three freely moving rats in mic README.md:48: - **Epochs**: one interval per original epoch name (`BL1`, `ES`, `ES2`, `MC`, ...), with its task type: Escher Staircase, Magic Carpet, or Baseline (rest in the README.md:49: - **Position**: video tracker pixel coordinates at spike times, within the rat's ES and MC epochs (flight sessions only) conversion_notes.md:28: - Fields: id, t_px, t_py, t_pa, t_pb (spike waveform peak amplitudes on 4 channels), t_mx, t_my, t_ma, t_mb (integrated waveform areas), t_maxwd, t_maxht (wavef conversion_notes.md:53: - **ES**: Escher Staircase — 3D track with 90° yaw/pitch turns conversion_notes.md:73: - Camera view and pixel size on the apparatus. OSD-968 says two synchronized cameras recorded position; the spike files hold one x, y pair.
src/mcnaughton2000_to_nwb/convert_session.py:4: with interleaved task epochs (BL=Baseline, ES=Escher Staircase, MC=Magic Carpet).
src/mcnaughton2000_to_nwb/convert_session.py:234: epochs.add_column("session_type", "Task type: BL (Baseline), ES (Escher Staircase), or MC (Magic Carpet)")
src/mcnaughton2000_to_nwb/convert_session.py:362: reference_frame=(
src/mcnaughton2000_to_nwb/convert_session.py:363: "Pixel coordinates (0 to 255) of the video tracker image, as written in the original spike files. "
src/mcnaughton2000_to_nwb/convert_session.py:364: "The camera view and the pixel size on the apparatus are not recorded in the source files."
src/mcnaughton2000_to_nwb/metadata.py:29: "to traverse a three-dimensional track ('Escher staircase') and a flat "
src/mcnaughton2000_to_nwb/metadata.py:35: "turns: while one rat ran the Escher staircase or Magic Carpet, the other "
src/mcnaughton2000_to_nwb/metadata.py:49: "ES": "Escher Staircase — three-dimensional track with 90° yaw and pitch turns",
tests/test_convert_session.py:57: # Rat1 rests (BL1) from 10 to 20 s and runs the staircase (ES) from 30 to 40 s;
si-comparable-to-preflight¶
In flight, spatial information content of place cells is comparable to preflight levels. (Dandiset description)
Criterion: two-sided Mann-Whitney U test at alpha 0.05 of the spatial information of the dandiset's rate maps, FD9 ES against preflight ES, Rats 1 and 2 pooled, every distinct map with a positive mean rate. Reproduced when the test finds no difference.
In [26]:
m = maps[(~maps.duplicate) & (maps.mean_rate > 0)].copy()
print("Maps dropped: duplicates", int(maps.duplicate.sum()), "; zero mean rate", int(((~maps.duplicate) & (maps.mean_rate <= 0)).sum()),
"; kept", len(m), "of", len(maps))
print("Spatial information of rate maps, range:", round(m["info"].min(), 3), "to", round(m["info"].max(), 3), "bits/spike")
print(m.groupby(["type", "rat", "day"])["info"].describe()[["count", "50%", "25%", "75%", "min", "max"]].round(3).to_string())
def compare(a, b, alternative="two-sided"):
if len(a) == 0 or len(b) == 0:
return {"n_a": len(a), "n_b": len(b)}
u = stats.mannwhitneyu(a, b, alternative=alternative)
return {"n_a": len(a), "n_b": len(b), "median_a": round(float(np.median(a)), 3), "median_b": round(float(np.median(b)), 3),
"diff": round(float(np.median(a) - np.median(b)), 3), "ratio": round(float(np.median(a) / np.median(b)), 2),
"U": float(u.statistic), "p": float(u.pvalue), "rank_biserial": round(1 - 2 * u.statistic / (len(a) * len(b)), 3)}
sel = lambda d, typ, days, rats: d[(d.type == typ) & d.day.isin(days) & d.rat.isin(rats)]["info"].values
main = compare(sel(m, "ES", ["FD9"], ["Rat1", "Rat2"]), sel(m, "ES", ["preflight"], ["Rat1", "Rat2"]))
print("\nMain test, ES FD9 (a) against ES preflight (b), Rats 1 and 2:", main)
Maps dropped: duplicates 0 ; zero mean rate 6 ; kept 528 of 534
Spatial information of rate maps, range: 0.019 to 7.583 bits/spike
count 50% 25% 75% min max
type rat day
ES Rat1 FD4 33.0 1.222 0.755 2.203 0.256 3.431
FD9 26.0 1.749 1.094 1.978 0.142 2.627
preflight 27.0 1.869 1.122 2.515 0.709 3.070
Rat2 FD4 31.0 0.872 0.614 1.975 0.200 7.583
FD9 38.0 1.343 1.059 1.841 0.487 2.821
preflight 26.0 1.321 0.954 1.556 0.111 3.583
Rat3 FD4 40.0 1.245 0.844 1.783 0.164 3.670
preflight 45.0 1.413 0.844 2.156 0.065 3.350
MC Rat1 FD4 31.0 0.704 0.496 1.236 0.140 3.930
FD9 38.0 1.216 0.916 1.393 0.019 2.589
preflight 35.0 0.918 0.692 1.682 0.387 4.051
Rat2 FD4 25.0 1.178 0.672 1.707 0.498 3.411
FD9 36.0 0.945 0.801 1.377 0.420 4.773
preflight 33.0 1.263 0.849 1.571 0.081 2.756
Rat3 FD4 23.0 1.558 0.656 1.809 0.024 6.497
preflight 41.0 1.216 0.763 1.718 0.063 3.068
Main test, ES FD9 (a) against ES preflight (b), Rats 1 and 2: {'n_a': 64, 'n_b': 53, 'median_a': 1.385, 'median_b': 1.448, 'diff': -0.064, 'ratio': 0.96, 'U': 1559.0, 'p': 0.4548208416095775, 'rank_biserial': np.float64(0.081)}
In [27]:
rows = []
variants = {
"per rat Rat1": (m, "ES", ["FD9"], ["Rat1"], ["preflight"], ["Rat1"]),
"per rat Rat2": (m, "ES", ["FD9"], ["Rat2"], ["preflight"], ["Rat2"]),
"CEL names only": (m[m.cel_name], "ES", ["FD9"], ["Rat1", "Rat2"], ["preflight"], ["Rat1", "Rat2"]),
"peak >= 1 Hz": (m[m.peak_rate >= 1], "ES", ["FD9"], ["Rat1", "Rat2"], ["preflight"], ["Rat1", "Rat2"]),
"MC maps": (m, "MC", ["FD9"], ["Rat1", "Rat2"], ["preflight"], ["Rat1", "Rat2"]),
"all flight ES (FD4+FD9, 3 rats)": (m, "ES", ["FD4", "FD9"], ["Rat1", "Rat2", "Rat3"], ["preflight"], ["Rat1", "Rat2", "Rat3"]),
}
for name, (d, typ, da, ra, db, rb) in variants.items():
a = d[(d.type == typ) & d.day.isin(da) & d.rat.isin(ra)]["info"].values
b = d[(d.type == typ) & d.day.isin(db) & d.rat.isin(rb)]["info"].values
rows.append({"variant": name} | compare(a, b))
print("Variants (a = flight, b = preflight):")
print(pd.DataFrame(rows).to_string(index=False))
fig, ax = plt.subplots(figsize=(7, 3.5))
groups = [(f"{rat}\n{day}", m[(m.type == "ES") & (m.rat == rat) & (m.day == day)]["info"].values)
for rat in ["Rat1", "Rat2", "Rat3"] for day in ["preflight", "FD4", "FD9"]]
groups = [(g, v) for g, v in groups if len(v)]
ax.boxplot([v for _, v in groups], showfliers=True)
ax.set_xticks(range(1, len(groups) + 1), [g for g, _ in groups], fontsize=8)
ax.set_ylabel("spatial information of ES rate maps (bits/spike)")
plt.tight_layout()
plt.show()
Variants (a = flight, b = preflight):
variant n_a n_b median_a median_b diff ratio U p rank_biserial
per rat Rat1 26 27 1.749 1.869 -0.120 0.94 286.0 0.251140 0.185
per rat Rat2 38 26 1.343 1.321 0.021 1.02 506.0 0.875087 -0.024
CEL names only 45 42 1.371 1.394 -0.022 0.98 851.0 0.427079 0.099
peak >= 1 Hz 38 48 1.481 1.463 0.019 1.01 895.0 0.885908 0.019
MC maps 74 68 1.082 1.153 -0.071 0.94 2490.0 0.917063 0.010
all flight ES (FD4+FD9, 3 rats) 168 98 1.314 1.441 -0.127 0.91 7299.0 0.123391 0.113
In [28]:
# Sampling bias of spatial information: epoch durations and spike counts per ES epoch (all files,
# from the units), the rate maps' total occupancy and rate x occupancy (proportional to the spike
# count where the frame rate is the same; see the data check), and two variants that limit the bias.
rows = []
for s in sessions:
for a, b, typ in s["epochs"][["start", "stop", "type"]].itertuples(index=False):
if typ != "ES":
continue
n = np.array([int(((sp >= a) & (sp <= b)).sum()) for sp in s["units"].spikes])
rows.append({"session": nl.label(s), "epoch": f"ES {a:.0f}-{b:.0f}", "duration_s": b - a, "units_with_spikes": int((n > 0).sum()),
"units_>=50_spikes": int((n >= 50).sum()), "median_spikes": float(np.median(n[n > 0])),
"median_rate_hz": round(float(np.median(n[n > 0]) / (b - a)), 2)})
print("Unit spike counts in each ES epoch (all units with a spike in the epoch):")
print(pd.DataFrame(rows).to_string(index=False))
mes = m[m.type == "ES"].copy()
print("\nES rate maps per session:")
print(mes.groupby(["rat", "day"]).agg(maps=("info", "size"), median_occupancy_total=("occupancy_total", "median"),
median_rate_x_occupancy=("rate_x_occupancy", "median"), median_info=("info", "median"),
median_info_occ5=("info_occ5", "median")).round(3).to_string())
x = np.log(mes.rate_x_occupancy)
print(f"\nSpearman rho of info with log(rate x occupancy), all ES maps: {stats.spearmanr(x, mes['info']).statistic:.3f} (n {len(mes)})")
print("Within each session:", {f"{r} {d}": round(float(stats.spearmanr(np.log(g.rate_x_occupancy), g['info']).statistic), 3) for (r, d), g in mes.groupby(["rat", "day"])})
fit = np.polyfit(x, mes["info"], 1)
print(f"Linear fit over all ES maps: info = {fit[1]:.3f} + {fit[0]:.3f} x log(rate x occupancy)")
m["info_adj"] = m["info"] - np.polyval(fit, np.log(m.rate_x_occupancy))
rows = []
rows.append({"variant": "info over bins with occupancy >= 5"} | compare(sel(m.assign(info=m.info_occ5), "ES", ["FD9"], ["Rat1", "Rat2"]), sel(m.assign(info=m.info_occ5), "ES", ["preflight"], ["Rat1", "Rat2"])))
rows.append({"variant": "count-adjusted info (residual of the fit)"} | compare(sel(m.assign(info=m.info_adj), "ES", ["FD9"], ["Rat1", "Rat2"]), sel(m.assign(info=m.info_adj), "ES", ["preflight"], ["Rat1", "Rat2"])))
print("\nBias variants, FD9 ES (a) against preflight ES (b), Rats 1 and 2, two-sided:")
print(pd.DataFrame(rows).to_string(index=False))
Unit spike counts in each ES epoch (all units with a spike in the epoch):
session epoch duration_s units_with_spikes units_>=50_spikes median_spikes median_rate_hz
Rat1 preflight ES 276-741 465.0 27 26 359.0 0.77
Rat1 FD4 ES 5935-6420 485.0 35 20 74.0 0.15
Rat1 FD9 ES 3178-4789 1611.0 26 25 436.5 0.27
Rat2 preflight ES 1013-1573 560.0 26 21 351.5 0.63
Rat2 FD4 ES 3850-4800 950.0 39 19 44.0 0.05
Rat2 FD9 ES 4871-7535 2664.0 41 28 205.0 0.08
Rat3 preflight ES 148-762 614.0 45 38 228.0 0.37
Rat3 FD4 ES 32-1345 1313.0 19 16 238.0 0.18
Rat3 FD4 ES 2397-3729 1332.0 23 20 449.0 0.34
ES rate maps per session:
maps median_occupancy_total median_rate_x_occupancy median_info median_info_occ5
rat day
Rat1 FD4 33 24849.0 1891.158 1.222 1.208
FD9 26 81138.0 9104.746 1.749 1.749
preflight 27 27733.0 7861.768 1.869 1.865
Rat2 FD4 31 43977.0 1818.619 0.872 0.864
FD9 38 117841.0 5226.288 1.343 1.339
preflight 26 33326.0 7627.513 1.321 1.314
Rat3 FD4 40 62327.0 7232.501 1.245 1.236
preflight 45 36660.0 5165.623 1.413 1.416
Spearman rho of info with log(rate x occupancy), all ES maps: -0.226 (n 266)
Within each session: {'Rat1 FD4': -0.274, 'Rat1 FD9': -0.283, 'Rat1 preflight': -0.349, 'Rat2 FD4': -0.792, 'Rat2 FD9': -0.26, 'Rat2 preflight': -0.229, 'Rat3 FD4': -0.118, 'Rat3 preflight': 0.042}
Linear fit over all ES maps: info = 4.172 + -0.315 x log(rate x occupancy)
Bias variants, FD9 ES (a) against preflight ES (b), Rats 1 and 2, two-sided:
variant n_a n_b median_a median_b diff ratio U p rank_biserial
info over bins with occupancy >= 5 64 53 1.383 1.449 -0.065 0.95 1556.0 0.444969 0.083
count-adjusted info (residual of the fit) 64 53 0.101 0.150 -0.049 0.67 1532.0 0.370659 0.097
In [29]:
# Does the rate maps' spatial information separate the units that pass the spike-based shuffle test
# from those that do not? Each tested flight ES unit joined to its rate map (same rat, day, epoch
# name, tetrode and cluster). Then a bootstrap 95% interval for the main test's median difference.
name_of = {(s["subject"], s["day"], f"{t} {a:.0f}-{b:.0f}"): n for s in sessions
for a, b, t, n in s["epochs"][["start", "stop", "type", "name"]].itertuples(index=False)}
j = tested[tested.type == "ES"].copy()
j["epoch_name"] = [name_of[k] for k in zip(j.rat, j.day, j.epoch)]
j = j.merge(m[["rat", "day", "epoch_name", "tetrode", "cluster", "info"]].rename(columns={"info": "map_info"}),
on=["rat", "day", "epoch_name", "tetrode", "cluster"], how="left")
print(f"Tested flight ES units: {len(j)}; with a rate map: {int(j.map_info.notna().sum())}")
rows = []
for key, g in list(j.groupby(["rat", "day", "epoch_name"])) + [(("all", "", ""), j)]:
g = g.dropna(subset=["map_info"])
a, b = g[g.sig].map_info.values, g[~g.sig].map_info.values
rows.append({"session": " ".join(key).strip(), "significant": len(a), "not_significant": len(b),
"median_map_info_sig": round(float(np.median(a)), 3) if len(a) else np.nan,
"median_map_info_not_sig": round(float(np.median(b)), 3) if len(b) else np.nan,
"mann_whitney_p": stats.mannwhitneyu(a, b).pvalue if len(a) and len(b) else np.nan})
print(pd.DataFrame(rows).to_string(index=False))
jj = j.dropna(subset=["map_info"])
print(f"Spearman rho of rate-map information with spike-based information: {stats.spearmanr(jj.map_info, jj['info']).statistic:.3f}; "
f"with spike-based information minus the shuffle mean: {stats.spearmanr(jj.map_info, jj['info'] - jj.null_mean).statistic:.3f} (n {len(jj)})")
rng_bs = np.random.default_rng(20000304)
for rats in (["Rat1", "Rat2"], ["Rat1"]):
a, b = sel(m, "ES", ["FD9"], rats), sel(m, "ES", ["preflight"], rats)
diffs = np.array([np.median(rng_bs.choice(a, len(a))) - np.median(rng_bs.choice(b, len(b))) for _ in range(10000)])
lo, hi = np.percentile(diffs, [2.5, 97.5])
print(f"FD9 minus preflight median, {'+'.join(rats)}: {np.median(a) - np.median(b):.3f} bits/spike, bootstrap 95% interval "
f"{lo:.3f} to {hi:.3f} ({lo / np.median(b):+.0%} to {hi / np.median(b):+.0%} of the preflight median)")
Tested flight ES units: 124; with a rate map: 124
session significant not_significant median_map_info_sig median_map_info_not_sig mann_whitney_p
Rat1 FD4 ES 11 9 1.924 0.752 4.938356e-03
Rat1 FD9 ES 10 15 1.919 1.100 9.898474e-03
Rat2 FD4 ES 1 15 1.829 0.614 2.500000e-01
Rat2 FD9 ES 2 26 2.178 1.283 4.761905e-02
Rat3 FD4 ES1 11 5 1.677 0.555 1.451465e-01
Rat3 FD4 ES2 1 18 1.248 1.228 9.473684e-01
all 36 88 1.855 1.076 8.967003e-07
Spearman rho of rate-map information with spike-based information: 0.791; with spike-based information minus the shuffle mean: 0.744 (n 124)
FD9 minus preflight median, Rat1+Rat2: -0.064 bits/spike, bootstrap 95% interval -0.425 to 0.312 (-29% to +22% of the preflight median)
FD9 minus preflight median, Rat1: -0.120 bits/spike, bootstrap 95% interval -0.750 to 0.543 (-40% to +29% of the preflight median)
Result: reproduced. Flight day 9 ES maps (64, Rats 1 and 2) have a median of 1.385 bits/spike against 1.448 for preflight ES maps (53); Mann-Whitney p = 0.45, a difference of -0.064 bits/spike (ratio 0.96, rank-biserial 0.081). The main test's maps and numbers are the same as in 0.260728.1352. Every variant gives p of 0.12 or more: per rat (0.25 and 0.88), maps with archived names containing CEL (0.43), maps with a peak of 1 Hz or more (0.89), MC maps (0.92), all flight maps (0.12), bins with 5 or more occupancy counts (0.44), and spatial information adjusted for rate x occupancy (0.37), which is proportional to the spike count in the flight maps (about 22 per spike, data check) and is assumed to be so in the preflight maps, where it cannot be checked. Spatial information falls as the spike count rises, as sampling bias predicts (Spearman rho -0.226 over all ES maps; -0.118 to -0.792 within seven files, 0.042 in Rat3's preflight). Rat1's compared sessions have similar spike counts and Rat2's do not: units in Rat1's and Rat2's preflight ES epochs have medians of 359 and 352 spikes in 465 and 560 s, and on FD9 437 and 205 in 1611 and 2664 s; the maps' median rate x occupancy is 7862 and 7628 preflight against 9105 and 5226 on FD9. Rat2's FD9 maps thus rest on fewer spikes than its preflight maps, which by sampling bias raises their spatial information. The rate maps are smoothed with an unrecorded kernel, so no shuffle-corrected value can be computed for them. The dandiset's description says the preflight recordings were made on a flat rectangular track, so this may be a comparison between different apparatus, and Rat2's preflight file shares no tetrode with its flight files (see the data check). Rat2's FD9 maps (median 1.343 bits/spike) rest on a trace that could not be attributed to Rat2, while only 2 of its 28 units are significant against shuffles, so the rate maps' spatial information does not show that Rat2's cells were selective. Over all 124 tested flight ES units, the rate maps of units significant against their shuffles have a higher median than the others (1.855 against 1.076 bits/spike, Mann-Whitney p = 9.0e-7; Spearman rho 0.791 with spike-based information), so the measure does track selectivity, but Rat2's 26 FD9 units that are not significant have a median of 1.283, close to the preflight median of 1.448. The bootstrap 95% interval of the FD9 minus preflight median difference is -0.425 to 0.312 bits/spike (-29% to +22% of the preflight median); for Rat1 alone, whose trace is attributed, it is -0.750 to 0.543 (p = 0.25). The criterion's test has no equivalence margin, and a drop of up to 29% of the preflight median (40% for Rat1 alone) is not excluded, so the result shows that no difference was detected.
fd4-abnormal-selectivity¶
On the first in-flight exposures (flight day 4) spatial selectivity was abnormal. (Dandiset description)
Criterion: abnormal is read as lower spatial information. FD4 ES rate maps (Rats 1 to 3) against preflight ES rate maps (Rats 1 to 3); reproduced when the FD4 median is lower and a one-sided Mann-Whitney U test gives p < 0.05.
In [30]:
fd4 = sel(m, "ES", ["FD4"], ["Rat1", "Rat2", "Rat3"])
pre = sel(m, "ES", ["preflight"], ["Rat1", "Rat2", "Rat3"])
main4 = compare(fd4, pre, alternative="less")
print("Main test, ES FD4 (a) below ES preflight (b), Rats 1 to 3, one-sided:", main4)
rows = []
for rat in ["Rat1", "Rat2", "Rat3"]:
rows.append({"comparison": f"{rat} FD4 < preflight"} | compare(sel(m, "ES", ["FD4"], [rat]), sel(m, "ES", ["preflight"], [rat]), "less"))
for rat in ["Rat1", "Rat2"]:
rows.append({"comparison": f"{rat} FD4 < FD9"} | compare(sel(m, "ES", ["FD4"], [rat]), sel(m, "ES", ["FD9"], [rat]), "less"))
rows.append({"comparison": "Rats 1-2 FD4 < FD9"} | compare(sel(m, "ES", ["FD4"], ["Rat1", "Rat2"]), sel(m, "ES", ["FD9"], ["Rat1", "Rat2"]), "less"))
rows.append({"comparison": "CEL names, FD4 < preflight"} | compare(sel(m[m.cel_name], "ES", ["FD4"], ["Rat1", "Rat2", "Rat3"]), sel(m[m.cel_name], "ES", ["preflight"], ["Rat1", "Rat2", "Rat3"]), "less"))
rows.append({"comparison": "peak >= 1 Hz, FD4 < preflight"} | compare(sel(m[m.peak_rate >= 1], "ES", ["FD4"], ["Rat1", "Rat2", "Rat3"]), sel(m[m.peak_rate >= 1], "ES", ["preflight"], ["Rat1", "Rat2", "Rat3"]), "less"))
rows.append({"comparison": "MC maps, FD4 < preflight"} | compare(sel(m, "MC", ["FD4"], ["Rat1", "Rat2", "Rat3"]), sel(m, "MC", ["preflight"], ["Rat1", "Rat2", "Rat3"]), "less"))
print(pd.DataFrame(rows).to_string(index=False))
print("\nSpike-based fraction of significant ES units by day (from flight-spatial-code):")
print(es.groupby(["day", "rat"]).agg(tested=("sig", "size"), significant=("sig", "sum"), fraction=("sig", "mean"),
median_info=("info", "median")).round(3).to_string())
k4, n4 = int(es[es.day == "FD4"].sig.sum()), int((es.day == "FD4").sum())
k9, n9 = int(es[es.day == "FD9"].sig.sum()), int((es.day == "FD9").sum())
print(f"FD4 {k4}/{n4} against FD9 {k9}/{n9}; Fisher exact one-sided (FD4 lower) p = "
f"{stats.fisher_exact([[k4, n4 - k4], [k9, n9 - k9]], alternative='less')[1]:.3g}")
print("\nVariants for the spatial-information bias (see si-comparable-to-preflight):")
rows = []
rows.append({"comparison": "info over bins with occupancy >= 5, FD4 < preflight"} | compare(sel(m.assign(info=m.info_occ5), "ES", ["FD4"], ["Rat1", "Rat2", "Rat3"]), sel(m.assign(info=m.info_occ5), "ES", ["preflight"], ["Rat1", "Rat2", "Rat3"]), "less"))
rows.append({"comparison": "count-adjusted info, FD4 < preflight"} | compare(sel(m.assign(info=m.info_adj), "ES", ["FD4"], ["Rat1", "Rat2", "Rat3"]), sel(m.assign(info=m.info_adj), "ES", ["preflight"], ["Rat1", "Rat2", "Rat3"]), "less"))
corr_es = es.assign(info=es["info"] - es.null_mean)
rows.append({"comparison": "spike-based info minus shuffle mean, Rats 1-2 FD4 < FD9"} | compare(corr_es[(corr_es.day == "FD4") & corr_es.rat.isin(["Rat1", "Rat2"])]["info"].values, corr_es[corr_es.day == "FD9"]["info"].values, "less"))
rows.append({"comparison": "spike-based info minus shuffle mean, all FD4 < FD9"} | compare(corr_es[corr_es.day == "FD4"]["info"].values, corr_es[corr_es.day == "FD9"]["info"].values, "less"))
print(pd.DataFrame(rows).to_string(index=False))
Main test, ES FD4 (a) below ES preflight (b), Rats 1 to 3, one-sided: {'n_a': 104, 'n_b': 98, 'median_a': 1.177, 'median_b': 1.441, 'diff': -0.264, 'ratio': 0.82, 'U': 4271.0, 'p': 0.023535842534057817, 'rank_biserial': np.float64(0.162)}
comparison n_a n_b median_a median_b diff ratio U p rank_biserial
Rat1 FD4 < preflight 33 27 1.222 1.869 -0.647 0.65 320.0 0.031630 0.282
Rat2 FD4 < preflight 31 26 0.872 1.321 -0.449 0.66 329.0 0.119475 0.184
Rat3 FD4 < preflight 40 45 1.245 1.413 -0.168 0.88 796.0 0.181077 0.116
Rat1 FD4 < FD9 33 26 1.222 1.749 -0.527 0.70 397.0 0.315283 0.075
Rat2 FD4 < FD9 31 38 0.872 1.343 -0.470 0.65 453.0 0.051066 0.231
Rats 1-2 FD4 < FD9 64 64 1.081 1.385 -0.304 0.78 1773.0 0.095411 0.134
CEL names, FD4 < preflight 72 76 1.177 1.288 -0.111 0.91 2357.0 0.073240 0.139
peak >= 1 Hz, FD4 < preflight 62 81 1.201 1.518 -0.317 0.79 1935.0 0.009531 0.229
MC maps, FD4 < preflight 79 109 0.866 1.179 -0.313 0.73 3731.0 0.059541 0.133
Spike-based fraction of significant ES units by day (from flight-spatial-code):
tested significant fraction median_info
day rat
FD4 Rat1 20 11 0.550 2.854
Rat2 16 1 0.062 1.954
Rat3 35 12 0.343 1.833
FD9 Rat1 25 10 0.400 2.555
Rat2 28 2 0.071 2.250
FD4 24/71 against FD9 12/53; Fisher exact one-sided (FD4 lower) p = 0.941
Variants for the spatial-information bias (see si-comparable-to-preflight):
comparison n_a n_b median_a median_b diff ratio U p rank_biserial
info over bins with occupancy >= 5, FD4 < preflight 104 98 1.179 1.442 -0.263 0.82 4264.0 0.022615 0.163
count-adjusted info, FD4 < preflight 104 98 -0.275 0.105 -0.380 -2.61 3819.0 0.001055 0.251
spike-based info minus shuffle mean, Rats 1-2 FD4 < FD9 36 53 0.025 0.211 -0.186 0.12 816.0 0.125189 0.145
spike-based info minus shuffle mean, all FD4 < FD9 71 53 0.232 0.211 0.021 1.10 1783.0 0.310304 0.052
In [31]:
# The same test on the rate maps of the earlier version (0.260728.1352), each distinct map once, as
# the earlier attempt counted them: the draft leaves out every copy in the six groups of identical maps.
rows = []
for s in sessions:
o = old_files[s["path"]]
in_draft = {np.where(np.isnan(r), -999, r).astype(np.float32).tobytes() for r in s["rate_maps"]}
seen = set()
for r_map, occ, typ, tt, f in zip(o["maps"], o["occ"], o["map_types"], o["map_tetrodes"], o["map_files"]):
if r_map.tobytes() in seen:
continue
seen.add(r_map.tobytes())
info, mean, _ = nl.rma_info(np.where(r_map == -999, np.nan, r_map), occ)
rows.append({"rat": s["subject"], "day": s["day"], "type": typ, "map": f"{tt}/{f}", "info": info, "mean_rate": mean,
"in_draft": r_map.tobytes() in in_draft})
old_maps = pd.DataFrame(rows)
old_maps = old_maps[old_maps.mean_rate > 0]
print(f"Distinct maps with a positive mean rate: {len(old_maps)} in {OLD_VERSION}, {len(m)} in the draft")
print("Distinct maps of the earlier version not in the draft:")
print(old_maps[~old_maps.in_draft].round(3).to_string(index=False))
print("ES maps per rat and day, earlier version and draft:")
print(pd.concat({OLD_VERSION: old_maps[old_maps.type == "ES"].groupby(["rat", "day"]).size(),
"draft": m[m.type == "ES"].groupby(["rat", "day"]).size()}, axis=1).to_string())
rats3 = ["Rat1", "Rat2", "Rat3"]
rows = [{"maps": f"{OLD_VERSION}, all distinct"} | compare(sel(old_maps, "ES", ["FD4"], rats3), sel(old_maps, "ES", ["preflight"], rats3), "less"),
{"maps": f"{OLD_VERSION}, those in the draft"} | compare(sel(old_maps[old_maps.in_draft], "ES", ["FD4"], rats3), sel(old_maps[old_maps.in_draft], "ES", ["preflight"], rats3), "less"),
{"maps": "draft"} | main4]
print("\nMain test (FD4 ES below preflight ES, Rats 1 to 3, one-sided) on each set of maps:")
print(pd.DataFrame(rows).to_string(index=False))
Distinct maps with a positive mean rate: 532 in 0.260728.1352, 528 in the draft
Distinct maps of the earlier version not in the draft:
rat day type map info mean_rate in_draft
Rat2 FD4 ES TT13/ES2BC4~1.RMA 7.583 0.0 False
Rat2 FD4 ES TT13/ES3BCA~1.RMA 7.583 0.0 False
Rat2 FD4 ES TT19/ES3BC8~1.RMA 7.583 0.0 False
Rat2 FD9 MC TT13/MC2DC6~1.RMA 4.773 0.0 False
ES maps per rat and day, earlier version and draft:
0.260728.1352 draft
rat day
Rat1 FD4 33 33
FD9 26 26
preflight 27 27
Rat2 FD4 34 31
FD9 38 38
preflight 26 26
Rat3 FD4 40 40
preflight 45 45
Main test (FD4 ES below preflight ES, Rats 1 to 3, one-sided) on each set of maps:
maps n_a n_b median_a median_b diff ratio U p rank_biserial
0.260728.1352, all distinct 107 98 1.214 1.441 -0.227 0.84 4565.0 0.055150 0.129
0.260728.1352, those in the draft 104 98 1.177 1.441 -0.264 0.82 4271.0 0.023536 0.162
draft 104 98 1.177 1.441 -0.264 0.82 4271.0 0.023536 0.162
In [32]:
# The main test without Rat2, whose position trace could not be attributed (data check).
rows = [{"comparison": "Rats 1 and 3, FD4 < preflight"} | compare(sel(m, "ES", ["FD4"], ["Rat1", "Rat3"]), sel(m, "ES", ["preflight"], ["Rat1", "Rat3"]), "less"),
{"comparison": "Rats 1 and 3, FD4 < preflight, old version, all distinct maps"} | compare(sel(old_maps, "ES", ["FD4"], ["Rat1", "Rat3"]), sel(old_maps, "ES", ["preflight"], ["Rat1", "Rat3"]), "less"),
{"comparison": "Rat1, FD9 against preflight (two-sided)"} | compare(sel(m, "ES", ["FD9"], ["Rat1"]), sel(m, "ES", ["preflight"], ["Rat1"]))]
print(pd.DataFrame(rows).to_string(index=False))
comparison n_a n_b median_a median_b diff ratio U p rank_biserial
Rats 1 and 3, FD4 < preflight 73 72 1.242 1.510 -0.268 0.82 2168.0 0.034603 0.175
Rats 1 and 3, FD4 < preflight, old version, all distinct maps 73 72 1.242 1.510 -0.268 0.82 2168.0 0.034603 0.175
Rat1, FD9 against preflight (two-sided) 26 27 1.749 1.869 -0.120 0.94 286.0 0.251140 0.185
Result: reproduced. Flight day 4 ES maps (104, Rats 1 to 3) have a median of 1.177 bits/spike against 1.441 for preflight (98); the one-sided Mann-Whitney p is 0.024, below the criterion's 0.05 (ratio 0.82, rank-biserial 0.162). The outcome differs from the first attempt's (p = 0.055) because of three maps: the draft leaves out every copy in the groups of identical maps, and three Rat2 FD4 ES maps the first attempt counted once each (7.583 bits/spike each, mean rate below 0.001 Hz) are no longer in the dandiset. On the earlier version's maps the same test gives p = 0.055 with all 107 distinct maps and p = 0.024 with the 104 that are in the draft, so with Rat2 included the result rests on whether those three maps are counted. Without Rat2, whose trace could not be attributed, Rats 1 and 3 give 1.242 against 1.510 (p = 0.035), the same on the earlier version, as the three maps are Rat2's. The FD4 median is lower in each rat (ratios 0.65 to 0.88), significantly only for Rat1 (p = 0.032; Rat2 p = 0.12, Rat3 p = 0.18). The test also passes with maps of peak rate 1 Hz or more (p = 0.0095), over bins with 5 or more occupancy counts (p = 0.023) and with spatial information adjusted for rate x occupancy (p = 0.0011), and not with maps whose archived names contain CEL (p = 0.073) or with MC maps (p = 0.060). Rat1's and Rat2's FD4 maps have less rate x occupancy, and so fewer spikes, than their preflight maps (medians 1891 and 1819 against 7862 and 7628), and their spatial information is lower even so: sampling bias works against the claim here. For Rats 1 and 2, FD4 is lower than FD9 without significance on the rate maps (1.081 against 1.385, p = 0.095) and on spike-based information minus the shuffle mean (0.025 against 0.211, p = 0.13). With Rat3's FD4 units added, that measure is higher on FD4 than on FD9 (0.232 against 0.211, p = 0.31). The spike-based fraction of significant units is not lower on FD4 (24 of 71 against 12 of 53 on FD9, p = 0.94). The claim's "abnormal" may refer to something other than lower spatial information, which the description does not say, and the preflight maps may come from another apparatus (see the data check).
fig3-stability-across-sessions¶
Place fields are stable across sessions. (Figure 3 title; which sessions the figure compares is not known)
Criterion: Rat3 FD4, its two ES epochs. Same tetrode and cluster id in both, at least 50 spikes at a known position in each. Smoothed 32 x 32 rate maps; median same-unit correlation against 1000 random re-pairings. Reproduced when above the 95th percentile of the re-paired medians and above 0.5. A unit is the same sorted cluster in both epochs, as the units' description states; that it is the same neuron is assumed (see the data check).
In [33]:
rng_c4 = np.random.default_rng(72910)
s3 = by_label["Rat3 FD4"]
es_epochs = s3["epochs"][s3["epochs"].type == "ES"].sort_values("start")
lookup3 = nl.PositionLookup(s3)
maps3, occs = {}, []
for k, (a, b, name) in enumerate(zip(es_epochs.start, es_epochs.stop, es_epochs.name)):
occ = nl.epoch_occupancy(s3, name)
occs.append(nl.merge_bins(nl.clean_occupancy(occ)))
print(f"ES epoch {k + 1}: {name} {a:.0f}-{b:.0f} s, occupancy map match {nl.occupancy_match(s3, a, b)}")
for _, u in s3["units"].iterrows():
sp = u.spikes[(u.spikes >= a) & (u.spikes <= b)]
found, r, c = lookup3(sp)
if found.sum() >= nl.MIN_SPIKES:
maps3[(u.tetrode, u.cluster, k)] = nl.smoothed_rate_map(r, c, occ)[0]
both = sorted({(t, c) for t, c, k in maps3 if (t, c, 0) in maps3 and (t, c, 1) in maps3})
only = sorted({(t, c) for t, c, k in maps3} - set(both))
print(f"Units with >= {nl.MIN_SPIKES} spikes in both epochs: {len(both)}; in one only: {len(only)}; units in file: {len(s3['units'])}")
shared = (occs[0] >= nl.MIN_OCC) & (occs[1] >= nl.MIN_OCC)
print("Bins with occupancy >= 5 in both epochs:", int(shared.sum()))
def corr(m1, m2):
ok = shared & np.isfinite(m1) & np.isfinite(m2)
return np.corrcoef(m1[ok], m2[ok])[0, 1]
A = [maps3[(t, c, 0)] for t, c in both]
B = [maps3[(t, c, 1)] for t, c in both]
C = np.array([[corr(a, b) for b in B] for a in A])
same = np.diag(C)
for (t, c), r in zip(both, same):
print(f" {t} cluster {c}: r = {r:.3f}")
null = np.array([np.median(C[np.arange(len(both)), rng_c4.permutation(len(both))]) for _ in range(1000)])
print(f"\nMedian same-unit r = {np.median(same):.3f} (range {same.min():.3f} to {same.max():.3f}); "
f"re-paired medians: mean {null.mean():.3f}, 95th percentile {np.percentile(null, 95):.3f}; "
f"p = {(np.sum(null >= np.median(same)) + 1) / 1001:.4f}")
off = C[~np.eye(len(both), dtype=bool)]
print(f"Different-unit pairs: median r {np.median(off):.3f}, range {off.min():.3f} to {off.max():.3f}")
print("Same-unit r above 0.5:", int((same > 0.5).sum()), "of", len(same), "; above 0.3:", int((same > 0.3).sum()))
print("\nVariants of the smoothing and the spike threshold (same re-pairing test, 1000 re-pairings each):")
rows = []
for sigma in [0.5, 1.0, 2.0]:
for min_spikes in [50, 100]:
mm = {}
for k, (a, b, name) in enumerate(zip(es_epochs.start, es_epochs.stop, es_epochs.name)):
occ = nl.epoch_occupancy(s3, name)
for _, u in s3["units"].iterrows():
sp = u.spikes[(u.spikes >= a) & (u.spikes <= b)]
found, r, c = lookup3(sp)
if found.sum() >= min_spikes:
mm[(u.tetrode, u.cluster, k)] = nl.smoothed_rate_map(r, c, occ, sigma)[0]
bb = sorted({(t, c) for t, c, k in mm if (t, c, 0) in mm and (t, c, 1) in mm})
CC = np.array([[corr(mm[(t1, c1, 0)], mm[(t2, c2, 1)]) for t2, c2 in bb] for t1, c1 in bb])
nn = np.array([np.median(CC[np.arange(len(bb)), rng_c4.permutation(len(bb))]) for _ in range(1000)])
rows.append({"sigma_bins": sigma, "min_spikes": min_spikes, "units": len(bb), "median_r": round(float(np.median(np.diag(CC))), 3),
"null_95th": round(float(np.percentile(nn, 95)), 3), "p": round((np.sum(nn >= np.median(np.diag(CC))) + 1) / 1001, 4)})
print(pd.DataFrame(rows).to_string(index=False))
ES epoch 1: ES1 32-1345 s, occupancy map match {'ES1': 0.8050601767055183, 'ES2': 0.6658463876929288, 'MC': 0.3174668304823492}
ES epoch 2: ES2 2397-3729 s, occupancy map match {'ES1': 0.6331878244119306, 'ES2': 0.8346318447648601, 'MC': 0.32440097449378746}
Units with >= 50 spikes in both epochs: 15; in one only: 5; units in file: 47
Bins with occupancy >= 5 in both epochs: 267
TT0 cluster 1: r = 0.277 TT10 cluster 2: r = 0.791 TT10 cluster 3: r = 0.284 TT3 cluster 1: r = 0.702 TT3 cluster 2: r = 0.849 TT3 cluster 3: r = 0.232 TT4 cluster 1: r = 0.820 TT4 cluster 3: r = 0.639 TT5 cluster 1: r = 0.838 TT5 cluster 2: r = -0.107 TT8 cluster 1: r = 0.010 TT8 cluster 3: r = 0.537 TT8 cluster 4: r = 0.076 TT8 cluster 5: r = 0.436 TT8 cluster 6: r = 0.762 Median same-unit r = 0.537 (range -0.107 to 0.849); re-paired medians: mean -0.003, 95th percentile 0.102; p = 0.0010 Different-unit pairs: median r -0.019, range -0.354 to 0.831 Same-unit r above 0.5: 8 of 15 ; above 0.3: 9 Variants of the smoothing and the spike threshold (same re-pairing test, 1000 re-pairings each):
sigma_bins min_spikes units median_r null_95th p
0.5 50 15 0.315 0.029 0.001
0.5 100 13 0.341 0.040 0.001
1.0 50 15 0.537 0.106 0.001
1.0 100 13 0.639 0.102 0.001
2.0 50 15 0.548 0.233 0.001
2.0 100 13 0.707 0.260 0.001
In [34]:
# Per-session results, and the variant on the dandiset's own rate maps: the ES1 and ES2 maps that
# the rate-map table links to the same unit.
per_epoch = tested[(tested.rat == "Rat3") & (tested.day == "FD4") & (tested.type == "ES")]
print(per_epoch.groupby("epoch").agg(tested=("sig", "size"), significant=("sig", "sum"), median_info=("info", "median")).round(3).to_string())
mp = s3["maps"]
rows = []
for unit, g in mp[mp.type == "ES"].groupby("unit"):
if {"ES1", "ES2"} <= set(g.epoch_name):
i, j = g.index[g.epoch_name == "ES1"][0], g.index[g.epoch_name == "ES2"][0]
r1, r2 = s3["rate_maps"][i], s3["rate_maps"][j]
ok = np.isfinite(r1) & np.isfinite(r2) & (s3["occupancy_maps"][i] > 0) & (s3["occupancy_maps"][j] > 0)
ok[nl.SENTINEL_BIN] = False
rows.append({"tetrode": mp.tetrode[i], "cluster": int(mp.cluster[i]), "in_main_test": (mp.tetrode[i], mp.cluster[i]) in both,
"r": np.corrcoef(r1[ok], r2[ok])[0, 1]})
rv = pd.DataFrame(rows)
print("\nDandiset rate maps, ES1 against ES2 map of the same unit (linked by the table's unit column):")
print(rv.round(3).to_string(index=False))
print(f"median r {rv.r.median():.3f} over {len(rv)} units with both maps; over the {int(rv.in_main_test.sum())} units of the main test "
f"{rv[rv.in_main_test].r.median():.3f}; units of the main test without both maps: {len(both) - int(rv.in_main_test.sum())}")
tested significant median_info
epoch
ES 2397-3729 19 1 1.629
ES 32-1345 16 11 2.155
Dandiset rate maps, ES1 against ES2 map of the same unit (linked by the table's unit column):
tetrode cluster in_main_test r
TT0 1 True 0.353
TT0 2 False -0.098
TT3 1 True 0.506
TT3 2 True 0.646
TT3 3 True 0.139
TT4 1 True 0.757
TT4 2 False 0.498
TT4 3 True 0.157
TT5 1 True 0.500
TT5 2 True -0.027
TT8 1 True -0.030
TT8 3 True 0.189
TT8 4 True 0.019
TT8 5 True 0.340
TT8 6 True 0.669
TT10 1 False NaN
TT10 2 True 0.378
TT10 3 True 0.024
median r 0.340 over 18 units with both maps; over the 15 units of the main test 0.340; units of the main test without both maps: 0
In [35]:
# Per unit in Rat3's two FD4 ES epochs: spikes, spike-based information and shuffle drops in each,
# and the same-unit correlation; then the stability of mean rates as a check on cluster identity.
ep_names = list(es_epochs.apply(lambda e: f"ES {e.start:.0f}-{e.stop:.0f}", axis=1))
u3 = units[(units.rat == "Rat3") & (units.day == "FD4") & (units.type == "ES")].copy()
u3["shuffle_spikes_dropped"] = 1 - u3.null_spikes_used / u3.spikes_in_epoch
wide = u3.pivot_table(index=["tetrode", "cluster"], columns="epoch", values=["spikes_in_epoch", "spikes_used", "info", "null95", "p"], aggfunc="first")
wide.columns = [f"{v}_{'e1' if e == ep_names[0] else 'e2'}" for v, e in wide.columns]
wide = wide[[f"{v}_{k}" for v in ["spikes_in_epoch", "spikes_used", "info", "null95", "p"] for k in ["e1", "e2"]]]
wide["r"] = pd.Series({(t, c): r for (t, c), r in zip(both, same)})
wide = wide.reset_index()
print(f"Units (e1 = {ep_names[0]}, e2 = {ep_names[1]}) with a spike in either epoch; "
f"{int(((wide.spikes_in_epoch_e1 == 0) & (wide.spikes_in_epoch_e2 == 0)).sum())} units of the file have none:")
print(wide[(wide.spikes_in_epoch_e1 > 0) | (wide.spikes_in_epoch_e2 > 0)].round(3).to_string(index=False))
print("\nPer epoch:")
rows = []
for k, (a, b) in enumerate(zip(es_epochs.start, es_epochs.stop)):
d = u3[u3.epoch == ep_names[k]]
t = d.dropna(subset=["info"])
rows.append({"epoch": ep_names[k], "duration_s": b - a, "occupancy_counts": int(occs[k].sum()),
"bins_occ>=5": int((occs[k] >= nl.MIN_OCC).sum()), "units_with_spikes": int((d.spikes_in_epoch > 0).sum()),
"tested": len(t), "significant": int((t.p < 0.05).sum()), "median_spikes_used": float(t.spikes_used.median()),
"median_info": round(float(t["info"].median()), 3), "median_null_mean": round(float(t.null_mean.median()), 3),
"median_null95": round(float(t.null95.median()), 3),
"median_shuffle_spikes_dropped": round(float(t.shuffle_spikes_dropped.median()), 3)})
print(pd.DataFrame(rows).to_string(index=False))
in_both = wide.dropna(subset=["spikes_in_epoch_e1", "spikes_in_epoch_e2"])
in_both = in_both[(in_both.spikes_in_epoch_e1 > 0) | (in_both.spikes_in_epoch_e2 > 0)]
d1 = np.diff(es_epochs[["start", "stop"]].values[0])[0]
d2 = np.diff(es_epochs[["start", "stop"]].values[1])[0]
r1, r2 = in_both.spikes_in_epoch_e1.values / d1, in_both.spikes_in_epoch_e2.values / d2
rho = stats.spearmanr(r1, r2).statistic
null_rho = np.array([stats.spearmanr(r1, rng_c4.permutation(r2)).statistic for _ in range(1000)])
print(f"\nUnits with spikes in either epoch: {len(in_both)}; with spikes in both: {int(((r1 > 0) & (r2 > 0)).sum())}; "
f"in the first only: {int(((r1 > 0) & (r2 == 0)).sum())}; in the second only: {int(((r1 == 0) & (r2 > 0)).sum())}")
print(f"Mean rate, epoch 1 against epoch 2, Spearman rho = {rho:.3f}; re-paired 95th percentile {np.percentile(null_rho, 95):.3f}; "
f"p = {(np.sum(null_rho >= rho) + 1) / 1001:.4f}")
Units (e1 = ES 32-1345, e2 = ES 2397-3729) with a spike in either epoch; 23 units of the file have none:
tetrode cluster spikes_in_epoch_e1 spikes_in_epoch_e2 spikes_used_e1 spikes_used_e2 info_e1 info_e2 null95_e1 null95_e2 p_e1 p_e2 r
TT0 1 95 233 92 208 3.957 3.255 4.482 6.624 0.150 0.192 0.277
TT0 2 22 491 20 469 NaN 2.099 NaN 5.120 NaN 0.518 NaN
TT10 1 46 190 46 188 NaN 5.385 NaN 5.566 NaN 0.064 NaN
TT10 2 331 749 320 694 3.124 1.648 2.526 1.676 0.004 0.082 0.791
TT10 3 978 1045 939 983 1.053 1.621 1.034 1.595 0.044 0.040 0.284
TT10 4 0 50 0 49 NaN NaN NaN NaN NaN NaN NaN
TT10 5 0 46 0 45 NaN NaN NaN NaN NaN NaN NaN
TT10 6 0 7 0 7 NaN NaN NaN NaN NaN NaN NaN
TT11 1 0 986 0 927 NaN 1.478 NaN 3.260 NaN 0.406 NaN
TT3 1 294 474 282 455 1.687 1.629 1.599 1.980 0.032 0.580 0.702
TT3 2 120 163 118 157 2.947 2.127 2.573 2.777 0.000 0.634 0.849
TT3 3 85 175 82 164 1.860 1.395 2.307 2.180 0.690 0.918 0.232
TT4 1 1603 449 1582 423 2.449 3.592 1.326 4.305 0.000 0.160 0.820
TT4 2 212 30 201 29 2.597 NaN 2.387 NaN 0.010 NaN NaN
TT4 3 638 449 594 432 1.839 1.602 1.639 2.014 0.006 0.344 0.639
TT5 1 9102 2322 8844 2221 0.910 1.472 0.504 2.393 0.000 0.338 0.838
TT5 2 422 816 376 802 1.519 1.255 1.976 1.490 0.132 0.238 -0.107
TT5 3 0 368 0 349 NaN 1.415 NaN 1.681 NaN 0.378 NaN
TT8 1 479 461 459 432 0.528 0.955 0.804 1.396 1.000 0.820 0.010
TT8 2 3 0 3 0 NaN NaN NaN NaN NaN NaN NaN
TT8 3 124 500 120 461 3.222 2.039 3.064 2.162 0.038 0.126 0.537
TT8 4 238 624 200 588 3.197 1.833 4.550 2.608 0.184 0.304 0.076
TT8 5 1851 847 1766 817 1.101 1.293 0.785 2.183 0.000 0.406 0.436
TT8 6 176 361 140 300 4.099 2.283 2.900 2.343 0.000 0.080 0.762
Per epoch:
epoch duration_s occupancy_counts bins_occ>=5 units_with_spikes tested significant median_spikes_used median_info median_null_mean median_null95 median_shuffle_spikes_dropped
ES 32-1345 1313.0 71496 348 19 16 11 301.0 2.155 1.667 2.141 0.205
ES 2397-3729 1332.0 62327 299 23 19 1 455.0 1.629 1.646 2.183 0.215
Units with spikes in either epoch: 24; with spikes in both: 18; in the first only: 1; in the second only: 5 Mean rate, epoch 1 against epoch 2, Spearman rho = 0.598; re-paired 95th percentile 0.339; p = 0.0010
Result: reproduced. In Rat3's two flight day 4 ES epochs, the 15 units with 50 or more spikes in each have a median rate-map correlation of 0.537 (range -0.107 to 0.849), above the 95th percentile of re-paired medians (0.102, p = 0.001) and above 0.5; 8 of 15 exceed 0.5. These are the same as in 0.260728.1352. With other smoothing the median is 0.315 to 0.707, always above its re-paired 95th percentile (p = 0.001), but at 0.5 bins of smoothing it is below 0.5 (0.315 and 0.341), which would make the claim partially reproduced. The dandiset's own ES1 and ES2 maps, now linked to the same unit, give a median of 0.340 over 18 units and over the 15 units of the main test. A unit is the same sorted cluster in both epochs, as the units' description now states (each tetrode was cut once); whether a cluster is one neuron throughout is not shown, and the files hold no waveforms. The units' mean rates are related across the two epochs (Spearman rho 0.598 over the 24 units with a spike in either, re-paired 95th percentile 0.339). Of the file's 47 units, 23 have no spike in either ES epoch, 18 have spikes in both, 1 only in the first and 5 only in the second. In the first epoch 11 of 16 units are significant against shuffles and in the second 1 of 19. The second epoch has similar occupancy (62327 against 71496 counts), more spikes per tested unit (median 455 against 301) and similar shuffle drops (0.215 against 0.205), but its units' information is at the level of their shuffles (median 1.629 against a shuffle mean of 1.646; in the first epoch 2.155 against 1.667). Several units whose fields correlate across the epochs are significant in the first epoch only, and for them the 95th percentile of the shuffles rises more in the second epoch than their information does: TT4 cluster 1 (r = 0.820) has 2.449 against 1.326 in the first epoch and 3.592 against 4.305 in the second, TT5 cluster 1 (r = 0.838) 0.910 against 0.504 and 1.472 against 2.393, TT8 cluster 5 (r = 0.436) 1.101 against 0.785 and 1.293 against 2.183. The epochs use 348 and 299 bins with 5 or more occupancy counts. With Poisson trains in place of the spikes the test gives 0 to 2 false positives of 16 to 20 units in either epoch (flight-spatial-code), so it is not too lenient there. The 95th percentile of a unit's shuffles follows how few 10 s periods hold its spikes and how few spikes it has (flight-spatial-code): TT4 cluster 1 has 1582 spikes filling an effective 0.267 of the 10 s bins in the first epoch and 423 filling 0.096 in the second, and TT5 cluster 1 8844 spikes at 0.575 and 2221 at 0.307. In the second epoch these units fire fewer spikes in fewer periods, so the shuffles reach higher information, and the test loses power for them; the stability measure does not depend on the shuffle test. Which sessions Figure 3 compares is not known.
theta-normal and ripples-normal¶
Hippocampal EEG showed normal theta rhythm during locomotion and normal sharp waves and ripples during quiet rest in flight. (Dandiset description)
Blocked. No file holds EEG or any continuous voltage: the unfiltered listing in the data check
shows an empty acquisition/ group and no ElectricalSeries or LFP in any file, and the draft's
experiment description says no EEG or other continuous voltage was archived. The NASA record of the
study (OSD-968) describes its files as spike files, rate-map files, legacy metadata files and NWB
copies, and its metadata does not mention EEG, theta, ripples or LFP (printed under
three-orthogonal-surfaces). Its archives were not opened, so an EEG file inside one cannot be ruled
out.
In [36]:
print(f"Run time of the notebook: {(time.time() - started) / 60:.1f} minutes")
Run time of the notebook: 6.6 minutes
Summary¶
| Claim | Outcome |
|---|---|
flight-spatial-code |
Not reproduced: Rat1 FD9 0.40 of units significant (0.40 to 0.44 over the repeats), Rat2 FD9 0.071 |
three-orthogonal-surfaces |
Blocked: metadata insufficient (no mapping of pixels to the track's surfaces) |
si-comparable-to-preflight |
Reproduced: FD9 1.385 against preflight 1.448 bits/spike, p = 0.45 |
fd4-abnormal-selectivity |
Reproduced: FD4 1.177 against preflight 1.441 bits/spike, one-sided p = 0.024 (0.055 with the three maps the draft leaves out) |
fig3-stability-across-sessions |
Reproduced: median r 0.537 between Rat3's two FD4 ES epochs (null 95th percentile 0.102) |
theta-normal |
Blocked: data missing (no EEG) |
ripples-normal |
Blocked: data missing (no EEG) |
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.
Against the first attempt on 0.260728.1352, fd4-abnormal-selectivity changes from not reproduced
to reproduced, because the draft leaves out three maps that moved the test across 0.05; the other
outcomes are the same.
How the analysis differs from the paper's. The paper's full text was not available, so its Methods, tests and thresholds are unknown; every method here is a common choice fixed before the analysis. Four claims come from the dandiset's description of the paper's findings, and one from a figure title alone. Spatial information is computed from the spikes for the spike-based claims, with the occupancy maps of the original analysis, as the dandiset holds position only at spike times; the comparisons with preflight use the original rate maps, as the preflight files have no position. The rate maps are smoothed with an unrecorded kernel and their spatial information is not corrected for the number of spikes; a count-adjusted variant is shown. The stability claim takes a unit to be the same neuron in two epochs because it is the same sorted cluster. With Poisson trains moved to position sample times, which reproduce the real drop rates, the shuffle test gives 0.096 false positives over all flight ES epochs and up to 0.214 in one session (0.074 and 0.160 with plain Poisson trains), against the 0.05 the criterion assumes; judged against these fractions, the outcome of flight-spatial-code is the same. Real trains reversed in time, which keep their bursts and changes in rate, give 2 of 122: with real temporal structure the test is conservative, and it has less power for units whose spikes fall in a few periods, such as most of Rat2's.
What was found about the dandiset (draft).
- Each flight file now holds the shared tracker's trace only in its own rat's task epochs, and says the trace there is taken to be that rat's. Rat2's units are not selective with it in any window, and Rat1's are selective in Rat2's FD4 ES window, so whose position the trace is in Rat2's epochs is still not established.
- The dandiset description still says the preflight recordings were on a flat rectangular track, while the preflight epochs and rate maps are labelled ES and MC, and still reports EEG, which no file holds.
- Rat2's preflight file shares no tetrode with its flight files, and the assignment of the preflight recordings to rats is an open question in the conversion notes.
- Position is given only at spike times, so occupancy cannot be computed from the files' position.
- The maps are smoothed with a kernel that is not recorded, and the rats' age, the camera view and the pixel size are not recorded.
- Fixed since the first attempt: the rate maps are linked to their units; NaN replaces -999 and the lost position, and the [63, 0] occupancy is documented; identical maps are left out (all copies, so six maps are gone); units list their spike files and windows and are documented as one sorted cluster across epochs; the spike-time resolution is documented and explains the zero intervals; the epochs no longer overlap.