Reproducing Mu et al. (2019) from dandiset 000350¶
Paper. Mu Y, Bennett DV, Rubinov M, Narayan S, Yang C-T, Tanimoto M, Mensh BD, Looger LL, Ahrens MB (2019). Glia accumulate evidence that actions are futile and suppress unsuccessful behavior. Cell 178, 27-43. doi:10.1016/j.cell.2019.05.050. The paper is open access (CC BY-NC-ND). This attempt used its full text (main text and STAR Methods); the Supplemental Information was not read.
Dandiset. DANDI:000350, version
0.240822.1759. The paper cites no dandiset version (it predates the dandiset). Two versions are
published, 0.221219.1506 and 0.240822.1759. The later one is pinned: it is the most recent, and
it corrects the imaging plane of the three files with radial astrocytes only, which the earlier
version names NeuronVolume ("for the neuron filter"). The data the analysis uses (swim bouts,
cell traces) are the same in both; the data check shows this.
What the dandiset holds. 12 light-sheet sessions of paralysed larval zebrafish swimming fictively in virtual reality, each with the swim and stimulus signals of the rig, the detected swim bouts and a table of behavioural states. Eight files also hold segmented cell traces: 4 of neurons and radial astrocytes imaged together, 4 of neurons only. The other 4 hold raw images only. None of the perturbation experiments of Figures 4 to 7 is in it.
Code. The paper links its cell segmentation code, https://github.com/mikarubi/segmentation
(now voluseg), read at commit 54d1907b1d911ad6445428f8a79ca44875ed42b1 (27 June 2019, the last
before publication). It holds motion correction, brain masking, cell detection, baseline
estimation and the first NMF of cell traces (z6_components.py), with censor_tau = [30, 30]
and baseline_tau = 300 s. It holds no behavioural analysis, no passivity detection, no
clustering of components across fish and no onset analysis, so every choice below on those comes
from the Methods or is made here and stated.
| Claim | Tested here | Data |
|---|---|---|
| fig1h-epoch-durations: passive 24.4 s, active 36.1 s, 29 bouts per active epoch | tested | 12 fish |
| methods-isi: inter-swim interval 1.2 s (closed loop), 1.19 s (active epochs) | tested | 12 fish |
| fig1j-open-loop-passivity: more passivity in open than closed loop | tested | 12 fish |
| fig1c-vigor-open-loop-onset: higher vigor after the switch to open loop | tested | 12 fish |
| fig1i-isi-before-passivity: intervals and vigor rise before passivity | tested | 12 fish |
| fig2bc-neurons-lower-in-passivity | tested | 8 fish |
| fig2c-glia-rise-at-passivity | tested | 4 fish |
| figs2a-open-loop-higher | tested | 8 and 4 fish |
| fig2h-glia-weaker-closed-loop-pauses | tested | 4 fish |
| fig2dg-nmf-components | blocked: methods underspecified | |
| fig2kl-lmo-glia-most-consistent | blocked: metadata insufficient | |
| fig1k-passive-unresponsive, fig1f-replay | blocked: metadata insufficient | |
| fig4-glia-necessary, fig5-glia-sufficient | blocked: data missing | |
| fig6a-ne-ramp | blocked: metadata insufficient | |
| fig6-ne-drives-glia | blocked: data missing | |
| fig7-ne-mismatch, fig7-glia-accumulate | blocked: metadata insufficient |
This is the second version of the notebook. After the first verification round it adds the protocol-code check, a variant that keeps the epochs cut by the edges of a condition period, the closed-loop passivity by gain level before and after the first open loop, the SD within periods, the first onsets after each switch and a stretch of the raw traces in two fish, and it limits the bouts counted after a switch to open loop to that open-loop period (fig1c). After the second round it adds the passivity onsets of the two-colour fish by protocol code and the SD within 10 s windows (figs2a), and corrects text the outputs did not support.
In [1]:
import os
os.environ.setdefault("OMP_NUM_THREADS", "1")
import time
import warnings
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from scipy import stats
import fictive as fv
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()
sessions = fv.load_all()
S = {s["fish"]: s for s in sessions}
print(f"loaded {len(sessions)} files from dandiset {fv.DANDISET_ID} version {fv.DANDISET_VERSION} "
f"in {time.time() - started:.0f} s")
# Whole-brain signal: mean DfOverF over an evenly spaced subset of ROI blocks, for every processed
# series; and, for one fish, a disjoint subset to check that the subset does not matter.
t0 = time.time()
jobs = [(s, name, 0.0) for s in sessions for name in s["imaging"]]
jobs += [(S["20170113-4"], "NeuronDfOverF", 0.5), (S["20170228-4"], "GliaDfOverF", 0.5)]
WB = fv.whole_brain_means(jobs)
print(f"whole-brain signals of {len(WB)} series in {time.time() - t0:.0f} s")
loaded 12 files from dandiset 000350 version 0.240822.1759 in 70 s
whole-brain signals of 14 series in 432 s
Data check¶
Each check of the metadata against the values is stated under its output.
In [2]:
rows = []
for s in sessions:
rows.append({
"fish": s["fish"], "subject_id": s["subject"]["subject_id"], "age": s["subject"]["age"],
"species": s["subject"]["species"], "start": s["session_start"][:19],
"minutes": round(s["duration_s"] / 60, 1),
"planes": "; ".join(f"{n}: {p.get('indicator')} ex {p.get('excitation_lambda')} em {p.get('emission_lambda')}"
for n, p in s["planes"].items()),
"raw_series": "; ".join(f"{k} {v['shape'][0]} vols ({v['neurodata_type']})" for k, v in s["raw_series"].items()),
"dfof_series": "; ".join(f"{k} {v['shape']}" for k, v in s["imaging"].items()) or "none",
"bouts": len(s["swims"]), "state_rows": len(s["states"]), "burst_events": s["n_burst_events"]})
info = pd.DataFrame(rows)
print(info.to_string(index=False))
print("\nsession_description (same in all files):", {s["session_description"] for s in sessions})
print("Distinct session start times:", info.start.nunique(), "of", len(info))
kinds = info.dfof_series.map(lambda x: "neurons and glia" if "Glia" in x and "Neuron" in x
else "neurons" if "Neuron" in x else "glia" if "Glia" in x else "raw images only")
print("\nFiles by processed traces:", kinds.value_counts().to_dict())
print("Raw-only files and their imaging plane:", {f: p for f, p, k in zip(info.fish, info.planes, kinds) if k == "raw images only"})
fish subject_id age species start minutes planes raw_series dfof_series bouts state_rows burst_events
20161022-1 20161022_1 P6D Danio rerio 2016-10-22T15:10:03 58.3 GliaVolume: GCaMP6f ex 488.0 em 537.5 GliaOnePhotonSeries 8650 vols (TwoPhotonSeries) none 446 3 5729
20161109-2 20161109_2 P6D Danio rerio 2016-11-09T21:19:50 74.1 GliaVolume: GCaMP6f ex 561.0 em 660.0; NeuronVolume: GCaMP6f ex 488.0 em 537.5 GliaOnePhotonSeries 12500 vols (TwoPhotonSeries); NeuronOnePhotonSeries 12500 vols (TwoPhotonSeries) GliaDfOverF (12500, 100840); NeuronDfOverF (12500, 109918) 2550 80 14122
20161206-1 20161206_1 P6D Danio rerio 2016-12-06T12:35:01 62.3 GliaVolume: GCaMP6f ex 488.0 em 537.5 GliaOnePhotonSeries 9000 vols (TwoPhotonSeries) none 1295 3 6261
20161206-2 20161206_2 P6D Danio rerio 2016-12-06T16:11:11 54.7 GliaVolume: GCaMP6f ex 488.0 em 537.5 GliaOnePhotonSeries 8499 vols (TwoPhotonSeries) none 1437 3 9094
20170111-1 20170111_1 P6D Danio rerio 2017-01-11T12:33:55 64.5 NeuronVolume: GCaMP6f ex 488.0 em 537.5 NeuronOnePhotonSeries 10050 vols (TwoPhotonSeries) none 1014 108 7073
20170113-2 20170113_2 P7D Danio rerio 2017-01-13T12:49:07 80.9 NeuronVolume: GCaMP6f ex 488.0 em 537.5 NeuronOnePhotonSeries 7616 vols (TwoPhotonSeries) NeuronDfOverF (7616, 85452) 1126 196 22226
20170113-3 20170113_3 P7D Danio rerio 2017-01-13T15:08:34 60.2 NeuronVolume: GCaMP6f ex 488.0 em 537.5 NeuronOnePhotonSeries 9813 vols (TwoPhotonSeries) NeuronDfOverF (9813, 90705) 1177 135 9781
20170113-4 20170113_4 P7D Danio rerio 2017-01-13T17:12:41 76.1 NeuronVolume: GCaMP6f ex 488.0 em 537.5 NeuronOnePhotonSeries 3065 vols (TwoPhotonSeries) NeuronDfOverF (12050, 93344) 1453 136 14574
20170127-1 20170127_1 P7D Danio rerio 2017-01-27T11:16:52 66.8 NeuronVolume: GCaMP6f ex 488.0 em 537.5 NeuronOnePhotonSeries 10822 vols (TwoPhotonSeries) NeuronDfOverF (10822, 99377) 1656 68 10138
20170203-1 20170203_1 P8D Danio rerio 2017-02-03T12:20:38 55.7 GliaVolume: GCaMP6f ex 561.0 em 660.0; NeuronVolume: GCaMP6f ex 488.0 em 537.5 GliaOnePhotonSeries 8499 vols (TwoPhotonSeries); NeuronOnePhotonSeries 8499 vols (TwoPhotonSeries) GliaDfOverF (8500, 100159); NeuronDfOverF (8500, 111235) 1525 434 15837
20170228-3 20170228_3 P6D Danio rerio 2017-02-28T16:57:30 59.7 GliaVolume: GCaMP6f ex 561.0 em 660.0; NeuronVolume: GCaMP6f ex 488.0 em 537.5 GliaOnePhotonSeries 6067 vols (TwoPhotonSeries); NeuronOnePhotonSeries 6067 vols (TwoPhotonSeries) GliaDfOverF (6067, 126467); NeuronDfOverF (6067, 111378) 559 246 13773
20170228-4 20170228_4 P6D Danio rerio 2017-02-28T18:50:02 96.5 GliaVolume: GCaMP6f ex 561.0 em 660.0; NeuronVolume: GCaMP6f ex 488.0 em 537.5 GliaOnePhotonSeries 5776 vols (TwoPhotonSeries); NeuronOnePhotonSeries 5776 vols (TwoPhotonSeries) GliaDfOverF (7200, 122296); NeuronDfOverF (7200, 124078) 2793 602 37212
session_description (same in all files): {'See the experiment_description.', 'A dual-color optic channel recording of both neuron and glia populations.', 'A single-color optic channel recording of either a neuron or a glia population.'}
Distinct session start times: 12 of 12
Files by processed traces: {'raw images only': 4, 'neurons and glia': 4, 'neurons': 4}
Raw-only files and their imaging plane: {'20161022-1': 'GliaVolume: GCaMP6f ex 488.0 em 537.5', '20161206-1': 'GliaVolume: GCaMP6f ex 488.0 em 537.5', '20161206-2': 'GliaVolume: GCaMP6f ex 488.0 em 537.5', '20170111-1': 'NeuronVolume: GCaMP6f ex 488.0 em 537.5'}
The files hold 12 fish aged 6 to 8 days, as in the Methods (5 to 7 dpf for behaviour, one at P8D).
The Methods (Figure 2D and 2E) used 5 fish with neuronal GCaMP6f and 6 with glial GCaMP6f in
separate lines. Here the neuron-only GCaMP6f files are 5 (4 with traces, 1 with raw images only),
and the glia-only GCaMP6f files are 3, all with raw images only. The 4 two-colour files label both
planes GCaMP6f, but their glial plane is excited at 561 nm and detected at 660 nm, which is the
jRGECO1b channel of the Methods (Tg(gfap:jRGECO1b); Tg(elavl3:GCaMP6f), the line of Figure 2B
and 2C). The data disagree with the indicator stated for the glial plane of these files. The raw
series are stored as TwoPhotonSeries and the archive lists "two-photon microscopy technique",
while the Methods describe a one-photon light-sheet microscope.
In [3]:
older = fv.plane_summary("0.221219.1506")
newer = fv.plane_summary(fv.DANDISET_VERSION)
both = older.merge(newer, on="fish", suffixes=("_0.221219", "_0.240822"))
for col in ["planes", "raw_series", "swim_bouts", "sum_swim_starts", "dfof_series"]:
same = (both[f"{col}_0.221219"] == both[f"{col}_0.240822"])
print(f"{col}: same in {same.sum()} of {len(both)} files")
if not same.all():
print(both.loc[~same, ["fish", f"{col}_0.221219", f"{col}_0.240822"]].to_string(index=False))
planes: same in 9 of 12 files
fish planes_0.221219 planes_0.240822
20161022-1 NeuronVolume (GCaMP6f, 488.0 nm) GliaVolume (GCaMP6f, 488.0 nm)
20161206-1 NeuronVolume (GCaMP6f, 488.0 nm) GliaVolume (GCaMP6f, 488.0 nm)
20161206-2 NeuronVolume (GCaMP6f, 488.0 nm) GliaVolume (GCaMP6f, 488.0 nm)
raw_series: same in 0 of 12 files
fish raw_series_0.221219 raw_series_0.240822
20161022-1 TwoPhotonSeries GliaOnePhotonSeries
20161109-2 GliaTwoPhotonSeries, NeuronTwoPhotonSeries GliaOnePhotonSeries, NeuronOnePhotonSeries
20161206-1 TwoPhotonSeries GliaOnePhotonSeries
20161206-2 TwoPhotonSeries GliaOnePhotonSeries
20170111-1 TwoPhotonSeries NeuronOnePhotonSeries
20170113-2 TwoPhotonSeries NeuronOnePhotonSeries
20170113-3 TwoPhotonSeries NeuronOnePhotonSeries
20170113-4 TwoPhotonSeries NeuronOnePhotonSeries
20170127-1 TwoPhotonSeries NeuronOnePhotonSeries
20170203-1 GliaTwoPhotonSeries, NeuronTwoPhotonSeries GliaOnePhotonSeries, NeuronOnePhotonSeries
20170228-3 GliaTwoPhotonSeries, NeuronTwoPhotonSeries GliaOnePhotonSeries, NeuronOnePhotonSeries
20170228-4 GliaTwoPhotonSeries, NeuronTwoPhotonSeries GliaOnePhotonSeries, NeuronOnePhotonSeries
swim_bouts: same in 12 of 12 files
sum_swim_starts: same in 12 of 12 files
dfof_series: same in 12 of 12 files
Between the two published versions, only the names change: the raw series (TwoPhotonSeries
to NeuronOnePhotonSeries or GliaOnePhotonSeries) and, in the three glia-only files, the imaging
plane (NeuronVolume to GliaVolume). Swim bouts and processed traces are the same.
In [4]:
# Duplicate recordings: compare bout times between every pair of files.
dup = []
fish_list = list(S)
for i, a in enumerate(fish_list):
for b in fish_list[i + 1:]:
x, y = S[a]["swims"].start_time.to_numpy(), S[b]["swims"].start_time.to_numpy()
n = min(len(x), len(y), 200)
dup.append({"a": a, "b": b, "same_first_200_starts": bool(np.allclose(x[:n], y[:n])),
"same_n_samples": S[a]["n_samples"] == S[b]["n_samples"]})
dup = pd.DataFrame(dup)
print("Pairs of files with the same first bout times:", int(dup.same_first_200_starts.sum()), "of", len(dup))
print("Pairs with the same number of rig samples:", int(dup.same_n_samples.sum()))
# Units and ranges of the behavioural data.
rows = []
for s in sessions:
sw = s["swims"]
rows.append({"fish": s["fish"], "rate_hz": s["rate"], "bouts": len(sw),
"dur_med_s": sw.eval("stop_time - start_time").median(), "width_minus_dur_max": (sw.width - sw.eval("stop_time - start_time")).abs().max(),
"power_p5": sw.power.quantile(0.05), "power_med": sw.power.median(), "power_p95": sw.power.quantile(0.95),
"power_zero": int((sw.power == 0).sum()),
"gain_values": np.unique(np.round(s["gain"], 3)).size, "gain_max": round(s["gain"].max(), 3),
"velocity_values": sorted(set(np.round(s["velocity"], 2))),
"protocol_codes": sorted(set(s["protocol_code"].astype(int)))})
beh = pd.DataFrame(rows)
print("\n", beh.round(3).to_string(index=False))
print("\nStated units:", sessions[0]["behavior_units"])
Pairs of files with the same first bout times: 0 of 66
Pairs with the same number of rig samples: 0
fish rate_hz bouts dur_med_s width_minus_dur_max power_p5 power_med power_p95 power_zero gain_values gain_max velocity_values protocol_codes
20161022-1 5989.6 446 0.366 0.008 0.089 1.514 3.528 0 5 0.036 [0.0, 0.5, 1.0] [22, 23, 32, 33]
20161109-2 5989.6 2550 0.204 0.004 0.118 0.344 0.850 0 4 0.030 [0.0, 0.5, 1.0, 2.0] [22, 23, 32, 33]
20161206-1 5989.6 1295 0.190 0.002 0.059 0.219 0.476 0 3 0.018 [0.0, 0.5, 1.0] [22, 23, 32, 33]
20161206-2 5989.6 1437 0.181 0.005 0.119 0.583 2.357 0 2 0.020 [0.0, 0.5, 1.0, 1.5] [22, 23, 32, 33]
20170111-1 5989.6 1014 0.248 0.003 0.125 0.338 0.691 0 5 0.036 [0.5, 1.0, 1.5] [22, 23, 32, 33, 51, 53]
20170113-2 5989.6 1126 0.214 0.016 0.021 0.121 3.743 0 4 0.030 [0.5, 1.0, 1.5] [22, 23, 32, 33, 51, 53]
20170113-3 5989.6 1177 0.255 0.006 0.022 0.333 1.251 0 4 0.036 [0.5, 1.0, 1.5] [22, 23, 32, 33, 51, 53]
20170113-4 5989.6 1453 0.233 0.022 0.041 0.300 1.430 0 18 0.030 [0.0, 0.5, 1.0, 1.5] [22, 23, 32, 33, 51, 52, 53]
20170127-1 5989.6 1656 0.156 0.057 0.039 0.180 0.344 0 23 0.030 [0.0, 0.5, 1.0, 1.5] [22, 23, 32, 33, 51, 52, 53]
20170203-1 5989.6 1525 0.185 0.002 0.006 0.146 0.320 70 2 0.015 [0.0, 0.5, 1.0] [22, 23, 32, 33, 51, 53]
20170228-3 5989.6 559 0.503 0.008 0.046 0.808 7.923 0 3 0.020 [0.0, 0.5] [12, 22, 23, 51, 53]
20170228-4 5989.6 2793 0.277 0.008 0.000 0.206 0.768 279 2 0.020 [0.0, 1.0] [22, 23, 51, 53]
Stated units: {'SwimSignals': 'a.u.', 'StimulusGainSignal': 'a.u.'}
No two files share bout times or length, so each recording is counted once. Swim bout width
equals stop minus start (in seconds), though the table describes it as "Estimated width spanned by
the swim event". Swim power and the stimulus signals are in arbitrary units, as stated. In two
files (20170203-1 and 20170228-4) some bouts have a power of exactly 0. The gain signal takes
values from 0 to 0.036 (arbitrary units; the Methods give none for the gain), and the grating's
baseline velocity takes 0, 0.5, 1, 1.5 or 2 (arbitrary units, against 2.86 mm/s in the Methods).
The first column of StimulusParameterSignals holds protocol codes (12, 22, 23, 32, 33, 51, 52,
53); no file describes them, and the analysis uses the gain and velocity signals instead. What the
codes hold is checked below.
In [5]:
rows = []
for s in sessions:
m = fv.condition_masks(s)
dt = 1 / fv.COND_HZ
rows.append({"fish": s["fish"], "closed_loop_s": m["closed"].sum() * dt, "open_loop_s": m["open"].sum() * dt,
"stationary_gain_on_s": (m["gain_on"] & ~m["closed"]).sum() * dt,
"stationary_gain_0_s": (~m["gain_on"] & ~m["open"]).sum() * dt,
"closed_periods": len(fv.segments(s["cond_t"], m["closed"])),
"open_periods": len(fv.segments(s["cond_t"], m["open"])),
"first_open_s": fv.segments(s["cond_t"], m["open"])[0][0]})
cond = pd.DataFrame(rows)
print("Time per condition (closed loop: gain above 0, grating moving; open loop: gain 0, grating moving):")
print(cond.round(1).to_string(index=False))
Time per condition (closed loop: gain above 0, grating moving; open loop: gain 0, grating moving):
fish closed_loop_s open_loop_s stationary_gain_on_s stationary_gain_0_s closed_periods open_periods first_open_s
20161022-1 530.5 2779.4 190.8 0.0 8 2 360.7
20161109-2 666.8 3352.4 144.7 282.4 10 29 565.6
20161206-1 828.4 2296.1 612.3 0.0 50 2 601.1
20161206-2 753.4 1983.4 544.9 0.0 48 2 601.1
20170111-1 1360.0 2511.1 0.0 0.0 20 19 180.4
20170113-2 864.8 3990.7 0.0 0.0 14 14 300.6
20170113-3 1040.5 2571.3 0.0 0.0 19 19 300.6
20170113-4 1944.5 2619.7 0.0 0.1 30 32 193.1
20170127-1 1775.8 1501.8 727.7 0.1 68 3 1217.6
20170203-1 1006.9 2121.8 213.6 0.0 17 17 601.1
20170228-3 192.3 3306.4 82.2 0.0 15 2 250.0
20170228-4 291.7 5173.1 324.8 0.0 13 3 240.5
Every fish has closed-loop time (gain above 0) and open-loop time (gain 0), with the open loop first starting 3 to 20 minutes into the session. Several sessions also alternate closed and open loop later on, and several closed-loop blocks alternate with stationary grating.
In [6]:
# Protocol codes: what each undocumented code of StimulusParameterSignals (first column) holds.
rows = []
for s in sessions:
c, t = s["protocol_code"].astype(int), s["cond_t"]
edges = np.r_[0, np.flatnonzero(np.diff(c)) + 1, len(c)]
for code_ in np.unique(c):
m = c == code_
blocks = [t[min(j, len(t) - 1)] - t[i] for i, j in zip(edges[:-1], edges[1:]) if c[i] == code_]
g = np.round(s["gain"][m], 3)
rows.append({"fish": s["fish"], "code": code_, "seconds": round(m.sum() / fv.COND_HZ), "blocks": len(blocks),
"median_block_s": round(float(np.median(blocks)), 1), "first_s": round(t[m][0]), "last_s": round(t[m][-1]),
"gain_levels": len(np.unique(g)), "gain_min": g.min(), "gain_max": g.max(),
"velocities": sorted(set(np.round(s["velocity"][m], 2)))})
codes = pd.DataFrame(rows)
print(codes.to_string(index=False))
# Code 33: the baseline velocity just before, during and just after each block.
rows = []
for s in sessions:
c, v = s["protocol_code"].astype(int), s["velocity"]
starts = np.flatnonzero((c[1:] == 33) & (c[:-1] != 33)) + 1
for i in starts:
j = i + np.argmax(c[i:] != 33)
rows.append({"fish": s["fish"], "before": v[i - 2], "during": v[i], "after": v[j + 1], "dur_s": (j - i) / fv.COND_HZ,
"gain_during": s["gain"][i]})
c33 = pd.DataFrame(rows)
print(f"\nCode 33 blocks: {len(c33)} in {c33.fish.nunique()} files; duration {c33.dur_s.min():.1f} to {c33.dur_s.max():.1f} s; "
f"gain during {c33.gain_during.max():.3f} at most; velocity during above before in {(c33.during > c33.before).sum()}, "
f"equal in {(c33.during == c33.before).sum()}, below in {(c33.during < c33.before).sum()}, zero during in {(c33.during == 0).sum()}")
# TotalStimulusVelocitySignal, described as "(drift speed - motosensory gain) x swim power".
rows = []
for s in sessions:
tv = s["total_velocity"]
rows.append({"fish": s["fish"], "total_velocity_min": tv.min(), "total_velocity_max": tv.max(),
"baseline_velocity_max": s["velocity"].max(),
"corr_with_baseline_velocity": np.corrcoef(tv, s["velocity"][:len(tv)])[0, 1]})
print("\n", pd.DataFrame(rows).round(4).to_string(index=False))
# One session of each kind of protocol: gain, baseline velocity and code over time.
show = ["20161022-1", "20170113-4", "20170228-3", "20170228-4"]
fig, axes = plt.subplots(len(show), 1, figsize=(11, 2.2 * len(show)), sharex=False)
for ax, f in zip(axes, show):
s = S[f]
ax.plot(s["cond_t"], s["gain"] / max(s["gain"].max(), 1e-9), lw=0.6, label="gain (scaled to its max)")
ax.plot(s["cond_t"], s["velocity"] / s["velocity"].max(), lw=0.6, label="baseline velocity (scaled)")
ax2 = ax.twinx(); ax2.plot(s["cond_t"], s["protocol_code"], color="k", lw=0.6); ax2.set_ylabel("code")
ax.set_title(f"{f}: gain max {s['gain'].max():.3f}, velocity max {s['velocity'].max():.1f}, codes {sorted(set(s['protocol_code'].astype(int)))}", fontsize=8)
ax.set_ylim(-0.05, 1.05)
axes[0].legend(fontsize=7, loc="lower right"); axes[-1].set_xlabel("time (s)")
plt.tight_layout(); plt.show()
# Switches from closed to open loop in the two-colour files, with the open period's length and bouts.
rows = []
for s in sessions:
if "GliaDfOverF" not in s["imaging"]:
continue
m, t = fv.condition_masks(s), s["cond_t"]
segs = fv.segments(t, m["open"], 0)
for t0 in t[1:][m["closed"][:-1] & m["open"][1:]]:
end = min((b for a, b in segs if abs(a - t0) < 0.2), default=np.nan)
rows.append({"fish": s["fish"], "switch_s": t0, "open_s": end - t0,
"bouts": int(((s["swims"].start_time >= t0) & (s["swims"].start_time < end)).sum())})
sw2 = pd.DataFrame(rows)
print("\nClosed-to-open switches in the two-colour files:")
print(sw2.groupby("fish").agg(switches=("switch_s", "size"), open_s_min=("open_s", "min"), open_s_max=("open_s", "max"),
bouts_min=("bouts", "min"), bouts_max=("bouts", "max"),
switches_with_16_bouts=("bouts", lambda b: int((b >= 16).sum()))).round(1).to_string())
fish code seconds blocks median_block_s first_s last_s gain_levels gain_min gain_max velocities
20161022-1 22 360 15 24.3 0 2500 3 0.012 0.018 [0.0, 1.0]
20161022-1 23 361 15 24.3 24 2524 3 0.015 0.036 [0.0, 1.0]
20161022-1 32 2752 71 40.0 361 3501 1 0.000 0.000 [0.5, 1.0]
20161022-1 33 28 69 0.4 364 3480 1 0.000 0.000 [0.5, 1.0]
20161109-2 22 698 56 10.5 0 2943 3 0.000 0.020 [0.0, 0.5]
20161109-2 23 689 57 10.5 0 2946 3 0.000 0.030 [0.0, 0.5]
20161109-2 32 3031 86 34.6 845 4446 1 0.000 0.000 [0.5, 1.0, 2.0]
20161109-2 33 29 84 0.3 879 4342 1 0.000 0.000 [0.5, 1.0, 2.0]
20161206-1 22 828 51 16.6 0 3736 2 0.015 0.018 [0.0, 1.0]
20161206-1 23 612 51 12.5 21 3737 2 0.015 0.018 [0.0, 1.0]
20161206-1 32 2273 58 41.1 601 3498 1 0.000 0.000 [0.5, 1.0]
20161206-1 33 23 56 0.4 622 3487 1 0.000 0.000 [0.5, 1.0]
20161206-2 22 753 50 15.4 0 3282 1 0.020 0.020 [0.0, 1.0]
20161206-2 23 545 49 11.6 16 3262 1 0.020 0.020 [0.0, 0.5, 1.0]
20161206-2 32 1963 55 38.1 0 3126 1 0.000 0.000 [0.5, 1.0, 1.5]
20161206-2 33 20 52 0.4 616 3119 1 0.000 0.000 [0.5, 1.0, 1.5]
20170111-1 22 92 6 15.4 0 169 1 0.012 0.012 [0.5]
20170111-1 23 88 7 15.4 0 180 1 0.036 0.036 [0.5]
20170111-1 32 1011 28 38.1 180 1202 1 0.000 0.000 [1.0, 1.5]
20170111-1 33 10 27 0.4 192 1194 1 0.000 0.000 [1.0, 1.5]
20170111-1 51 1180 19 69.3 1202 3871 3 0.012 0.015 [1.0]
20170111-1 53 1489 18 84.7 1225 3851 1 0.000 0.000 [1.0]
20170113-2 22 286 20 14.7 0 3561 1 0.010 0.010 [0.5]
20170113-2 23 315 21 14.6 0 3606 1 0.030 0.030 [0.5]
20170113-2 32 2392 35 36.3 301 4808 1 0.000 0.000 [1.0, 1.5]
20170113-2 33 12 33 0.4 329 1502 1 0.000 0.000 [1.0, 1.5]
20170113-2 51 264 12 22.0 1612 3246 1 0.020 0.020 [1.0]
20170113-2 53 1587 14 124.6 1503 4856 1 0.000 0.000 [1.0]
20170113-3 22 212 15 14.6 0 3239 1 0.012 0.012 [0.5]
20170113-3 23 208 16 14.7 0 3246 1 0.036 0.036 [0.5]
20170113-3 32 1374 40 36.3 301 3612 1 0.000 0.000 [1.0, 1.5]
20170113-3 33 14 38 0.4 329 3591 1 0.000 0.000 [1.0, 1.5]
20170113-3 51 620 17 36.6 1322 3114 1 0.020 0.020 [1.0]
20170113-3 53 1183 17 73.3 1356 3126 1 0.000 0.000 [1.0]
20170113-4 22 191 13 14.7 0 4382 1 0.010 0.010 [0.5]
20170113-4 23 183 14 14.6 15 4388 1 0.030 0.030 [0.5]
20170113-4 32 1056 27 36.3 193 4565 1 0.000 0.000 [0.0, 0.5, 1.0, 1.5]
20170113-4 33 9 25 0.4 229 4397 1 0.000 0.000 [0.5, 1.0, 1.5]
20170113-4 51 1057 29 36.6 1099 4207 3 0.010 0.020 [0.5]
20170113-4 52 514 15 36.6 2600 4140 16 0.000 0.015 [0.5]
20170113-4 53 1554 30 36.7 1082 4177 1 0.000 0.000 [0.5]
20170127-1 22 1194 67 18.5 0 4006 2 0.010 0.012 [0.0, 1.0]
20170127-1 23 728 68 11.1 0 3988 2 0.010 0.012 [0.0, 1.0]
20170127-1 32 1468 42 36.6 1218 3772 1 0.000 0.000 [0.0, 0.5, 1.0, 1.5]
20170127-1 33 15 40 0.4 1254 3767 1 0.000 0.000 [1.0, 1.5]
20170127-1 51 309 6 59.1 2437 3005 3 0.010 0.030 [1.0]
20170127-1 52 272 6 51.7 2404 2972 22 0.009 0.030 [1.0]
20170127-1 53 19 1 19.4 2972 2991 1 0.000 0.000 [1.0]
20170203-1 22 388 76 3.9 0 601 1 0.015 0.015 [0.0, 1.0]
20170203-1 23 213 76 2.0 0 597 1 0.015 0.015 [0.0, 1.0]
20170203-1 32 1488 39 38.9 601 2104 1 0.000 0.000 [0.5, 1.0]
20170203-1 33 15 38 0.4 640 2095 1 0.000 0.000 [0.5, 1.0]
20170203-1 51 619 17 39.3 2104 3341 1 0.015 0.015 [0.5]
20170203-1 53 619 17 39.3 2123 3342 1 0.000 0.000 [0.5]
20170228-3 12 3 1 3.4 247 251 2 0.000 0.020 [0.5]
20170228-3 22 165 15 11.7 0 247 1 0.020 0.020 [0.0, 0.5]
20170228-3 23 82 15 5.9 0 247 1 0.020 0.020 [0.0, 0.5]
20170228-3 51 1740 147 11.7 251 3581 2 0.000 0.010 [0.5]
20170228-3 53 1590 146 11.7 262 3559 1 0.000 0.000 [0.5]
20170228-4 22 413 34 12.7 0 4141 1 0.020 0.020 [0.0, 1.0]
20170228-4 23 204 33 6.3 13 4147 1 0.020 0.020 [0.0, 1.0]
20170228-4 51 2995 142 12.7 241 5790 1 0.000 0.000 [1.0]
20170228-4 53 2178 142 12.7 241 4566 1 0.000 0.000 [1.0]
Code 33 blocks: 462 in 10 files; duration 0.3 to 0.5 s; gain during 0.000 at most; velocity during above before in 390, equal in 72, below in 0, zero during in 0
fish total_velocity_min total_velocity_max baseline_velocity_max corr_with_baseline_velocity
20161022-1 -0.0008 0.0075 1.0 -0.3237
20161109-2 -0.0017 0.0040 2.0 -0.8899
20161206-1 -0.0008 0.0040 1.0 -0.4220
20161206-2 -0.0012 0.0079 1.5 -0.3418
20170111-1 -0.0012 0.0062 1.5 -0.4908
20170113-2 -0.0012 0.0079 1.5 -0.7190
20170113-3 -0.0012 0.0051 1.5 -0.7039
20170113-4 -0.0012 0.0079 1.5 -0.1377
20170127-1 -0.0012 0.0034 1.5 -0.6480
20170203-1 -0.0008 0.0030 1.0 0.0003
20170228-3 -0.0008 0.0063 0.5 -0.1432
20170228-4 -0.0008 0.0031 1.0 -0.2031
Closed-to-open switches in the two-colour files:
switches open_s_min open_s_max bouts_min bouts_max switches_with_16_bouts
fish
20161109-2 3 0.1 1.6 0 1 0
20170203-1 18 1.1 1502.6 1 616 7
20170228-3 4 5.9 2871.4 1 337 2
20170228-4 3 1642.6 1803.1 570 1218 3
The protocol codes are not described in any file; what they hold, from the signals:
- 22 and 23 alternate in blocks with median lengths of 2 to 24 s, with the gain above 0 and the grating moving or stationary (closed loop, with stationary periods). In 20161109-2 these blocks also run with the gain at 0, so its short open-loop periods follow stationary grating, not closed loop.
- 32 is prolonged open loop, in blocks of 35 to 41 s (medians) separated by code 33. Code 33 lasts 0.3 to 0.5 s with the gain at 0; the grating moves faster during it than before it in 390 of the 462 blocks and at the same speed in 72, and it is never stationary. It is a brief speed-up of the grating, not the static then forward stimulus of Figure 1K, and nothing says what it is for.
- 51 and 53 alternate in blocks with median lengths of 12 to 125 s. In six files 51 has the gain above 0 and 53 has it at 0 (alternating closed and open loop), but in 20170228-4 both have the gain at 0, and in 20170228-3 code 51 has the gain at 0 for part of the time, so the two codes cannot be read as closed and open loop there.
- 52 appears in 20170113-4 and 20170127-1, with 16 and 22 gain levels (from 0 and from 0.009), consistent with a randomised gain such as that of Figure 7C.
- 12 appears once, for 3.4 s in 20170228-3.
TotalStimulusVelocitySignal is described as the drift speed minus the gain, times the swim power,
but its values lie between -0.002 and 0.008 in every file, against baseline velocities of 0.5 to
2, so it is not that velocity. The visual feedback each bout produced, and so a replayed stimulus
or a per-bout feedback probability, cannot be recovered from the files; only the baseline
velocity and the gain are usable. Among the two-colour files, only 20170203-1 has more than 4
switches from closed to open loop (18; 7 of them followed by 16 or more bouts in the open-loop
period). In 20170228-3 and 20170228-4, 2 of 4 and 3 of 3 switches are followed by 16 or more bouts,
and the 3 open-loop periods of 20170228-4 hold 570 to 1218 bouts each.
In [7]:
rows = []
for s in sessions:
st = s["states"]
dur = st.stop - st.start
act = st[st.type == "active"]
trans = st[st.type == "transient"]
same_as_active = trans.merge(act, on=["start", "stop"]).shape[0]
rows.append({"fish": s["fish"], "rows": len(st), "types": dict(st.type.value_counts()),
"first_s": round(st.start.min(), 1), "last_s": round(st.stop.max(), 1),
"first_open_s": round(cond.set_index("fish").first_open_s[s["fish"]], 1),
"active_under_5s": int((act.stop - act.start <= 5).sum()), "active_rows": len(act),
"transient_rows_also_active": f"{same_as_active} of {len(trans)}",
"passive_under_5s": int(((st.type == "passive") & (dur < 5)).sum())})
states = pd.DataFrame(rows)
print(states.to_string(index=False))
fish rows types first_s last_s first_open_s active_under_5s active_rows transient_rows_also_active passive_under_5s
20161022-1 3 {'passive': 1, 'active': 1, 'transient': 1} 360.7 492.7 360.7 0 1 0 of 1 0
20161109-2 80 {'active': 37, 'passive': 35, 'transient': 8} 844.7 4446.6 565.6 7 37 7 of 8 2
20161206-1 3 {'passive': 1, 'active': 1, 'transient': 1} 601.1 1220.0 601.1 0 1 0 of 1 1
20161206-2 3 {'active': 1, 'passive': 1, 'transient': 1} 601.1 618.1 601.1 0 1 0 of 1 0
20170111-1 108 {'active': 54, 'passive': 54} 182.8 1208.6 180.4 6 54 0 of 0 0
20170113-2 196 {'active': 97, 'passive': 95, 'transient': 4} 601.0 3305.8 300.6 16 97 4 of 4 0
20170113-3 135 {'active': 60, 'passive': 57, 'transient': 18} 304.7 3596.5 300.6 18 60 18 of 18 0
20170113-4 136 {'active': 65, 'passive': 56, 'transient': 15} 601.0 4416.1 193.1 16 65 15 of 15 0
20170127-1 68 {'active': 32, 'passive': 32, 'transient': 4} 1217.6 3799.5 1217.6 4 32 4 of 4 1
20170203-1 434 {'active': 191, 'passive': 186, 'transient': 57} 601.1 3342.4 601.1 95 191 57 of 57 11
20170228-3 246 {'active': 131, 'passive': 67, 'transient': 48} 261.9 3581.1 250.0 67 131 48 of 48 0
20170228-4 602 {'active': 292, 'passive': 248, 'transient': 62} 256.5 5789.9 240.5 66 292 62 of 62 30
The ActivityStates table is described only as "Classified periods of activity (passive, active,
or transient)". In the three glia-only files it holds 3 rows; in the others it starts at or after
the first open loop, and in sessions that later alternate closed and open loop it runs through both
without saying which. Its "transient" rows are not defined in the paper, and in the files with more
than 3 rows all of them, or all but one, are also listed with the same times as active rows. Many active rows
last 5 s or less, which the Methods' definition of an active epoch excludes. The table is therefore
not the Methods' classification as written; the claims use epochs built from the swim bouts by
the Methods' rule, and the table as a variant.
In [8]:
rows = []
for s in sessions:
for name, im in s["imaging"].items():
t = im["t"]
ttl = s["ttl_volume_starts"]
raw_name = name.replace("DfOverF", "OnePhotonSeries")
raw_t = s["raw_series"][raw_name]["t"]
y = WB[(s["fish"], name, 0.0)][0]
cens = fv.censored(y)
rows.append({"fish": s["fish"], "series": name, "volumes": len(t), "rois": im["shape"][1],
"rois_used": WB[(s["fish"], name, 0.0)][1], "dtype_unit": im["unit"],
"vol_hz": round(1 / np.median(np.diff(t)), 2), "ttl_starts": len(ttl),
"t_matches_ttl": bool(np.allclose(t, ttl[:len(t)], atol=1e-3)),
"raw_vols": len(raw_t), "raw_first_s": round(raw_t[0], 1), "dfof_last_s": round(t[-1], 1),
"session_s": round(s["duration_s"], 1),
"trace_mean": float(np.mean(y[~cens])), "censored_start_s": round(t[np.argmin(cens)] - t[0], 1) if cens[0] else 0.0,
"censored_end_s": round(t[-1] - t[len(cens) - 1 - np.argmin(cens[::-1])], 1) if cens[-1] else 0.0})
img = pd.DataFrame(rows)
print(img.to_string(index=False))
samples = {k: fv.dfof_sample(S[f]["asset_id"], k) for f, k in [("20170113-4", "NeuronDfOverF"), ("20170228-4", "GliaDfOverF")]}
for k, d in samples.items():
x, F = d["dfof"], d["fluorescence"]
print(f"\n{k} (first 200 ROIs, stored as {d['dtype']}, unit {d['unit']}): DfOverF median {np.median(x):.4f}, "
f"range {x.min():.3f} to {x.max():.3f}, per-ROI mean {x.mean(0).min():.4f} to {x.mean(0).max():.4f}")
print(f" Fluorescence median {np.median(F):.4f}, range {F.min():.3f} to {F.max():.3f}; "
f"median absolute change between frames {np.median(np.abs(np.diff(F, axis=0))):.2e} "
f"(DfOverF {np.median(np.abs(np.diff(x, axis=0))):.2e})")
fish series volumes rois rois_used dtype_unit vol_hz ttl_starts t_matches_ttl raw_vols raw_first_s dfof_last_s session_s trace_mean censored_start_s censored_end_s 20161109-2 GliaDfOverF 12500 100840 11264 a.u. 2.86 12501 True 12500 0.0 4363.2 4446.6 -0.000026 60.0 0.0 20161109-2 NeuronDfOverF 12500 109918 11264 a.u. 2.86 12501 True 12500 0.0 4363.2 4446.6 -0.000134 60.0 0.0 20170113-2 NeuronDfOverF 7616 85452 9328 a.u. 2.73 9751 True 7616 0.0 2790.1 4855.8 1.000000 0.0 0.0 20170113-3 NeuronDfOverF 9813 90705 8424 a.u. 2.73 9814 True 9813 0.0 3595.0 3612.0 1.000000 0.0 0.0 20170113-4 NeuronDfOverF 12050 93344 7448 a.u. 2.73 12051 True 3065 3292.0 4414.7 4564.6 1.000000 0.0 0.0 20170127-1 NeuronDfOverF 10822 99377 8424 a.u. 2.71 10823 True 10822 0.0 3996.1 4005.6 1.000000 0.0 0.0 20170203-1 GliaDfOverF 8500 100159 13200 a.u. 2.54 8501 True 8499 0.0 3340.9 3342.4 -0.000040 60.1 0.0 20170203-1 NeuronDfOverF 8500 111235 14352 a.u. 2.54 8501 True 8499 0.0 3340.9 3342.4 -0.000011 60.1 0.0 20170228-3 GliaDfOverF 6067 126467 17760 a.u. 1.70 6068 True 6067 0.0 3561.6 3581.1 -0.000116 60.5 0.0 20170228-3 NeuronDfOverF 6067 111378 16848 a.u. 1.70 6068 True 6067 0.0 3561.6 3581.1 -0.000048 60.5 0.0 20170228-4 GliaDfOverF 7200 122296 16000 a.u. 1.58 7201 True 5776 0.0 4565.3 5789.9 -0.000048 60.2 0.0 20170228-4 NeuronDfOverF 7200 124078 16456 a.u. 1.58 7201 True 5776 0.0 4565.3 5789.9 -0.000010 60.2 0.0
NeuronDfOverF (first 200 ROIs, stored as float64, unit a.u.): DfOverF median 0.9930, range 0.728 to 2.833, per-ROI mean 1.0000 to 1.0000 Fluorescence median 0.9176, range 0.669 to 1.008; median absolute change between frames 2.44e-06 (DfOverF 1.02e-02) GliaDfOverF (first 200 ROIs, stored as float32, unit a.u.): DfOverF median -0.0009, range -0.230 to 1.428, per-ROI mean -0.0005 to 0.0009 Fluorescence median -0.0100, range -0.167 to 0.001; median absolute change between frames 0.00e+00 (DfOverF 4.05e-03)
- The imaging timestamps equal the rising edges of the camera trigger (
FrameCaptureTTLSignalabove 3.8), so imaging and behaviour share the rig's clock. In 20170113-2 the trigger has 9751 volume starts but the imaging series only 7616: the last 2135 volumes are not in the file. - The volume rate is 1.6 to 2.9 Hz, against "3 volumes/second" in the Methods; the two-colour sessions are the slowest.
- In 20170113-4 the raw image series holds 3065 volumes, starting at 3292 s, while the traces cover 12050 volumes from the start; in 20170228-4 the raw series holds 5776 volumes against 7200 traces. Raw images are missing for part of these sessions.
- The series called
DfOverFis not a ΔF/F as the Methods define it (baseline-subtracted and divided by the baseline, so near 0 at rest). In the single-colour files it is each ROI's trace divided by its mean (every ROI mean is 1, float64); in the two-colour files it is centred near 0 (float32) with its first 60 s held constant, as the code's censoring (censor_tau) does. The two groups were processed differently. TheFluorescenceseries changes by about 2e-6 per frame in the single-colour file and not at all from frame to frame in the median ROI of the two-colour file, which is consistent with a slowly varying baseline, not raw fluorescence; its description is the same as that ofDfOverF. The analysis therefore compares the whole-brain signal only within a fish, against a baseline window or between conditions, and drops the censored frames.
In [9]:
# Does the ROI subset matter? Correlate the whole-brain signal from two disjoint subsets.
for fish, name in [("20170113-4", "NeuronDfOverF"), ("20170228-4", "GliaDfOverF")]:
a, na = WB[(fish, name, 0.0)]
b, nb = WB[(fish, name, 0.5)]
ok = ~fv.censored(a)
print(f"{fish} {name}: subsets of {na} and {nb} ROIs, correlation {np.corrcoef(a[ok], b[ok])[0, 1]:.3f}")
# Clock check: whole-brain neuronal signal against the swim rate, with the true bout times and
# with the bouts shifted by 60 to 600 s (circularly), for every file with neuronal traces.
rng = np.random.default_rng(0)
rows = []
for s in sessions:
if "NeuronDfOverF" not in s["imaging"]:
continue
t = s["imaging"]["NeuronDfOverF"]["t"]
y = WB[(s["fish"], "NeuronDfOverF", 0.0)][0]
ok = ~fv.censored(y)
t, y = t[ok], y[ok]
rate = fv.swim_rate(s, t)
dt = np.median(np.diff(t))
lags = np.arange(-int(10 / dt), int(10 / dt) + 1)
def best(r):
cc = [np.corrcoef(r[max(0, -l):len(r) - max(0, l)], y[max(0, l):len(y) - max(0, -l)])[0, 1] for l in lags]
return np.max(cc), lags[np.argmax(cc)] * dt
true_cc, true_lag = best(rate)
null = [best(np.roll(rate, int(sh / dt)))[0] for sh in rng.uniform(60, 600, 20) * rng.choice([-1, 1], 20)]
rows.append({"fish": s["fish"], "peak_corr": true_cc, "peak_lag_s": true_lag,
"shifted_max": np.max(null), "shifted_median": np.median(null)})
clock = pd.DataFrame(rows)
print("\nNeuronal whole-brain signal against swim rate (positive lag: signal after swims):")
print(clock.round(3).to_string(index=False))
print("Files where the true peak correlation exceeds every shifted one:", int((clock.peak_corr > clock.shifted_max).sum()), "of", len(clock))
20170113-4 NeuronDfOverF: subsets of 7448 and 7448 ROIs, correlation 0.981 20170228-4 GliaDfOverF: subsets of 16000 and 16000 ROIs, correlation 0.998
Neuronal whole-brain signal against swim rate (positive lag: signal after swims):
fish peak_corr peak_lag_s shifted_max shifted_median
20161109-2 0.162 0.349 0.087 0.008
20170113-2 0.219 4.764 0.095 0.033
20170113-3 0.339 2.199 0.133 0.025
20170113-4 0.294 4.398 0.056 0.026
20170127-1 0.421 0.369 0.188 0.126
20170203-1 0.190 1.966 0.056 0.011
20170228-3 0.285 1.762 0.175 0.137
20170228-4 0.384 1.902 0.209 0.082
Files where the true peak correlation exceeds every shifted one: 8 of 8
The whole-brain signal from two disjoint ROI subsets is nearly the same, so the subset stands for the whole brain. In every file with neuronal traces, the neuronal signal follows swimming more closely with the true bout times than with any of 20 shifted copies, with its peak 0.3 to 5 s after the swims, which is consistent with the stated shared clock. The four raw-only files are used for behaviour only, which is on the rig's clock by construction.
Epochs used by the claims¶
Passive and active epochs are built from the swim bouts by the Methods' rule (see fictive.py),
within each closed- or open-loop period; epochs cut by a period's edge are left out, as the claims
file fixes. This drops the chain of bouts that runs into a switch to open loop, and with it the
first passivity onset of an open-loop period when the fish swims across the switch. A variant
keeps these epochs (keep_edges); the synthetic check below shows the difference, and the table
gives the onsets each fish loses under the main rule.
In [10]:
EP = {s["fish"]: {c: fv.session_epochs(s, c) for c in ["closed", "open"]} for s in sessions}
ON = {s["fish"]: fv.passivity_onsets(s) for s in sessions}
rows = []
for f, e in EP.items():
rows.append({"fish": f, **{f"{c}_{k}": sum(len(x[k]) for x in e[c]) for c in ["closed", "open"] for k in ["passive", "active"]},
"passivity_onsets": len(ON[f])})
print(pd.DataFrame(rows).to_string(index=False))
# Variant: keep the epochs cut by a period's edges (fictive.epochs_in_segment, keep_edges). The
# main rule drops the chain of bouts that runs into a switch to open loop, and with it the first
# passivity onset of each open-loop period. Synthetic check: bouts every 0.5 s from 100.3 to 130 s,
# then a pause, in an open-loop period starting at 100 s.
syn = pd.DataFrame({"start_time": np.r_[np.arange(100.3, 129.9, 0.5), 140.0], "power": 1.0})
syn["stop_time"] = syn.start_time + 0.2
for keep in [False, True]:
e = fv.epochs_in_segment(syn, (100.0, 160.0), keep_edges=keep)
print(f"synthetic, keep_edges={keep}: active epochs {[(round(x['start'], 1), round(x['stop'], 1), x['followed_by_passive']) for x in e['active']]}")
EPX = {s["fish"]: {c: fv.session_epochs(s, c, keep_edges=True) for c in ["closed", "open"]} for s in sessions}
ONX = {s["fish"]: fv.passivity_onsets(s, keep_edges=True) for s in sessions}
dropped = pd.DataFrame({"fish": list(ON), "onsets_main": [len(ON[f]) for f in ON], "onsets_edges_kept": [len(ONX[f]) for f in ON],
"open_periods": [len(EPX[f]["open"]) for f in ON]})
dropped["dropped_by_main_rule"] = dropped.onsets_edges_kept - dropped.onsets_main
print("\nPassivity onsets under the main rule and with the edge epochs kept:")
print(dropped.to_string(index=False))
fish closed_passive closed_active open_passive open_active passivity_onsets 20161022-1 11 8 71 32 31 20161109-2 18 2 32 26 26 20161206-1 1 0 41 34 33 20161206-2 4 0 60 38 38 20170111-1 19 0 44 44 37 20170113-2 13 10 147 104 97 20170113-3 20 3 43 32 23 20170113-4 45 7 49 35 27 20170127-1 12 1 29 25 24 20170203-1 24 2 96 50 46 20170228-3 0 0 129 68 67 20170228-4 0 0 142 108 108 synthetic, keep_edges=False: active epochs [] synthetic, keep_edges=True: active epochs [(np.float64(100.3), np.float64(130.0), np.True_)]
Passivity onsets under the main rule and with the edge epochs kept:
fish onsets_main onsets_edges_kept open_periods dropped_by_main_rule
20161022-1 31 31 2 0
20161109-2 26 28 29 2
20161206-1 33 35 2 2
20161206-2 38 40 2 2
20170111-1 37 40 19 3
20170113-2 97 102 14 5
20170113-3 23 29 19 6
20170113-4 27 35 32 8
20170127-1 24 25 3 1
20170203-1 46 48 17 2
20170228-3 67 68 2 1
20170228-4 108 109 3 1
fig1h-epoch-durations¶
In prolonged open loop, fish cycle between passive epochs lasting 24.4 +/- 2.5 s and active epochs lasting 36.1 +/- 4.2 s (mean +/- SEM, n = 74 fish), swimming 29 +/- 6 times per active epoch.
Readings fixed in the claims file: (a) epochs from the swim bouts in open loop, mean over fish of
each fish's mean; (b) the dandiset's ActivityStates rows; (c) as (a), pooled over epochs.
Tolerances 1.96 SEM: 24.4 +/- 4.9 s, 36.1 +/- 8.2 s, 29 +/- 11.8 bouts.
In [11]:
EXPECTED = {"passive_s": (24.4, 4.9), "active_s": (36.1, 8.2), "bouts": (29, 11.8)}
rows, pooled = [], {"passive_s": [], "active_s": [], "bouts": []}
for s in sessions:
f = s["fish"]
pas = [b - a for x in EP[f]["open"] for a, b in x["passive"]]
act = [e["stop"] - e["start"] for x in EP[f]["open"] for e in x["active"]]
nb = [e["n_bouts"] for x in EP[f]["open"] for e in x["active"]]
st = s["states"]
tab = st.assign(d=st.stop - st.start)
rows.append({"fish": f, "n_passive": len(pas), "passive_s": np.mean(pas), "n_active": len(act), "active_s": np.mean(act),
"bouts": np.mean(nb), "table_passive_s": tab[tab.type == "passive"].d.mean(),
"table_active_s": tab[tab.type == "active"].d.mean(),
"table_active_over_5s": tab[(tab.type == "active") & (tab.d > 5)].d.mean()})
pooled["passive_s"] += pas; pooled["active_s"] += act; pooled["bouts"] += nb
fig1h = pd.DataFrame(rows)
print(fig1h.round(2).to_string(index=False))
readings = {
"(a) bouts, mean of fish means": {k: fig1h[k].mean() for k in EXPECTED},
"(b) ActivityStates table, mean of fish means": {"passive_s": fig1h.table_passive_s.mean(), "active_s": fig1h.table_active_s.mean(), "bouts": np.nan},
"(c) bouts, pooled epochs": {k: np.mean(v) for k, v in pooled.items()},
}
print()
for name, vals in readings.items():
ok = {k: abs(vals[k] - EXPECTED[k][0]) <= EXPECTED[k][1] for k in EXPECTED}
print(f"{name}: " + "; ".join(f"{k} {vals[k]:.2f} ({'within' if ok[k] else 'outside'} {EXPECTED[k][0]} +/- {EXPECTED[k][1]})" for k in EXPECTED)
+ f" -> {sum(ok.values())} of 3 within")
print(f"\nSEM over fish under (a): passive {fig1h.passive_s.sem():.2f} s, active {fig1h.active_s.sem():.2f} s, bouts {fig1h.bouts.sem():.2f}")
print(f"Table variant with active rows over 5 s only: active {fig1h.table_active_over_5s.mean():.2f} s")
print(f"Range over fish (a): passive {fig1h.passive_s.min():.1f} to {fig1h.passive_s.max():.1f} s, active {fig1h.active_s.min():.1f} to {fig1h.active_s.max():.1f} s")
# Variant: the threshold of 5 s is the paper's; show 4 and 6 s.
for gap in [4.0, 6.0]:
p, a = [], []
for s in sessions:
ep = [fv.epochs_in_segment(s["swims"], seg, gap) for seg in fv.segments(s["cond_t"], fv.condition_masks(s)["open"])]
p.append(np.mean([y - x for e in ep for x, y in e["passive"]]))
a.append(np.mean([e["stop"] - e["start"] for x in ep for e in x["active"]]))
print(f"Pause threshold {gap:.0f} s: passive {np.nanmean(p):.2f} s, active {np.nanmean(a):.2f} s (mean of fish means)")
# Variant: epochs cut by the edges of an open-loop period kept (reading (a) otherwise).
rows = []
for s in sessions:
f = s["fish"]
pas = [b - a for x in EPX[f]["open"] for a, b in x["passive"]]
act = [e["stop"] - e["start"] for x in EPX[f]["open"] for e in x["active"]]
nb = [e["n_bouts"] for x in EPX[f]["open"] for e in x["active"]]
rows.append({"fish": f, "n_passive": len(pas), "passive_s": np.mean(pas), "n_active": len(act), "active_s": np.mean(act), "bouts": np.mean(nb)})
fig1hx = pd.DataFrame(rows)
print("\nVariant with edge epochs kept:")
print(fig1hx.round(2).to_string(index=False))
vals = {k: fig1hx[k].mean() for k in EXPECTED}
ok = {k: abs(vals[k] - EXPECTED[k][0]) <= EXPECTED[k][1] for k in EXPECTED}
print("Edge epochs kept, mean of fish means: " + "; ".join(f"{k} {vals[k]:.2f} ({'within' if ok[k] else 'outside'} {EXPECTED[k][0]} +/- {EXPECTED[k][1]})" for k in EXPECTED)
+ f" -> {sum(ok.values())} of 3 within")
fish n_passive passive_s n_active active_s bouts table_passive_s table_active_s table_active_over_5s 20161022-1 71 32.54 32 8.52 6.31 12.55 5.98 5.98 20161109-2 32 30.73 26 64.15 64.77 28.62 60.19 74.00 20161206-1 41 36.65 34 19.64 17.50 3.51 9.91 9.91 20161206-2 60 9.68 38 29.28 19.79 6.38 9.08 9.08 20170111-1 44 21.58 44 19.29 16.75 22.55 15.44 17.35 20170113-2 147 18.30 104 7.31 4.10 19.01 6.46 6.90 20170113-3 43 18.14 32 26.69 20.84 27.63 19.38 27.26 20170113-4 49 16.65 35 20.03 16.49 22.73 14.34 18.52 20170127-1 29 16.92 25 25.91 23.64 25.81 21.64 24.70 20170203-1 96 11.96 50 11.14 12.24 13.86 5.76 9.97 20170228-3 129 19.50 68 7.32 3.90 22.38 5.11 7.23 20170228-4 142 12.98 108 28.99 22.41 13.21 15.96 20.09 (a) bouts, mean of fish means: passive_s 20.47 (within 24.4 +/- 4.9); active_s 22.36 (outside 36.1 +/- 8.2); bouts 19.06 (within 29 +/- 11.8) -> 2 of 3 within (b) ActivityStates table, mean of fish means: passive_s 18.19 (outside 24.4 +/- 4.9); active_s 15.77 (outside 36.1 +/- 8.2); bouts nan (outside 29 +/- 11.8) -> 0 of 3 within (c) bouts, pooled epochs: passive_s 18.81 (outside 24.4 +/- 4.9); active_s 19.66 (outside 36.1 +/- 8.2); bouts 15.99 (outside 29 +/- 11.8) -> 0 of 3 within SEM over fish under (a): passive 2.46 s, active 4.49 s, bouts 4.61 Table variant with active rows over 5 s only: active 19.25 s Range over fish (a): passive 9.7 to 36.6 s, active 7.3 to 64.1 s Pause threshold 4 s: passive 19.27 s, active 19.99 s (mean of fish means)
Pause threshold 6 s: passive 21.46 s, active 25.31 s (mean of fish means)
Variant with edge epochs kept:
fish n_passive passive_s n_active active_s bouts
20161022-1 74 32.26 32 8.52 6.31
20161109-2 42 25.57 47 47.02 48.53
20161206-1 43 37.09 36 18.97 16.97
20161206-2 60 9.68 42 32.20 22.55
20170111-1 68 20.87 57 18.76 16.63
20170113-2 164 18.03 111 7.28 4.06
20170113-3 67 19.86 43 28.05 22.56
20170113-4 81 16.72 58 19.98 16.33
20170127-1 32 25.24 27 25.60 23.89
20170203-1 112 12.47 55 11.17 12.33
20170228-3 131 19.61 69 7.44 3.97
20170228-4 144 13.04 110 29.39 22.55
Edge epochs kept, mean of fish means: passive_s 20.87 (within 24.4 +/- 4.9); active_s 21.20 (outside 36.1 +/- 8.2); bouts 18.06 (within 29 +/- 11.8) -> 2 of 3 within
Result: partially reproduced. Under reading (a), epochs built from the swim bouts in open
loop, passive epochs last 20.47 s on average over the 12 fish (within 24.4 +/- 4.9) and active
epochs hold 19.06 bouts (within 29 +/- 11.8), but active epochs last 22.36 s, outside 36.1 +/- 8.2.
Readings (b), the ActivityStates table (18.19 and 15.77 s), and (c), pooled epochs (18.81 s,
19.66 s, 15.99 bouts), put none of the values within tolerance. Fish differ widely (active epochs
7.3 to 64.1 s). With a pause threshold of 6 s the means are 21.46 and 25.31 s. Keeping the epochs
cut by the edges of an open-loop period (a variant) gives 20.87 s, 21.20 s and 18.06 bouts, again
2 of 3 within tolerance, so leaving them out does not explain the short active epochs.
methods-isi¶
In closed loop fish swim in discrete bouts about once per second (inter-swim interval 1.2 +/- 0.05 s, n = 74), and within active epochs of open loop the interval is 1.19 +/- 0.047 s.
Readings: (a) start to start, (b) end to start; in closed loop (i) intervals under 5 s, (ii) all. Tolerances 1.96 SEM: 1.2 +/- 0.098 s and 1.19 +/- 0.092 s.
In [12]:
rows = []
for s in sessions:
f = s["fish"]
r = {"fish": f}
for c in ["closed"]:
ss, ee = [], []
for seg in fv.segments(s["cond_t"], fv.condition_masks(s)[c]):
w = s["swims"][(s["swims"].start_time >= seg[0]) & (s["swims"].stop_time <= seg[1])].sort_values("start_time")
ss += list(np.diff(w.start_time.to_numpy()))
ee += list(w.start_time.to_numpy()[1:] - w.stop_time.to_numpy()[:-1])
ss, ee = np.array(ss), np.array(ee)
r.update({"closed_a_i": ss[ee < 5].mean(), "closed_a_ii": ss.mean(), "closed_b_i": ee[ee < 5].mean(), "closed_b_ii": ee.mean(),
"closed_n": len(ss)})
act = [e for x in EP[f]["open"] for e in x["active"]]
r["active_a"] = np.mean(np.concatenate([np.diff(e["bout_starts"]) for e in act]))
sw = s["swims"].set_index("start_time")
gaps = []
for e in act:
b = s["swims"][(s["swims"].start_time >= e["start"]) & (s["swims"].stop_time <= e["stop"])].sort_values("start_time")
gaps += list(b.start_time.to_numpy()[1:] - b.stop_time.to_numpy()[:-1])
r["active_b"] = np.mean(gaps)
rows.append(r)
isi = pd.DataFrame(rows)
print(isi.round(3).to_string(index=False))
print()
for cl, op in [("closed_a_i", "active_a"), ("closed_a_ii", "active_a"), ("closed_b_i", "active_b"), ("closed_b_ii", "active_b")]:
c, o = isi[cl].mean(), isi[op].mean()
ok_c, ok_o = abs(c - 1.2) <= 0.098, abs(o - 1.19) <= 0.092
print(f"{cl} / {op}: closed loop {c:.3f} s ({'within' if ok_c else 'outside'} 1.2 +/- 0.098), "
f"active epochs {o:.3f} s ({'within' if ok_o else 'outside'} 1.19 +/- 0.092); SEM {isi[cl].sem():.3f} and {isi[op].sem():.3f}")
print(f"\nRange over fish, closed_a_i {isi.closed_a_i.min():.2f} to {isi.closed_a_i.max():.2f} s; active_a {isi.active_a.min():.2f} to {isi.active_a.max():.2f} s")
fish closed_a_i closed_a_ii closed_b_i closed_b_ii closed_n active_a active_b 20161022-1 1.281 3.236 1.080 3.033 111 1.501 1.021 20161109-2 0.955 2.445 0.840 2.332 162 1.002 0.784 20161206-1 0.916 0.927 0.745 0.755 568 1.177 0.946 20161206-2 1.351 1.420 1.225 1.291 370 1.528 1.270 20170111-1 1.743 13.695 1.581 13.538 35 1.207 0.923 20170113-2 0.909 1.237 0.707 1.037 539 1.243 0.464 20170113-3 1.191 2.693 1.088 2.588 163 1.331 0.982 20170113-4 1.063 2.858 0.786 2.590 412 1.247 0.816 20170127-1 0.992 1.312 0.838 1.160 899 1.133 0.950 20170203-1 1.314 2.816 1.197 2.700 148 0.976 0.790 20170228-3 0.843 0.843 0.673 0.673 146 2.299 0.700 20170228-4 1.156 1.156 1.009 1.009 94 1.336 0.932 closed_a_i / active_a: closed loop 1.143 s (within 1.2 +/- 0.098), active epochs 1.332 s (outside 1.19 +/- 0.092); SEM 0.074 and 0.100 closed_a_ii / active_a: closed loop 2.886 s (outside 1.2 +/- 0.098), active epochs 1.332 s (outside 1.19 +/- 0.092); SEM 1.013 and 0.100 closed_b_i / active_b: closed loop 0.981 s (outside 1.2 +/- 0.098), active epochs 0.881 s (outside 1.19 +/- 0.092); SEM 0.077 and 0.057 closed_b_ii / active_b: closed loop 2.725 s (outside 1.2 +/- 0.098), active epochs 0.881 s (outside 1.19 +/- 0.092); SEM 1.013 and 0.057 Range over fish, closed_a_i 0.84 to 1.74 s; active_a 0.98 to 2.30 s
Result: partially reproduced. Read as start to start, over intervals under 5 s, the closed-loop inter-swim interval is 1.143 s (within 1.2 +/- 0.098), but within active epochs of open loop it is 1.332 s, outside 1.19 +/- 0.092. No other reading puts both values within tolerance; end-to-start intervals give 0.981 and 0.881 s, and all closed-loop intervals 2.886 s.
fig1j-open-loop-passivity¶
Fish spend much more time in the passive state when swimming produces no visual feedback (open loop) than in closed loop with a naturalistic gain, where passivity is rare.
Criterion: fraction of time in pauses of 5 s or more, open against closed loop, one-sided Wilcoxon over fish at p < 0.05, and median closed-loop fraction below 10%.
In [13]:
rows = []
for s in sessions:
f = s["fish"]
r = {"fish": f}
for c in ["closed", "open"]:
tot = sum(x["seg"][1] - x["seg"][0] for x in EP[f][c])
r[f"{c}_s"] = tot
r[f"{c}_passive_frac"] = sum(b - a for x in EP[f][c] for a, b in x["pauses"]) / tot
# Variant: closed loop before the first open loop only.
first_open = cond.set_index("fish").first_open_s[f]
pre = [x for x in EP[f]["closed"] if x["seg"][1] <= first_open]
tot = sum(x["seg"][1] - x["seg"][0] for x in pre)
r["closed_before_first_open_frac"] = sum(b - a for x in pre for a, b in x["pauses"]) / tot if tot else np.nan
r["closed_before_first_open_s"] = tot
rows.append(r)
frac = pd.DataFrame(rows)
print(frac.round(3).to_string(index=False))
w = stats.wilcoxon(frac.open_passive_frac, frac.closed_passive_frac, alternative="greater")
print(f"\nOpen loop: median {frac.open_passive_frac.median():.3f}; closed loop: median {frac.closed_passive_frac.median():.3f}; "
f"open higher in {(frac.open_passive_frac > frac.closed_passive_frac).sum()} of {len(frac)} fish; Wilcoxon p = {w.pvalue:.3g}")
v = frac.dropna(subset=["closed_before_first_open_frac"])
w2 = stats.wilcoxon(v.open_passive_frac, v.closed_before_first_open_frac, alternative="greater")
print(f"Variant, closed loop before the first open loop: median {v.closed_before_first_open_frac.median():.3f}; "
f"open higher in {(v.open_passive_frac > v.closed_before_first_open_frac).sum()} of {len(v)} fish; Wilcoxon p = {w2.pvalue:.3g}")
# By gain level in closed loop.
rows = []
for s in sessions:
m = fv.condition_masks(s)
for g in np.unique(np.round(s["gain"][m["closed"]], 3)):
mask = m["closed"] & (np.round(s["gain"], 3) == g)
segs = fv.segments(s["cond_t"], mask)
tot = sum(b - a for a, b in segs)
if tot < 60:
continue
ep = [fv.epochs_in_segment(s["swims"], seg) for seg in segs]
rows.append({"fish": s["fish"], "gain": g, "seconds": tot, "passive_frac": sum(b - a for e in ep for a, b in e["pauses"]) / tot})
print("\nClosed-loop passive fraction by gain level (levels with 60 s or more):")
print(pd.DataFrame(rows).round(3).to_string(index=False))
# By gain level, split into closed loop before and after the first open loop, with the protocol
# codes of the blocks (levels with 30 s or more in a part).
rows = []
for s in sessions:
m = fv.condition_masks(s)
first_open = cond.set_index("fish").first_open_s[s["fish"]]
g_r = np.round(s["gain"], 3)
for part, pm in [("before", s["cond_t"] < first_open), ("after", s["cond_t"] >= first_open)]:
for g in np.unique(g_r[m["closed"] & pm]):
mask = m["closed"] & pm & (g_r == g)
segs = fv.segments(s["cond_t"], mask)
tot = sum(b - a for a, b in segs)
if tot < 30:
continue
ep = [fv.epochs_in_segment(s["swims"], seg) for seg in segs]
rows.append({"fish": s["fish"], "part": part, "gain": g, "codes": sorted(set(s["protocol_code"][mask].astype(int))),
"seconds": tot, "passive_frac": sum(b - a for e in ep for a, b in e["pauses"]) / tot})
bygain = pd.DataFrame(rows)
print("\nClosed-loop passive fraction by gain level, before and after the first open loop:")
print(bygain.round(3).to_string(index=False))
for part in ["before", "after"]:
d = bygain[bygain.part == part]
ends = d.sort_values("gain").groupby("fish").filter(lambda x: len(x) > 1).groupby("fish").passive_frac.agg(["first", "last"])
print(f"{part} the first open loop: {d.fish.nunique()} fish, {len(d)} gain levels, passive fraction median {d.passive_frac.median():.3f} "
f"(range {d.passive_frac.min():.3f} to {d.passive_frac.max():.3f}); of {len(ends)} fish with 2 or more levels, "
f"passivity lower at the highest gain than at the lowest in {(ends['last'] < ends['first']).sum()}")
fish closed_s closed_passive_frac open_s open_passive_frac closed_before_first_open_frac closed_before_first_open_s 20161022-1 530.535 0.740 2779.586 0.859 0.797 360.724 20161109-2 637.043 0.736 3349.924 0.321 0.589 401.627 20161206-1 823.355 0.197 2296.253 0.695 0.135 347.023 20161206-2 744.250 0.142 1983.432 0.293 0.063 338.923 20170111-1 1360.091 0.968 2511.268 0.565 NaN 0.000 20170113-2 864.858 0.409 3990.967 0.741 0.017 300.620 20170113-3 1040.569 0.808 2571.472 0.517 0.615 300.620 20170113-4 1943.030 0.766 2619.875 0.517 0.000 193.113 20170127-1 1769.718 0.397 1501.900 0.538 0.099 746.250 20170203-1 619.441 0.683 2120.842 0.658 NaN 0.000 20170228-3 176.212 0.169 3294.820 0.780 0.110 164.511 20170228-4 164.411 0.055 5173.445 0.363 0.059 152.310 Open loop: median 0.551; closed loop: median 0.546; open higher in 7 of 12 fish; Wilcoxon p = 0.285 Variant, closed loop before the first open loop: median 0.104; open higher in 8 of 10 fish; Wilcoxon p = 0.0137
Closed-loop passive fraction by gain level (levels with 60 s or more):
fish gain seconds passive_frac
20161022-1 0.012 190.613 0.594
20161022-1 0.015 71.905 0.740
20161022-1 0.018 97.106 0.528
20161022-1 0.036 170.111 0.946
20161109-2 0.010 366.324 0.524
20161109-2 0.030 244.216 0.937
20161206-1 0.015 98.007 0.256
20161206-1 0.018 725.348 0.189
20161206-2 0.020 744.250 0.142
20170111-1 0.012 1186.779 0.968
20170111-1 0.015 82.105 0.903
20170111-1 0.036 87.906 0.999
20170113-2 0.010 278.619 0.167
20170113-2 0.020 263.718 0.625
20170113-2 0.030 308.721 0.399
20170113-3 0.012 205.014 0.668
20170113-3 0.020 619.741 0.873
20170113-3 0.036 201.413 0.668
20170113-4 0.010 190.813 0.359
20170113-4 0.015 959.964 0.882
20170113-4 0.020 109.807 0.921
20170113-4 0.030 174.412 0.441
20170127-1 0.010 141.609 0.348
20170127-1 0.012 1041.670 0.100
20170127-1 0.030 302.720 0.937
20170203-1 0.015 619.441 0.683
20170228-3 0.020 164.511 0.110
20170228-4 0.020 164.411 0.055
Closed-loop passive fraction by gain level, before and after the first open loop:
fish part gain codes seconds passive_frac
20161022-1 before 0.012 [22] 190.613 0.594
20161022-1 before 0.036 [23] 170.111 0.946
20161022-1 after 0.015 [22] 71.905 0.740
20161022-1 after 0.018 [22, 23] 97.106 0.528
20161109-2 before 0.010 [22, 23] 261.817 0.334
20161109-2 before 0.030 [23] 139.709 0.910
20161109-2 after 0.010 [22, 23] 104.507 1.000
20161109-2 after 0.030 [23] 104.507 0.972
20161206-1 before 0.015 [22, 23] 98.007 0.256
20161206-1 before 0.018 [22, 23] 249.017 0.087
20161206-1 after 0.018 [22, 23] 476.332 0.243
20161206-2 before 0.020 [22, 23] 338.923 0.063
20161206-2 after 0.020 [22, 23] 405.327 0.209
20170111-1 before 0.012 [22] 92.406 0.959
20170111-1 before 0.036 [23] 87.906 0.999
20170111-1 after 0.012 [51] 1094.373 0.968
20170111-1 after 0.015 [51] 82.105 0.903
20170113-2 before 0.010 [22] 146.610 0.000
20170113-2 before 0.030 [23] 146.510 0.000
20170113-2 after 0.010 [22] 132.009 0.353
20170113-2 after 0.020 [51] 263.718 0.625
20170113-2 after 0.030 [23] 162.211 0.760
20170113-3 before 0.012 [22] 146.610 0.576
20170113-3 before 0.036 [23] 146.510 0.570
20170113-3 after 0.012 [22] 58.404 0.898
20170113-3 after 0.020 [51] 619.741 0.873
20170113-3 after 0.036 [23] 54.904 0.929
20170113-4 before 0.010 [22] 102.707 0.000
20170113-4 before 0.030 [23] 87.906 0.000
20170113-4 after 0.010 [22, 51, 52] 88.106 0.776
20170113-4 after 0.015 [51, 52] 959.964 0.882
20170113-4 after 0.020 [51] 109.807 0.921
20170113-4 after 0.030 [23] 86.506 0.888
20170127-1 before 0.012 [22, 23] 746.250 0.099
20170127-1 after 0.010 [22, 23, 51, 52] 141.609 0.348
20170127-1 after 0.012 [22, 23, 52] 295.420 0.102
20170127-1 after 0.030 [51, 52] 302.720 0.937
20170203-1 after 0.015 [51] 619.441 0.683
20170228-3 before 0.020 [12, 22, 23] 164.511 0.110
20170228-4 before 0.020 [22, 23] 164.411 0.055
before the first open loop: 11 fish, 18 gain levels, passive fraction median 0.183 (range 0.000 to 0.999); of 7 fish with 2 or more levels, passivity lower at the highest gain than at the lowest in 2
after the first open loop: 10 fish, 22 gain levels, passive fraction median 0.825 (range 0.102 to 1.000); of 7 fish with 2 or more levels, passivity lower at the highest gain than at the lowest in 3
Result: not reproduced. Over all closed-loop time the fish spend a median 0.546 of the time in pauses of 5 s or more, against 0.551 in open loop; open loop is higher in 7 of 12 fish (Wilcoxon p = 0.285), and the closed-loop median is far above 10%. Restricted to closed loop before the first open loop (a variant), the median is 0.104 and open loop is higher in 8 of 10 fish (p = 0.0137). By gain level, closed-loop passivity has a median of 0.183 before the first open loop and 0.825 after it. Passivity shows no consistent fall with the gain: it is lower at the highest gain than at the lowest in 2 of 7 fish before the first open loop and 3 of 7 after it, though before the first open loop 2 fish are tied at 0.000 and one differs by 0.006. The gain is in arbitrary units and no level can be identified as the paper's naturalistic gain.
The outcome rests on the choice of closed-loop time. The closed-loop blocks after the first open loop, where passivity is high, decide it. The variant restricted to closed loop before the first open loop is closer to the paper's comparison and meets the criterion for partially reproduced (test passes, median 0.104 above 10%). The main reading was fixed before the analysis and is kept.
fig1c-vigor-open-loop-onset¶
When switched from closed to open loop, fish first swim with higher vigor, then stop swimming.
Criterion: per switch, mean power of the bouts from the switch to the first passive epoch over the mean power of the closed-loop bouts in the 30 s before; per fish the median log ratio; one-sided Wilcoxon over fish. Switches with fewer than 3 bouts on either side are left out.
In [14]:
def switch_ratios(s, strict=False):
m = fv.condition_masks(s)
t = s["cond_t"]
sw = s["swims"].sort_values("start_time")
out = []
switches = t[1:][m["closed"][:-1] & m["open"][1:]]
for t0 in switches:
before = sw[(sw.start_time >= t0 - 30) & (sw.stop_time <= t0)]
before = before[fv.label_at(s, before.start_time.to_numpy()) == 1]
if strict and not m["closed"][(t >= t0 - 30) & (t < t0)].all():
continue
t_end = t[np.flatnonzero((t > t0) & ~m["open"])[0]] if ((t > t0) & ~m["open"]).any() else t[-1]
after = sw[(sw.start_time >= t0) & (sw.start_time < t_end)]
st, sp = after.start_time.to_numpy(), after.stop_time.to_numpy()
if len(st) == 0 or st[0] - t0 >= fv.PASSIVE_GAP_S:
out.append({"t": t0, "n_before": len(before), "n_after": 0, "log_ratio": np.nan})
continue
g = np.flatnonzero(st[1:] - sp[:-1] >= fv.PASSIVE_GAP_S)
after = after.iloc[:(g[0] + 1) if len(g) else len(st)]
after = after[fv.label_at(s, after.start_time.to_numpy()) == 2]
ok = len(before) >= 3 and len(after) >= 3 and before.power.mean() > 0
out.append({"t": t0, "n_before": len(before), "n_after": len(after),
"log_ratio": np.log(after.power.mean() / before.power.mean()) if ok else np.nan})
return pd.DataFrame(out)
# Synthetic check: closed loop 0 to 60 s, open loop 60 to 70 s, then closed loop. Bouts every 2 s
# before the switch (power 1), every 1 s in the open loop (power 3), and every 1 s afterwards
# (power 100) with no pause of 5 s. Only the open-loop bouts may count: the ratio must be 3.
syn_t = np.arange(0, 200, 0.1)
syn = {"cond_t": syn_t, "gain": np.where((syn_t >= 60) & (syn_t < 70), 0.0, 0.02), "velocity": np.ones_like(syn_t),
"swims": pd.DataFrame({"start_time": np.r_[np.arange(30, 60, 2.0), np.arange(60.5, 70, 1.0), np.arange(70.5, 150, 1.0)]})}
syn["swims"]["stop_time"] = syn["swims"].start_time + 0.2
syn["swims"]["power"] = np.where(syn["swims"].start_time < 60, 1.0, np.where(syn["swims"].start_time < 70, 3.0, 100.0))
print("synthetic switch, power ratio:", np.round(np.exp(switch_ratios(syn).log_ratio.to_numpy()), 3))
rows = []
for s in sessions:
r = switch_ratios(s)
rs = switch_ratios(s, strict=True)
rows.append({"fish": s["fish"], "switches": len(r), "used": int(r.log_ratio.notna().sum()) if len(r) else 0,
"median_log_ratio": r.log_ratio.median() if len(r) else np.nan,
"strict_used": int(rs.log_ratio.notna().sum()) if len(rs) else 0,
"strict_median_log_ratio": rs.log_ratio.median() if len(rs) else np.nan})
vig = pd.DataFrame(rows)
vig["median_ratio"] = np.exp(vig.median_log_ratio)
print(vig.round(3).to_string(index=False))
d = vig.median_log_ratio.dropna()
w = stats.wilcoxon(d, alternative="greater")
print(f"\n{len(d)} fish with a usable switch; ratio above 1 in {(d > 0).sum()}; median ratio {np.exp(d.median()):.2f} "
f"(range {np.exp(d.min()):.2f} to {np.exp(d.max()):.2f}); Wilcoxon p = {w.pvalue:.3g}")
ds = vig.strict_median_log_ratio.dropna()
if len(ds) > 1:
print(f"Variant requiring 30 s of closed loop before the switch: {len(ds)} fish, ratio above 1 in {(ds > 0).sum()}, "
f"median ratio {np.exp(ds.median()):.2f}, Wilcoxon p = {stats.wilcoxon(ds, alternative='greater').pvalue:.3g}")
synthetic switch, power ratio: [3.]
fish switches used median_log_ratio strict_used strict_median_log_ratio median_ratio 20161022-1 1 0 NaN 0 NaN NaN 20161109-2 3 0 NaN 0 NaN NaN 20161206-1 1 1 0.623 0 NaN 1.865 20161206-2 2 2 1.682 0 NaN 5.375 20170111-1 19 0 NaN 0 NaN NaN 20170113-2 14 3 2.249 2 2.556 9.483 20170113-3 19 4 1.918 4 1.918 6.805 20170113-4 31 10 1.249 10 1.249 3.487 20170127-1 2 1 0.097 0 NaN 1.102 20170203-1 18 8 -0.193 7 -0.167 0.824 20170228-3 4 1 2.332 0 NaN 10.297 20170228-4 3 1 -0.251 0 NaN 0.778 9 fish with a usable switch; ratio above 1 in 7; median ratio 3.49 (range 0.78 to 10.30); Wilcoxon p = 0.0195 Variant requiring 30 s of closed loop before the switch: 4 fish, ratio above 1 in 3, median ratio 4.87, Wilcoxon p = 0.125
Result: reproduced. In the 9 fish with a usable switch from closed to open loop, bouts after the switch are more powerful than closed-loop bouts in the 30 s before in 7 fish, with a median ratio of 3.49 (0.78 to 10.30), Wilcoxon p = 0.0195. Requiring 30 s of uninterrupted closed loop before the switch leaves 4 fish (median ratio 4.87, p = 0.125).
fig1i-isi-before-passivity¶
Before the switch to passivity the inter-swim interval increases progressively, and fish swim more vigorously and less frequently in the seconds before passivity.
Criterion: in active epochs of 10 or more bouts that end in passivity (open loop), intervals k = 1 (last) to 7 back and bout powers; interval at k = 1 longer than at k = 7 (one-sided Wilcoxon over fish), negative Spearman correlation of k with the group mean, and the last 3 bouts more powerful than bouts 8 to 10 back (one-sided Wilcoxon).
In [15]:
def before_passivity(min_bouts, epochs=EP):
rows = []
for s in sessions:
isis, pws = [], []
for x in epochs[s["fish"]]["open"]:
for e in x["active"]:
if e["followed_by_passive"] and e["n_bouts"] >= min_bouts:
isis.append(np.diff(e["bout_starts"])[::-1][:7])
pw = e["bout_powers"][::-1][:10]
pws.append(np.r_[pw, np.full(10 - len(pw), np.nan)])
if isis:
I, P = np.mean(isis, 0), np.nanmean(pws, 0)
rows.append({"fish": s["fish"], "epochs": len(isis), **{f"isi_k{k + 1}": I[k] for k in range(7)},
"power_last3": P[:3].mean(), "power_8to10": P[7:10].mean()})
return pd.DataFrame(rows)
def judge(df, label):
w_isi = stats.wilcoxon(df.isi_k1, df.isi_k7, alternative="greater")
means = df[[f"isi_k{k}" for k in range(1, 8)]].mean()
rho = stats.spearmanr(np.arange(1, 8), means.to_numpy())
w_pow = stats.wilcoxon(df.power_last3, df.power_8to10, alternative="greater", nan_policy="omit")
print(f"{label}: {len(df)} fish; mean interval k=1 {df.isi_k1.mean():.3f} s, k=7 {df.isi_k7.mean():.3f} s, "
f"k=1 longer in {(df.isi_k1 > df.isi_k7).sum()} fish, Wilcoxon p = {w_isi.pvalue:.3g}; "
f"Spearman rho (k, mean) = {rho.statistic:.3f}; power last 3 {df.power_last3.mean():.3f} against bouts 8-10 back "
f"{df.power_8to10.mean():.3f}, higher in {(df.power_last3 > df.power_8to10).sum()} fish, Wilcoxon p = {w_pow.pvalue:.3g}")
return means
bp = before_passivity(10)
print(bp.round(3).to_string(index=False))
print()
means = judge(bp, "Main (10 or more bouts)")
print("Mean interval by k (1 = last):", means.round(3).to_dict())
judge(before_passivity(8), "Variant (8 or more bouts)")
judge(before_passivity(15), "Variant (15 or more bouts)")
_ = judge(before_passivity(10, EPX), "Variant (10 or more bouts, edge epochs kept)")
fish epochs isi_k1 isi_k2 isi_k3 isi_k4 isi_k5 isi_k6 isi_k7 power_last3 power_8to10
20161022-1 2 1.495 1.891 1.394 1.890 1.285 1.539 2.031 2.082 2.014
20161109-2 25 1.283 1.272 1.343 1.031 1.010 1.114 1.171 0.505 0.451
20161206-1 25 1.538 1.198 1.372 1.330 1.353 1.256 1.289 0.290 0.304
20161206-2 17 1.977 1.891 1.882 1.661 1.999 1.756 1.716 1.621 1.037
20170111-1 31 1.614 1.440 1.366 1.375 1.347 1.408 1.185 0.401 0.443
20170113-2 7 0.869 0.687 0.608 0.940 0.840 1.178 1.202 0.218 0.220
20170113-3 16 2.311 1.331 1.612 1.570 1.460 1.460 1.092 0.573 0.468
20170113-4 20 1.522 1.145 1.802 1.290 1.776 1.020 1.610 0.487 0.651
20170127-1 22 1.237 1.176 1.169 1.138 1.169 1.132 1.061 0.274 0.220
20170203-1 36 1.316 1.005 0.978 0.997 0.874 0.909 0.747 0.161 0.167
20170228-4 71 1.399 1.299 1.401 1.351 1.480 1.244 1.322 0.382 0.349
Main (10 or more bouts): 11 fish; mean interval k=1 1.506 s, k=7 1.312 s, k=1 longer in 8 fish, Wilcoxon p = 0.103; Spearman rho (k, mean) = -0.536; power last 3 0.636 against bouts 8-10 back 0.575, higher in 6 fish, Wilcoxon p = 0.16
Mean interval by k (1 = last): {'isi_k1': 1.506, 'isi_k2': 1.303, 'isi_k3': 1.357, 'isi_k4': 1.325, 'isi_k5': 1.327, 'isi_k6': 1.274, 'isi_k7': 1.312}
Variant (8 or more bouts): 12 fish; mean interval k=1 1.502 s, k=7 1.249 s, k=1 longer in 10 fish, Wilcoxon p = 0.00806; Spearman rho (k, mean) = -0.286; power last 3 0.692 against bouts 8-10 back 0.539, higher in 7 fish, Wilcoxon p = 0.0269
Variant (15 or more bouts): 10 fish; mean interval k=1 1.534 s, k=7 1.116 s, k=1 longer in 9 fish, Wilcoxon p = 0.00195; Spearman rho (k, mean) = -0.607; power last 3 0.429 against bouts 8-10 back 0.406, higher in 4 fish, Wilcoxon p = 0.423
Variant (10 or more bouts, edge epochs kept): 11 fish; mean interval k=1 1.468 s, k=7 1.322 s, k=1 longer in 8 fish, Wilcoxon p = 0.0874; Spearman rho (k, mean) = -0.357; power last 3 0.643 against bouts 8-10 back 0.574, higher in 6 fish, Wilcoxon p = 0.139
Result: not reproduced. In active epochs of 10 or more bouts that end in passivity (11 fish), the last interval (1.506 s) is longer than the seventh from last (1.312 s) in 8 fish, but not significantly (Wilcoxon p = 0.103). The Spearman rho of k against the mean interval is -0.536, but only k = 1 stands out: the means for k = 2 to 7 lie between 1.274 and 1.357 s with no trend. The last 3 bouts are not more powerful than bouts 8 to 10 back (0.636 against 0.575, higher in 6 fish, p = 0.16). The outcome depends on the minimum epoch length, which the paper does not state: with epochs of 8 or more bouts both parts hold (intervals p = 0.00806, rho = -0.286; power p = 0.0269), and with 15 or more the interval part holds (p = 0.00195) and the power part does not (p = 0.423). Keeping the epochs cut by either edge of an open-loop period (a variant) gives p = 0.0874 for the intervals and p = 0.139 for the power.
Signals around passivity onset¶
The whole-brain signal of each cell type around each passivity onset, against the baseline 10 to 5 s before onset (in the units of the stored series; see the data check).
In [16]:
LAGS = np.arange(-15, 20.01, 0.5)
TRACES, rows = {}, []
for s in sessions:
for name in s["imaging"]:
f = s["fish"]
t = s["imaging"][name]["t"]
y = WB[(f, name, 0.0)][0].astype(float)
ok = ~fv.censored(y)
t, y = t[ok], y[ok]
on = ON[f]
base = fv.window_mean(t, y, on, fv.BASELINE)
E = fv.event_windows(t, y, on, LAGS) - base[:, None]
TRACES[(f, name)] = (t, y, on, base, E)
rows.append({"fish": f, "series": name, "onsets": len(on), "trace_sd": y.std()})
print(pd.DataFrame(rows).round(5).to_string(index=False))
fig, axes = plt.subplots(1, 2, figsize=(11, 3.5))
for ax, name in zip(axes, ["NeuronDfOverF", "GliaDfOverF"]):
for (f, n), (t, y, on, base, E) in TRACES.items():
if n == name:
avg = np.nanmean(E, 0) / y.std()
ax.plot(LAGS, avg, label=f)
print(f"{name} {f}: onset-triggered mean in SD units ranges {np.nanmin(avg):.2f} to {np.nanmax(avg):.2f}")
ax.axvline(0, color="k", lw=0.5); ax.axhline(0, color="k", lw=0.5)
ax.set_xlabel("time from passivity onset (s)"); ax.set_ylabel("change from baseline (SD of trace)")
ax.set_title(name.replace("DfOverF", "") + " whole-brain signal"); ax.legend(fontsize=7)
plt.tight_layout(); plt.show()
fish series onsets trace_sd 20161109-2 GliaDfOverF 26 0.00160 20161109-2 NeuronDfOverF 26 0.00842 20170113-2 NeuronDfOverF 97 0.01457 20170113-3 NeuronDfOverF 23 0.01800 20170113-4 NeuronDfOverF 27 0.01105 20170127-1 NeuronDfOverF 24 0.01149 20170203-1 GliaDfOverF 46 0.00109 20170203-1 NeuronDfOverF 46 0.00161 20170228-3 GliaDfOverF 67 0.00970 20170228-3 NeuronDfOverF 67 0.01104 20170228-4 GliaDfOverF 108 0.00723 20170228-4 NeuronDfOverF 108 0.00690 NeuronDfOverF 20161109-2: onset-triggered mean in SD units ranges -0.33 to 0.26 NeuronDfOverF 20170113-2: onset-triggered mean in SD units ranges -0.25 to 2.85 NeuronDfOverF 20170113-3: onset-triggered mean in SD units ranges -1.12 to 0.25 NeuronDfOverF 20170113-4: onset-triggered mean in SD units ranges -1.62 to 0.25 NeuronDfOverF 20170127-1: onset-triggered mean in SD units ranges -1.30 to 0.96 NeuronDfOverF 20170203-1: onset-triggered mean in SD units ranges -0.46 to 1.07 NeuronDfOverF 20170228-3: onset-triggered mean in SD units ranges -1.05 to 1.62 NeuronDfOverF 20170228-4: onset-triggered mean in SD units ranges -1.80 to 0.25 GliaDfOverF 20161109-2: onset-triggered mean in SD units ranges -0.96 to 0.55 GliaDfOverF 20170203-1: onset-triggered mean in SD units ranges -1.07 to 0.72 GliaDfOverF 20170228-3: onset-triggered mean in SD units ranges -0.06 to 0.68 GliaDfOverF 20170228-4: onset-triggered mean in SD units ranges -0.24 to 0.88
fig2bc-neurons-lower-in-passivity¶
Whole-brain average neuronal activity is higher during active behavior than during passivity, and decays after passivity onset.
Criterion: per fish, mean neuronal signal over open-loop passive epochs minus that over open-loop active epochs; below 0 by a one-sided Wilcoxon over the 8 fish.
In [17]:
rows = []
for s in sessions:
if "NeuronDfOverF" not in s["imaging"]:
continue
f = s["fish"]
t, y, on, base, E = TRACES[(f, "NeuronDfOverF")]
pas = [p for x in EP[f]["open"] for p in x["passive"]]
act = [(e["start"], e["stop"]) for x in EP[f]["open"] for e in x["active"]]
diff = fv.mean_in_intervals(t, y, pas) - fv.mean_in_intervals(t, y, act)
late = fv.window_mean(t, y, on, (5, 10)) - base
good = np.isfinite(late)
rows.append({"fish": f, "passive_minus_active": diff, "in_trace_sd": diff / y.std(), "onsets": int(good.sum()),
"after_5to10_minus_base": np.nanmean(late), "p_decrease": stats.wilcoxon(late[good], alternative="less").pvalue})
neu = pd.DataFrame(rows)
print(neu.to_string(index=False, float_format=lambda v: f"{v:.3g}"))
w = stats.wilcoxon(neu.passive_minus_active, alternative="less")
print(f"\nPassive below active in {(neu.passive_minus_active < 0).sum()} of {len(neu)} fish; Wilcoxon p = {w.pvalue:.3g}; "
f"effect {neu.in_trace_sd.min():.2f} to {neu.in_trace_sd.max():.2f} SD of the trace")
print(f"Decrease 5 to 10 s after onset (p < 0.05 within fish) in {(neu.p_decrease < 0.05).sum()} of {len(neu)} fish")
# Variant: ActivityStates rows in place of the Methods' epochs.
v = []
for s in sessions:
if "NeuronDfOverF" not in s["imaging"]:
continue
t, y = TRACES[(s["fish"], "NeuronDfOverF")][:2]
st = s["states"]
v.append(fv.mean_in_intervals(t, y, st[st.type == "passive"][["start", "stop"]].to_numpy())
- fv.mean_in_intervals(t, y, st[st.type == "active"][["start", "stop"]].to_numpy()))
v = np.array(v)
print(f"Variant with ActivityStates rows: passive below active in {(v < 0).sum()} of {len(v)} fish, Wilcoxon p = {stats.wilcoxon(v, alternative='less').pvalue:.3g}")
fish passive_minus_active in_trace_sd onsets after_5to10_minus_base p_decrease 20161109-2 -0.00541 -0.642 26 -0.00163 0.0102 20170113-2 -0.0214 -1.47 64 -0.00294 0.0123 20170113-3 -0.0269 -1.5 23 -0.0187 5.13e-06 20170113-4 -0.0135 -1.22 21 -0.0159 3.34e-06 20170127-1 -0.013 -1.13 24 -0.0145 5.96e-08 20170203-1 -0.00127 -0.784 46 -0.000241 0.0737 20170228-3 -0.0216 -1.96 67 -0.0112 2.34e-12 20170228-4 -0.00892 -1.29 97 -0.0118 1.28e-17 Passive below active in 8 of 8 fish; Wilcoxon p = 0.00391; effect -1.96 to -0.64 SD of the trace Decrease 5 to 10 s after onset (p < 0.05 within fish) in 7 of 8 fish Variant with ActivityStates rows: passive below active in 8 of 8 fish, Wilcoxon p = 0.00391
Result: reproduced. The whole-brain neuronal signal is lower in open-loop passive epochs
than in active epochs in all 8 fish (Wilcoxon p = 0.00391), by 0.64 to 1.96 SD of each fish's
trace. It falls 5 to 10 s after passivity onset, against the baseline, with p < 0.05 in 7 of 8
fish. The ActivityStates rows give the same (8 of 8).
fig2c-glia-rise-at-passivity¶
Whole-brain average radial astrocyte calcium increases before passivity onset and peaks soon after it.
Criterion: per fish, 1 to 4 s after onset and 2 to 0 s before onset, each against 10 to 5 s before (one-tailed Wilcoxon over onsets), and the peak of the triggered average between 0 and 10 s; reproduced when all hold in at least 3 of 4 fish, partially when the after-onset increase holds in at least 3.
In [18]:
def onset_tests(f, name, onsets=None):
t, y, on, base, E = TRACES[(f, name)]
if onsets is not None:
on = onsets
base = fv.window_mean(t, y, on, fv.BASELINE)
E = fv.event_windows(t, y, on, LAGS) - base[:, None]
resp = fv.window_mean(t, y, on, fv.RESPONSE) - base
pre = fv.window_mean(t, y, on, (-2, 0)) - base
good = np.isfinite(resp) & np.isfinite(pre)
avg = np.nanmean(E, 0)
return {"fish": f, "onsets": int(good.sum()), "after_1to4": np.median(resp[good]) / y.std(),
"p_after": stats.wilcoxon(resp[good], alternative="greater").pvalue,
"before_2to0": np.median(pre[good]) / y.std(),
"p_before": stats.wilcoxon(pre[good], alternative="greater").pvalue,
"peak_lag_s": LAGS[np.nanargmax(avg)]}
glia = pd.DataFrame([onset_tests(s["fish"], "GliaDfOverF") for s in sessions if "GliaDfOverF" in s["imaging"]])
glia["after_ok"] = glia.p_after < 0.05
glia["all_ok"] = glia.after_ok & (glia.p_before < 0.05) & glia.peak_lag_s.between(0, 10)
print("Medians in SD units of each fish's trace:")
print(glia.to_string(index=False, float_format=lambda v: f"{v:.3g}"))
print(f"\nAfter-onset increase in {glia.after_ok.sum()} of 4 fish; all three parts in {glia.all_ok.sum()} of 4 fish")
print("Onsets per fish:", glia.set_index("fish").onsets.to_dict(), "(paper's example fish: 64 events)")
# Variant: ActivityStates passive rows that start in open loop, as onsets.
rows = []
for s in sessions:
if "GliaDfOverF" not in s["imaging"]:
continue
st = s["states"]
on = st[st.type == "passive"].start.to_numpy()
on = on[fv.label_at(s, on) == 2]
rows.append(onset_tests(s["fish"], "GliaDfOverF", on))
print("\nVariant with ActivityStates passive rows in open loop as onsets:")
print(pd.DataFrame(rows).to_string(index=False, float_format=lambda v: f"{v:.3g}"))
# Neurons in the same two-colour fish, for comparison.
print("\nNeurons in the two-colour fish, same tests:")
print(pd.DataFrame([onset_tests(s["fish"], "NeuronDfOverF") for s in sessions if "GliaDfOverF" in s["imaging"]]).to_string(index=False, float_format=lambda v: f"{v:.3g}"))
# Variant: onsets with the edge epochs kept, which adds the first passivity of each open-loop period.
gx = pd.DataFrame([onset_tests(s["fish"], "GliaDfOverF", ONX[s["fish"]]) for s in sessions if "GliaDfOverF" in s["imaging"]])
gx["all_ok"] = (gx.p_after < 0.05) & (gx.p_before < 0.05) & gx.peak_lag_s.between(0, 10)
print("\nVariant with edge epochs kept (onsets added per fish: "
f"{ {f: len(ONX[f]) - len(ON[f]) for f in gx.fish} }):")
print(gx.to_string(index=False, float_format=lambda v: f"{v:.3g}"))
print(f"After-onset increase in {(gx.p_after < 0.05).sum()} of 4 fish; all three parts in {gx.all_ok.sum()} of 4 fish")
# The first five onsets after the start of each open-loop period (edge epochs kept).
rows = []
for s in sessions:
if "GliaDfOverF" not in s["imaging"]:
continue
f = s["fish"]
early = np.concatenate([ONX[f][(ONX[f] >= a) & (ONX[f] < b)][:5] for a, b in [x["seg"] for x in EPX[f]["open"]]])
for name in ["GliaDfOverF", "NeuronDfOverF"]:
r = onset_tests(f, name, early)
r["series"] = name
rows.append(r)
print("\nFirst five onsets after the start of each open-loop period (edge epochs kept):")
print(pd.DataFrame(rows).to_string(index=False, float_format=lambda v: f"{v:.3g}"))
# A stretch of the whole-brain traces with behaviour, in the two fish whose glia do not rise.
fig, axes = plt.subplots(2, 1, figsize=(11, 5))
for ax, f in zip(axes, ["20161109-2", "20170203-1"]):
t0 = ON[f][len(ON[f]) // 2] - 150
for name, col in [("GliaDfOverF", "C3"), ("NeuronDfOverF", "C2")]:
t, y = TRACES[(f, name)][:2]
sel = (t >= t0) & (t <= t0 + 300)
z = (y[sel] - y.mean()) / y.std()
ax.plot(t[sel], z, color=col, lw=0.8, label=f"{name.replace('DfOverF', '')} (z, range {z.min():.2f} to {z.max():.2f})")
sw = S[f]["swims"]
sw = sw[(sw.start_time >= t0) & (sw.start_time <= t0 + 300)]
ax.plot(sw.start_time, np.full(len(sw), -3.2), "|k", ms=6)
lab = fv.label_at(S[f], np.arange(t0, t0 + 300, 0.1))
ax.fill_between(np.arange(t0, t0 + 300, 0.1), -3.6, -3.4, where=lab == 2, color="0.6", step="mid")
for o in ON[f][(ON[f] >= t0) & (ON[f] <= t0 + 300)]:
ax.axvline(o, color="k", lw=0.4, ls=":")
ax.set_ylim(-3.8, 4); ax.set_title(f"{f}: bouts (ticks), open loop (grey bar), passivity onsets (dotted)", fontsize=8)
ax.legend(fontsize=7, loc="upper right")
print(f"{f}: {len(sw)} bouts and {((ON[f] >= t0) & (ON[f] <= t0 + 300)).sum()} onsets in the 300 s shown from {t0:.0f} s; "
f"glia trace SD {TRACES[(f, 'GliaDfOverF')][1].std():.4f}, neurons {TRACES[(f, 'NeuronDfOverF')][1].std():.4f}")
axes[-1].set_xlabel("time (s)")
plt.tight_layout(); plt.show()
for f in sorted({s["fish"] for s in sessions if "GliaDfOverF" in s["imaging"]}):
print(f"{f}: glia trace SD {TRACES[(f, 'GliaDfOverF')][1].std():.4f}, neuron trace SD {TRACES[(f, 'NeuronDfOverF')][1].std():.4f}, "
f"correlation of glia and neurons {np.corrcoef(*[np.interp(TRACES[(f, 'GliaDfOverF')][0], TRACES[(f, n)][0], TRACES[(f, n)][1]) for n in ['GliaDfOverF', 'NeuronDfOverF']])[0, 1]:.2f}")
# Passivity onsets by the protocol code at onset, in the two-colour fish (main rule).
rows = []
for s in sessions:
if "GliaDfOverF" not in s["imaging"]:
continue
f = s["fish"]
i = np.clip(np.searchsorted(s["cond_t"], ON[f]), 0, len(s["cond_t"]) - 1)
rows.append({"fish": f, "onsets": len(ON[f]), **pd.Series(s["protocol_code"][i].astype(int)).value_counts().sort_index()
.rename(lambda c: f"code_{c}").to_dict()})
print("\nPassivity onsets by protocol code at onset:")
print(pd.DataFrame(rows).fillna(0).convert_dtypes().to_string(index=False))
Medians in SD units of each fish's trace:
fish onsets after_1to4 p_after before_2to0 p_before peak_lag_s after_ok all_ok
20161109-2 26 -0.906 1 -0.117 0.48 0 False False
20170203-1 46 -0.639 1 0.617 0.000551 -4 False False
20170228-3 67 0.612 8.03e-13 0.405 1.7e-11 3.5 True True
20170228-4 98 0.525 1.1e-14 0.486 3.61e-17 2.5 True True
After-onset increase in 2 of 4 fish; all three parts in 2 of 4 fish
Onsets per fish: {'20161109-2': 26, '20170203-1': 46, '20170228-3': 67, '20170228-4': 98} (paper's example fish: 64 events)
Variant with ActivityStates passive rows in open loop as onsets:
fish onsets after_1to4 p_after before_2to0 p_before peak_lag_s
20161109-2 35 -0.772 1 0.0718 0.294 0
20170203-1 186 -0.23 1 0.509 5.41e-13 0
20170228-3 67 0.523 3.34e-11 0.357 3.02e-09 5
20170228-4 226 0.495 5.86e-19 0.429 1.04e-25 2.5
Neurons in the two-colour fish, same tests:
fish onsets after_1to4 p_after before_2to0 p_before peak_lag_s
20161109-2 26 -0.179 0.986 -0.0118 0.146 0
20170203-1 46 0.277 0.0023 0.856 2.72e-10 0
20170228-3 67 -0.354 1 1.32 2.08e-09 -3.5
20170228-4 98 -1.51 1 -0.14 0.428 -3
Variant with edge epochs kept (onsets added per fish: {'20161109-2': 2, '20170203-1': 2, '20170228-3': 1, '20170228-4': 1}):
fish onsets after_1to4 p_after before_2to0 p_before peak_lag_s all_ok
20161109-2 28 -0.928 1 -0.117 0.496 0 False
20170203-1 48 -0.639 1 0.504 0.00081 -4 False
20170228-3 68 0.619 5.45e-13 0.417 1.11e-11 3 True
20170228-4 99 0.525 7.01e-15 0.487 2.42e-17 2.5 True
After-onset increase in 2 of 4 fish; all three parts in 2 of 4 fish
First five onsets after the start of each open-loop period (edge epochs kept):
fish onsets after_1to4 p_after before_2to0 p_before peak_lag_s series
20161109-2 10 -1.13 0.998 0.0268 0.278 0 GliaDfOverF
20161109-2 10 -0.0286 0.722 0.00822 0.278 0 NeuronDfOverF
20170203-1 11 -0.624 0.993 1.02 0.00488 -3 GliaDfOverF
20170203-1 11 0.481 0.00342 1.24 0.000488 0 NeuronDfOverF
20170228-3 8 0.826 0.0117 0.467 0.0977 2 GliaDfOverF
20170228-3 8 -0.0842 0.809 0.645 0.23 -3.5 NeuronDfOverF
20170228-4 15 1.36 0.00336 0.708 0.00058 3 GliaDfOverF
20170228-4 15 -1.77 0.999 -0.388 0.262 -2 NeuronDfOverF
20161109-2: 197 bouts and 7 onsets in the 300 s shown from 3427 s; glia trace SD 0.0016, neurons 0.0084
20170203-1: 129 bouts and 10 onsets in the 300 s shown from 1307 s; glia trace SD 0.0011, neurons 0.0016
20161109-2: glia trace SD 0.0016, neuron trace SD 0.0084, correlation of glia and neurons 0.58
20170203-1: glia trace SD 0.0011, neuron trace SD 0.0016, correlation of glia and neurons 0.59
20170228-3: glia trace SD 0.0097, neuron trace SD 0.0110, correlation of glia and neurons 0.06
20170228-4: glia trace SD 0.0072, neuron trace SD 0.0069, correlation of glia and neurons 0.15
Passivity onsets by protocol code at onset:
fish onsets code_32 code_53 code_51
20161109-2 26 26 0 0
20170203-1 46 42 4 0
20170228-3 67 0 33 34
20170228-4 108 0 65 43
Result: not reproduced. The whole-brain glial signal rises before passivity onset and peaks after it in 2 of the 4 two-colour fish (20170228-3 and 20170228-4: 1 to 4 s after onset 0.61 and 0.53 SD above baseline, p = 8.0e-13 and 1.1e-14; 2 to 0 s before, p = 1.7e-11 and 3.6e-17; peaks 3.5 and 2.5 s after onset). In the other two it falls after onset (-0.91 and -0.64 SD, p = 1 for an increase). The criterion needed 3 of 4. One of the matching fish has 67 onsets, close to the 64 events of the paper's example fish. All 26 onsets of 20161109-2 fall in code 32, prolonged open loop with the grating moving, and 42 of the 46 onsets of 20170203-1, so the short alternating blocks do not explain the fall in these two fish.
The split does not depend on the onsets the main rule leaves out: with the edge epochs kept, 1 or 2 onsets are added per fish and the same 2 fish pass. Over the first five onsets after the start of each open-loop period, the glial response 1 to 4 s after onset is again negative in 20161109-2 and 20170203-1 (-1.13 and -0.62 SD) and positive in the other two (0.83 and 1.36 SD). In 20161109-2 and 20170203-1 the glial trace has an SD of 0.0016 and 0.0011, against 0.0097 and 0.0072, and it correlates with the neuronal trace at 0.58 and 0.59, against 0.06 and 0.15 in the other two fish. In 20170203-1 the neuronal trace also has a low SD (0.0016, against 0.0069 to 0.0180 in the other fish) and, unlike the other three two-colour fish, rises after onset (+0.28 SD, p = 0.0023), so the low SD there affects both channels and does not point to the glial channel alone. The data here do not show why these two files differ.
figs2a-open-loop-higher¶
Whole-brain average neuronal and glial calcium are higher and fluctuate more strongly in open loop than in closed loop.
Criterion: per fish, mean and SD of the whole-brain signal in open and closed loop, leaving out 30 s after each switch; neurons by one-sided Wilcoxon over 8 fish, glia in all 4 fish.
In [19]:
rows = []
for (f, name), (t, y, on, base, E) in TRACES.items():
s = S[f]
lab = fv.label_at(s, t)
keep = fv.seconds_since_switch(s, t) >= 30
r = {"fish": f, "series": name}
for c, code_ in [("closed", 1), ("open", 2)]:
sel = keep & (lab == code_)
r[f"{c}_frames"] = int(sel.sum())
r[f"{c}_mean"] = (y[sel].mean() - y.mean()) / y.std()
r[f"{c}_sd"] = y[sel].std() / y.std()
rows.append(r)
s2a = pd.DataFrame(rows)
s2a["mean_higher_open"] = s2a.open_mean > s2a.closed_mean
s2a["sd_higher_open"] = s2a.open_sd > s2a.closed_sd
print("Means and SDs in units of each fish's whole-trace SD:")
print(s2a.round(3).to_string(index=False))
for name in ["NeuronDfOverF", "GliaDfOverF"]:
d = s2a[s2a.series == name]
d = d.dropna(subset=["closed_mean"])
pm = stats.wilcoxon(d.open_mean, d.closed_mean, alternative="greater").pvalue
ps = stats.wilcoxon(d.open_sd, d.closed_sd, alternative="greater").pvalue
print(f"{name}: mean higher in open loop in {d.mean_higher_open.sum()} of {len(d)} fish (Wilcoxon p = {pm:.3g}, "
f"median difference {np.median(d.open_mean - d.closed_mean):.2f} SD); SD higher in {d.sd_higher_open.sum()} of {len(d)} "
f"(p = {ps:.3g}, median ratio {np.median(d.open_sd / d.closed_sd):.2f})")
# Variant for the SD: computed within each continuous closed- or open-loop period (after the first
# 30 s), raw or after removing a linear trend, then averaged over periods weighted by frames.
rows = []
for (f, name), (t, y, on, base, E) in TRACES.items():
s = S[f]
lab = fv.label_at(s, t)
keep = fv.seconds_since_switch(s, t) >= 30
r = {"fish": f, "series": name}
for c, code_ in [("closed", 1), ("open", 2)]:
sel = keep & (lab == code_)
runs = np.split(np.flatnonzero(sel), np.flatnonzero(np.diff(np.flatnonzero(sel)) > 1) + 1)
runs = [i for i in runs if len(i) >= 10]
n = np.array([len(i) for i in runs])
sd_raw = np.array([y[i].std() for i in runs]) / y.std()
sd_det = np.array([(y[i] - np.polyval(np.polyfit(t[i], y[i], 1), t[i])).std() for i in runs]) / y.std()
r[f"{c}_periods"] = len(runs)
r[f"{c}_sd_within"] = np.average(sd_raw, weights=n) if len(runs) else np.nan
r[f"{c}_sd_detrended"] = np.average(sd_det, weights=n) if len(runs) else np.nan
rows.append(r)
s2v = pd.DataFrame(rows)
print("\nSD within periods (periods of 10 frames or more), in units of each fish's whole-trace SD:")
print(s2v.round(3).to_string(index=False))
for name in ["NeuronDfOverF", "GliaDfOverF"]:
d = s2v[s2v.series == name].dropna()
for k in ["sd_within", "sd_detrended"]:
p = stats.wilcoxon(d[f"open_{k}"], d[f"closed_{k}"], alternative="greater").pvalue
print(f"{name} {k}: higher in open loop in {(d[f'open_{k}'] > d[f'closed_{k}']).sum()} of {len(d)} fish "
f"(Wilcoxon p = {p:.3g}, median ratio {np.median(d[f'open_{k}'] / d[f'closed_{k}']):.2f})")
# Variant for the SD: windows of equal length (10 s, after the first 30 s of each period), the SD
# within each window, averaged with equal weight over the windows of each condition. This compares
# the same time scale in both conditions, whatever the length of the periods.
WIN_S = 10.0
rows = []
for (f, name), (t, y, on, base, E) in TRACES.items():
s = S[f]
lab = fv.label_at(s, t)
keep = fv.seconds_since_switch(s, t) >= 30
n_win = int(round(WIN_S / np.median(np.diff(t))))
r = {"fish": f, "series": name, "frames_per_window": n_win}
for c, code_ in [("closed", 1), ("open", 2)]:
sel = keep & (lab == code_)
runs = np.split(np.flatnonzero(sel), np.flatnonzero(np.diff(np.flatnonzero(sel)) > 1) + 1)
sds = [y[i[k:k + n_win]].std() / y.std() for i in runs for k in range(0, len(i) - n_win + 1, n_win)]
r[f"{c}_windows"] = len(sds)
r[f"{c}_sd_10s"] = np.mean(sds) if sds else np.nan
rows.append(r)
s2w = pd.DataFrame(rows)
print("\nSD within 10 s windows, in units of each fish's whole-trace SD:")
print(s2w.round(3).to_string(index=False))
for name in ["NeuronDfOverF", "GliaDfOverF"]:
d = s2w[s2w.series == name].dropna()
p = stats.wilcoxon(d.open_sd_10s, d.closed_sd_10s, alternative="greater").pvalue
print(f"{name} SD in 10 s windows: higher in open loop in {(d.open_sd_10s > d.closed_sd_10s).sum()} of {len(d)} fish "
f"(Wilcoxon p = {p:.3g}, median ratio {np.median(d.open_sd_10s / d.closed_sd_10s):.2f})")
Means and SDs in units of each fish's whole-trace SD:
fish series closed_frames closed_mean closed_sd open_frames open_mean open_sd mean_higher_open sd_higher_open
20161109-2 GliaDfOverF 1643 -0.265 0.842 9262 0.106 0.954 True True
20161109-2 NeuronDfOverF 1643 -0.403 0.393 9262 0.057 0.889 True True
20170113-2 NeuronDfOverF 746 -0.101 0.291 5557 0.050 1.075 True True
20170113-3 NeuronDfOverF 1290 -0.694 0.244 5473 0.242 1.024 True True
20170113-4 NeuronDfOverF 2857 -0.494 0.608 4216 0.261 0.934 True True
20170127-1 NeuronDfOverF 4542 0.102 0.798 3853 0.072 1.176 False True
20170203-1 GliaDfOverF 1244 0.092 1.051 4106 -0.004 0.933 False False
20170203-1 NeuronDfOverF 1244 -0.073 1.135 4106 -0.003 0.900 True False
20170228-3 GliaDfOverF 213 -1.190 0.117 5477 0.028 0.670 True True
20170228-3 NeuronDfOverF 213 0.011 0.251 5477 0.000 1.025 False True
20170228-4 GliaDfOverF 356 -0.965 0.181 6085 0.122 0.907 True True
20170228-4 NeuronDfOverF 356 -0.141 0.266 6085 0.048 1.023 True True
NeuronDfOverF: mean higher in open loop in 6 of 8 fish (Wilcoxon p = 0.0195, median difference 0.17 SD); SD higher in 7 of 8 (p = 0.00781, median ratio 2.98)
GliaDfOverF: mean higher in open loop in 3 of 4 fish (Wilcoxon p = 0.125, median difference 0.73 SD); SD higher in 3 of 4 (p = 0.188, median ratio 3.06)
SD within periods (periods of 10 frames or more), in units of each fish's whole-trace SD:
fish series closed_periods closed_sd_within closed_sd_detrended open_periods open_sd_within open_sd_detrended
20161109-2 GliaDfOverF 13 0.563 0.494 26 0.842 0.822
20161109-2 NeuronDfOverF 13 0.202 0.094 26 0.717 0.694
20170113-2 NeuronDfOverF 1 0.273 0.262 9 1.065 1.043
20170113-3 NeuronDfOverF 18 0.127 0.122 18 0.881 0.743
20170113-4 NeuronDfOverF 30 0.300 0.245 30 0.844 0.765
20170127-1 NeuronDfOverF 65 0.551 0.348 2 0.836 0.827
20170203-1 GliaDfOverF 67 0.703 0.451 16 0.907 0.900
20170203-1 NeuronDfOverF 67 0.354 0.205 16 0.868 0.861
20170228-3 GliaDfOverF 10 0.065 0.054 2 0.557 0.548
20170228-3 NeuronDfOverF 10 0.219 0.127 2 1.025 1.024
20170228-4 GliaDfOverF 22 0.057 0.040 3 0.906 0.886
20170228-4 NeuronDfOverF 22 0.234 0.167 3 1.020 1.018
NeuronDfOverF sd_within: higher in open loop in 8 of 8 fish (Wilcoxon p = 0.00391, median ratio 3.72)
NeuronDfOverF sd_detrended: higher in open loop in 8 of 8 fish (Wilcoxon p = 0.00391, median ratio 5.15)
GliaDfOverF sd_within: higher in open loop in 4 of 4 fish (Wilcoxon p = 0.0625, median ratio 5.04)
GliaDfOverF sd_detrended: higher in open loop in 4 of 4 fish (Wilcoxon p = 0.0625, median ratio 6.11)
SD within 10 s windows, in units of each fish's whole-trace SD:
fish series frames_per_window closed_windows closed_sd_10s open_windows open_sd_10s
20161109-2 GliaDfOverF 29 51 0.383 317 0.463
20161109-2 NeuronDfOverF 29 51 0.074 317 0.136
20170113-2 NeuronDfOverF 27 27 0.154 201 0.625
20170113-3 NeuronDfOverF 27 36 0.090 196 0.339
20170113-4 NeuronDfOverF 27 90 0.177 140 0.371
20170127-1 NeuronDfOverF 27 115 0.349 142 0.408
20170203-1 GliaDfOverF 25 12 0.940 149 0.593
20170203-1 NeuronDfOverF 25 12 0.510 149 0.369
20170228-3 GliaDfOverF 17 10 0.067 321 0.196
20170228-3 NeuronDfOverF 17 10 0.212 321 0.627
20170228-4 GliaDfOverF 16 9 0.061 379 0.239
20170228-4 NeuronDfOverF 16 9 0.203 379 0.611
NeuronDfOverF SD in 10 s windows: higher in open loop in 7 of 8 fish (Wilcoxon p = 0.0195, median ratio 2.52)
GliaDfOverF SD in 10 s windows: higher in open loop in 3 of 4 fish (Wilcoxon p = 0.438, median ratio 2.07)
Result: partially reproduced. For neurons, the whole-brain signal is higher in open than closed loop in 6 of 8 fish (Wilcoxon p = 0.0195) and fluctuates more in 7 of 8 (p = 0.0078, median SD ratio 2.98). For glia, both hold in 3 of the 4 fish, not in all 4 as the criterion required; 20170203-1 is the exception. With the SD taken within each continuous period, raw or after removing a linear trend (a variant), it is higher in open loop in 8 of 8 neuron fish (p = 0.0039, median ratios 3.72 and 5.15) and in 4 of 4 glia fish (median ratios 5.04 and 6.11). These periods differ in length (in 20170228-3, 10 closed-loop periods over 213 frames against 2 open-loop periods over 5477), and a long period holds slow fluctuations that a short one cannot. With the SD taken in windows of 10 s in both conditions (a variant), it is higher in open loop in 7 of 8 neuron fish (p = 0.0195, median ratio 2.52) and 3 of 4 glia fish (median ratio 2.07), with 20170203-1 again the exception, where the closed loop gives 12 windows.
fig2h-glia-weaker-closed-loop-pauses¶
When fish in closed loop stop swimming because the stimulus stops moving, glial calcium signals are weaker than during futility-induced passivity.
Criterion: closed-loop pauses (5 s or more, starting while the gain is above 0) against open-loop passivity onsets, response 1 to 4 s after minus 10 to 5 s before; one-sided Mann-Whitney per fish with 3 or more pauses.
In [20]:
rows = []
for s in sessions:
if "GliaDfOverF" not in s["imaging"]:
continue
f = s["fish"]
t, y, on, base, E = TRACES[(f, "GliaDfOverF")]
r_open = fv.window_mean(t, y, on, fv.RESPONSE) - base
cp = fv.closed_loop_pause_onsets(s)
r_cl = fv.window_mean(t, y, cp, fv.RESPONSE) - fv.window_mean(t, y, cp, fv.BASELINE)
r_open, r_cl = r_open[np.isfinite(r_open)], r_cl[np.isfinite(r_cl)]
stationary = np.mean(s["velocity"][np.clip(np.searchsorted(s["cond_t"], cp), 0, len(s["cond_t"]) - 1)] == 0) if len(cp) else np.nan
rows.append({"fish": f, "closed_pauses": len(r_cl), "share_in_stationary": stationary, "open_onsets": len(r_open),
"closed_median": np.median(r_cl) / y.std() if len(r_cl) else np.nan, "open_median": np.median(r_open) / y.std(),
"p": stats.mannwhitneyu(r_cl, r_open, alternative="less").pvalue if len(r_cl) >= 3 else np.nan})
h = pd.DataFrame(rows)
h["smaller_median"] = h.closed_median < h.open_median
print("Medians in units of each fish's trace SD:")
print(h.to_string(index=False, float_format=lambda v: f"{v:.3g}"))
kept = h[h.closed_pauses >= 3]
print(f"\n{len(kept)} fish with 3 or more closed-loop pauses; pause response smaller with p < 0.05 in {(kept.p < 0.05).sum()}, "
f"smaller in median in {kept.smaller_median.sum()}")
# Variants: (1) the gain above 0 over the whole window, 10 s before to 4 s after the pause onset;
# (2) open-loop onsets with the edge epochs kept; (3) both.
rows = []
for s in sessions:
if "GliaDfOverF" not in s["imaging"]:
continue
f = s["fish"]
t, y = TRACES[(f, "GliaDfOverF")][:2]
for label, whole, onsets in [("whole window", True, ON[f]), ("edge epochs kept", False, ONX[f]), ("both", True, ONX[f])]:
r_open = fv.window_mean(t, y, onsets, fv.RESPONSE) - fv.window_mean(t, y, onsets, fv.BASELINE)
cp = fv.closed_loop_pause_onsets(s, whole_window=whole)
r_cl = fv.window_mean(t, y, cp, fv.RESPONSE) - fv.window_mean(t, y, cp, fv.BASELINE)
r_open, r_cl = r_open[np.isfinite(r_open)], r_cl[np.isfinite(r_cl)]
rows.append({"variant": label, "fish": f, "closed_pauses": len(r_cl), "open_onsets": len(r_open),
"closed_median": np.median(r_cl) / y.std() if len(r_cl) else np.nan, "open_median": np.median(r_open) / y.std(),
"p": stats.mannwhitneyu(r_cl, r_open, alternative="less").pvalue if len(r_cl) >= 3 else np.nan})
hv = pd.DataFrame(rows)
print(hv.to_string(index=False, float_format=lambda v: f"{v:.3g}"))
for label, d in hv.groupby("variant", sort=False):
d = d[d.closed_pauses >= 3]
print(f"{label}: {len(d)} fish with 3 or more pauses; smaller with p < 0.05 in {(d.p < 0.05).sum()}, smaller in median in {(d.closed_median < d.open_median).sum()}")
Medians in units of each fish's trace SD:
fish closed_pauses share_in_stationary open_onsets closed_median open_median p smaller_median
20161109-2 26 0.267 26 -0.134 -0.906 0.782 False
20170203-1 43 0.256 46 -0.0894 -0.639 0.982 False
20170228-3 10 0.538 67 -0.093 0.612 1.07e-05 True
20170228-4 30 0.424 98 -0.000607 0.525 5.03e-11 True
4 fish with 3 or more closed-loop pauses; pause response smaller with p < 0.05 in 2, smaller in median in 2
variant fish closed_pauses open_onsets closed_median open_median p
whole window 20161109-2 25 26 -0.273 -0.906 0.717
edge epochs kept 20161109-2 26 28 -0.134 -0.928 0.836
both 20161109-2 25 28 -0.273 -0.928 0.781
whole window 20170203-1 34 46 -0.205 -0.639 0.914
edge epochs kept 20170203-1 43 48 -0.0894 -0.639 0.986
both 20170203-1 34 48 -0.205 -0.639 0.926
whole window 20170228-3 8 67 -0.093 0.612 3.97e-09
edge epochs kept 20170228-3 10 68 -0.093 0.619 9.9e-06
both 20170228-3 8 68 -0.093 0.619 3.55e-09
whole window 20170228-4 28 98 0.00987 0.525 2.73e-10
edge epochs kept 20170228-4 30 99 -0.000607 0.525 4.23e-11
both 20170228-4 28 99 0.00987 0.525 2.31e-10
whole window: 4 fish with 3 or more pauses; smaller with p < 0.05 in 2, smaller in median in 2
edge epochs kept: 4 fish with 3 or more pauses; smaller with p < 0.05 in 2, smaller in median in 2
both: 4 fish with 3 or more pauses; smaller with p < 0.05 in 2, smaller in median in 2
Result: not reproduced. The glial response to closed-loop pauses is smaller than to open-loop passivity onsets in 2 of 4 fish (20170228-3 and 20170228-4, Mann-Whitney p = 1.1e-05 and 5.0e-11), the two fish whose glia rise at passivity. In the other two, glia fall after passivity onset, and the pause response is not smaller (p = 0.78 and 0.98). Requiring the gain above 0 over the whole window from 10 s before to 4 s after the pause, keeping the edge epochs, or both (variants) give the same 2 of 4.
Summary¶
A failure to reproduce here may be a failure of this attempt, not of the paper: see the reasoning and the differences from the paper's method.
| Claim | Data | Outcome |
|---|---|---|
| fig1h-epoch-durations | dandiset, 12 fish | partially reproduced |
| methods-isi | dandiset, 12 fish | partially reproduced |
| fig1j-open-loop-passivity | dandiset, 12 fish | not reproduced |
| fig1c-vigor-open-loop-onset | dandiset, 9 fish | reproduced |
| fig1i-isi-before-passivity | dandiset, 11 fish | not reproduced |
| fig2bc-neurons-lower-in-passivity | dandiset, 8 fish | reproduced |
| fig2c-glia-rise-at-passivity | dandiset, 4 fish | not reproduced |
| figs2a-open-loop-higher | dandiset, 8 and 4 fish | partially reproduced |
| fig2h-glia-weaker-closed-loop-pauses | dandiset, 4 fish | not reproduced |
| fig2dg-nmf-components | blocked: methods underspecified | |
| fig2kl-lmo-glia-most-consistent | blocked: metadata insufficient | |
| fig1k-passive-unresponsive | blocked: metadata insufficient | |
| fig1f-replay | blocked: metadata insufficient | |
| fig4-glia-necessary | blocked: data missing | |
| fig5-glia-sufficient | blocked: data missing | |
| fig6a-ne-ramp | blocked: metadata insufficient | |
| fig6-ne-drives-glia | blocked: data missing | |
| fig7-ne-mismatch | blocked: metadata insufficient | |
| fig7-glia-accumulate | blocked: metadata insufficient |
How the analysis differs from the paper's.
- The behavioural claims use the 12 imaged fish; the paper's behavioural numbers come from 74 (61 for Figure 1I). The Methods' epoch rule is applied to the dandiset's swim bouts, with pauses measured from the end of one bout to the start of the next; the dandiset's own state table is shown as a variant. Epochs cut by the edges of a condition period are left out, which drops the first passivity onset of some open-loop periods (0 to 8 per fish); a variant keeps them for fig1h, fig1i, fig2c and fig2h and changes none of these outcomes. It was not run for methods-isi, fig1j or fig2bc.
- Closed and open loop are read from the gain signal, and stationary-grating time is left out; the protocol codes in the files have no description.
- The cell traces are used as stored. They are not ΔF/F as the Methods define it, and the two groups of files were processed differently, so every imaging comparison is made within a fish and scaled by the SD of that fish's trace.
- The whole-brain signal is the mean over about 8 to 15% of ROIs (8 evenly spaced blocks); two disjoint subsets agree (correlation 0.981 and 0.998).
- The paper's glial claims rest on NMF components clustered across fish and registered to a template, mostly in a lateral hindbrain region (L-MO); here only the whole-brain average is tested, as the paper shows for its two-colour example fish (Figure 2C).
- Glia were tested in the 4 two-colour fish, whose glia carry jRGECO1b; the paper's NMF analysis used 6 fish with glial GCaMP6f, of which the dandiset holds 3 without processed traces.
What was found about the dandiset.
- It holds the 12 imaging sessions of the paper's Figure 2 and related panels, with behaviour, but none of the perturbation or NE imaging experiments of Figures 4 to 6. Undescribed protocols (randomised gain, brief grating speed-ups) may belong to Figures 1K and 7 but cannot be identified.
- Four files have raw images only; three of them are the only glia-only GCaMP6f recordings.
- The glial plane of the two-colour files is labelled GCaMP6f but excited at 561 nm and detected
at 660 nm (jRGECO1b). The raw images are typed
TwoPhotonSeriesfrom a light-sheet microscope. - The
DfOverFseries are not ΔF/F (mean-normalised in 4 files, zero-centred and partly censored in the other 4), and theFluorescenceseries look like baselines, with the same description. - Imaging is missing for part of three sessions (20170113-2: 2135 of 9751 triggered volumes; raw images of 20170113-4 and 20170228-4).
- The
ActivityStatestable is undocumented, does not follow the Methods' definition (active rows of 5 s or less, transient rows duplicating active rows) and has 3 rows in the glia-only files. - The protocol codes, gain and velocity signals carry no units or descriptions. Code 33 is a 0.3 to 0.5 s speed-up of the grating in open loop, code 52 a randomised gain, and codes 51 and 53 do not consistently mark closed and open loop.
TotalStimulusVelocitySignal, described as the drift speed minus the gain times swim power, stays between -0.002 and 0.008, so the visual feedback of each bout is not recorded.- In two of the four two-colour files (20161109-2, 20170203-1) the glial whole-brain trace has a SD of 0.0016 and 0.0011 (0.0097 and 0.0072 in the others) and correlates 0.58 and 0.59 with the neuronal trace (0.06 and 0.15 in the others). In 20170203-1 the neuronal trace also has a low SD (0.0016).