Reproducing Carcea et al. (2021) from dandiset 000114¶
Paper. Carcea I, Caraballo NL, Marlin BJ, et al. (2021). Oxytocin neurons enable social transmission of maternal behaviour. Nature 596, 553-557. doi:10.1038/s41586-021-03814-7. This attempt used the full text, which is open access (CC BY, PMC8387235), with its Source Data spreadsheets for the reported numbers.
Dandiset. DANDI:000114, version 0.230602.1643. The paper cites the dandiset without a version (it was published in 2021, before this version). This is the only published version. The draft holds the same 30 files; each is 240 bytes or 6 KB larger, and the check below compares the spike times, the events and the shape and rate of the voltage series (the photometry trace itself) of every draft file with the pinned version.
What the dandiset holds. PVN recordings from five virgins co-housed with a dam and litter (V85, V110, V111, V116, V118): raw voltage from a 16-wire bundle, sorted units for four of them in 8 of 14 files, and video annotations. Fibre photometry from six non-co-housed virgins during observation of dam retrievals (ROV40 to ROV49). It does not hold the optical tagging that identifies oxytocin (OT-PVN) units, the behaviour of Figure 1, the unit recordings during observation, or the auditory cortex data of Figure 4.
Outside input. The files do not name their behaviour labels. The paper's code (GitHub
ic283/Carcea-et-al.-2021, commit 2db2a75) uses "Mom agressing" for shepherding, "Virgin in rest"
for nest entry and "Mom retrieval" for dam retrieval; claims that need these names are tested
through the claim label-link. The analysis follows the Methods and, where they leave a choice
open, the code's calc_spike_freq_behavior_significance.m; the choices are listed in the claims
file and in pvn.py.
| Claim | Tested here |
|---|---|
| label-link: the labels can be named through the code | tested |
| ed7a-nest-entry-fraction: units activated during nest entry | tested |
| ed7a-shepherding-fraction: units activated during shepherding | tested |
| ed7a-dam-retrieval-fraction: units activated during dam retrieval | tested |
| fig2j-shepherded-entry: higher firing for shepherded than voluntary nest entries | tested |
| fig1j-shepherding-rate: shepherding above 0.2 per hour | tested |
| ed8a-dam-retrieval-rate: dams retrieve about once every two hours | tested |
| ot-pvn-units: responses and correlations of OT-PVN units | blocked, data missing (no tagging) |
| fig2f-virgin-chasing: less firing when the virgin chases the dam | blocked, metadata insufficient |
| ed6-dam-vs-pup: dam against pup interactions, pup calls | blocked, data missing |
| fig1cd-retrieval-onset, fig2a-dreadd, fig3fg-observational-learning, fig4jk-oxytocin-in-ac, ed9-colliculus: behavioural cohorts | blocked, data missing |
| fig3ij-observation-units: PVN units during observation | blocked, data missing |
| fig4c-ac-photometry, fig4f-pvn-ac-correlation: auditory cortex | blocked, data missing |
| fig4i-projection-photometry: learners against non-learners | blocked, metadata insufficient |
import time
import warnings
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from scipy import stats
import pvn
pd.set_option("display.width", 220)
pd.set_option("display.max_columns", 40)
pd.set_option("display.max_rows", 200)
warnings.filterwarnings("ignore", category=RuntimeWarning)
started = time.time()
ASSETS = pvn.assets()
eph = [pvn.load_ecephys(p, u) for p, size, u in ASSETS if p.endswith("_ecephys.nwb")]
oph = [pvn.load_ophys(p, u) for p, size, u in ASSETS if p.endswith("_ophys.nwb")]
with_units = [s for s in eph if s["units"]]
print(f"dandiset {pvn.DANDISET_ID} version {pvn.DANDISET_VERSION}: {len(ASSETS)} files, "
f"{len(eph)} ecephys and {len(oph)} ophys, loaded in {time.time() - started:.0f} s")
dandiset 000114 version 0.230602.1643: 30 files, 14 ecephys and 16 ophys, loaded in 73 s
Data check¶
Each check of the metadata against the values is stated under its output.
rows = []
for s in eph:
ser = s["series"]
rows.append({"session": s["session"], "subject": s["subject"]["subject_id"], "session_id": s["session_id"],
"start_time": s["session_start_time"][:19], "series": len(ser),
"rate_hz": sorted({v["rate"] for v in ser.values()}),
"rec_h": round(max(v["duration_s"] for v in ser.values()) / 3600, 2),
"units": len(s["units"]), "events": len(s["events"]),
"last_event_h": round(s["events"]["stop"].max() / 3600, 2) if len(s["events"]) else None,
"last_spike_h": round(max(u.max() for u in s["units"]) / 3600, 2) if s["units"] else None})
files = pd.DataFrame(rows)
print(files.to_string(index=False))
print("\nSubjects (ecephys):", sorted(files.subject.unique()), "; (ophys):", sorted({s['subject']['subject_id'] for s in oph}))
print("Distinct subjects in all files:", len(set(files.subject) | {s['subject']['subject_id'] for s in oph}),
"(archive summary: 11)")
print("\nSubject fields (ecephys, all files):",
{k: sorted({s['subject'][k] for s in eph}) for k in ['species', 'sex', 'age', 'genotype', 'strain', 'description']})
print("Subject fields (ophys, all files):",
{k: sorted({s['subject'][k] for s in oph}) for k in ['genotype', 'description']})
print("Electrode locations:", sorted({l for s in eph for l in s['electrode_locations']}),
"; channels per series:", sorted({len(s['channels']) for s in eph}))
print("Units table columns:", sorted({tuple(s['units_columns']) for s in with_units}))
print("Stated voltage unit and conversion:", sorted({(v['unit'], v['conversion']) for s in eph for v in s['series'].values()}))
print("unit_name values per session are 0 to n-1:",
{s["session"]: s["unit_names"] == [str(i) for i in range(len(s["units"]))] for s in with_units})
session subject session_id start_time series rate_hz rec_h units events last_event_h last_spike_h
sub-V110_ses-03032018 V110 03032018 2018-03-04T14:00:40 1 [30000.0] 12.04 31 123 12.04 12.04
sub-V111_ses-03032018 V111 03032018 2018-03-04T09:15:48 1 [30000.0] 4.62 31 105 4.62 4.62
sub-V111_ses-03052018 V111 03052018 2018-03-06T13:55:38 1 [30000.0] 10.77 31 236 10.77 10.77
sub-V111_ses-03062018 V111 03062018 2018-03-07T03:47:12 1 [30000.0] 12.41 31 279 12.41 12.41
sub-V116_ses-03082018 V116 03082018 2018-05-07T13:47:05 1 [30000.0] 12.22 23 92 12.21 12.22
sub-V116_ses-03092018 V116 03092018 2018-03-11T23:50:56 1 [30000.0] 15.82 21 164 15.80 15.82
sub-V116_ses-03102018 V116 03102018 2018-05-08T12:51:47 1 [30000.0] 4.46 26 86 13.45 13.46
sub-V116_ses-03122018 V116 03122018 2018-03-12T04:10:23 1 [30000.0] 12.53 0 53 12.52 NaN
sub-V118_ses-03182018 V118 03182018 2018-03-23T18:42:22 1 [30000.0] 12.95 34 112 12.96 12.95
sub-V118_ses-03192018 V118 03192018 2018-03-20T04:03:13 1 [30000.0] 10.50 0 69 9.89 NaN
sub-V118_ses-03212018 V118 03212018 2018-03-21T20:30:06 4 [30000.0] 16.96 0 0 NaN NaN
sub-V118_ses-03222018 V118 03222018 2018-03-23T01:13:10 4 [30000.0] 2.27 0 16 2.52 NaN
sub-V85_ses-12102017-2 V85 12102017_2 2017-12-12T07:47:55 1 [30000.0] 1.88 0 0 NaN NaN
sub-V85_ses-12102017 V85 12102017 2017-12-11T19:41:23 1 [30000.0] 12.01 0 497 12.01 NaN
Subjects (ecephys): ['V110', 'V111', 'V116', 'V118', 'V85'] ; (ophys): ['ROV40', 'ROV41', 'ROV43', 'ROV44', 'ROV45', 'ROV49']
Distinct subjects in all files: 11 (archive summary: 11)
Subject fields (ecephys, all files): {'species': ['Mus musculus'], 'sex': ['F'], 'age': ['P8W/P10W'], 'genotype': ['Oxy-IRES-Cre het'], 'strain': ['mixed C57BL/6 x 129S'], 'description': ['Virgin female co-housed with dam and pups']}
Subject fields (ophys, all files): {'genotype': ['Oxy-IRES-Cre het'], 'description': ['Virgin female co-housed with dam and pups']}
Electrode locations: ['unknown'] ; channels per series: [16]
Units table columns: [('spike_times', 'unit_name')]
Stated voltage unit and conversion: [('volts', 2.5e-07)]
unit_name values per session are 0 to n-1: {'sub-V110_ses-03032018': True, 'sub-V111_ses-03032018': True, 'sub-V111_ses-03052018': True, 'sub-V111_ses-03062018': True, 'sub-V116_ses-03082018': True, 'sub-V116_ses-03092018': True, 'sub-V116_ses-03102018': True, 'sub-V118_ses-03182018': True}
The files hold 11 subjects, as the archive's summary says. The paper reports 541 PVN units from 17 virgins; the dandiset has units from 4 virgins. Every electrode location is "unknown", and every ecephys and ophys subject is described as "Virgin female co-housed with dam and pups" with genotype Oxy-IRES-Cre het, although the photometry animals were not co-housed (they watched dam retrievals in a separate arena). The units table holds spike times and a unit name only: no waveform, channel, quality or OT-PVN label.
The session_start_time does not match the date in the session id in several files: V116
03082018 and 03102018 start in May 2018, V116 03092018 starts on 11 March, and V118 03182018
starts on 23 March, after V118 03192018 and on the same day as V118 03222018. For V116 03102018
the raw voltage lasts 4.46 h but spikes and events run to 13.4 h, so its spike times and
annotations do not share the extent of the raw series; V118 03222018 likewise has events up to
2.52 h against 2.27 h of raw voltage. The other units files end their spikes within a minute of
the raw series.
q = []
for s in with_units:
sp = s["units"]
rate = np.array([len(u) / (u.max() - u.min()) for u in sp])
q.append({"session": s["session"], "units": len(sp), "spikes": sum(len(u) for u in sp),
"sorted": all(np.all(np.diff(u) >= 0) for u in sp),
"rate_median_hz": round(np.median(rate), 2), "rate_min": round(rate.min(), 3), "rate_max": round(rate.max(), 1),
"isi_lt2ms_median": round(np.median([pvn.isi_violation(u) for u in sp]), 4),
"isi_lt2ms_max": round(max(pvn.isi_violation(u) for u in sp), 4),
"coincident_0.1ms": round(pvn.coincident_fraction(sp), 3),
"units_over_excel_rows": int(sum(len(u) > pvn.EXCEL_ROWS for u in sp))})
quality = pd.DataFrame(q)
print(quality.to_string(index=False))
# Recordings deposited twice: identical spike trains or annotations across files.
sig = {}
for s in with_units:
for i, u in enumerate(s["units"]):
sig.setdefault((len(u), round(float(u[0]), 4), round(float(u[-1]), 4)), []).append((s["session"], i))
print("\nUnits with identical spike count, first and last spike in two files:", [v for v in sig.values() if len(v) > 1])
ev_sig = {}
for s in eph:
if len(s["events"]):
ev_sig.setdefault(tuple(np.round(s["events"].sort_values("start").start.to_numpy()[:20], 2)), []).append(s["session"])
print("Files with identical first 20 annotation starts:", [v for v in ev_sig.values() if len(v) > 1])
# Units that look the same across days: correlation of log spike counts between sessions with equal unit counts.
names = [s["session"] for s in with_units]
for i in range(len(with_units)):
for j in range(i + 1, len(with_units)):
a, b = with_units[i]["units"], with_units[j]["units"]
if len(a) == len(b) and with_units[i]["subject"]["subject_id"] == with_units[j]["subject"]["subject_id"]:
r = np.corrcoef(np.log([len(u) for u in a]), np.log([len(u) for u in b]))[0, 1]
print(f"{names[i]} and {names[j]}: {len(a)} units each, correlation of log spike counts unit by unit {r:.3f}")
session units spikes sorted rate_median_hz rate_min rate_max isi_lt2ms_median isi_lt2ms_max coincident_0.1ms units_over_excel_rows sub-V110_ses-03032018 31 30618470 True 24.57 3.576 87.5 0.0922 0.1187 0.788 17 sub-V111_ses-03032018 31 6111694 True 15.42 0.227 28.8 0.0725 0.0986 0.864 0 sub-V111_ses-03052018 31 14299265 True 16.16 0.027 30.9 0.0784 0.1003 0.901 2 sub-V111_ses-03062018 31 10454175 True 8.34 0.010 41.3 0.0698 0.0865 0.786 1 sub-V116_ses-03082018 23 405978 True 0.23 0.006 3.2 0.0165 0.0342 0.280 0 sub-V116_ses-03092018 21 993626 True 0.37 0.083 5.3 0.0413 0.0594 0.152 0 sub-V116_ses-03102018 26 2262635 True 0.27 0.008 15.7 0.0155 0.0582 0.272 0 sub-V118_ses-03182018 34 19750322 True 9.31 0.073 51.1 0.0153 0.0595 0.446 6 Units with identical spike count, first and last spike in two files: [] Files with identical first 20 annotation starts: [] sub-V111_ses-03032018 and sub-V111_ses-03052018: 31 units each, correlation of log spike counts unit by unit 0.974 sub-V111_ses-03032018 and sub-V111_ses-03062018: 31 units each, correlation of log spike counts unit by unit 0.934 sub-V111_ses-03052018 and sub-V111_ses-03062018: 31 units each, correlation of log spike counts unit by unit 0.963
Spike times are sorted, and no spike train or annotation table appears in two files, so no recording is deposited twice. Unit quality differs between mice. In V110 and V111 a median 7.0 to 9.2% of inter-spike intervals are under 2 ms, and 79 to 90% of all spikes lie within 0.1 ms of a spike of another unit; the Methods say spikes violating a 2 to 3 ms refractory period, and spikes recorded at once on 10 or more channels, were removed. These units are consistent with multi-unit clusters or one neuron split over several wires; the claims below are also given without V110 and V111. In V116 the values are 1.6 to 4.1% and 15 to 28%, and in V118 1.5% and 45%.
V111's three sessions each hold 31 units whose spike counts are correlated unit by unit (0.93 to 0.97), so they are probably the same clusters on consecutive days. The paper counts units per day (541 units from 17 virgins); here too each session's units are counted separately, and results are also given per session and per mouse.
lab_rows = []
for s in eph:
ev = s["events"]
for lab, g in ev.groupby("label"):
d = g.stop - g.start
lab_rows.append({"session": s["session"], "label": lab, "n": len(g), "dur_median_s": round(d.median(), 1),
"dur_min": round(d.min(), 1), "dur_max": round(d.max(), 1)})
labels = pd.DataFrame(lab_rows)
print(labels.pivot_table(index="label", columns="session", values="n", aggfunc="sum").fillna(0).astype(int).to_string())
print("\nNo stop before its start:", bool(all((s["events"].stop >= s["events"].start).all() for s in eph)))
print("Annotation table description:", "'events and states of pups'", "(the same in every file; it describes no label)")
session sub-V110_ses-03032018 sub-V111_ses-03032018 sub-V111_ses-03052018 sub-V111_ses-03062018 sub-V116_ses-03082018 sub-V116_ses-03092018 sub-V116_ses-03102018 sub-V116_ses-03122018 sub-V118_ses-03182018 sub-V118_ses-03192018 sub-V118_ses-03222018 sub-V85_ses-12102017 label Mom Building 2 5 4 2 0 2 1 1 0 0 0 0 Mom aggressing 0 0 0 0 0 0 0 0 0 0 0 19 Mom agressing 4 4 2 1 6 22 13 7 0 0 0 0 Mom retrieval 0 1 0 2 4 11 2 2 0 0 1 0 Mother in nest 41 73 86 69 58 94 51 0 112 69 11 338 Pup out 0 0 0 1 2 5 2 2 0 0 1 0 Self gen retrieval 0 0 0 0 0 0 0 2 0 0 0 0 Shepherding 0 0 0 0 0 0 0 7 0 0 0 0 Virgin Agressing 4 1 11 15 0 1 0 0 0 0 0 0 Virgin in nest 0 0 0 0 0 0 0 19 0 0 0 140 Virgin in rest 72 21 133 189 22 29 17 13 0 0 3 0 No stop before its start: True Annotation table description: 'events and states of pups' (the same in every file; it describes no label)
The labels have no description. "Mom agressing" is spelled "Mom aggressing" in V85. V118's first two files carry only "Mother in nest"; V118 03212018 and the short V85 file have no annotations. The annotations describe the virgin and the dam, not the pups, as the table's description says.
# Annotation coverage: first and last annotation of any label, and the longest stretch without one.
cov = []
for s in eph:
e = s["events"].sort_values("start")
if not len(e):
continue
gaps = e.start.to_numpy()[1:] - np.maximum.accumulate(e.stop.to_numpy())[:-1]
cov.append({"session": s["session"], "rec_h": round(max(v["duration_s"] for v in s["series"].values()) / 3600, 2),
"annotations": len(e), "first_h": round(e.start.min() / 3600, 2), "last_h": round(e.stop.max() / 3600, 2),
"longest_gap_h": round(gaps.max() / 3600, 2) if len(gaps) else None})
print(pd.DataFrame(cov).to_string(index=False))
# Coverage of the labels the rate claims count: virgin labels, "Mom agressing" and "Mom retrieval".
lab_cov = []
for s in eph:
row = {"session": s["session"]}
for name, labs in (("virgin", pvn.VIRGIN_LABELS), ("shepherding", pvn.SHEPHERDING), ("dam_retrieval", pvn.DAM_RETRIEVAL)):
e = s["events"][s["events"].label.isin(labs)].sort_values("start")
row[f"{name}_n"] = len(e)
if len(e):
g = e.start.to_numpy()[1:] - np.maximum.accumulate(e.stop.to_numpy())[:-1]
row[f"{name}_first_h"] = round(e.start.min() / 3600, 2)
row[f"{name}_last_h"] = round(e.stop.max() / 3600, 2)
row[f"{name}_gap_h"] = round(g.max() / 3600, 2) if len(g) else None
if row["virgin_n"]:
row["rec_h"] = round(pvn.session_length_s(s) / 3600, 2)
lab_cov.append(row)
print("\nCoverage by label group, sessions with a virgin label (first and last annotation, longest gap between two, in hours):")
print(pd.DataFrame(lab_cov).to_string(index=False))
session rec_h annotations first_h last_h longest_gap_h
sub-V110_ses-03032018 12.04 123 0.00 12.04 0.27
sub-V111_ses-03032018 4.62 105 0.00 4.62 0.08
sub-V111_ses-03052018 10.77 236 0.03 10.77 0.12
sub-V111_ses-03062018 12.41 279 0.00 12.41 0.32
sub-V116_ses-03082018 12.22 92 0.00 12.21 0.09
sub-V116_ses-03092018 15.82 164 0.00 15.80 0.05
sub-V116_ses-03102018 4.46 86 0.00 13.45 0.03
sub-V116_ses-03122018 12.53 53 0.00 12.52 0.10
sub-V118_ses-03182018 12.95 112 0.00 12.96 0.41
sub-V118_ses-03192018 10.50 69 0.00 9.89 0.81
sub-V118_ses-03222018 2.27 16 0.00 2.52 0.01
sub-V85_ses-12102017 12.01 497 0.01 12.01 0.08
Coverage by label group, sessions with a virgin label (first and last annotation, longest gap between two, in hours):
session virgin_n virgin_first_h virgin_last_h virgin_gap_h shepherding_n shepherding_first_h shepherding_last_h shepherding_gap_h dam_retrieval_n rec_h dam_retrieval_first_h dam_retrieval_last_h dam_retrieval_gap_h
sub-V110_ses-03032018 76 0.00 12.04 1.93 4 3.88 11.65 3.84 0 12.04 NaN NaN NaN
sub-V111_ses-03032018 22 0.61 4.62 0.14 4 0.48 1.17 0.58 1 4.62 0.15 0.15 NaN
sub-V111_ses-03052018 144 0.03 10.69 0.34 2 3.95 10.52 6.56 0 10.77 NaN NaN NaN
sub-V111_ses-03062018 204 0.15 12.35 0.52 1 12.30 12.31 NaN 2 12.41 4.90 4.91 0.00
sub-V116_ses-03082018 22 0.00 12.21 0.18 6 4.54 9.30 4.51 4 12.22 1.48 11.91 9.89
sub-V116_ses-03092018 30 0.00 15.80 0.10 22 2.00 14.49 2.37 11 15.82 0.38 14.53 9.03
sub-V116_ses-03102018 17 0.00 13.45 0.21 13 1.56 13.10 4.18 2 13.45 2.62 8.80 6.19
sub-V116_ses-03122018 32 0.00 12.52 0.10 7 8.02 10.15 2.07 2 12.53 9.76 9.99 0.23
sub-V118_ses-03222018 3 0.00 2.52 0.07 0 NaN NaN NaN 1 2.52 1.91 1.91 NaN
sub-V85_ses-12102017 140 0.02 11.86 1.81 19 2.23 11.86 5.55 0 12.01 NaN NaN NaN
Annotations of any label start within 0.03 h of the start of every file. They end at or after the end of the recording except in V118 03192018 (9.89 against 10.50 h, 0.61 h short) and in three V116 files (0.01 to 0.02 h short), and no stretch without an annotation lasts more than 0.81 h (V118 03192018). Most of these annotations are "Mother in nest".
The labels the rate claims count are sparser. The virgin labels run from within 0.61 h of each session's start (V111 03032018; within 0.15 h in the others) to within 0.15 h of its end, with gaps of up to 1.93 h (V110 03032018) and 1.81 h (V85). Shepherding and dam retrieval are rare events, with gaps between them of up to 6.56 h and 9.89 h, so their spacing cannot tell unscored stretches from stretches without the behaviour. This is consistent with each session being scored over its whole length, which the rate claims assume; it does not show that every event in a session was annotated.
# Shared clock: firing change at annotation onsets (5 s after minus 5 s before), with the onsets
# shifted in time. If units respond to the onsets and share the annotations' clock, the contrast
# is largest near shift 0.
shifts = np.arange(-120, 121, 2.0)
clock = []
for s in with_units:
for lab in ("Virgin in rest", "Mom agressing", "Mother in nest"):
st, _ = pvn.episodes(s["events"], (lab,))
if len(st) < 4:
continue
c = pvn.onset_contrast(s["units"], st, shifts)
c0 = c[shifts == 0][0]
clock.append({"session": s["session"], "label": lab, "onsets": len(st), "contrast_at_0": round(c0, 2),
"median_over_shifts": round(np.median(c), 2), "best_shift_s": shifts[np.argmax(c)],
"rank_of_0": f"{int(np.sum(c > c0)) + 1}/{len(c)}"})
print(pd.DataFrame(clock).to_string(index=False))
session label onsets contrast_at_0 median_over_shifts best_shift_s rank_of_0 sub-V110_ses-03032018 Virgin in rest 72 0.57 0.74 88.0 76/121 sub-V110_ses-03032018 Mom agressing 4 0.54 0.99 82.0 108/121 sub-V110_ses-03032018 Mother in nest 41 0.85 0.81 -60.0 56/121 sub-V111_ses-03032018 Virgin in rest 21 0.72 0.71 -62.0 60/121 sub-V111_ses-03032018 Mom agressing 4 0.75 0.99 -6.0 96/121 sub-V111_ses-03032018 Mother in nest 73 1.26 0.78 4.0 17/121 sub-V111_ses-03052018 Virgin in rest 133 1.08 0.71 -94.0 25/121 sub-V111_ses-03052018 Mother in nest 86 0.41 0.75 68.0 111/121 sub-V111_ses-03062018 Virgin in rest 189 0.58 0.67 -84.0 73/121 sub-V111_ses-03062018 Mother in nest 69 0.37 0.77 24.0 117/121 sub-V116_ses-03082018 Virgin in rest 22 0.65 0.85 -118.0 108/121 sub-V116_ses-03082018 Mom agressing 6 0.89 0.81 -74.0 36/121 sub-V116_ses-03082018 Mother in nest 58 0.57 0.87 -14.0 109/121 sub-V116_ses-03092018 Virgin in rest 29 0.53 0.85 8.0 113/121 sub-V116_ses-03092018 Mom agressing 22 0.49 0.89 -6.0 117/121 sub-V116_ses-03092018 Mother in nest 94 0.73 0.73 30.0 60/121 sub-V116_ses-03102018 Virgin in rest 17 0.84 0.84 74.0 60/121 sub-V116_ses-03102018 Mom agressing 13 0.56 0.86 -32.0 115/121 sub-V116_ses-03102018 Mother in nest 51 0.48 0.79 70.0 119/121 sub-V118_ses-03182018 Mother in nest 112 0.72 0.78 -58.0 80/121
The test does not establish a shared clock between spikes and annotations in any session: the contrast at shift 0 is not above the contrast at most other shifts, and its best rank is 17 of 121 ("Mother in nest" in V111 03032018). This can mean that few units respond to these onsets within 5 s, or that the clocks differ; the outputs cannot tell which. The clock is taken as the files give it, and every session is named here as not established. The shuffle test of the claims below uses the same comparison of real against shifted episodes. In V116 03102018 spikes and events run to 13.46 h against 4.46 h of raw voltage, so their relation to the raw series is not known either; each unit claim is also given without that session.
ph = []
for s in oph:
ev = s["events"]
t = s["t"]
ph.append({"session": s["session"], "fs_hz": round(1 / np.median(np.diff(t)), 2), "dur_s": round(t[-1] - t[0], 1),
"unit": s["unit"], "f_min": round(s["f"].min(), 2), "f_median": round(np.median(s["f"]), 1),
"f_max": round(s["f"].max(), 1), "second_series": s["commanded_voltage"] is not None,
"second_min": None if s["commanded_voltage"] is None else round(s["commanded_voltage"].min(), 2),
"second_max": None if s["commanded_voltage"] is None else round(s["commanded_voltage"].max(), 1),
"pup_out": int((ev.label == "Pup out").sum()), "retrieval": int((ev.label == "retrieval").sum()),
"fiber_location": s["fiber_location"], "fluorophore": s["fluorophore"]})
print(pd.DataFrame(ph).to_string(index=False))
session fs_hz dur_s unit f_min f_median f_max second_series second_min second_max pup_out retrieval fiber_location fluorophore
sub-ROV40_ses-Day 1-obs 1017.25 1134.7 W 1.71 998.2 1253.9 True 1.96 337.3 11 11 unknown unknown
sub-ROV40_ses-Day 1-retrieval 1017.25 1284.8 W 1.53 889.4 1055.0 True 1.27 322.2 10 6 unknown unknown
sub-ROV40_ses-Day 2-obs 1017.25 1160.4 W 0.85 983.1 1398.1 True 3.03 422.2 10 11 unknown unknown
sub-ROV40_ses-Day 2-retrieval 1017.25 1531.5 W 3.62 898.3 1142.6 True 2.92 386.6 10 6 unknown unknown
sub-ROV40_ses-Day 3-obs 1017.25 891.1 W 3.53 971.2 1420.6 True 0.29 471.0 10 10 unknown unknown
sub-ROV40_ses-Day 3-retrieval 1017.25 1634.4 W 2.08 868.7 1160.6 True 4.59 335.8 10 1 unknown unknown
sub-ROV41_ses-Day 1-obs 1017.25 816.5 W 4.38 1154.8 1563.3 True 6.16 357.0 10 10 unknown unknown
sub-ROV41_ses-Day 1-retrieval 1017.25 1676.7 W 3.30 1056.6 1355.8 True 2.92 420.9 10 0 unknown unknown
sub-ROV41_ses-Day 2-obs 1017.25 1408.0 W 0.11 1092.3 1521.2 True 0.07 482.1 11 11 unknown unknown
sub-ROV41_ses-Day 2-retrieval 1017.25 1275.9 W 1.08 1011.5 1366.5 True 0.99 418.5 10 6 unknown unknown
sub-ROV41_ses-Day 3-obs 1017.25 886.7 W 2.56 1027.4 1343.0 True 7.35 391.8 10 10 unknown unknown
sub-ROV41_ses-Day 3-retrieval 1017.25 1714.9 W 2.78 960.5 1181.3 True 3.02 317.1 10 0 unknown unknown
sub-ROV43_ses-Day 1-obs 1017.25 1098.2 W 0.54 895.3 1126.5 False NaN NaN 10 11 unknown unknown
sub-ROV44_ses-Day 1-obs 1017.25 1347.1 W 0.59 1010.7 1284.4 False NaN NaN 10 10 unknown unknown
sub-ROV45_ses-Day 1-obs 1017.25 764.3 W 1.71 859.0 1250.8 False NaN NaN 10 10 unknown unknown
sub-ROV49_ses-Day 1-obs 1017.25 981.2 W 1.44 1145.4 1464.3 False NaN NaN 10 10 unknown unknown
The photometry trace is stored in watts (unit "W") with values of about 1000, which cannot be
watts from a femtowatt photoreceiver; the unit is not the one the values have. The
commanded_voltage series, described as the LED command with power in microwatts, varies from
about 0 to 480 within a session, which is not a constant 400 Hz command at 30 µW; what it holds is
not stated. Fibre location and fluorophore are "unknown". The four animals with only a Day 1
observation file (ROV43, ROV44, ROV45, ROV49) have no retrieval test, so whether they learned is
not in the dandiset; ROV40 retrieved on Day 1 and ROV41 on Day 2. Fig. 4i compares three learners
with three non-learners, so the claim fig4i-projection-photometry is blocked.
# Draft against the pinned version: sizes of all files, and the content the analysis uses in every file.
from dandi.dandiapi import DandiAPIClient
draft = {a.path: a for a in DandiAPIClient().get_dandiset(pvn.DANDISET_ID, "draft").get_assets()}
sizes = pd.DataFrame([{"path": p, "pinned": size, "draft": draft[p].size} for p, size, u in ASSETS])
print("Files in both:", len(sizes), "of", len(ASSETS), "pinned and", len(draft), "draft")
print("Draft minus pinned size (bytes):", sorted(set((sizes.draft - sizes.pinned).tolist())))
same = []
for s in eph:
d = pvn.load_ecephys(s["path"], draft[s["path"]].get_content_url(follow_redirects=1, strip_query=True))
same.append({"path": s["path"], "spikes": len(d["units"]) == len(s["units"]) and all(np.array_equal(a, b) for a, b in zip(d["units"], s["units"])),
"events": d["events"].equals(s["events"]), "series": d["series"] == s["series"]})
for s in oph:
d = pvn.load_ophys(s["path"], draft[s["path"]].get_content_url(follow_redirects=1, strip_query=True))
same.append({"path": s["path"], "spikes": None, "events": d["events"].equals(s["events"]),
"series": np.array_equal(d["t"], s["t"]) and np.array_equal(d["f"], s["f"])})
same = pd.DataFrame(same)
print("Files whose spike times, events and series (shape and rate, or photometry trace) are the same in the draft:",
int((same.spikes.fillna(True) & same.events & same.series).sum()), "of", len(same))
print(same[~(same.spikes.fillna(True) & same.events & same.series)].to_string(index=False))
Files in both: 30 of 30 pinned and 30 draft Draft minus pinned size (bytes): [240, 6144]
Files whose spike times, events and series (shape and rate, or photometry trace) are the same in the draft: 30 of 30 Empty DataFrame Columns: [path, spikes, events, series] Index: []
label-link¶
The dandiset's behaviour labels can be named through the paper's code; "Mom agressing" is shepherding, "Virgin in rest" is nest entry (a nest dwell) and "Mom retrieval" is spontaneous pup retrieval by the dam.
The file V116 03122018 carries both the code's labels and labels with the paper's names. Extended Data Fig. 7b gives the nest dwell times: 70.5% of them last 40 s or less.
s12 = next(s for s in eph if s["session"] == "sub-V116_ses-03122018")
ev = s12["events"]
def coincide(a_lab, b_lab, tol=0.5):
a = ev[ev.label == a_lab].sort_values("start")[["start", "stop"]].to_numpy()
b = ev[ev.label == b_lab].sort_values("start")[["start", "stop"]].to_numpy()
hit = sum(any(np.all(np.abs(x - y) <= tol) for y in b) for x in a)
return len(a), len(b), hit
for a_lab, b_lab in [("Mom agressing", "Shepherding"), ("Mom retrieval", "Self gen retrieval"), ("Virgin in rest", "Virgin in nest")]:
na, nb, hit = coincide(a_lab, b_lab)
print(f"{a_lab!r} ({na}) coinciding with {b_lab!r} ({nb}) within 0.5 s: {hit} of {na}")
rest = ev[ev.label == "Virgin in rest"]; nest = ev[ev.label == "Virgin in nest"]
inside = sum(((nest.start <= r.start + 0.5) & (nest.stop >= r.stop - 0.5)).any() for r in rest.itertuples())
print(f"'Virgin in rest' intervals lying inside a 'Virgin in nest' interval: {inside} of {len(rest)}")
dw = []
for s in eph:
for lab in ("Virgin in rest", "Virgin in nest"):
g = s["events"][s["events"].label == lab]
if len(g):
d = (g.stop - g.start).to_numpy()
dw.append({"session": s["session"], "label": lab, "n": len(d), "le40": int(np.sum(d <= 40)),
"pct_le40": round(100 * np.mean(d <= 40), 1), "median_s": round(np.median(d), 1)})
dw = pd.DataFrame(dw)
print("\n", dw.to_string(index=False))
main = dw[dw.label == "Virgin in rest"]
pct = 100 * main.le40.sum() / main.n.sum()
print(f"\nPooled 'Virgin in rest': {main.le40.sum()} of {main.n.sum()} = {pct:.1f}% last 40 s or less "
f"(paper 70.5%, tolerance 63.5 to 77.6)")
alt = dw[(dw.label == "Virgin in rest") | ~dw.session.isin(main.session)]
print(f"Variant, adding 'Virgin in nest' from files without 'Virgin in rest' (V85): "
f"{100 * alt.le40.sum() / alt.n.sum():.1f}% of {alt.n.sum()}")
by_mouse = main.assign(mouse=main.session.str.split("_").str[0]).groupby("mouse")[["le40", "n"]].sum()
print("By mouse (% of 'Virgin in rest' 40 s or less):", (100 * by_mouse.le40 / by_mouse.n).round(1).to_dict())
'Mom agressing' (7) coinciding with 'Shepherding' (7) within 0.5 s: 7 of 7
'Mom retrieval' (2) coinciding with 'Self gen retrieval' (2) within 0.5 s: 2 of 2
'Virgin in rest' (13) coinciding with 'Virgin in nest' (19) within 0.5 s: 13 of 13
'Virgin in rest' intervals lying inside a 'Virgin in nest' interval: 13 of 13
session label n le40 pct_le40 median_s
sub-V110_ses-03032018 Virgin in rest 72 52 72.2 12.6
sub-V111_ses-03032018 Virgin in rest 21 10 47.6 52.0
sub-V111_ses-03052018 Virgin in rest 133 115 86.5 8.2
sub-V111_ses-03062018 Virgin in rest 189 181 95.8 5.0
sub-V116_ses-03082018 Virgin in rest 22 1 4.5 1318.4
sub-V116_ses-03092018 Virgin in rest 29 5 17.2 1103.7
sub-V116_ses-03102018 Virgin in rest 17 1 5.9 2492.7
sub-V116_ses-03122018 Virgin in rest 13 1 7.7 1970.6
sub-V116_ses-03122018 Virgin in nest 19 2 10.5 796.6
sub-V118_ses-03222018 Virgin in rest 3 0 0.0 64.1
sub-V85_ses-12102017 Virgin in nest 140 128 91.4 4.8
Pooled 'Virgin in rest': 366 of 499 = 73.3% last 40 s or less (paper 70.5%, tolerance 63.5 to 77.6)
Variant, adding 'Virgin in nest' from files without 'Virgin in rest' (V85): 77.3% of 639
By mouse (% of 'Virgin in rest' 40 s or less): {'sub-V110': 72.2, 'sub-V111': 89.2, 'sub-V116': 9.9, 'sub-V118': 0.0}
Result: reproduced. In V116 03122018 all 7 "Mom agressing" intervals coincide with "Shepherding" intervals, both "Mom retrieval" intervals with "Self gen retrieval" intervals, and all 13 "Virgin in rest" intervals lie inside "Virgin in nest" intervals. Pooled over files, 366 of 499 "Virgin in rest" intervals (73.3%) last 40 s or less, within the tolerance of the paper's 70.5%. The pooled value hides a split by mouse: 72.2% in V110 and 89.2% in V111, but 9.9% in V116 and 0 of 3 in V118. In V116 the label marks long dwells (medians 1103.7 to 2492.7 s), so the 40 s cap decides which part of a dwell is measured.
ed7a-nest-entry-fraction¶
A significant fraction of virgin PVN units were activated during nest entry (6 of 21 OT-PVN and 33 of 126 PVN units; one-sided Fisher's exact test p = 0.01 and 0.0001).
Nest entry episodes are "Virgin in rest" intervals, cut at 40 s; the baseline has the episode's duration and ends 0.5 s before it. Each unit's mean change is compared with 1000 shifts of all the session's episodes, and units pass a Benjamini-Hochberg correction over the session's units.
def summarise(tab, name):
k, n = int(tab.activated.sum()), len(tab)
print(f"{name}: {k} of {n} units activated = {100 * k / max(n, 1):.1f}%, suppressed {int(tab.suppressed.sum())}, "
f"Fisher one-sided p = {pvn.fisher_vs_zero(k, n):.2g}; units without a value {int(tab.mi.isna().sum())}")
return k, n
def breakdown(tab):
g = tab.groupby(["mouse", "session"]).agg(units=("unit", "size"), episodes=("episodes", "first"),
activated=("activated", "sum"), suppressed=("suppressed", "sum"),
p_lt_05=("p", lambda x: int(np.sum(x < 0.05))),
mi_median=("mi", "median"))
g["pct_activated"] = (100 * g.activated / g.units).round(1)
print(g.round(2).to_string())
m = tab.groupby("mouse").agg(units=("unit", "size"), activated=("activated", "sum"))
m["pct"] = (100 * m.activated / m.units).round(1)
print(m.to_string())
rng_nest = np.random.default_rng(114001)
nest = pvn.unit_table(with_units, pvn.NEST_ENTRY, "equal", pvn.NEST_CAP_S, rng_nest)
k, n = summarise(nest, "Nest entry (Methods' 40 s cap)")
print("Readings: (a) pooled 26.5% (23.9 to 29.2); (b) PVN only 26.2% (23.6 to 28.8)")
breakdown(nest)
print("\nChange in firing, all units: median %.1f, range %.1f to %.1f" % (nest.mi.median(), nest.mi.min(), nest.mi.max()))
Nest entry (Methods' 40 s cap): 42 of 194 units activated = 21.6%, suppressed 7, Fisher one-sided p = 1.9e-14; units without a value 0
Readings: (a) pooled 26.5% (23.9 to 29.2); (b) PVN only 26.2% (23.6 to 28.8)
units episodes activated suppressed p_lt_05 mi_median pct_activated
mouse session
V110 sub-V110_ses-03032018 31 72 1 0 8 4.20 3.2
V111 sub-V111_ses-03032018 31 21 0 0 0 -6.11 0.0
sub-V111_ses-03052018 31 133 2 1 3 -0.23 6.5
sub-V111_ses-03062018 31 189 2 6 8 -0.55 6.5
V116 sub-V116_ses-03082018 23 22 8 0 12 17.51 34.8
sub-V116_ses-03092018 21 29 12 0 12 44.80 57.1
sub-V116_ses-03102018 26 17 17 0 20 22.36 65.4
units activated pct
mouse
V110 31 1 3.2
V111 93 4 4.3
V116 70 37 52.9
Change in firing, all units: median 1.6, range -27.8 to 75.1
# Variants: the code's uncapped episodes; without V110 and V111; without units above the spreadsheet row limit;
# other thresholds.
rng_nest_v = np.random.default_rng(114002)
nest_nocap = pvn.unit_table(with_units, pvn.NEST_ENTRY, "equal", None, rng_nest_v)
summarise(nest_nocap, "Variant, no 40 s cap (code)")
clean = nest[~nest.mouse.isin(["V110", "V111"])]
summarise(clean, "Variant, V116 only (V110 and V111 left out for unit quality)")
summarise(nest[nest.spikes <= pvn.EXCEL_ROWS], "Variant, units with at most 1,048,500 spikes")
for qq in (0.01, 0.1):
act = np.concatenate([pvn.bh(g.p.to_numpy(), qq) & (g.mi.to_numpy() > 0) for _, g in nest.groupby("session", sort=False)])
print(f"Variant, q = {qq}: {act.sum()} of {len(nest)} = {100 * act.mean():.1f}%")
raw = (nest.p < 0.05) & (nest.mi > 0)
print(f"Variant, p < 0.05 without correction: {raw.sum()} of {len(nest)} = {100 * raw.mean():.1f}%")
def common_variants(tab, labels, baseline, cap, seed):
"""Variants shared by the three fraction claims: empty episodes, false-positive level, V111 once, the clock."""
rng = np.random.default_rng(seed)
summarise(pvn.unit_table(with_units, labels, baseline, cap, rng, real_empty_as_zero=True),
"Variant, real episodes without spikes counted as 0, as in the shuffles")
print(f"Units with at least one real episode left out for lack of spikes: {int((tab.n_empty > 0).sum())} of {len(tab)}; "
f"median share of episodes left out among them {(tab.n_empty / tab.episodes)[tab.n_empty > 0].median():.2f}")
null = pvn.unit_table(pvn.poisson_sessions(with_units, rng), labels, baseline, cap, rng)
print(f"False-positive level, Poisson units with each unit's spike count and no response, same episodes and test: "
f"{int(null.activated.sum())} of {len(null)} activated ({100 * null.activated.mean():.1f}%), "
f"{int(null.suppressed.sum())} suppressed")
for keep in sorted(tab.session[tab.mouse == "V111"].unique()):
summarise(tab[(tab.mouse != "V111") | (tab.session == keep)], f"Variant, V111 counted once ({keep} only)")
summarise(tab[tab.session != "sub-V116_ses-03102018"], "Variant, without V116 03102018 (clock unknown)")
common_variants(nest, pvn.NEST_ENTRY, "equal", pvn.NEST_CAP_S, 114007)
# Baselines of capped nest entries that overlap an earlier nest dwell or a shepherding episode.
ov = []
for s in with_units:
st, en = pvn.episodes(s["events"], pvn.NEST_ENTRY, pvn.NEST_CAP_S)
if not len(st):
continue
b0, b1 = pvn.baseline_windows(st, en, "equal")
ns, ne = pvn.episodes(s["events"], pvn.NEST_ENTRY)
ss, se = pvn.episodes(s["events"], pvn.SHEPHERDING)
prev = np.array([np.any((ns < b1[i]) & (ne > b0[i]) & (ns < st[i])) for i in range(len(st))])
shp = np.array([np.any((ss < b1[i]) & (se > b0[i])) for i in range(len(st))])
ov.append({"session": s["session"], "entries": len(st), "baseline_in_earlier_dwell": int(prev.sum()),
"baseline_in_shepherding": int(shp.sum()), "either": int((prev | shp).sum())})
print("\nBaselines overlapping earlier episodes:\n", pd.DataFrame(ov).to_string(index=False))
fig, ax = plt.subplots(1, 2, figsize=(10, 3.2))
ax[0].hist(nest.mi.dropna(), bins=40, color="0.6")
ax[0].hist(nest.mi[nest.activated], bins=40, color="C3", label="activated")
ax[0].set_xlabel("change in firing during nest entry (%)"); ax[0].set_ylabel("units"); ax[0].legend()
ax[1].hist(nest.p.dropna(), bins=20, color="0.5"); ax[1].set_xlabel("shuffle p"); ax[1].set_ylabel("units")
plt.tight_layout(); plt.show()
print("Plotted: all units with a value; change from %.1f to %.1f, p from %.3f to %.3f" %
(nest.mi.min(), nest.mi.max(), nest.p.min(), nest.p.max()))
Variant, no 40 s cap (code): 2 of 194 units activated = 1.0%, suppressed 7, Fisher one-sided p = 0.25; units without a value 0 Variant, V116 only (V110 and V111 left out for unit quality): 37 of 70 units activated = 52.9%, suppressed 0, Fisher one-sided p = 1e-14; units without a value 0 Variant, units with at most 1,048,500 spikes: 42 of 174 units activated = 24.1%, suppressed 7, Fisher one-sided p = 1.4e-14; units without a value 0 Variant, q = 0.01: 28 of 194 = 14.4% Variant, q = 0.1: 55 of 194 = 28.4% Variant, p < 0.05 without correction: 55 of 194 = 28.4%
Variant, real episodes without spikes counted as 0, as in the shuffles: 31 of 194 units activated = 16.0%, suppressed 1, Fisher one-sided p = 1.3e-10; units without a value 0 Units with at least one real episode left out for lack of spikes: 116 of 194; median share of episodes left out among them 0.17
False-positive level, Poisson units with each unit's spike count and no response, same episodes and test: 2 of 194 activated (1.0%), 1 suppressed
Variant, V111 counted once (sub-V111_ses-03032018 only): 38 of 132 units activated = 28.8%, suppressed 0, Fisher one-sided p = 1.6e-13; units without a value 0
Variant, V111 counted once (sub-V111_ses-03052018 only): 40 of 132 units activated = 30.3%, suppressed 1, Fisher one-sided p = 2.8e-14; units without a value 0
Variant, V111 counted once (sub-V111_ses-03062018 only): 40 of 132 units activated = 30.3%, suppressed 6, Fisher one-sided p = 2.8e-14; units without a value 0
Variant, without V116 03102018 (clock unknown): 25 of 168 units activated = 14.9%, suppressed 7, Fisher one-sided p = 1.1e-08; units without a value 0
Baselines overlapping earlier episodes:
session entries baseline_in_earlier_dwell baseline_in_shepherding either
sub-V110_ses-03032018 72 20 1 20
sub-V111_ses-03032018 21 4 1 5
sub-V111_ses-03052018 133 19 0 19
sub-V111_ses-03062018 189 19 0 19
sub-V116_ses-03082018 22 2 1 3
sub-V116_ses-03092018 29 3 16 16
sub-V116_ses-03102018 17 2 8 8
Plotted: all units with a value; change from -27.8 to 75.1, p from 0.000 to 0.504
Result: partially reproduced. 42 of 194 units (21.6%) are activated during nest entry, below the tolerance of both readings (23.9 to 29.2% pooled, 23.6 to 28.8% PVN only), and the Fisher test against zero gives p = 1.9e-14. Almost all activated units are from V116 (37 of 70, 52.9%); V110 has 1 of 31 and V111 4 of 93. The result depends on the 40 s cap: without it, as in the code, 2 of 194 units (1.0%) are activated (p = 0.25). Leaving out the 20 units above the spreadsheet row limit gives 42 of 174 (24.1%), inside the tolerance of both readings, and q = 0.01 and 0.1 give 14.4% and 28.4%.
Poisson units with each unit's spike count and no response, on the same episodes and test, give 2 of 194 activated (1.0%), so the 42 are far above the test's false-positive level. Real episodes without spikes are left out in 116 of 194 units (a median 17% of episodes among them); counting them as 0, as the shuffles do, gives 31 of 194 (16.0%, p = 1.3e-10). Counting V111 once, with any one of its three sessions, gives 28.8 to 30.3% of 132 units; 28.8% (with V111 03032018) is inside the pooled reading's tolerance. Without V116 03102018, whose clock is not known, 25 of 168 (14.9%, p = 1.1e-08). The baselines of 3 to 20 entries per session overlap an earlier nest dwell or a shepherding episode; the largest count is in V110 03032018 (20 of 72, all an earlier dwell) and the largest share in V116 03092018 (16 of 29, all overlapping shepherding).
ed7a-shepherding-fraction¶
A significant fraction of virgin PVN units were activated during shepherding by the dam (5 of 21 OT-PVN and 9 of 97 PVN units; p = 0.02 and 0.001).
Shepherding episodes are "Mom agressing" intervals; the baseline is the 10 s before each, as the Methods state. The code's equal-duration baseline is a variant.
rng_shep = np.random.default_rng(114003)
shep = pvn.unit_table(with_units, pvn.SHEPHERDING, "pre10", None, rng_shep)
summarise(shep, "Shepherding (10 s baseline, Methods)")
print("Readings: (a) pooled 11.9% (10.7 to 13.1); (b) PVN only 9.3% (8.3 to 10.2)")
breakdown(shep)
rng_shep_v = np.random.default_rng(114004)
shep_code = pvn.unit_table(with_units, pvn.SHEPHERDING, "equal", None, rng_shep_v)
summarise(shep_code, "Variant, equal-duration baseline ending 0.5 s before (code)")
summarise(shep[~shep.mouse.isin(["V110", "V111"])], "Variant, V116 only")
raw = (shep.p < 0.05) & (shep.mi > 0)
print(f"Variant, p < 0.05 without correction: {raw.sum()} of {len(shep)} = {100 * raw.mean():.1f}%")
print("Change in firing: median %.1f, range %.1f to %.1f" % (shep.mi.median(), shep.mi.min(), shep.mi.max()))
common_variants(shep, pvn.SHEPHERDING, "pre10", None, 114008)
Shepherding (10 s baseline, Methods): 3 of 194 units activated = 1.5%, suppressed 10, Fisher one-sided p = 0.12; units without a value 15
Readings: (a) pooled 11.9% (10.7 to 13.1); (b) PVN only 9.3% (8.3 to 10.2)
units episodes activated suppressed p_lt_05 mi_median pct_activated
mouse session
V110 sub-V110_ses-03032018 31 4 0 0 0 -11.69 0.0
V111 sub-V111_ses-03032018 31 4 0 0 0 13.20 0.0
sub-V111_ses-03052018 31 2 1 1 2 28.39 3.2
sub-V111_ses-03062018 31 1 0 4 5 -65.89 0.0
V116 sub-V116_ses-03082018 23 6 0 0 2 -7.79 0.0
sub-V116_ses-03092018 21 22 1 0 5 -3.12 4.8
sub-V116_ses-03102018 26 13 1 5 8 -9.06 3.8
units activated pct
mouse
V110 31 0 0.0
V111 93 1 1.1
V116 70 2 2.9
Variant, equal-duration baseline ending 0.5 s before (code): 9 of 194 units activated = 4.6%, suppressed 14, Fisher one-sided p = 0.0018; units without a value 12 Variant, V116 only: 2 of 70 units activated = 2.9%, suppressed 5, Fisher one-sided p = 0.25; units without a value 2 Variant, p < 0.05 without correction: 5 of 194 = 2.6% Change in firing: median -8.1, range -100.0 to 100.0
Variant, real episodes without spikes counted as 0, as in the shuffles: 1 of 194 units activated = 0.5%, suppressed 5, Fisher one-sided p = 0.5; units without a value 0 Units with at least one real episode left out for lack of spikes: 70 of 194; median share of episodes left out among them 0.50
False-positive level, Poisson units with each unit's spike count and no response, same episodes and test: 3 of 194 activated (1.5%), 4 suppressed Variant, V111 counted once (sub-V111_ses-03032018 only): 2 of 132 units activated = 1.5%, suppressed 5, Fisher one-sided p = 0.25; units without a value 2 Variant, V111 counted once (sub-V111_ses-03052018 only): 3 of 132 units activated = 2.3%, suppressed 6, Fisher one-sided p = 0.12; units without a value 7 Variant, V111 counted once (sub-V111_ses-03062018 only): 2 of 132 units activated = 1.5%, suppressed 9, Fisher one-sided p = 0.25; units without a value 10 Variant, without V116 03102018 (clock unknown): 2 of 168 units activated = 1.2%, suppressed 5, Fisher one-sided p = 0.25; units without a value 15
Result: not reproduced. With the Methods' 10 s baseline, 3 of 194 units (1.5%) are activated during shepherding (Fisher p = 0.12), against 11.9% pooled or 9.3% PVN only in the paper; 10 units are suppressed and 15 have no value. With the code's equal-duration baseline, 9 of 194 (4.6%) are activated (p = 0.0018), still below both readings. V116 alone gives 2 of 70 (p = 0.25). Sessions hold 1 to 22 shepherding episodes. Poisson units with no response give 3 of 194 activated on the same episodes, as many as the real units. Counting empty episodes as 0 gives 1 of 194; counting V111 once gives 2 or 3 of 132; leaving out V116 03102018 gives 2 of 168; none has p < 0.05.
ed7a-dam-retrieval-fraction¶
Spontaneous pup retrieval by the dam, in which the virgin was not involved, activated a significant fraction of virgin PVN units (5 of 17 OT-PVN and 10 of 69 PVN units; p = 0.02 and 0.0007).
Dam retrieval episodes are "Mom retrieval" intervals, with the 10 s baseline of the Methods.
rng_ret = np.random.default_rng(114005)
ret = pvn.unit_table(with_units, pvn.DAM_RETRIEVAL, "pre10", None, rng_ret)
summarise(ret, "Dam retrieval (10 s baseline, Methods)")
print("Readings: (a) pooled 17.4% (15.7 to 19.1); (b) PVN only 14.5% (13.1 to 16.0)")
breakdown(ret)
rng_ret_v = np.random.default_rng(114006)
ret_code = pvn.unit_table(with_units, pvn.DAM_RETRIEVAL, "equal", None, rng_ret_v)
summarise(ret_code, "Variant, equal-duration baseline ending 0.5 s before (code)")
summarise(ret[ret.episodes >= 2], "Variant, sessions with at least 2 episodes")
raw = (ret.p < 0.05) & (ret.mi > 0)
print(f"Variant, p < 0.05 without correction: {raw.sum()} of {len(ret)} = {100 * raw.mean():.1f}%")
print("Change in firing: median %.1f, range %.1f to %.1f" % (ret.mi.median(), ret.mi.min(), ret.mi.max()))
common_variants(ret, pvn.DAM_RETRIEVAL, "pre10", None, 114009)
# Are the activated units of V111 03032018 separate neurons? Fraction of each unit's spikes within 0.1 ms of the other's.
s3 = next(s for s in with_units if s["session"] == "sub-V111_ses-03032018")
ids = ret.unit[(ret.session == s3["session"]) & ret.activated].tolist()
rest_ids = [i for i in range(len(s3["units"])) if i not in ids]
pair = lambda group: np.array([pvn.coincident_pair(s3["units"][i], s3["units"][j]) for i in group for j in group if i != j])
pa, pr = pair(ids), pair(rest_ids)
print(f"V111 03032018, {len(ids)} activated units: shared spikes per ordered pair median {np.median(pa):.2f} "
f"(range {pa.min():.2f} to {pa.max():.2f}); the other {len(rest_ids)} units: median {np.median(pr):.3f}")
Dam retrieval (10 s baseline, Methods): 13 of 132 units activated = 9.8%, suppressed 20, Fisher one-sided p = 9e-05; units without a value 4
Readings: (a) pooled 17.4% (15.7 to 19.1); (b) PVN only 14.5% (13.1 to 16.0)
units episodes activated suppressed p_lt_05 mi_median pct_activated
mouse session
V111 sub-V111_ses-03032018 31 1 9 1 15 98.10 29.0
sub-V111_ses-03062018 31 2 2 1 3 -29.41 6.5
V116 sub-V116_ses-03082018 23 4 1 1 6 -33.33 4.3
sub-V116_ses-03092018 21 11 0 10 10 -40.65 0.0
sub-V116_ses-03102018 26 2 1 7 9 -35.73 3.8
units activated pct
mouse
V111 62 11 17.7
V116 70 2 2.9
Variant, equal-duration baseline ending 0.5 s before (code): 16 of 132 units activated = 12.1%, suppressed 16, Fisher one-sided p = 9.4e-06; units without a value 11 Variant, sessions with at least 2 episodes: 4 of 101 units activated = 4.0%, suppressed 19, Fisher one-sided p = 0.061; units without a value 1 Variant, p < 0.05 without correction: 18 of 132 = 13.6% Change in firing: median -22.4, range -100.0 to 100.0
Variant, real episodes without spikes counted as 0, as in the shuffles: 12 of 132 units activated = 9.1%, suppressed 8, Fisher one-sided p = 0.00019; units without a value 0 Units with at least one real episode left out for lack of spikes: 50 of 132; median share of episodes left out among them 0.50
False-positive level, Poisson units with each unit's spike count and no response, same episodes and test: 3 of 132 activated (2.3%), 12 suppressed Variant, V111 counted once (sub-V111_ses-03032018 only): 11 of 101 units activated = 10.9%, suppressed 19, Fisher one-sided p = 0.00037; units without a value 3 Variant, V111 counted once (sub-V111_ses-03062018 only): 4 of 101 units activated = 4.0%, suppressed 19, Fisher one-sided p = 0.061; units without a value 1 Variant, without V116 03102018 (clock unknown): 12 of 106 units activated = 11.3%, suppressed 13, Fisher one-sided p = 0.00018; units without a value 4
V111 03032018, 9 activated units: shared spikes per ordered pair median 0.71 (range 0.33 to 0.81); the other 22 units: median 0.017
# V111 03032018 has one dam retrieval. Spike counts of each unit in the 10 s baseline and in the episode,
# and the population's spike count in 5 s bins over the 60 s around the episode's start.
r0, r1 = pvn.episodes(s3["events"], pvn.DAM_RETRIEVAL)
print(f"V111 03032018 retrieval: {r0[0]:.1f} to {r1[0]:.1f} s ({r1[0] - r0[0]:.1f} s); baseline {r0[0] - pvn.PRE_S:.1f} to {r0[0]:.1f} s")
cnt = pd.DataFrame([{"unit": i, "activated": i in ids,
"baseline_spikes": int(np.sum((u >= r0[0] - pvn.PRE_S) & (u < r0[0]))),
"episode_spikes": int(np.sum((u >= r0[0]) & (u <= r1[0]))),
"session_rate_hz": round(len(u) / (u.max() - u.min()), 2)} for i, u in enumerate(s3["units"])])
cnt["baseline_hz"] = (cnt.baseline_spikes / pvn.PRE_S).round(2)
cnt["episode_hz"] = (cnt.episode_spikes / (r1[0] - r0[0])).round(2)
cnt["change"] = ret.set_index(["session", "unit"]).loc[s3["session"]].mi.round(1).to_numpy()
print(cnt.to_string(index=False))
print(f"Units with no spike in the baseline: {int((cnt.baseline_spikes == 0).sum())} of {len(cnt)} "
f"(activated {int(((cnt.baseline_spikes == 0) & cnt.activated).sum())} of {len(ids)})")
allsp = np.sort(np.concatenate(s3["units"]))
edges = np.arange(r0[0] - 30, r0[0] + 30.1, 5.0)
pop = np.histogram(allsp, edges)[0]
print("Population spikes per 5 s bin, from 30 s before to 30 s after the episode's start:", pop.tolist())
print(f"Session population rate {len(allsp) / (allsp[-1] - allsp[0]) * 5:.0f} spikes per 5 s; "
f"longest interval without any spike within 60 s of the start {np.diff(allsp[(allsp >= r0[0] - 60) & (allsp <= r0[0] + 60)]).max():.2f} s")
# The same counts for every 10 s window of the session: how often is the population this quiet?
w = np.histogram(allsp, np.arange(allsp[0], allsp[-1], pvn.PRE_S))[0]
b = int(np.sum((allsp >= r0[0] - pvn.PRE_S) & (allsp < r0[0])))
print(f"Population spikes in the 10 s baseline: {b}; share of the session's 10 s windows with at most as many: {np.mean(w <= b):.4f} "
f"(median window {np.median(w):.0f})")
w5 = np.histogram(allsp, np.arange(allsp[0], allsp[-1], 5.0))[0]
print(f"Session's 5 s windows with at least {pop[6]} population spikes, as in the episode's first 5 s: {int(np.sum(w5 >= pop[6]))} of {len(w5)}; "
f"with at least 5,000: {int(np.sum(w5 >= 5000))}")
V111 03032018 retrieval: 527.4 to 532.6 s (5.2 s); baseline 517.4 to 527.4 s
unit activated baseline_spikes episode_spikes session_rate_hz baseline_hz episode_hz change
0 False 275 139 24.04 27.5 26.65 -1.6
1 False 2 7 1.97 0.2 1.34 74.1
2 True 0 800 15.96 0.0 153.41 100.0
3 True 0 772 15.33 0.0 148.04 100.0
4 False 11 9 2.10 1.1 1.73 22.1
5 False 5 4 1.09 0.5 0.77 21.1
6 False 2 752 15.74 0.2 144.20 99.7
7 False 4 831 18.16 0.4 159.35 99.5
8 False 10 17 2.73 1.0 3.26 53.1
9 False 1 4 0.66 0.1 0.77 76.9
10 False 4 916 19.80 0.4 175.65 99.5
11 False 2 819 15.75 0.2 157.05 99.7
12 False 85 60 6.45 8.5 11.51 15.0
13 False 10 4 1.49 1.0 0.77 -13.2
14 True 0 781 15.56 0.0 149.76 100.0
15 True 0 768 15.42 0.0 147.27 100.0
16 False 0 0 0.30 0.0 0.00 NaN
17 False 3 2 0.37 0.3 0.38 12.2
18 False 3 863 17.84 0.3 165.49 99.6
19 False 10 948 20.87 1.0 181.79 98.9
20 False 0 0 0.24 0.0 0.00 NaN
21 False 0 0 0.23 0.0 0.00 NaN
22 True 19 1022 23.70 1.9 195.98 98.1
23 False 21 1041 24.31 2.1 199.62 97.9
24 False 4 0 0.27 0.4 0.00 -100.0
25 False 4 13 3.48 0.4 2.49 72.3
26 True 10 975 22.14 1.0 186.97 98.9
27 True 11 1016 23.06 1.1 194.83 98.9
28 False 1 31 1.11 0.1 5.94 96.7
29 True 21 1154 28.83 2.1 221.29 98.1
30 True 19 1150 28.54 1.9 220.52 98.3
Units with no spike in the baseline: 7 of 31 (activated 4 of 9)
Population spikes per 5 s bin, from 30 s before to 30 s after the episode's start: [269, 255, 226, 228, 246, 291, 14346, 1185, 220, 218, 243, 3864]
Session population rate 1838 spikes per 5 s; longest interval without any spike within 60 s of the start 0.22 s
Population spikes in the 10 s baseline: 537; share of the session's 10 s windows with at most as many: 0.4344 (median window 730) Session's 5 s windows with at least 14346 population spikes, as in the episode's first 5 s: 12 of 3325; with at least 5,000: 453
Result: partially reproduced. 13 of 132 units (9.8%) are activated during dam retrieval (Fisher p = 9e-05), below the tolerance of both readings (15.7 to 19.1% pooled, 13.1 to 16.0% PVN only). 9 of the 13 come from V111 03032018, which has a single retrieval episode, and these 9 share a median 71% of their spikes pairwise within 0.1 ms (33 to 81%), against 1.7% among the session's other 22 units, so they are probably far fewer than 9 neurons. Over sessions with at least 2 episodes, 4 of 101 units are activated (p = 0.061). Counting V111 once gives 11 of 101 with V111 03032018 (p = 0.00037) and 4 of 101 with V111 03062018 (p = 0.061). More units are suppressed (20) than activated, and the median change is -22.4; Poisson units with no response give 3 activated and 12 suppressed of 132, so the suppressed count is not far above the test's false-positive level. Counting empty episodes as 0 gives 12 of 132 (p = 0.00019), and leaving out V116 03102018 gives 12 of 106 (p = 0.00018). With the code's baseline, 16 of 132 (12.1%) are activated.
The baseline of the single retrieval in V111 03032018 is not a gap in the recording: the population fires 537 spikes in it, more than in 43% of the session's 10 s windows, and no interval within 60 s of the episode is longer than 0.22 s. 7 of the 31 units have no spike in the baseline (4 of the 9 activated ones), and the 9 fire 0 to 2.1 Hz in it. During the 5.2 s episode about 16 units fire at 144 to 221 Hz together, and the population fires 14,346 spikes in the first 5 s against 1,838 per 5 s over the session; 12 of the session's 3,325 windows of 5 s hold as many, and 453 hold at least 5,000. The retrieval coincides with one of the session's largest synchronous bursts, among units that share most of their spikes. These outputs do not tell a response from a burst that happened to fall on the episode, so the outcome follows the criterion, which counts the units as given.
fig2j-shepherded-entry¶
PVN units fired more when the virgin entered the nest after being shepherded than when she entered it voluntarily (18 units, Wilcoxon matched-pairs signed rank test, p = 0.0056).
A nest entry is shepherded when it starts between the start of a "Mom agressing" episode and 10 s after its end.
def split_entries(s, window):
st, en = pvn.episodes(s["events"], pvn.NEST_ENTRY, pvn.NEST_CAP_S)
ss, se = pvn.episodes(s["events"], pvn.SHEPHERDING)
shep_mask = np.array([np.any((x >= ss) & (x <= se + window)) for x in st], bool)
return st, en, shep_mask
def shepherded_vs_voluntary(window, units_filter=None):
rows = []
for s in with_units:
st, en, m = split_entries(s, window)
if m.sum() == 0 or (~m).sum() == 0:
continue
for i, u in enumerate(s["units"]):
if units_filter is not None and (s["session"], i) not in units_filter:
continue
a = pvn.modulation(u, st[m], en[m], "equal")[0]
b = pvn.modulation(u, st[~m], en[~m], "equal")[0]
rows.append({"session": s["session"], "mouse": s["subject"]["subject_id"], "unit": i,
"shepherded_entries": int(m.sum()), "voluntary_entries": int((~m).sum()), "shep": a, "vol": b})
t = pd.DataFrame(rows)
ok = t.dropna(subset=["shep", "vol"])
w = stats.wilcoxon(ok.shep, ok.vol) if len(ok) > 1 else None
return t, ok, w
t10, ok10, w10 = shepherded_vs_voluntary(10)
print("Sessions and entries:\n", t10.groupby("session")[["shepherded_entries", "voluntary_entries"]].first().to_string())
print(f"\nWindow 10 s: {len(ok10)} units with both values ({len(t10) - len(ok10)} dropped without spikes); "
f"median shepherded {ok10.shep.median():.1f}, voluntary {ok10.vol.median():.1f}, median difference "
f"{(ok10.shep - ok10.vol).median():.1f}; higher for shepherded in {int((ok10.shep > ok10.vol).sum())}; "
f"Wilcoxon p = {w10.pvalue:.3g}")
for sess, g in ok10.groupby("session"):
ww = stats.wilcoxon(g.shep, g.vol) if len(g) > 5 else None
print(f" {sess}: {len(g)} units, median difference {(g.shep - g.vol).median():.1f}, "
f"shepherded higher in {int((g.shep > g.vol).sum())}" + (f", p = {ww.pvalue:.3g}" if ww else ""))
for win in (5, 30):
t, ok, w = shepherded_vs_voluntary(win)
print(f"Variant, window {win} s: {len(ok)} units, median difference {(ok.shep - ok.vol).median():.1f}, p = {w.pvalue:.3g}")
act = set(zip(nest.session[nest.activated], nest.unit[nest.activated]))
t, ok, w = shepherded_vs_voluntary(10, act)
if len(ok) > 1:
print(f"Variant, units activated during nest entry: {len(ok)} units, median difference {(ok.shep - ok.vol).median():.1f}, p = {w.pvalue:.3g}")
else:
print(f"Variant, units activated during nest entry: {len(ok)} units with both values; no test")
t, ok, w = shepherded_vs_voluntary(10)
ok = ok[~ok.mouse.isin(["V110", "V111"])]
print(f"Variant, V116 only: {len(ok)} units, median difference {(ok.shep - ok.vol).median():.1f}, p = {stats.wilcoxon(ok.shep, ok.vol).pvalue:.3g}")
for keep in sorted(ok10.session[ok10.mouse == "V111"].unique()):
g = ok10[(ok10.mouse != "V111") | (ok10.session == keep)]
print(f"Variant, V111 counted once ({keep} only): {len(g)} units, median difference {(g.shep - g.vol).median():.1f}, "
f"p = {stats.wilcoxon(g.shep, g.vol).pvalue:.3g}")
g = ok10[ok10.session != "sub-V116_ses-03102018"]
print(f"Variant, without V116 03102018: {len(g)} units, median difference {(g.shep - g.vol).median():.1f}, "
f"p = {stats.wilcoxon(g.shep, g.vol).pvalue:.3g}")
print("Plotted range: shepherded %.1f to %.1f, voluntary %.1f to %.1f" % (ok10.shep.min(), ok10.shep.max(), ok10.vol.min(), ok10.vol.max()))
fig, ax = plt.subplots(figsize=(3.5, 3.5))
ax.scatter(ok10.vol, ok10.shep, s=10, c="0.3"); ax.plot([-100, 100], [-100, 100], c="C3", lw=0.8)
ax.set_xlabel("voluntary entry change (%)"); ax.set_ylabel("shepherded entry change (%)"); plt.tight_layout(); plt.show()
Sessions and entries:
shepherded_entries voluntary_entries
session
sub-V110_ses-03032018 1 71
sub-V111_ses-03032018 1 20
sub-V116_ses-03082018 1 21
sub-V116_ses-03092018 15 14
sub-V116_ses-03102018 8 9
Window 10 s: 128 units with both values (4 dropped without spikes); median shepherded 24.8, voluntary 7.8, median difference 3.9; higher for shepherded in 70; Wilcoxon p = 0.182
sub-V110_ses-03032018: 31 units, median difference 22.2, shepherded higher in 30, p = 4.66e-09
sub-V111_ses-03032018: 31 units, median difference 35.8, shepherded higher in 16, p = 0.247
sub-V116_ses-03082018: 19 units, median difference 3.2, shepherded higher in 11, p = 0.651
sub-V116_ses-03092018: 21 units, median difference -4.9, shepherded higher in 7, p = 0.374
sub-V116_ses-03102018: 26 units, median difference -12.0, shepherded higher in 6, p = 0.00669
Variant, window 5 s: 128 units, median difference 3.1, p = 0.37
Variant, window 30 s: 128 units, median difference 4.2, p = 0.166
Variant, units activated during nest entry: 36 units, median difference -7.3, p = 0.148 Variant, V116 only: 66 units, median difference -6.0, p = 0.0732 Variant, V111 counted once (sub-V111_ses-03032018 only): 128 units, median difference 3.9, p = 0.182 Variant, without V116 03102018: 102 units, median difference 10.8, p = 0.0237 Plotted range: shepherded -100.0 to 100.0, voluntary -17.0 to 94.9
Result: not reproduced. Over 128 units, the change is higher for shepherded than for voluntary entries in 70, with a median difference of 3.9 and Wilcoxon p = 0.182. Sessions disagree: in V110 03032018, with one shepherded entry, 30 of 31 units are higher (p = 4.66e-09), while in V116 03102018 6 of 26 are (p = 0.00669). Windows of 5 and 30 s (p = 0.37 and 0.166), the 36 units activated during nest entry (median difference -7.3, p = 0.148) and V116 alone (-6.0, p = 0.0732) do not give a significant difference either. Only one V111 session enters the test, so counting V111 once changes nothing. Leaving out V116 03102018, whose clock is not known, gives a median difference of 10.8 over 102 units and p = 0.0237; this variant does not decide the outcome, and 31 of its units come from V110 03032018, compared on a single shepherded entry.
fig1j-shepherding-rate¶
Dams shepherded virgins more often than the 0.2 events per hour at which dams chased virgins without pups, at 3.5 +/- 0.9 events per hour with pups present.
Rates are per session, over the sessions whose annotations include a label for the virgin.
rate_rows = []
for s in eph:
ev = s["events"]
if not ev.label.isin(pvn.VIRGIN_LABELS).any():
continue
hours = pvn.session_length_s(s) / 3600
rate_rows.append({"session": s["session"], "mouse": s["subject"]["subject_id"], "hours": round(hours, 2),
"shepherding": int(ev.label.isin(pvn.SHEPHERDING).sum()),
"dam_retrieval": int(ev.label.isin(pvn.DAM_RETRIEVAL).sum())})
rates = pd.DataFrame(rate_rows)
rates["shep_per_h"] = rates.shepherding / rates.hours
rates["ret_per_h"] = rates.dam_retrieval / rates.hours
print(rates.round(3).to_string(index=False))
print("Left out, no virgin label:", [s["session"] for s in eph if not s["events"].label.isin(pvn.VIRGIN_LABELS).any()])
t = stats.ttest_1samp(rates.shep_per_h, 0.2)
print(f"\nShepherding: mean {rates.shep_per_h.mean():.2f} per h (s.e.m. {rates.shep_per_h.sem():.2f}, n = {len(rates)} sessions), "
f"against 0.2: t = {t.statistic:.2f}, two-sided p = {t.pvalue:.3g}; paper 3.5 +/- 0.9")
pm = rates.groupby("mouse").shep_per_h.mean()
tm = stats.ttest_1samp(pm, 0.2)
print("Per mouse:", pm.round(2).to_dict(), f"; mean {pm.mean():.2f}, t-test against 0.2 p = {tm.pvalue:.3g}")
print(f"Variant, sessions with units only: mean {rates[rates.session.isin([s['session'] for s in with_units])].shep_per_h.mean():.2f} per h")
session mouse hours shepherding dam_retrieval shep_per_h ret_per_h
sub-V110_ses-03032018 V110 12.04 4 0 0.332 0.000
sub-V111_ses-03032018 V111 4.62 4 1 0.866 0.216
sub-V111_ses-03052018 V111 10.77 2 0 0.186 0.000
sub-V111_ses-03062018 V111 12.41 1 2 0.081 0.161
sub-V116_ses-03082018 V116 12.22 6 4 0.491 0.327
sub-V116_ses-03092018 V116 15.82 22 11 1.391 0.695
sub-V116_ses-03102018 V116 13.45 13 2 0.967 0.149
sub-V116_ses-03122018 V116 12.53 7 2 0.559 0.160
sub-V118_ses-03222018 V118 2.52 0 1 0.000 0.397
sub-V85_ses-12102017 V85 12.01 19 0 1.582 0.000
Left out, no virgin label: ['sub-V118_ses-03182018', 'sub-V118_ses-03192018', 'sub-V118_ses-03212018', 'sub-V85_ses-12102017-2']
Shepherding: mean 0.65 per h (s.e.m. 0.17, n = 10 sessions), against 0.2: t = 2.59, two-sided p = 0.0292; paper 3.5 +/- 0.9
Per mouse: {'V110': 0.33, 'V111': 0.38, 'V116': 0.85, 'V118': 0.0, 'V85': 1.58} ; mean 0.63, t-test against 0.2 p = 0.193
Variant, sessions with units only: mean 0.62 per h
Result: partially reproduced. Over 10 sessions dams shepherd at a mean of 0.65 events per hour (s.e.m. 0.17), above 0.2 (t = 2.59, p = 0.0292), but far below the 3.5 +/- 0.9 of Extended Data Fig. 3b. Averaged per mouse first, the mean is 0.63 per hour and the test against 0.2 gives p = 0.193. Session rates range from 0 (V118 03222018, 2.52 h) to 1.58 per hour (V85). The virgin labels span each session with gaps of up to 1.93 h (data check), which is consistent with sessions scored in full but does not show that every shepherding episode was annotated.
ed8a-dam-retrieval-rate¶
About once every two hours, dams dropped and then retrieved pups during co-housing (daily means 0.34 to 0.60 retrievals per hour, N = 14 dams).
print(rates[["session", "mouse", "hours", "dam_retrieval", "ret_per_h"]].round(3).to_string(index=False))
print(f"\nMean {rates.ret_per_h.mean():.2f} per h over {len(rates)} sessions (median {rates.ret_per_h.median():.2f}, "
f"range {rates.ret_per_h.min():.2f} to {rates.ret_per_h.max():.2f}); total {rates.dam_retrieval.sum()} in {rates.hours.sum():.1f} h "
f"= {rates.dam_retrieval.sum() / rates.hours.sum():.2f} per h")
pm = rates.groupby("mouse").ret_per_h.mean()
print("Per mouse:", pm.round(2).to_dict())
print("Paper: 0.5 per h (0.45 to 0.55), or the daily range 0.34 to 0.60 widened to 0.31 to 0.66; singly housed dams 0.2 +/- 0.1 per h (Extended Data Fig. 8a legend, N = 5)")
session mouse hours dam_retrieval ret_per_h
sub-V110_ses-03032018 V110 12.04 0 0.000
sub-V111_ses-03032018 V111 4.62 1 0.216
sub-V111_ses-03052018 V111 10.77 0 0.000
sub-V111_ses-03062018 V111 12.41 2 0.161
sub-V116_ses-03082018 V116 12.22 4 0.327
sub-V116_ses-03092018 V116 15.82 11 0.695
sub-V116_ses-03102018 V116 13.45 2 0.149
sub-V116_ses-03122018 V116 12.53 2 0.160
sub-V118_ses-03222018 V118 2.52 1 0.397
sub-V85_ses-12102017 V85 12.01 0 0.000
Mean 0.21 per h over 10 sessions (median 0.16, range 0.00 to 0.70); total 23 in 108.4 h = 0.21 per h
Per mouse: {'V110': 0.0, 'V111': 0.13, 'V116': 0.33, 'V118': 0.4, 'V85': 0.0}
Paper: 0.5 per h (0.45 to 0.55), or the daily range 0.34 to 0.60 widened to 0.31 to 0.66; singly housed dams 0.2 +/- 0.1 per h (Extended Data Fig. 8a legend, N = 5)
Result: not reproduced. Dams retrieve 23 times in 108.4 h, 0.21 per hour (session mean 0.21, range 0 to 0.70), below the 0.31 to 0.66 of the criterion and close to the 0.2 +/- 0.1 per hour the legend of Extended Data Fig. 8a gives for singly housed dams (N = 5). Two of 10 sessions fall inside the range, V116 03082018 (0.33) and V118 03222018 (0.40). The annotations span each session (data check), but the dandiset does not say whether every retrieval in a session was annotated.
Outputs for blocked claims¶
These outputs were asked for in the first verification round. They do not test the blocked claims; they show how far the dandiset goes towards them.
fig2f-virgin-chasing: how "Virgin Agressing" lies against "Mom agressing", and the change in firing of the shepherding-activated units during it (10 s baseline, as for shepherding).
va_rows = []
for s in eph:
va, vae = pvn.episodes(s["events"], ("Virgin Agressing",))
ma, mae = pvn.episodes(s["events"], pvn.SHEPHERDING)
if not len(va):
continue
overlap = sum(np.any((ma < vae[i]) & (mae > va[i])) for i in range(len(va)))
near = np.array([np.min(np.abs(ma - x)) for x in va]) if len(ma) else np.array([np.nan])
va_rows.append({"session": s["session"], "virgin_agressing": len(va), "mom_agressing": len(ma), "overlapping": int(overlap),
"nearest_mom_onset_s_min": round(np.nanmin(near), 1), "nearest_mom_onset_s_median": round(np.nanmedian(near), 1)})
print(pd.DataFrame(va_rows).to_string(index=False))
rng_f2f = np.random.default_rng(114010)
for r in shep[shep.activated].itertuples():
s = next(x for x in with_units if x["session"] == r.session)
st, en = pvn.episodes(s["events"], ("Virgin Agressing",))
if not len(st):
print(f"{r.session} unit {r.unit}: shepherding change {r.mi:.1f}; no 'Virgin Agressing' episode")
continue
real, null, p = pvn.shuffle_test(s["units"][r.unit], st, en, "pre10", pvn.draw_shifts(rng_f2f))
print(f"{r.session} unit {r.unit}: shepherding change {r.mi:.1f}; 'Virgin Agressing' change {real:.1f} over {len(st)} episodes, shuffle p = {p:.3f}")
session virgin_agressing mom_agressing overlapping nearest_mom_onset_s_min nearest_mom_onset_s_median sub-V110_ses-03032018 4 4 0 828.5 3763.1 sub-V111_ses-03032018 1 4 0 1169.3 1169.3 sub-V111_ses-03052018 11 2 0 162.3 5319.6 sub-V111_ses-03062018 15 1 0 1652.5 17872.1 sub-V116_ses-03092018 1 22 0 212.4 212.4
sub-V111_ses-03052018 unit 5: shepherding change 100.0; 'Virgin Agressing' change -25.3 over 11 episodes, shuffle p = 0.054 sub-V116_ses-03092018 unit 5: shepherding change 39.5; 'Virgin Agressing' change -100.0 over 1 episodes, shuffle p = 0.000 sub-V116_ses-03102018 unit 22: shepherding change 33.3; no 'Virgin Agressing' episode
"Virgin Agressing" never overlaps a "Mom agressing" episode, and its onsets lie at least 162 s from the nearest one (median per session 212 to 17,872 s), so it is a separate behaviour, which may be the virgin chasing the dam. Of the 3 shepherding-activated units, 2 have "Virgin Agressing" episodes: one fires less over 11 episodes (change -25.3, shuffle p = 0.054), the other stops firing in the single episode of its session. Two units do not test the paper's 6, and the label is still unnamed, so the claim stays blocked.
fig4i-projection-photometry: the change in the photometry trace in the 5 s after each dam retrieval against the 5 s before, per observation session, and the virgin's own retrievals in the retrieval-test files.
ph_rows = []
for s in oph:
tr = s["events"][s["events"].label == "retrieval"].start.to_numpy()
c = pvn.triggered_change(s["t"], s["f"], tr) if len(tr) else np.array([np.nan])
ph_rows.append({"session": s["session"], "retrievals": len(tr), "mean_change_pct": round(np.nanmean(c), 2),
"min": round(np.nanmin(c), 2), "max": round(np.nanmax(c), 2)})
ph_rows = pd.DataFrame(ph_rows)
print(ph_rows.to_string(index=False))
day1 = ph_rows[ph_rows.session.str.endswith("Day 1-obs")]
print("Day 1 observation, mean change per animal:", dict(zip(day1.session.str.split("_").str[0], day1.mean_change_pct)))
session retrievals mean_change_pct min max
sub-ROV40_ses-Day 1-obs 11 0.84 -1.53 2.86
sub-ROV40_ses-Day 1-retrieval 6 0.89 0.02 1.80
sub-ROV40_ses-Day 2-obs 11 0.57 -1.20 3.25
sub-ROV40_ses-Day 2-retrieval 6 -0.07 -0.96 0.68
sub-ROV40_ses-Day 3-obs 10 -0.18 -2.36 1.42
sub-ROV40_ses-Day 3-retrieval 1 -0.37 -0.37 -0.37
sub-ROV41_ses-Day 1-obs 10 1.11 -3.27 4.24
sub-ROV41_ses-Day 1-retrieval 0 NaN NaN NaN
sub-ROV41_ses-Day 2-obs 11 0.49 -2.02 2.94
sub-ROV41_ses-Day 2-retrieval 6 0.45 -1.31 2.44
sub-ROV41_ses-Day 3-obs 10 -0.19 -1.39 1.98
sub-ROV41_ses-Day 3-retrieval 0 NaN NaN NaN
sub-ROV43_ses-Day 1-obs 11 0.86 -1.84 2.44
sub-ROV44_ses-Day 1-obs 10 1.85 -7.47 7.99
sub-ROV45_ses-Day 1-obs 10 0.45 -6.06 10.72
sub-ROV49_ses-Day 1-obs 10 0.22 -1.19 1.38
Day 1 observation, mean change per animal: {'sub-ROV40': 0.84, 'sub-ROV41': 1.11, 'sub-ROV43': 0.86, 'sub-ROV44': 1.85, 'sub-ROV45': 0.45, 'sub-ROV49': 0.22}
Every animal's mean change after the dam's retrievals on Day 1 is positive (0.22 to 1.85%). ROV40 retrieves in its Day 1 test and ROV41 in its Day 2 test; the other four have no test file, so the learner and non-learner groups cannot be formed, and the Source Data values (-12.1 to 65.7) are on a scale this measure does not match. The claim stays blocked.
Summary¶
| Claim | Outcome |
|---|---|
| label-link | reproduced |
| ed7a-nest-entry-fraction | partially reproduced |
| ed7a-shepherding-fraction | not reproduced |
| ed7a-dam-retrieval-fraction | partially reproduced |
| fig2j-shepherded-entry | not reproduced |
| fig1j-shepherding-rate | partially reproduced |
| ed8a-dam-retrieval-rate | not reproduced |
| ot-pvn-units | blocked: data missing |
| fig2f-virgin-chasing | blocked: metadata insufficient |
| ed6-dam-vs-pup | blocked: data missing |
| fig1cd-retrieval-onset | blocked: data missing |
| fig2a-dreadd | blocked: data missing |
| fig3fg-observational-learning | blocked: data missing |
| fig3ij-observation-units | blocked: data missing |
| fig4c-ac-photometry | blocked: data missing |
| fig4f-pvn-ac-correlation | blocked: data missing |
| fig4i-projection-photometry | blocked: metadata insufficient |
| fig4jk-oxytocin-in-ac | blocked: data missing |
| ed9-colliculus | blocked: data missing |
A failure to reproduce here may be a failure of this attempt, not of the paper: see the reasoning and the differences from the paper's method.
How the analysis differs from the paper's.
- The dandiset has units from 4 of the paper's 17 virgins (194 units in 7 sessions with nest entries, against the paper's 147 units), and no OT-PVN identity, so OT-PVN and PVN units are pooled and the false discovery rate is applied over each session's units together.
- The label names come from the paper's code. The definition of a shepherded entry, the session selection for the rates, and the readings of each reported fraction were fixed in the claims file before the analysis.
- Where the Methods and the code differ (the 40 s cap, the 10 s baseline), the Methods are followed and the code's choice is shown; for nest entry the cap changes the fraction from 21.6% to 1.0%.
- Units of V111's three days are counted as separate units, as the paper counts units per day, though their spike counts suggest the same clusters; counting V111 once is shown as a variant.
- As in the code, real episodes without spikes are left out of a unit's mean while shuffled ones count as 0; counting both as 0 is shown as a variant (nest entry 16.0% against 21.6%).
- Poisson units with no response give 1.0 to 2.3% activated on the same episodes and test, the Fisher test against zero, as in the paper, does not take this level into account.
What was found about the dandiset.
- The behaviour labels have no description; "Mom agressing" (shepherding) and "Virgin in rest" (nest entry) can be read only through the paper's code. One file carries both naming schemes.
- The OT-PVN identity and the optical tagging recordings are not in the dandiset, nor are the behavioural cohorts, the observation-session unit recordings, or the auditory cortex data.
- Units in V110 and V111 have a median 7.0 to 9.2% of inter-spike intervals under 2 ms, and 79 to 90% of their spikes coincide within 0.1 ms with another unit's, against the Methods' removal of such spikes.
- Session start times disagree with the dates in the session ids for V116 and V118, and the spikes and events of V116 03102018 run to 13.4 h while its raw voltage lasts 4.46 h.
- No test here could show that spikes and annotations share a clock. Without V116 03102018, whose spikes and events outlast its raw voltage, nest entry gives 14.9% and the shepherded-entry comparison p = 0.0237.
- In V111 03032018 the 9 units activated by the dam's single retrieval share a median 71% of their spikes pairwise within 0.1 ms, and the retrieval falls on one of the session's largest synchronous bursts (14,346 population spikes in 5 s against 1,838 on average).
- Subject descriptions say every animal was co-housed, including the photometry animals of the
observation experiment; electrode and fibre locations and the fluorophore are "unknown"; the
photometry trace is labelled in watts; the
commanded_voltageseries is not a command signal. - Which photometry animals learned to retrieve is known for 2 of 6.
- The draft holds the same spike times and events, the same shape and rate of the voltage series, and the same photometry traces as the pinned version in all 30 files; the voltage values were not compared.
print(f"Total run time {time.time() - started:.0f} s")
Total run time 183 s