Reproducing Brezovec et al. (2022) from dandiset 000727¶
Paper. Brezovec LE, Berger AB, Druckmann S, Clandinin TR (2022). Mapping the neural dynamics of locomotion across the Drosophila brain. bioRxiv. doi:10.1101/2022.03.20.485047. This attempt used the full text of the preprint, version 1 (21 March 2022), which is open access (CC BY-NC-ND) and the only version on bioRxiv. The paper later appeared in Current Biology (doi:10.1016/j.cub.2023.12.063), which is not open access and was not used.
Dandiset. DANDI:000727, version 0.240106.0043. The paper cites no DANDI version (its data statement points to Dryad and gives the DOIs in a key resources table that is not in the preprint). This is the only published version. The draft holds the same 18 assets (same blobs and sizes); it differs only in metadata added later (keywords, contributor details, and links to the conversion code and an example notebook), none of which the analysis uses.
What the dandiset holds. One NWB file per fly for nine flies, and for each a behaviour video (AVI) that the analysis does not use. Each NWB file holds raw GCaMP6f and tdTomato volumes (256 x 128 x 49 voxels, 3384 volumes, not motion corrected), the GCaMP6f volumes after motion correction, high-pass filtering and z-scoring, warped into the Functional Drosophila Atlas (FDA; 314 x 146 x 91 voxels), an anatomical scan, and the FicTrac output of the ball at 50 Hz. It does not hold the paper's mean brain, a neuropil atlas for any of these spaces, the supervoxels, the per-plane acquisition times, or the motion-corrected tdTomato channel.
How the data are used. The analysis runs on the processed (FDA) volumes, since motion
correction and registration with ANTs cannot be repeated here. Twelve of the 91 FDA planes, every
seventh from plane 10, are streamed for all nine flies (22 GB compressed); the files total
610 GB.
Two raw planes of both channels are also streamed for a check on the tdTomato channel. The
paper's supervoxels are replaced by 6 x 6 voxel blocks within each plane. Behaviour follows the
paper's code (ClandininLab/brezovec_volition_2023, commit 8e2b275). No data from outside the
dandiset are used.
| Claim | Tested |
|---|---|
| dataset-structure: nine flies, 30,456 volumes at 1.8 Hz, 1.6M voxels, 50 Hz tracking | yes |
| fig2c-fraction-correlated: 39% of the brain correlates with velocity | yes |
| fig2c-venn-categories: 15% forward only, 20% forward and turning, 4% turning only | yes |
| fig2-symmetry: forward symmetric, turning anti-symmetric | yes |
| fig2e-region-fractions, fig3-topographic-maps | blocked: no neuropil labels |
| fig4bcd-pc-structure: PC1 state, PC2 forward, PC3 turn sign | yes |
| fig4fg-ridge-r2: R2 about 0.3 with 100 PCs | yes |
| fig4h-r2-vs-movement: flies that move more are better predicted | yes |
| fig4i-weights-stereotyped | yes |
| fig5a-unique-variance: forward uniquely encoded, walking state not | yes |
| fig5c-global-below-threshold | yes |
| fig5f-ipsi-contra: ipsilateral increase, contralateral decrease | yes |
| fig6-peak-time-range: peaks from 300 ms before to over 1 s after | yes |
| fig6h-peak-vs-width | yes |
| fig6fg-ipsi-leads-contra | yes |
| fig6eg-region-sequence, fig7-within-region-gradients | blocked: no neuropil labels |
| figs2-tdtomato-control | blocked: no motion-corrected tdTomato (a check on the raw channels is shown) |
| figs4a-supervoxel-number | blocked: too large |
import time
import warnings
from concurrent.futures import ThreadPoolExecutor
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from dandi.dandiapi import DandiAPIClient
from scipy import signal, stats
import brainwide as bw
pd.set_option("display.width", 200)
pd.set_option("display.max_columns", 40)
pd.set_option("display.max_rows", 200)
warnings.filterwarnings("ignore", category=RuntimeWarning)
started = time.time()
sessions = bw.load_sessions()
flies = [s["subject_id"] for s in sessions]
by_fly = dict(zip(flies, sessions))
times = {f: bw.neural_times(s) for f, s in by_fly.items()}
beh = {f: bw.behaviour_at(s, times[f]) for f, s in by_fly.items()}
beh_start = {f: bw.behaviour_at(s, times[f] - 0.5 / s["processed"]["rate"]) for f, s in by_fly.items()}
# Null reference: each fly's behaviour rolled circularly by 800 to 2584 volumes (one shift per fly)
rng_roll = np.random.default_rng(2)
beh_rolled = {}
for f in flies:
shift = int(rng_roll.integers(800, 2584))
beh_rolled[f] = {k: np.roll(beh[f][k], shift) for k in ("fwd", "left", "right")}
print(f"metadata, FicTrac and chunk index of {len(sessions)} files in {time.time() - started:.0f} s")
planes = {}
with ThreadPoolExecutor(bw.N_THREADS) as pool:
for f, s in by_fly.items():
regressors = np.array([beh[f]["fwd"], beh[f]["left"], beh[f]["right"],
beh_start[f]["fwd"], beh_start[f]["left"], beh_start[f]["right"],
beh_rolled[f]["fwd"], beh_rolled[f]["left"], beh_rolled[f]["right"]])
planes[f] = bw.load_fly(s, regressors, pool)
print(f"streamed {len(bw.SLICES)} FDA planes of {len(flies)} flies from dandiset {bw.DANDISET_ID} "
f"version {bw.DANDISET_VERSION} in {time.time() - started:.0f} s")
t_raw = time.time()
raw_r = {}
with ThreadPoolExecutor(bw.N_THREADS) as pool:
for f, s in by_fly.items():
raw_r[f] = bw.raw_correlations(s, np.array([beh[f]["fwd"], beh[f]["left"], beh[f]["right"]]), pool)
print(f"streamed raw GCaMP6f and tdTomato planes {bw.RAW_SLICES} of {len(flies)} flies in {time.time() - t_raw:.0f} s")
metadata, FicTrac and chunk index of 9 files in 52 s
streamed 12 FDA planes of 9 flies from dandiset 000727 version 0.240106.0043 in 347 s
streamed raw GCaMP6f and tdTomato planes (18, 30) of 9 flies in 230 s
Data check¶
Each check of the metadata against the values is stated under its output.
with DandiAPIClient() as client:
meta = client.get_dandiset(bw.DANDISET_ID, bw.DANDISET_VERSION).get_raw_metadata()
sizes = {a.path: a.size for a in client.get_dandiset(bw.DANDISET_ID, bw.DANDISET_VERSION).get_assets()}
summary = meta["assetsSummary"]
print("Assets by type:", pd.Series([p.rsplit(".", 1)[-1] for p in sizes]).value_counts().to_dict())
print("Behaviour videos:", sorted(p for p in sizes if not p.endswith(".nwb"))[:2], "...")
print("Archive summary: subjects", summary["numberOfSubjects"], "; files", summary["numberOfFiles"],
"; species", [s["name"] for s in summary["species"]])
info = pd.DataFrame([{
"fly": f, "folder": s["path"].split("/")[0], "session_id": s["session_id"], "start": s["session_start"][:19],
"sex": s["subject"]["sex"], "age": s["subject"]["age"],
"volumes": s["processed"]["shape"][0], "raw_shape": s["raw"]["shape"][1:], "fda_shape": s["processed"]["shape"][1:],
"volume_hz": round(s["processed"]["rate"], 5), "raw_hz": round(s["raw"]["rate"], 5),
"fictrac_frames": s["fictrac"]["SpatialSeriesRotationDeltaLabFrame"]["data"].shape[0],
"fictrac_hz": s["fictrac"]["SpatialSeriesRotationDeltaLabFrame"]["rate"],
"anat_start_s": round(s["anatomical"]["start"], 1)} for f, s in by_fly.items()])
print("\nGenotypes:", sorted({s["subject"]["genotype"] for s in sessions}))
print("Species:", sorted({s["subject"]["species"] for s in sessions}))
print("\n", info.to_string(index=False))
print("\nFlies:", len(flies), "; distinct subject ids:", info.fly.nunique(), "; sessions per fly:", info.groupby("fly").size().max())
print("Volumes over all flies:", int(info.volumes.sum()))
nwb_gb = [sizes[s["path"]] / 1e9 for s in sessions]
streamed_gb = [sum(c[2] for z in s["chunk_index"] for c in s["chunk_index"][z]) / 1e9 for s in sessions]
print(f"NWB file sizes {min(nwb_gb):.1f} to {max(nwb_gb):.1f} GB; compressed bytes of the {len(bw.SLICES)} processed planes streamed per fly "
f"{min(streamed_gb):.2f} to {max(streamed_gb):.2f} GB, {sum(streamed_gb):.1f} GB in all; all files {sum(sizes.values()) / 1e9:.0f} GB")
dur_img = info.volumes / info.volume_hz
dur_fic = info.fictrac_frames / info.fictrac_hz
print(f"Imaging length {dur_img.min():.1f} to {dur_img.max():.1f} s; FicTrac length {dur_fic.min():.1f} to {dur_fic.max():.1f} s")
print("\nEvery array of two or more dimensions in", sessions[0]["path"].split("/")[-1])
print(pd.DataFrame(sessions[0]["datasets"]).to_string(index=False))
print("Same arrays, shapes and chunks in every file:", len({str(s["datasets"]) for s in sessions}) == 1)
print("Raw red channel description:", repr(sessions[0]["raw_red_description"]))
print("Start times of imaging and FicTrac series (s):",
sorted({s["processed"]["start"] for s in sessions} | {s["raw"]["start"] for s in sessions}),
sorted({v["start"] for s in sessions for v in s["fictrac"].values()}))
Assets by type: {'avi': 9, 'nwb': 9}
Behaviour videos: ['sub-fly-087/sub-fly-087_ses-20200228161226_behavior+image+ophys/715320fc-831c-4218-a723-32ca9e147cb2_external_file_0.avi', 'sub-fly-089/sub-fly-089_ses-20200618081653_behavior+image+ophys/5c033c49-ed7c-49c1-841a-10cc13faf1d0_external_file_0.avi'] ...
Archive summary: subjects 9 ; files 18 ; species ['Drosophila melanogaster - Fruit fly']
Genotypes: ['w+/w+;UAS-myr::tdTomato/UAS-GCaMP6f;nSyb-Gal4/+']
Species: ['Drosophila melanogaster']
fly folder session_id start sex age volumes raw_shape fda_shape volume_hz raw_hz fictrac_frames fictrac_hz anat_start_s
fly_087 sub-fly-087 20200228161226 2020-02-28T16:12:42 F P3D/P4D 3384 (256, 128, 49) (314, 146, 91) 1.87971 1.87971 90000 50.0 1926.6
fly_089 sub-fly-089 20200618081653 2020-06-18T08:17:27 F P3D/P4D 3384 (256, 128, 49) (314, 146, 91) 1.87971 1.87971 90000 50.0 2013.7
fly_094 sub-fly-094 20200620122048 2020-06-20T12:20:53 F P3D/P4D 3384 (256, 128, 49) (314, 146, 91) 1.87971 1.87971 90000 50.0 1964.4
fly_097 sub-fly-097 20200627113329 2020-06-27T11:33:37 F P3D/P4D 3384 (256, 128, 49) (314, 146, 91) 1.87971 1.87971 90000 50.0 1894.0
fly_098 sub-fly-098 20200627130208 2020-06-27T13:02:18 F P3D/P4D 3384 (256, 128, 49) (314, 146, 91) 1.87971 1.87971 90000 50.0 1884.4
fly_099 sub-fly-099 20200627142647 2020-06-27T14:27:04 F P3D/P4D 3384 (256, 128, 49) (314, 146, 91) 1.87971 1.87971 90000 50.0 2240.4
fly_100 sub-fly-100 20200702083849 2020-07-02T08:38:55 F P3D/P4D 3384 (256, 128, 49) (314, 146, 91) 1.87971 1.87971 90000 50.0 1949.2
fly_101 sub-fly-101 20200702100612 2020-07-02T10:06:21 F P3D/P4D 3384 (256, 128, 49) (314, 146, 91) 1.87971 1.87971 90000 50.0 1938.4
fly_105 sub-fly-105 20200704142706 2020-07-04T14:27:12 F P3D/P4D 3384 (256, 128, 49) (314, 146, 91) 1.87971 1.87971 90000 50.0 1900.8
Flies: 9 ; distinct subject ids: 9 ; sessions per fly: 1
Volumes over all flies: 30456
NWB file sizes 66.6 to 69.0 GB; compressed bytes of the 12 processed planes streamed per fly 2.29 to 2.61 GB, 22.1 GB in all; all files 610 GB
Imaging length 1800.3 to 1800.3 s; FicTrac length 1800.0 to 1800.0 s
Every array of two or more dimensions in sub-fly-087_ses-20200228161226_behavior+image+ophys.nwb
dataset shape dtype chunks compression
acquisition/TwoPhotonSeriesAnatomicalGreen/data (100, 1024, 512, 241) uint16 (9, 1024, 512, 1) gzip
acquisition/TwoPhotonSeriesAnatomicalRed/data (100, 1024, 512, 241) uint16 (9, 1024, 512, 1) gzip
acquisition/TwoPhotonSeriesFunctionalGreen/data (3384, 256, 128, 49) uint16 (152, 256, 128, 1) gzip
acquisition/TwoPhotonSeriesFunctionalRed/data (3384, 256, 128, 49) uint16 (152, 256, 128, 1) gzip
processing/behavior/FicTrac/SpatialSeriesAnimalHeading/data (90000, 1) float64 (1407, 1) gzip
processing/behavior/FicTrac/SpatialSeriesIntegratedMotion/data (90000, 2) float64 (2813, 1) gzip
processing/behavior/FicTrac/SpatialSeriesMovementDirection/data (90000, 1) float64 (1407, 1) gzip
processing/behavior/FicTrac/SpatialSeriesMovementSpeed/data (90000, 1) float64 (1407, 1) gzip
processing/behavior/FicTrac/SpatialSeriesPosition/data (90000, 2) float64 (2813, 1) gzip
processing/behavior/FicTrac/SpatialSeriesRotationCameraFrame/data (90000, 3) float64 (2813, 1) gzip
processing/behavior/FicTrac/SpatialSeriesRotationDeltaCameraFrame/data (90000, 3) float64 (2813, 1) gzip
processing/behavior/FicTrac/SpatialSeriesRotationDeltaError/data (90000, 1) float64 (1407, 1) gzip
processing/behavior/FicTrac/SpatialSeriesRotationDeltaLabFrame/data (90000, 3) float64 (2813, 1) gzip
processing/behavior/FicTrac/SpatialSeriesRotationLabFrame/data (90000, 3) float64 (2813, 1) gzip
processing/ophys/TwoPhotonSeriesFunctionalGreenProcessed/data (3384, 314, 146, 91) float32 (54, 314, 146, 1) gzip
Same arrays, shapes and chunks in every file: True
Raw red channel description: 'Functional imaging data (tdTomato)'
Start times of imaging and FicTrac series (s): [0.0] [0.0]
The archive's summary (9 subjects, 18 files, Drosophila melanogaster) matches the files: nine
NWB files and nine behaviour videos (AVI), which are not used here. Every NWB file holds the same
arrays with the same shapes and chunks: the raw GCaMP6f and tdTomato series (chunks of 152 volumes
x one plane), two anatomical scans, the processed series (chunks of 54 volumes x one plane) and
the FicTrac series. The flies are
females of the stated genotype (nSyb-Gal4 driving GCaMP6f and myr::tdTomato), one session each,
with the subject ids fly_087 to fly_105 that the paper's code lists (the code also drops a
fly_095 that failed its quality control; it is not in the dandiset). Each holds 3384 volumes at
1.88 Hz (1800.3 s) and 90,000 FicTrac frames at 50 Hz (1800.0 s), so the 30 minute sessions and
the 30,456 volumes of the paper are there. Imaging and FicTrac both start at 0 s; the conversion
code (catalystneuro/clandinin-lab-to-nwb, commit f083b09) wrote FicTrac with uniform
timestamps from 0 because synchronisation data were missing, and the paper's code makes the same
assumption (FicTrac frame 0 at the start of imaging). The clock check below tests it. The
anatomical scan starts 1884 to 2240 s after the functional scan, so it was taken after it, as
the Methods say. The age field (P3D/P4D) gives the 3 to 4 days after eclosion of the Methods.
s0 = sessions[0]
print("Raw functional:", s0["raw"]["shape"], s0["raw"]["dtype"], "grid spacing", s0["raw"]["grid_spacing"], s0["raw"]["grid_spacing_unit"])
print(" voxels per volume:", int(np.prod(s0["raw"]["shape"][1:])))
print("Processed:", s0["processed"]["shape"], s0["processed"]["dtype"], "chunks", s0["processed"]["chunks"],
"grid spacing", s0["processed"]["grid_spacing"], repr(s0["processed"]["grid_spacing_unit"]))
print(" description:", s0["processed"]["description"])
print("Every file has the same shapes:", len({(s["raw"]["shape"], s["processed"]["shape"]) for s in sessions}) == 1)
print("\nFicTrac series (unit, conversion, description):")
for name, v in s0["fictrac"].items():
print(f" {name:40s} {str(v['data'].shape):12s} {v['unit']:7s} {v['conversion']:<7} {v['description'][:110]}")
rows = []
for f, s in by_fly.items():
ft = s["fictrac"]
d = ft["SpatialSeriesRotationDeltaLabFrame"]["data"]
speed = ft["SpatialSeriesMovementSpeed"]["data"][:, 0]
err = ft["SpatialSeriesRotationDeltaError"]["data"][:, 0]
v = bw.smoothed_velocity(s)
rows.append({"fly": f, "dRotLab_x_1_99": np.round(np.percentile(d[:, 0], [1, 99]), 3).tolist(),
"dRotLab_y_1_99": np.round(np.percentile(d[:, 1], [1, 99]), 3).tolist(),
"dRotLab_z_1_99": np.round(np.percentile(d[:, 2], [1, 99]), 3).tolist(),
"speed_vs_xy_norm_r": round(np.corrcoef(speed, np.hypot(d[:, 0], d[:, 1]))[0, 1], 5),
"delta_error_median": round(float(np.median(err)), 0),
"fwd_mm_s_1_99": np.round(np.percentile(v["forward_mm_s"], [1, 99]), 1).tolist(),
"rot_deg_s_1_99": np.round(np.percentile(v["rotation_deg_s"], [1, 99]), 0).tolist(),
"walking_pct": round(100 * bw.walking_state(v["Y"], v["Z"]).mean(), 1),
"mean_VF_mm_s": round(float(np.clip(v["forward_mm_s"], 0, None).mean()), 2)})
beh_table = pd.DataFrame(rows)
print("\n", beh_table.to_string(index=False))
Raw functional: (3384, 256, 128, 49) uint16 grid spacing [2.61117299284587e-06, 2.61117299284587e-06, 5e-06] meters
voxels per volume: 1605632
Processed: (3384, 314, 146, 91) float32 chunks (54, 314, 146, 1) grid spacing [1.0, 1.0, 1.0] 'unknown'
description: motion corrected high-pass temporal filtered z-scored and registered into a common space, the Functional Drosophila Atlas
Every file has the same shapes: True
FicTrac series (unit, conversion, description):
SpatialSeriesAnimalHeading (90000, 1) meters 0.0045 Animal's heading direction in radians from the lab's perspective.
SpatialSeriesIntegratedMotion (90000, 2) meters 0.0045 Integrated x/y position of the sphere in laboratory coordinates, neglecting heading.
SpatialSeriesMovementDirection (90000, 1) meters 0.0045 Instantaneous running direction of the animal in the lab coordinates. Direction inferred by the ball's rotatio
SpatialSeriesMovementSpeed (90000, 1) meters 0.0045 Instantaneous running speed of the animal in radians per frame.
SpatialSeriesPosition (90000, 2) meters 0.0045 x and y positions in the lab frame in radians, inferred by integrating the rotation over time.
SpatialSeriesRotationCameraFrame (90000, 3) meters 0.0045 Orientation in radians from the camera's perspective. x: rotation to the sphere's right (pitch), y: rotation u
SpatialSeriesRotationDeltaCameraFrame (90000, 3) meters 0.0045 Change in orientation since last frame from the camera's perspective. x: rotation to the sphere's right (pitch
SpatialSeriesRotationDeltaError (90000, 1) meters 0.0045 Error in rotation delta in radians from the lab's perspective.
SpatialSeriesRotationDeltaLabFrame (90000, 3) meters 0.0045 Change in orientation since last frame from the lab's perspective. x: rotation in front of subject (roll), y:
SpatialSeriesRotationLabFrame (90000, 3) meters 0.0045 Orientation in radians from the lab's perspective. x: rotation in front of subject (roll), y: rotation to subj
fly dRotLab_x_1_99 dRotLab_y_1_99 dRotLab_z_1_99 speed_vs_xy_norm_r delta_error_median fwd_mm_s_1_99 rot_deg_s_1_99 walking_pct mean_VF_mm_s
fly_087 [-0.056, 0.039] [-0.031, 0.084] [-0.12, 0.097] 1.0 6929.0 [-3.1, 14.8] [-298.0, 228.0] 55.2 2.48
fly_089 [-0.027, 0.032] [-0.027, 0.047] [-0.056, 0.066] 1.0 6168.0 [-2.5, 5.8] [-138.0, 160.0] 48.5 0.91
fly_094 [-0.047, 0.02] [-0.018, 0.061] [-0.028, 0.102] 1.0 7109.0 [-0.8, 10.5] [-61.0, 280.0] 28.3 1.11
fly_097 [-0.028, 0.038] [-0.035, 0.031] [-0.077, 0.059] 1.0 8933.0 [-3.4, 2.5] [-208.0, 149.0] 54.3 0.21
fly_098 [-0.038, 0.044] [-0.036, 0.062] [-0.108, 0.097] 1.0 7614.0 [-3.3, 9.0] [-272.0, 246.0] 46.0 1.15
fly_099 [-0.043, 0.046] [-0.038, 0.055] [-0.137, 0.044] 1.0 10435.0 [-3.6, 7.5] [-361.0, 92.0] 44.8 0.85
fly_100 [-0.055, 0.043] [-0.028, 0.069] [-0.138, 0.115] 1.0 6179.0 [-2.3, 10.5] [-371.0, 286.0] 61.1 2.37
fly_101 [-0.057, 0.04] [-0.029, 0.072] [-0.098, 0.101] 1.0 6067.0 [-2.2, 11.1] [-238.0, 247.0] 53.2 1.78
fly_105 [-0.031, 0.04] [-0.023, 0.037] [-0.085, 0.073] 1.0 7654.0 [-2.0, 6.1] [-199.0, 170.0] 40.4 0.41
The raw volumes are 256 x 128 x 49 voxels of 2.61 x 2.61 x 5 um, as the Methods state (1,605,632 voxels per volume). The processed series is 314 x 146 x 91 voxels with a grid spacing of 1, 1, 1 in unit "unknown", so its voxel size is not given; its description names the FDA, a space the preprint does not use (the preprint's analyses are in its own mean brain, of the raw volume's 256 x 128 x 49 shape).
Every FicTrac series is stored with unit "meters" and conversion 0.0045 (the ball radius), while
the descriptions give radians for angles and positions and "radians per frame" for speed. The
values are radians (per frame for the deltas): the 1st to 99th percentiles of the lab-frame
deltas lie within 0.14 rad, and the movement speed equals the norm of the roll and pitch deltas
(r = 1.0). Multiplying a
heading or an orientation by the conversion gives no meaningful length, so the unit and
conversion are wrong for those series. SpatialSeriesRotationDeltaError, described as an error in
radians, has medians of 6067 to 10,435, which cannot be radians. With the 4.5 mm radius and
50 Hz, the 99th percentile of forward velocity is 2.5 to 14.8 mm/s and rotation reaches 92 to
371 deg/s (1st and 99th percentiles above), plausible for a tethered fly; fly_097 walks least
(mean VF 0.21 mm/s).
import hashlib
digest = {f: hashlib.md5(s["fictrac"]["SpatialSeriesRotationDeltaLabFrame"]["data"].tobytes()).hexdigest()[:12]
for f, s in by_fly.items()}
print("Distinct FicTrac series:", len(set(digest.values())), "of", len(digest))
fwd = np.array([bw.smoothed_velocity(s)["forward_mm_s"] for s in sessions])
cc = np.corrcoef(fwd)
print("Largest correlation of forward velocity between two flies:", round(float(np.max(cc[np.triu_indices(9, 1)])), 3))
img = {f: hashlib.md5(s["h5py_check"].tobytes()).hexdigest()[:12] for f, s in by_fly.items()}
print("Distinct processed blocks (first 54 volumes of plane", bw.SLICES[len(bw.SLICES) // 2], "):", len(set(img.values())), "of", len(img))
check_z = bw.SLICES[len(bw.SLICES) // 2]
print("Direct chunk reads equal the h5py values in every file:",
[planes[f][check_z].get("same_as_h5py") for f in flies])
Distinct FicTrac series: 9 of 9
Largest correlation of forward velocity between two flies: 0.133
Distinct processed blocks (first 54 volumes of plane 52 ): 9 of 9 Direct chunk reads equal the h5py values in every file: [True, True, True, True, True, True, True, True, True]
No recording is deposited twice: the nine FicTrac series and the nine processed volumes differ, and no two flies' forward velocities are correlated. The direct reads of compressed chunks used here give the same values as h5py.
Z = bw.SLICES
n_vol = {f: by_fly[f]["processed"]["shape"][0] for f in flies}
N = sum(n_vol.values())
inside = np.array([[planes[f][z]["inside"] for z in Z] for f in flies]) # fly, plane, x, y
brain = inside.all(0)
print("FDA planes used:", list(Z))
print("In-brain voxels per plane (in all 9 flies):", dict(zip(Z, brain.sum((1, 2)).tolist())))
print(f"In the brain of all 9 flies: {brain.sum()} voxels; of at least one: {inside.any(0).sum()}; "
f"per fly {inside.sum((1, 2, 3)).min()} to {inside.sum((1, 2, 3)).max()}")
rows = []
for i, f in enumerate(flies):
sx = np.array([planes[f][z]["sx"] for z in Z])[inside[i]]
sxx = np.array([planes[f][z]["sxx"] for z in Z])[inside[i]]
mu = sx / n_vol[f]
sd = np.sqrt(sxx / n_vol[f] - mu ** 2)
rows.append({"fly": f, "mean_abs_max": float(np.abs(mu).max()), "sd_1": np.percentile(sd, 1),
"sd_median": np.median(sd), "sd_99": np.percentile(sd, 99), "sd_max": sd.max()})
print("\nPer-voxel mean and s.d. over time of the processed series, in-brain voxels:")
print(pd.DataFrame(rows).round(4).to_string(index=False))
# Left-right axis: mirror symmetry of the brain mask on each in-plane axis
both_any = inside.any(0)
for axis, size in ((1, 314), (2, 146)):
best = []
for c2 in range(size - 12, size + 13): # mirror x -> c2 - x, c2 = twice the centre
idx = c2 - np.arange(size)
ok = (idx >= 0) & (idx < size)
a = np.take(both_any, np.arange(size)[ok], axis=axis)
b = np.take(both_any, idx[ok], axis=axis)
best.append((float((a & b).sum() / max((a | b).sum(), 1)), c2 / 2))
best.sort(reverse=True)
print(f"Axis of {size} voxels: best mirror centre {best[0][1]} (overlap {best[0][0]:.3f}); centre {(size - 1) / 2} gives "
f"{[o for o, c in best if c == (size - 1) / 2][0]:.3f}")
FDA planes used: [10, 17, 24, 31, 38, 45, 52, 59, 66, 73, 80, 87]
In-brain voxels per plane (in all 9 flies): {10: 181, 17: 3007, 24: 9552, 31: 16875, 38: 23010, 45: 26548, 52: 27272, 59: 26976, 66: 24671, 73: 19287, 80: 11457, 87: 2283}
In the brain of all 9 flies: 191119 voxels; of at least one: 309079; per fly 231265 to 265909
Per-voxel mean and s.d. over time of the processed series, in-brain voxels:
fly mean_abs_max sd_1 sd_median sd_99 sd_max
fly_087 0.0 0.0521 0.6439 0.8646 0.9883
fly_089 0.0 0.0571 0.6530 0.8716 0.9896
fly_094 0.0 0.0459 0.6467 0.8703 0.9823
fly_097 0.0 0.0533 0.6414 0.8638 0.9743
fly_098 0.0 0.0544 0.6455 0.8684 0.9806
fly_099 0.0 0.0489 0.6498 0.8735 0.9878
fly_100 0.0 0.0547 0.6508 0.8692 0.9848
fly_101 0.0 0.0475 0.6470 0.8700 0.9922
fly_105 0.0 0.0539 0.6488 0.8708 0.9844
Axis of 314 voxels: best mirror centre 156.5 (overlap 0.914); centre 156.5 gives 0.914
Axis of 146 voxels: best mirror centre 75.0 (overlap 0.665); centre 72.5 gives 0.662
The processed values are consistent with the description: every in-brain voxel has a mean of 0 over time (to the four decimals printed) and a s.d. below 1 (at most 0.992), as expected if each voxel was z-scored and then warped (interpolation lowers the spread; the Methods z-score before alignment). Voxels outside each fly's brain are constant 0. The in-brain masks differ between flies; the analysis uses voxels in the brain of all nine. On the 314-voxel axis the mask is close to mirror-symmetric about its centre, 156.5 (overlap 0.914), and less so on the 146-voxel axis (best 0.665), so the 314-voxel axis is left-right and its centre is the midline; the maps and blocks are mirrored about 156.5. Which half is the fly's left is not stated anywhere in the dandiset.
# Clock check. Brain-wide mean activity against walking speed, and the left-minus-right hemisphere
# difference against signed rotation, with the behaviour shifted in time. A shared clock gives a
# peak at a small positive lag (GCaMP6f follows behaviour); a broken one gives no peak near 0.
rng = np.random.default_rng(0)
nbx, nby = bw.block_layout()
counts = np.array([[planes[f][z]["counts"] for z in Z] for f in flies]) # fly, plane, block
keep = (counts >= bw.BLOCK ** 2 / 2).all(0) # plane, block
bx_of = np.tile(np.arange(nbx * nby) // nby, (len(Z), 1))
lags = np.arange(-10, 10.01, 0.2)
rows = []
for f, s in by_fly.items():
blocks = np.concatenate([planes[f][z]["blocks"][:, keep[i]] for i, z in enumerate(Z)], 1)
bxs = bx_of[keep]
glob = blocks.mean(1)
lr = blocks[:, bxs < nbx // 2].mean(1) - blocks[:, bxs >= nbx // 2].mean(1)
v = bw.smoothed_velocity(s)
speed = np.abs(v["forward_mm_s"])
out = {"fly": f}
for name, sig, src in (("global~speed", glob, speed), ("halfdiff~rotation", lr, v["rotation_deg_s"])):
def lag_r(shift):
M = np.interp((times[f][None, :] - shift - lags[:, None]) % 1800, v["t"], src) # lag, time
M = M - M.mean(1, keepdims=True)
sc = sig - sig.mean()
return (M @ sc) / np.sqrt((M ** 2).sum(1) * (sc ** 2).sum())
rs = lag_r(0.0)
k = np.argmax(np.abs(rs))
null = [np.abs(lag_r(rng.uniform(60, 1740))).max() for _ in range(100)] # same search over lags
out[name + "_peak_lag_s"] = round(lags[k], 1)
out[name + "_peak_r"] = round(rs[k], 3)
out[name + "_null_99pct"] = round(np.percentile(null, 99), 3)
rows.append(out)
clock = pd.DataFrame(rows)
print("Positive lag: neural activity follows behaviour. Null: the largest |r| over the same 101 lags, behaviour shifted by 60 to 1740 s.")
print(clock.to_string(index=False))
Positive lag: neural activity follows behaviour. Null: the largest |r| over the same 101 lags, behaviour shifted by 60 to 1740 s.
fly global~speed_peak_lag_s global~speed_peak_r global~speed_null_99pct halfdiff~rotation_peak_lag_s halfdiff~rotation_peak_r halfdiff~rotation_null_99pct
fly_087 0.4 0.412 0.103 -0.0 0.342 0.079
fly_089 0.8 0.577 0.087 0.4 0.184 0.066
fly_094 0.4 0.684 0.084 -0.2 0.208 0.072
fly_097 0.4 0.522 0.079 -0.2 0.176 0.069
fly_098 0.4 0.453 0.097 -0.0 0.159 0.067
fly_099 0.4 0.519 0.083 0.4 0.481 0.075
fly_100 0.4 0.455 0.088 -0.0 0.185 0.079
fly_101 0.4 0.439 0.069 -0.0 0.327 0.078
fly_105 0.6 0.499 0.073 0.8 -0.100 0.079
In every fly the brain-wide mean correlates with walking speed most strongly at a lag of 0.4 to 0.8 s (neural activity after behaviour), at r = 0.412 to 0.684. The null is the largest |r| over the same 101 lags with the behaviour shifted by a minute or more, and its 99th percentile is at most 0.103. In 8 of 9 flies the difference between the two image halves correlates positively with signed rotation (0.159 to 0.481, at lags of -0.2 to 0.4 s), against null 99th percentiles of at most 0.079; in fly_105 the largest correlation is -0.100, at 0.8 s, larger in size than its null (0.079) but of the other sign. The behaviour and imaging streams are consistent with one clock, to within the lag of GCaMP6f. A shared start at 0 s is an assumption of the conversion and of the paper's code; this check supports it to within a few hundred milliseconds per fly, the resolution set by the volume rate, and cannot test offsets smaller than that.
dataset-structure¶
Nine flies were imaged for 30 minutes each at a volume rate of about 1.8 Hz, with 256 x 128 x 49 voxels of 2.6 x 2.6 x 5 um (about 1.6 million voxels per volume), giving 30,456 volumes in all, while the ball was tracked at 50 Hz.
structure = pd.DataFrame([
("flies", 9, len(flies), 0),
("volumes over all flies", 30456, int(info.volumes.sum()), 0),
("volume rate (Hz)", 1.8, round(float(info.volume_hz.mean()), 4), 0.18),
("voxels per volume", 1600000, int(np.prod(s0["raw"]["shape"][1:])), 160000),
("ball tracking rate (Hz)", 50, float(info.fictrac_hz.mean()), 0),
], columns=["quantity", "expected", "measured", "tolerance"])
structure["within"] = (structure.measured - structure.expected).abs() <= structure.tolerance
print(structure.to_string(index=False))
print("Voxel size (um):", np.round(np.array(s0["raw"]["grid_spacing"]) * 1e6, 2).tolist(),
"; session length (min):", round(float(dur_img.mean()) / 60, 2))
quantity expected measured tolerance within
flies 9.0 9.0000 0.00 True
volumes over all flies 30456.0 30456.0000 0.00 True
volume rate (Hz) 1.8 1.8797 0.18 True
voxels per volume 1600000.0 1605632.0000 160000.00 True
ball tracking rate (Hz) 50.0 50.0000 0.00 True
Voxel size (um): [2.61, 2.61, 5.0] ; session length (min): 30.0
Result: reproduced. The dandiset holds 9 flies with 30,456 volumes in all, at 1.8797 Hz, of 1,605,632 voxels of 2.61 x 2.61 x 5 um, over 30.0 minutes, with ball tracking at 50 Hz; every value is within its tolerance.
Voxel correlations with velocity (Figure 2)¶
The correlation of each voxel with VF, VL and VR, pooled over the nine flies' 30,456 volumes, is computed from per-fly sums kept while streaming. The thresholds of the paper are |r| for two-sided p < 1e-4 (single voxels) and p < 1e-8 (supervoxels, Bonferroni over 2000 x 49). As in the paper, the p-values treat time points as independent.
def pooled_maps(k_offset=0):
sx = sum(np.array([planes[f][z]["sx"] for z in Z]) for f in flies)
sxx = sum(np.array([planes[f][z]["sxx"] for z in Z]) for f in flies)
sxb = sum(np.array([planes[f][z]["sxb"] for z in Z]) for f in flies) # plane, k, x, y
src = {0: beh, 3: beh_start, 6: beh_rolled}[k_offset]
out = []
for k, key in enumerate(("fwd", "left", "right")):
b = np.concatenate([src[f][key] for f in flies])
out.append(bw.pooled_r(sx, sxx, sxb[:, k + k_offset], b.sum(), (b ** 2).sum(), N))
return np.nan_to_num(np.array(out)) # k, plane, x, y
r_vox = pooled_maps(0)
r_vox_start = pooled_maps(3)
r_vox_rolled = pooled_maps(6)
thr_vox = bw.r_threshold(1e-4, N)
thr_sv = bw.r_threshold(1e-8, N)
print(f"N = {N}; |r| threshold for p < 1e-4: {thr_vox:.4f}; for p < 1e-8: {thr_sv:.4f}")
for k, key in enumerate(("VF", "VL", "VR")):
vals = r_vox[k][brain]
print(f"{key}: in-brain voxel r from {vals.min():.3f} to {vals.max():.3f}, median {np.median(vals):.3f}; "
f"{(vals > thr_vox).mean() * 100:.1f}% above +threshold, {(vals < -thr_vox).mean() * 100:.2f}% below -threshold")
# Block signals pooled over flies (the supervoxels of this attempt)
X = np.concatenate([np.concatenate([planes[f][z]["blocks"][:, keep[i]] for i, z in enumerate(Z)], 1) for f in flies], 0)
X -= X.mean(0)
block_plane = np.concatenate([np.full(keep[i].sum(), z) for i, z in enumerate(Z)])
block_index = np.concatenate([np.flatnonzero(keep[i]) for i in range(len(Z))])
block_bx = block_index // nby
brain_blocks = np.concatenate([brain[i, bw.X_CROP[0]:bw.X_CROP[1], bw.Y_CROP[0]:bw.Y_CROP[1]]
.reshape(nbx, bw.BLOCK, nby, bw.BLOCK).sum((1, 3)).ravel()[keep[i]] for i in range(len(Z))])
V = {key: np.concatenate([beh[f][key] for f in flies]) for key in ("fwd", "left", "right", "walk", "Y", "Z")}
def corr_cols(M, b):
bc = (b - b.mean()) / b.std()
return (M * bc[:, None]).mean(0) / M.std(0)
r_blk = np.array([corr_cols(X, V[key]) for key in ("fwd", "left", "right")])
print(f"\nBlocks kept: {X.shape[1]} of {keep.size} (at least half their voxels in the brain of every fly); "
f"they hold {brain_blocks.sum()} of {brain.sum()} in-brain voxels")
print("Block r ranges:", {key: (round(float(r_blk[k].min()), 3), round(float(r_blk[k].max()), 3)) for k, key in enumerate(("VF", "VL", "VR"))})
N = 30456; |r| threshold for p < 1e-4: 0.0223; for p < 1e-8: 0.0328 VF: in-brain voxel r from -0.028 to 0.123, median 0.013; 27.8% above +threshold, 0.00% below -threshold VL: in-brain voxel r from -0.025 to 0.091, median 0.009; 12.8% above +threshold, 0.01% below -threshold VR: in-brain voxel r from -0.032 to 0.100, median 0.007; 10.0% above +threshold, 0.01% below -threshold
Blocks kept: 5363 of 14976 (at least half their voxels in the brain of every fly); they hold 186298 of 191119 in-brain voxels
Block r ranges: {'VF': (-0.032, 0.252), 'VL': (-0.027, 0.162), 'VR': (-0.033, 0.167)}
fig2c-fraction-correlated¶
About 39% of the brain volume correlates with at least one of forward velocity, left and right angular velocity.
Result: partially reproduced. Under reading (a), 31.1% of in-brain voxels correlate with VF, VL or VR at p < 1e-4, below the expected 35.1 to 42.9% and above half of 39%; under reading (b), the blocks significant at p < 1e-8 cover 53.5% of the in-brain voxels, above it. No reading is within tolerance. With each fly's behaviour rolled in time, no voxel or block passes either threshold (0.0%). With behaviour sampled at the start of each volume in place of its midpoint, reading (a) gives 39.2%, so the value depends on neural times that the dandiset does not hold. The correlated fraction is 0.0 to 7.5% in the outermost planes (10, 80, 87) and 21.0 to 39.1% in the others. Within single flies, with the threshold for 3384 points, 2.4 to 19.7% of voxels pass.
fig2c-venn-categories¶
Of the brain volume, about 15% correlates only with forward velocity, 20% with forward velocity and turning, and 4% only with turning; only about 10% of correlated voxels correlate with all three components, so 90% show a tuning preference.
expected_venn = {"VF only %": (15, 1.5), "VF and turning %": (20, 2.0), "turning only %": (4, 0.4),
"all three, % of correlated": (10, 1.0)}
rows = []
for name, vals in readings.items():
row = {"reading": name}
for q, (e, tol) in expected_venn.items():
row[q] = round(vals[q], 1)
row[q + " ok"] = abs(vals[q] - e) <= tol
row["ranking holds"] = vals["VF and turning %"] > vals["VF only %"] > vals["turning only %"]
row["all three < 20% of correlated"] = vals["all three, % of correlated"] < 20
rows.append(row)
venn = pd.DataFrame(rows)
print(venn.to_string(index=False))
print("\nExpected:", expected_venn)
reading VF only % VF only % ok VF and turning % VF and turning % ok turning only % turning only % ok all three, % of correlated all three, % of correlated ok ranking holds all three < 20% of correlated
(a) voxels, p<1e-4 11.9 False 15.9 False 3.2 False 11.5 False True True
(b) blocks, p<1e-8 12.3 False 39.4 False 1.8 False 32.2 False True False
variant: voxels, behaviour at volume start 11.6 False 24.1 False 3.5 False 25.7 False True False
variant: blocks, p<1e-8, volume start 10.8 False 51.7 False 2.3 False 51.2 False True False
variant: voxels in the brain of any fly 7.9 False 10.0 False 2.4 False 11.0 False True True
variant: voxels, positive r only 11.9 False 15.9 False 3.2 False 11.5 False True True
variant: blocks, p<1e-4 10.2 False 55.0 False 3.4 False 46.2 False True False
null: voxels, behaviour rolled 0.0 False 0.0 False 0.0 False 0.0 False False True
null: blocks, p<1e-8, behaviour rolled 0.0 False 0.0 False 0.0 False 0.0 False False True
Expected: {'VF only %': (15, 1.5), 'VF and turning %': (20, 2.0), 'turning only %': (4, 0.4), 'all three, % of correlated': (10, 1.0)}
Result: partially reproduced. Under reading (a), 11.9% of in-brain voxels correlate with VF only (expected 15), 15.9% with VF and turning (20), 3.2% with turning only (4), and 11.5% of the correlated voxels with all three (10). None is within its tolerance, but the ranking holds (forward and turning above forward only, above turning only) and fewer than 20% of correlated voxels correlate with all three. Under reading (b) forward and turning is 39.4% and all three 32.2% of correlated volume. With behaviour at the start of each volume, reading (a) gives 11.6, 24.1, 3.5 and 25.7%.
fig2-symmetry¶
Forward velocity correlations are mirror-symmetric across the midline, while left and right turning correlations are anti-symmetric, with more correlated voxels in the hemisphere ipsilateral to the turn.
xs = np.arange(314)
xm = 313 - xs # mirror about 156.5
pair = brain & brain[:, xm, :]
def mirror(m):
return m[:, xm, :]
sym = {}
for k, key in enumerate(("VF", "VL", "VR")):
sym[key + " vs own mirror"] = np.corrcoef(r_vox[k][pair], mirror(r_vox[k])[pair])[0, 1]
sym["VL vs mirrored VR"] = np.corrcoef(r_vox[1][pair], mirror(r_vox[2])[pair])[0, 1]
print("Correlation of voxel r maps over in-brain mirror pairs:", {k: round(v, 3) for k, v in sym.items()})
halfA = (xs < 157)[None, :, None] & brain
halfB = (xs >= 157)[None, :, None] & brain
sig_counts = pd.DataFrame({key: {"half A (x < 157)": int(((np.abs(r_vox[k]) > thr_vox) & halfA).sum()),
"half B (x >= 157)": int(((np.abs(r_vox[k]) > thr_vox) & halfB).sum())}
for k, key in enumerate(("VF", "VL", "VR"))}).T
sig_counts["A / B"] = sig_counts.iloc[:, 0] / sig_counts.iloc[:, 1]
print("\nIn-brain voxels with p < 1e-4, by image half (in-brain voxels: A", int(halfA.sum()), ", B", int(halfB.sum()), "):")
print(sig_counts.round(2).to_string())
ipsi_half = {"left": "A" if sig_counts.loc["VL", "A / B"] > 1 else "B", "right": "A" if sig_counts.loc["VR", "A / B"] > 1 else "B"}
crit = {"(i) VF mirror r >= 0.7": sym["VF vs own mirror"] >= 0.7,
"(ii) VL~mirrored VR > VL~own mirror": sym["VL vs mirrored VR"] > sym["VL vs own mirror"],
"(iii) VL and VR >= 2x in opposite halves": (max(sig_counts.loc["VL", "A / B"], 1 / sig_counts.loc["VL", "A / B"]) >= 2
and max(sig_counts.loc["VR", "A / B"], 1 / sig_counts.loc["VR", "A / B"]) >= 2
and ipsi_half["left"] != ipsi_half["right"])}
print("\nCriteria:", crit)
print("Half with more VL voxels:", ipsi_half["left"], "; with more VR voxels:", ipsi_half["right"])
# Orientation check: each fly against the pooled maps of the other eight (leave-one-out), unflipped and mirrored
tot = {q: sum(np.array([planes[g][z][q] for z in Z]) for g in flies) for q in ("sx", "sxx", "sxb")}
rows = []
for f in flies:
sx = np.array([planes[f][z]["sx"] for z in Z]); sxx = np.array([planes[f][z]["sxx"] for z in Z])
sxb = np.array([planes[f][z]["sxb"] for z in Z])
rs = np.nan_to_num(np.array([bw.pooled_r(sx, sxx, sxb[:, k], beh[f][key].sum(), (beh[f][key] ** 2).sum(), n_vol[f])
for k, key in enumerate(("fwd", "left", "right"))]))
i = flies.index(f)
own = inside[i]
loo = []
for k, key in enumerate(("fwd", "left", "right")):
b = np.concatenate([beh[g][key] for g in flies if g != f])
loo.append(bw.pooled_r(tot["sx"] - sx, tot["sxx"] - sxx, tot["sxb"][:, k] - sxb[:, k], b.sum(), (b ** 2).sum(), N - n_vol[f]))
loo = np.nan_to_num(np.array(loo))
rows.append({"fly": f, "VF mirror r": np.corrcoef(rs[0][pair], mirror(rs[0])[pair])[0, 1],
"VL: sum r, A - B": rs[1][halfA].sum() - rs[1][halfB].sum(),
"VR: sum r, A - B": rs[2][halfA].sum() - rs[2][halfB].sum(),
"VF~others VF": np.corrcoef(rs[0][brain], loo[0][brain])[0, 1],
"VL~others VL": np.corrcoef(rs[1][pair], loo[1][pair])[0, 1],
"VL~mirrored others VL": np.corrcoef(rs[1][pair], mirror(loo[1])[pair])[0, 1],
"VR~others VR": np.corrcoef(rs[2][pair], loo[2][pair])[0, 1],
"VR~mirrored others VR": np.corrcoef(rs[2][pair], mirror(loo[2])[pair])[0, 1],
"mask overlap with all-fly mask": (own & brain).sum() / (own | brain).sum()})
per_fly_sym = pd.DataFrame(rows)
per_fly_sym["VL more in A"] = per_fly_sym["VL: sum r, A - B"] > 0
per_fly_sym["VR more in B"] = per_fly_sym["VR: sum r, A - B"] < 0
print("\nPer fly (each fly's own correlations; orientation checks against the pooled maps of the other eight flies, 'others'):")
print(per_fly_sym.round(3).to_string(index=False))
print("Flies with VL more in half A:", int(per_fly_sym["VL more in A"].sum()), "; VR more in half B:", int(per_fly_sym["VR more in B"].sum()),
"; both:", int((per_fly_sym["VL more in A"] & per_fly_sym["VR more in B"]).sum()), "of 9",
[f for f, a, b in zip(flies, per_fly_sym["VL more in A"], per_fly_sym["VR more in B"]) if a and b])
unflipped = (per_fly_sym["VL~others VL"] > per_fly_sym["VL~mirrored others VL"]) & (per_fly_sym["VR~others VR"] > per_fly_sym["VR~mirrored others VR"])
print("Flies whose VL and VR maps both match the others' maps better unflipped than mirrored:", int(unflipped.sum()), "of 9;",
"not:", [f for f, u in zip(flies, unflipped) if not u])
fig, axes = plt.subplots(1, 3, figsize=(12, 3.2))
zi = list(Z).index(59)
vals = r_vox[:, zi][:, brain[zi]]
cmax = float(np.abs(vals).max())
for k, (ax, key) in enumerate(zip(axes, ("VF", "VL", "VR"))):
im = ax.imshow(np.where(brain[zi], r_vox[k][zi], np.nan).T, cmap="RdBu_r", vmin=-cmax, vmax=cmax, origin="lower")
ax.set_title(f"{key}, FDA plane 59"); ax.set_xlabel("x (left-right, voxels)")
fig.colorbar(im, ax=axes, label="pooled r")
plt.show()
print(f"Plane 59 r range {vals.min():.3f} to {vals.max():.3f}; colour range +-{cmax:.3f}, so no voxel falls outside it")
Correlation of voxel r maps over in-brain mirror pairs: {'VF vs own mirror': np.float64(0.823), 'VL vs own mirror': np.float64(0.449), 'VR vs own mirror': np.float64(0.332), 'VL vs mirrored VR': np.float64(0.672)}
In-brain voxels with p < 1e-4, by image half (in-brain voxels: A 99154 , B 91965 ):
half A (x < 157) half B (x >= 157) A / B
VF 28062 25120 1.12
VL 17242 7234 2.38
VR 4552 14568 0.31
Criteria: {'(i) VF mirror r >= 0.7': np.True_, '(ii) VL~mirrored VR > VL~own mirror': np.True_, '(iii) VL and VR >= 2x in opposite halves': True}
Half with more VL voxels: A ; with more VR voxels: B
Per fly (each fly's own correlations; orientation checks against the pooled maps of the other eight flies, 'others'):
fly VF mirror r VL: sum r, A - B VR: sum r, A - B VF~others VF VL~others VL VL~mirrored others VL VR~others VR VR~mirrored others VR mask overlap with all-fly mask VL more in A VR more in B
fly_087 0.580 1284.577 -1054.030 0.590 0.365 0.148 0.429 0.158 0.758 True True
fly_089 0.289 289.967 -167.466 0.395 0.357 0.305 0.286 0.164 0.748 True True
fly_094 0.575 460.681 -238.552 0.630 0.438 0.314 0.369 0.262 0.784 True True
fly_097 0.171 23.935 -294.962 0.330 0.353 0.336 0.384 0.272 0.826 True True
fly_098 0.243 917.478 182.601 0.486 0.416 0.290 0.256 0.238 0.796 True False
fly_099 0.166 -278.132 -1250.241 0.423 0.158 0.190 0.439 0.199 0.766 False True
fly_100 0.358 377.192 -329.914 0.495 0.304 0.188 0.279 0.143 0.755 True True
fly_101 0.416 393.544 -705.686 0.525 0.273 0.123 0.386 0.178 0.788 True True
fly_105 0.168 287.815 169.434 0.356 0.274 0.234 0.212 0.214 0.719 True False
Flies with VL more in half A: 8 ; VR more in half B: 7 ; both: 6 of 9 ['fly_087', 'fly_089', 'fly_094', 'fly_097', 'fly_100', 'fly_101']
Flies whose VL and VR maps both match the others' maps better unflipped than mirrored: 7 of 9; not: ['fly_099', 'fly_105']
Plane 59 r range -0.028 to 0.100; colour range +-0.100, so no voxel falls outside it
Result: reproduced. The VF correlation map correlates with its mirror image at r = 0.823. The VL map correlates with the mirrored VR map at 0.672, more than with its own mirror image (0.449). Voxels significant for VL are 2.38 times as many in half A as in half B, and voxels significant for VR 0.31 times as many (3.2 times as many in half B). The pooled maps hide differences between flies: in single flies the VF map correlates with its mirror at 0.166 to 0.580; summed VL correlation is larger in half A in 8 of 9 flies and summed VR correlation larger in half B in 7 of 9, and both hold in 6 of 9 (fly_098 and fly_105 have VR on the same side as VL, fly_099 has VL on the other side). Against the pooled maps of the other eight flies, no fly matches better mirrored than unflipped for both turning directions, as a left-right flip would give. In 7 of 9 flies both the VL and the VR map match better unflipped. fly_099's VL map matches the mirrored map slightly better (0.158 against 0.190) while its VR map clearly does not (0.439 against 0.199). For fly_105 the check does not decide: VL 0.274 against 0.234, VR 0.212 against 0.214. A flip of fly_105 is therefore not excluded; if it were flipped, its maps would lower the pooled asymmetry. If turning correlates with the ipsilateral hemisphere, as the paper says, half A (x < 157) is the fly's left; the dandiset does not state which half is left, so the ipsilateral part is not tested. The symmetry threshold of 0.7 was chosen after seeing r = 0.82 in the prototype run (claims file header).
Principal components of the pooled blocks (Figure 4)¶
PCA of the block signals of all nine flies concatenated (the "superfly"), by eigendecomposition of the covariance matrix, as pca.py does on the paper's supervoxels.
t0 = time.time()
cov = (X.T @ X) / (X.shape[0] - 1)
evals, evecs = np.linalg.eigh(cov.astype(np.float64))
evals, evecs = evals[::-1], evecs[:, ::-1]
del cov
N_PCS = 500
pcs = X @ evecs[:, :N_PCS].astype(np.float32)
for i, key in ((0, "fwd"), (1, "fwd"), (2, "left")): # orientation fixed in the criterion
if np.corrcoef(pcs[:, i], V[key])[0, 1] < 0:
pcs[:, i] *= -1
evecs[:, i] *= -1
print(f"PCA of {X.shape[0]} x {X.shape[1]} in {time.time() - t0:.0f} s; variance explained by PCs 1-10 (%):",
np.round(100 * evals[:10] / evals.sum(), 2).tolist())
PCA of 30456 x 5363 in 19 s; variance explained by PCs 1-10 (%): [2.22, 0.41, 0.2, 0.18, 0.17, 0.15, 0.13, 0.1, 0.1, 0.09]
fig4bcd-pc-structure¶
The first principal components of brain-wide activity are structured by behaviour; PC1 is low when the fly is stopped and high otherwise, PC2 increases with forward velocity and not with angular velocity, and PC3 carries the sign of angular velocity.
def split_corr(sig, comp, reading):
"""Correlation of sig with comp below and above 0.75 s.d. of comp, with p-values."""
z = comp / comp.std() if reading == "a" else (comp - comp.mean()) / comp.std()
out = {}
for side, m in (("below", z <= 0.75), ("above", z > 0.75)):
r, p = stats.pearsonr(sig[m], comp[m])
out[side] = (r, p, int(m.sum()))
return out
pc_rows, pc_checks = [], {}
for reading in ("a", "b"):
res = {}
for i in range(3):
for key in ("fwd", "left", "right"):
sc = split_corr(pcs[:, i], V[key], reading)
res[(i, key)] = sc
pc_rows.append({"reading": reading, "PC": i + 1, "component": key,
"r below": sc["below"][0], "p below": sc["below"][1], "n below": sc["below"][2],
"r above": sc["above"][0], "p above": sc["above"][1], "n above": sc["above"][2]})
c1 = res[(0, "fwd")]["below"][0] > 0 and all(abs(res[(0, "fwd")]["below"][0]) > abs(res[(j, "fwd")]["below"][0]) for j in (1, 2))
c2 = (np.argmax([res[(j, "fwd")]["above"][0] for j in range(3)]) == 1
and abs(res[(1, "left")]["above"][0]) < abs(res[(1, "fwd")]["above"][0])
and abs(res[(1, "right")]["above"][0]) < abs(res[(1, "fwd")]["above"][0]))
l3, r3 = res[(2, "left")]["above"], res[(2, "right")]["above"]
c3 = np.sign(l3[0]) != np.sign(r3[0]) and l3[1] < 0.05 and r3[1] < 0.05
pc_checks[reading] = {"(i) PC1 below-threshold": c1, "(ii) PC2 forward above": c2, "(iii) PC3 turn sign": c3}
print(pd.DataFrame(pc_rows).to_string(index=False, float_format=lambda v: f"{v:.3g}"))
print("\nCriteria by reading:", pc_checks)
# Which of the first 10 PCs best match each role (all time points)
corr10 = pd.DataFrame({key: [np.corrcoef(pcs[:, i], V[key])[0, 1] for i in range(10)] for key in ("fwd", "left", "right", "walk")},
index=[f"PC{i + 1}" for i in range(10)])
corr10["left - right"] = corr10["left"] - corr10["right"]
print("\nCorrelations over all time points, PCs 1-10:")
print(corr10.round(3).to_string())
# PCs by fly: correlation of PC1-3 with the components within each fly
rows = []
for j, f in enumerate(flies):
sl = slice(j * 3384, (j + 1) * 3384)
rows.append({"fly": f, **{f"PC{i + 1}~{key}": np.corrcoef(pcs[sl, i], V[key][sl])[0, 1]
for i in range(3) for key in ("fwd", "left", "right")}})
print("\nPer fly:")
print(pd.DataFrame(rows).round(2).to_string(index=False))
fig, axes = plt.subplots(1, 3, figsize=(12, 3.2))
for i, ax in enumerate(axes):
ax.hexbin(V["Y"] * bw.BALL_RADIUS_MM * bw.FICTRAC_HZ, V["Z"] * 180 / np.pi * bw.FICTRAC_HZ, C=pcs[:, i],
gridsize=30, extent=(-5, 15, -300, 300), reduce_C_function=np.mean, cmap="viridis")
ax.set_title(f"mean PC{i + 1} in velocity space"); ax.set_xlabel("forward (mm/s)"); ax.set_ylabel("rotation (deg/s, + left)")
plt.show()
fv, rv = V["Y"] * bw.BALL_RADIUS_MM * bw.FICTRAC_HZ, V["Z"] * 180 / np.pi * bw.FICTRAC_HZ
print("Points outside the axes (forward -5 to 15 mm/s, rotation -300 to 300 deg/s):",
int(((fv < -5) | (fv > 15) | (np.abs(rv) > 300)).sum()), "of", len(fv))
reading PC component r below p below n below r above p above n above
a 1 fwd 0.202 3.41e-215 23614 0.162 1.13e-41 6842
a 1 left 0.128 6.5e-95 26049 0.0716 1.97e-06 4407
a 1 right 0.155 3.62e-138 25819 -0.00999 0.497 4637
a 2 fwd 0.065 1.58e-23 23614 0.272 2.29e-116 6842
a 2 left -0.0108 0.0823 26049 0.0909 1.47e-09 4407
a 2 right 0.0364 5.08e-09 25819 0.0872 2.67e-09 4637
a 3 fwd -0.0447 6.56e-12 23614 0.0864 7.97e-13 6842
a 3 left 0.0447 5.24e-13 26049 0.114 3.47e-14 4407
a 3 right -0.0665 9.87e-27 25819 -0.259 5.61e-72 4637
b 1 fwd 0.273 0 25500 0.166 8.07e-32 4956
b 1 left 0.165 2.89e-164 27054 0.0607 0.000397 3402
b 1 right 0.186 4.78e-208 26846 0.006 0.719 3610
b 2 fwd 0.0844 1.52e-41 25500 0.282 2.53e-91 4956
b 2 left 0.0279 4.48e-06 27054 0.0958 2.2e-08 3402
b 2 right 0.0598 1.04e-22 26846 0.0822 7.57e-07 3610
b 3 fwd -0.0563 2.28e-19 25500 0.0818 8.09e-09 4956
b 3 left 0.04 4.91e-11 27054 0.0907 1.17e-07 3402
b 3 right -0.096 5.33e-56 26846 -0.25 1.13e-52 3610
Criteria by reading: {'a': {'(i) PC1 below-threshold': True, '(ii) PC2 forward above': np.True_, '(iii) PC3 turn sign': np.True_}, 'b': {'(i) PC1 below-threshold': True, '(ii) PC2 forward above': np.True_, '(iii) PC3 turn sign': np.True_}}
Correlations over all time points, PCs 1-10:
fwd left right walk left - right
PC1 0.361 0.222 0.188 0.415 0.034
PC2 0.249 0.108 0.131 0.188 -0.023
PC3 0.002 0.110 -0.257 -0.087 0.367
PC4 -0.147 -0.204 -0.043 -0.166 -0.161
PC5 -0.015 -0.133 -0.121 -0.161 -0.012
PC6 0.090 -0.087 0.172 0.080 -0.259
PC7 0.120 0.009 0.026 0.034 -0.017
PC8 -0.091 -0.005 -0.030 -0.043 0.025
PC9 0.044 -0.023 -0.018 -0.019 -0.005
PC10 -0.067 -0.008 0.005 -0.004 -0.013
Per fly:
fly PC1~fwd PC1~left PC1~right PC2~fwd PC2~left PC2~right PC3~fwd PC3~left PC3~right
fly_087 0.36 0.13 0.16 0.30 0.09 0.13 0.18 0.30 -0.31
fly_089 0.41 0.35 0.20 0.12 0.07 0.07 -0.23 -0.07 -0.21
fly_094 0.57 0.36 0.24 0.35 0.05 0.12 -0.29 -0.02 -0.24
fly_097 0.22 0.35 0.38 0.23 0.30 0.23 -0.16 0.00 -0.27
fly_098 0.37 0.33 0.17 0.20 0.25 0.13 0.24 0.30 0.02
fly_099 0.42 0.16 0.39 0.31 0.11 0.36 -0.31 -0.10 -0.48
fly_100 0.41 0.20 0.21 0.18 0.12 0.02 0.07 0.08 -0.18
fly_101 0.29 0.10 0.17 0.33 0.06 0.15 -0.15 0.11 -0.38
fly_105 0.33 0.30 0.27 0.23 0.12 0.06 0.09 0.10 -0.01
Points outside the axes (forward -5 to 15 mm/s, rotation -300 to 300 deg/s): 392 of 30456
Result: reproduced. All three parts hold under both readings of the threshold. Under reading (a), PC1's below-threshold correlation with VF is 0.202, against 0.065 for PC2 and -0.045 for PC3. Above the threshold PC2 has the largest correlation with VF (0.272; PC1 0.162, PC3 0.086), with 0.091 for VL and 0.087 for VR. PC3's above-threshold correlations are 0.114 with VL and -0.259 with VR (p = 3.5e-14 and 5.6e-72). PC1 explains 2.22% of the block variance. Within each fly, over all time points, PC1 correlates positively with VF (0.22 to 0.57) and PC3 more with VL than with VR, in 9 of 9 flies.
fig4fg-ridge-r2¶
Linear models on the principal components predict forward, left and right velocity on held-out data with average R2 of about 0.3 using 100 PCs, performance improving over the first 100 PCs before declining.
ALPHAS = [1e-7, 1e-6, 1e-5, 1e-4, 1e-3, 1e-2, 1e-1, 1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7] # pca_fit_behavior_models.py
N_PC_LIST = [1, 10, 30, 100, 250, 500]
folds = bw.contiguous_folds(3384, 5)
t0 = time.time()
def fly_scores(scores_source):
rows, best_alpha = [], {}
for j, f in enumerate(flies):
P_f = scores_source(j)
Y = np.column_stack([beh[f][key] for key in ("fwd", "left", "right")])
for n in N_PC_LIST:
sc_all = bw.ridge_scores(P_f[:, :n], Y, folds, ALPHAS).mean(0) # alpha, component
for c, key in enumerate(("fwd", "left", "right")):
sc = sc_all[:, c]
rows.append({"fly": f, "component": key, "n_pcs": n, "R2": sc.max(), "alpha": ALPHAS[int(np.argmax(sc))]})
best_alpha[(f, key, n)] = ALPHAS[int(np.argmax(sc))]
return pd.DataFrame(rows), best_alpha
ridge, best_alpha = fly_scores(lambda j: pcs[j * 3384:(j + 1) * 3384])
r2_100 = ridge[ridge.n_pcs == 100].pivot(index="fly", columns="component", values="R2")[["fwd", "left", "right"]]
print(f"Fits done in {time.time() - t0:.0f} s; R2 with 100 superfly PCs, per fly:")
print(r2_100.round(3).to_string())
summary_r2 = pd.DataFrame({"(a) median": r2_100.median(), "(b) mean": r2_100.mean()}).T
print("\nAverage over flies (expected 0.3, 0.27 to 0.33):")
print(summary_r2.round(3).to_string())
print("\nMean R2 over flies by number of PCs:")
curve = ridge.groupby(["component", "n_pcs"]).R2.mean().unstack()
print(curve.round(3).to_string())
print("Chosen alphas at 100 PCs:", ridge[ridge.n_pcs == 100].alpha.value_counts().to_dict())
Fits done in 11 s; R2 with 100 superfly PCs, per fly:
component fwd left right
fly
fly_087 -0.307 0.196 0.339
fly_089 0.312 0.246 0.097
fly_094 0.534 0.362 0.026
fly_097 0.125 0.144 0.292
fly_098 -0.306 0.214 0.008
fly_099 0.156 0.016 0.373
fly_100 0.155 0.136 0.212
fly_101 0.174 0.231 0.266
fly_105 -0.503 -0.150 -0.014
Average over flies (expected 0.3, 0.27 to 0.33):
component fwd left right
(a) median 0.155 0.196 0.212
(b) mean 0.038 0.155 0.178
Mean R2 over flies by number of PCs:
n_pcs 1 10 30 100 250 500
component
fwd -0.105 -0.007 0.022 0.038 0.036 0.033
left 0.020 0.099 0.142 0.155 0.147 0.139
right 0.010 0.127 0.156 0.178 0.173 0.162
Chosen alphas at 100 PCs: {100.0: 14, 1000.0: 11, 10000.0: 2}
# Variant 1: alpha chosen within the training folds (inner contiguous 4-fold), not on the test fold.
t0 = time.time()
rows = []
for j, f in enumerate(flies):
P_f = pcs[j * 3384:(j + 1) * 3384, :100]
for key in ("fwd", "left", "right"):
y = beh[f][key]
sc = []
for test in folds:
train = np.setdiff1d(np.arange(3384), test)
inner = [train[ix] for ix in bw.contiguous_folds(len(train), 4)]
inner_sc = bw.ridge_scores(P_f[train], y[train], [np.searchsorted(train, ix) for ix in inner], ALPHAS)
a = ALPHAS[int(np.argmax(inner_sc.mean(0)))]
sc.append(bw.ridge_scores(P_f, y, [test], [a])[0, 0])
rows.append({"fly": f, "component": key, "R2": np.mean(sc)})
nested = pd.DataFrame(rows).pivot(index="fly", columns="component", values="R2")[["fwd", "left", "right"]]
print(f"Variant, nested alpha ({time.time() - t0:.0f} s): median {nested.median().round(3).to_dict()}, mean {nested.mean().round(3).to_dict()}")
# Variant 2: PCs of each fly alone (Figure 4G dashed lines)
t0 = time.time()
own = []
for j, f in enumerate(flies):
Xf = X[j * 3384:(j + 1) * 3384]
Xf = Xf - Xf.mean(0)
g_vals, g_vecs = np.linalg.eigh((Xf @ Xf.T).astype(np.float64)) # Gram matrix: same scores as the covariance PCs
order = np.argsort(g_vals)[::-1][:100]
own.append((g_vecs[:, order] * np.sqrt(np.clip(g_vals[order], 0, None))).astype(np.float32))
own_r2 = []
for j, f in enumerate(flies):
sc_all = bw.ridge_scores(own[j], np.column_stack([beh[f][k] for k in ("fwd", "left", "right")]), folds, ALPHAS).mean(0)
for c, key in enumerate(("fwd", "left", "right")):
own_r2.append({"fly": f, "component": key, "R2": sc_all[:, c].max()})
own_r2 = pd.DataFrame(own_r2).pivot(index="fly", columns="component", values="R2")[["fwd", "left", "right"]]
print(f"Variant, each fly's own 100 PCs ({time.time() - t0:.0f} s): median {own_r2.median().round(3).to_dict()}, mean {own_r2.mean().round(3).to_dict()}")
# Variant 3: interleaved folds (50 segments of the session, segment i in fold i mod 5)
seg = np.array_split(np.arange(3384), 50)
inter = [np.concatenate(seg[k::5]) for k in range(5)]
rows = []
for j, f in enumerate(flies):
sc_all = bw.ridge_scores(pcs[j * 3384:(j + 1) * 3384, :100], np.column_stack([beh[f][k] for k in ("fwd", "left", "right")]), inter, ALPHAS).mean(0)
rows.append({"fly": f, **{key: sc_all[:, c].max() for c, key in enumerate(("fwd", "left", "right"))}})
inter_r2 = pd.DataFrame(rows).set_index("fly")
print(f"Variant, interleaved folds: median {inter_r2.median().round(3).to_dict()}, mean {inter_r2.mean().round(3).to_dict()}")
# Reference: a model with no information (behaviour shifted circularly by half the session) gives
rng_ridge = np.random.default_rng(1)
null = []
for j, f in enumerate(flies):
for key in ("fwd", "left", "right"):
y = np.roll(beh[f][key], int(rng_ridge.integers(800, 2584)))
null.append(bw.ridge_scores(pcs[j * 3384:(j + 1) * 3384, :100], y, folds, ALPHAS).mean(0).max())
print(f"Reference, behaviour rolled by 800 to 2584 volumes: R2 median {np.median(null):.3f}, max {np.max(null):.3f} over 27 fits")
Variant, nested alpha (9 s): median {'fwd': 0.134, 'left': 0.18, 'right': 0.212}, mean {'fwd': -0.009, 'left': 0.153, 'right': 0.177}
Variant, each fly's own 100 PCs (33 s): median {'fwd': 0.151, 'left': 0.21, 'right': 0.217}, mean {'fwd': 0.04, 'left': 0.148, 'right': 0.173}
Variant, interleaved folds: median {'fwd': 0.318, 'left': 0.252, 'right': 0.271}, mean {'fwd': 0.329, 'left': 0.221, 'right': 0.243}
Reference, behaviour rolled by 800 to 2584 volumes: R2 median -0.023, max -0.002 over 27 fits
Result: partially reproduced. With 100 superfly PCs, the median R2 over flies is 0.155 for VF, 0.196 for VL and 0.212 for VR (reading a; means 0.038, 0.155 and 0.178), all below 0.27 to 0.33 and all at least 0.15. The mean R2 rises up to 100 PCs and declines after it for each component (VF 0.022, 0.038, 0.036 and 0.033 at 30, 100, 250 and 500 PCs), as the paper describes. Three flies have negative VF R2 (-0.307, -0.306 and -0.503). With interleaved folds in place of contiguous ones the medians are 0.318, 0.252 and 0.271; with alpha chosen within the training folds 0.134, 0.180 and 0.212; with each fly's own PCs 0.151, 0.210 and 0.217. With the behaviour rolled in time, R2 is at most -0.002.
fig4h-r2-vs-movement¶
Flies that moved less along a velocity component were poorly predicted (R2 near 0), while flies that explored it more were well predicted, with R2 up to 0.6.
amount = pd.DataFrame({f: {"fwd": beh[f]["fwd"].mean() * bw.BALL_RADIUS_MM * bw.FICTRAC_HZ,
"left": beh[f]["left"].mean() * 180 / np.pi * bw.FICTRAC_HZ,
"right": beh[f]["right"].mean() * 180 / np.pi * bw.FICTRAC_HZ} for f in flies}).T
tab = pd.concat({"mean amount": amount, "R2": r2_100}, axis=1)
print("Per fly: mean VF (mm/s), VL and VR (deg/s), and R2 with 100 PCs")
print(tab.round(3).to_string())
rho = {key: stats.spearmanr(amount[key], r2_100[key]) for key in ("fwd", "left", "right")}
print("\nSpearman per component:", {k: (round(v.statistic, 3), round(v.pvalue, 4)) for k, v in rho.items()})
ra = np.concatenate([stats.rankdata(amount[k]) for k in ("fwd", "left", "right")])
rr = np.concatenate([stats.rankdata(r2_100[k]) for k in ("fwd", "left", "right")])
pooled = stats.pearsonr(ra, rr)
print(f"Within-component ranks pooled over 27 points: r = {pooled.statistic:.3f}, p = {pooled.pvalue:.3g}")
print(f"Largest per-fly R2 {r2_100.values.max():.3f} ({r2_100.stack().idxmax()}); smallest {r2_100.values.min():.3f} ({r2_100.stack().idxmin()})")
print("Expected largest 0.6 (0.54 to 0.66), smallest below 0.1")
Per fly: mean VF (mm/s), VL and VR (deg/s), and R2 with 100 PCs
mean amount R2
fwd left right fwd left right
fly_087 2.311 21.477 47.205 -0.307 0.196 0.339
fly_089 0.889 20.157 9.499 0.312 0.246 0.097
fly_094 1.008 13.409 4.634 0.534 0.362 0.026
fly_097 0.212 7.913 14.909 0.125 0.144 0.292
fly_098 1.186 21.620 21.833 -0.306 0.214 0.008
fly_099 0.784 5.549 40.122 0.156 0.016 0.373
fly_100 2.166 24.644 44.208 0.155 0.136 0.212
fly_101 1.584 30.171 22.248 0.174 0.231 0.266
fly_105 0.435 11.075 7.893 -0.503 -0.150 -0.014
Spearman per component: {'fwd': (np.float64(0.017), np.float64(0.9661)), 'left': (np.float64(0.383), np.float64(0.3085)), 'right': (np.float64(0.683), np.float64(0.0424))}
Within-component ranks pooled over 27 points: r = 0.361, p = 0.0642
Largest per-fly R2 0.534 (('fly_094', 'fwd')); smallest -0.503 (('fly_105', 'fwd'))
Expected largest 0.6 (0.54 to 0.66), smallest below 0.1
Result: partially reproduced. The Spearman correlation between a fly's mean amount of a component and its R2 is positive for all three components (VF 0.017, VL 0.383, VR 0.683; p = 0.042 for VR only), but the pooled correlation of within-component ranks is 0.361 with p = 0.064, short of the p < 0.05 the criterion requires. The largest R2 is 0.534 (fly_094, VF), just below the band of 0.54 to 0.66, and the smallest is -0.503, below 0.1.
fig4i-weights-stereotyped¶
The model weights given to each PC are stereotyped across flies; PC1 is used for all components, PC2 mainly for forward velocity, and PC3 for left and right turning with opposite signs.
W = {}
for j, f in enumerate(flies):
P_f = pcs[j * 3384:(j + 1) * 3384, :100]
for key in ("fwd", "left", "right"):
y = beh[f][key] / beh[f][key].std()
W[(f, key)] = bw.ridge_fit(P_f, y, best_alpha[(f, key, 100)])[:20]
pair_r = {key: np.median([np.corrcoef(W[(a, key)], W[(b, key)])[0, 1] for i, a in enumerate(flies) for b in flies[i + 1:]])
for key in ("fwd", "left", "right")}
pc3_opposite = sum(np.sign(W[(f, "left")][2]) != np.sign(W[(f, "right")][2]) for f in flies)
pc2_fwd = sum(abs(W[(f, "fwd")][1]) > max(abs(W[(f, "left")][1]), abs(W[(f, "right")][1])) for f in flies)
pc1_same = sum(len({np.sign(W[(f, k)][0]) for k in ("fwd", "left", "right")}) == 1 for f in flies)
print("Median pairwise r of the 20 weights over 36 fly pairs:", {k: round(v, 3) for k, v in pair_r.items()})
print("PC3 weights of opposite sign for VL and VR:", pc3_opposite, "of 9 flies")
print("PC2 weight largest for VF:", pc2_fwd, "of 9 flies")
print("PC1 weights of one sign for all three:", pc1_same, "of 9 flies")
wt = pd.DataFrame({(f, k): W[(f, k)][:3] for f in flies for k in ("fwd", "left", "right")}, index=["PC1", "PC2", "PC3"]).T
print("\nWeights on PCs 1-3 (behaviour in s.d. units):")
print(wt.round(4).to_string())
fi_checks = {"(i) median pair r >= 0.5 for each": all(v >= 0.5 for v in pair_r.values()),
"(ii) PC3 opposite in >= 7": pc3_opposite >= 7, "(iii) PC2 forward in >= 7": pc2_fwd >= 7,
"(iv) PC1 same sign in >= 7": pc1_same >= 7}
print("\nCriteria:", fi_checks)
# Null: each fly's behaviour rolled by 800 to 2584 volumes, refitted at the same alpha; 20 rolls
rng_w = np.random.default_rng(3)
null_pair, null_pc3, null_pc1 = {k: [] for k in ("fwd", "left", "right")}, [], []
for _ in range(20):
Wn = {}
for j, f in enumerate(flies):
P_f = pcs[j * 3384:(j + 1) * 3384, :100]
shift = int(rng_w.integers(800, 2584))
for key in ("fwd", "left", "right"):
y = np.roll(beh[f][key], shift)
Wn[(f, key)] = bw.ridge_fit(P_f, y / y.std(), best_alpha[(f, key, 100)])[:20]
for key in ("fwd", "left", "right"):
null_pair[key].append(np.median([np.corrcoef(Wn[(a, key)], Wn[(b, key)])[0, 1] for i, a in enumerate(flies) for b in flies[i + 1:]]))
null_pc3.append(sum(np.sign(Wn[(f, "left")][2]) != np.sign(Wn[(f, "right")][2]) for f in flies))
null_pc1.append(sum(len({np.sign(Wn[(f, k)][0]) for k in ("fwd", "left", "right")}) == 1 for f in flies))
print("Null, behaviour rolled (20 rolls): median pair r, largest over rolls",
{k: round(float(np.max(v)), 3) for k, v in null_pair.items()}, "; mean", {k: round(float(np.mean(v)), 3) for k, v in null_pair.items()})
print("Null: PC3 opposite in", min(null_pc3), "to", max(null_pc3), "of 9 flies; PC1 one sign in", min(null_pc1), "to", max(null_pc1), "of 9")
Median pairwise r of the 20 weights over 36 fly pairs: {'fwd': np.float64(0.646), 'left': np.float64(0.732), 'right': np.float64(0.761)}
PC3 weights of opposite sign for VL and VR: 9 of 9 flies
PC2 weight largest for VF: 7 of 9 flies
PC1 weights of one sign for all three: 9 of 9 flies
Weights on PCs 1-3 (behaviour in s.d. units):
PC1 PC2 PC3
fly_087 fwd 0.0724 0.1296 0.0318
left 0.0165 0.0531 0.1378
right 0.0506 0.0657 -0.2708
fly_089 fwd 0.2108 0.1182 -0.0178
left 0.1213 0.0773 0.1345
right 0.0692 0.0852 -0.1170
fly_094 fwd 0.1408 0.2747 -0.1200
left 0.0885 0.0110 0.1589
right 0.0536 0.0932 -0.1425
fly_097 fwd 0.0790 0.1800 -0.1858
left 0.1238 0.1446 0.0225
right 0.1557 0.1157 -0.3300
fly_098 fwd 0.0956 0.0278 0.0176
left 0.0726 0.1148 0.1971
right 0.0692 0.0718 -0.1145
fly_099 fwd 0.1646 0.1369 -0.0196
left 0.0663 0.0542 0.0187
right 0.0578 0.1125 -0.3551
fly_100 fwd 0.1918 0.2550 0.1024
left 0.0824 0.1239 0.1872
right 0.0714 0.0794 -0.3953
fly_101 fwd 0.0973 0.2180 -0.0766
left 0.0476 0.0386 0.2469
right 0.0679 0.1894 -0.4436
fly_105 fwd 0.0895 0.0434 0.0075
left 0.1200 0.0869 0.0377
right 0.1137 0.0492 -0.0500
Criteria: {'(i) median pair r >= 0.5 for each': True, '(ii) PC3 opposite in >= 7': np.True_, '(iii) PC2 forward in >= 7': np.True_, '(iv) PC1 same sign in >= 7': True}
Null, behaviour rolled (20 rolls): median pair r, largest over rolls {'fwd': 0.151, 'left': 0.102, 'right': 0.089} ; mean {'fwd': 0.021, 'left': -0.001, 'right': -0.021}
Null: PC3 opposite in 3 to 8 of 9 flies; PC1 one sign in 0 to 6 of 9
Result: reproduced. The median correlation of the weights on PCs 1 to 20 between pairs of flies is 0.646 for VF, 0.732 for VL and 0.761 for VR, against at most 0.151 over 20 rolls of each fly's behaviour in time. PC3's weight is positive for VL and negative for VR in 9 of 9 flies (3 to 8 of 9 with rolled behaviour); PC2's weight is largest for VF in 7 of 9; PC1's weights have one sign for all three components in 9 of 9 (0 to 6 of 9 with rolled behaviour). Part (iv) follows in part from the orientation of PC1, which is signed to correlate positively with VF and correlates positively with VL and VR over all time points (0.222 and 0.188).
fig5a-unique-variance¶
Demixing with linear models shows that forward velocity is uniquely encoded by a spatially organized subset of voxels, while locomotor state (walking or stopped) is not uniquely encoded anywhere, although the single-variable maps of forward velocity and locomotor state are very similar.
Each behaviour variable is divided by its s.d. before fitting, so that the fixed ridge penalty (alpha 1) is small next to the 30,456 time points.
def design(walk_key="walk"):
cols = [V["fwd"], V["left"], V["right"], V[walk_key]]
return np.column_stack([c / c.std() for c in cols])
def demix(D, names, alpha=1.0):
"""Mean held-out R2 per block for every model; fold k is the k-th fifth of every fly."""
fifths = bw.contiguous_folds(3384, 5)
fold_idx = [np.concatenate([j * 3384 + fifths[k] for j in range(len(flies))]) for k in range(5)]
models = {"all": list(range(4))}
models.update({f"without {n}": [i for i in range(4) if i != k] for k, n in enumerate(names)})
models.update({f"only {n}": [k] for k, n in enumerate(names)})
r2 = {m: np.zeros(X.shape[1]) for m in models}
every = np.arange(X.shape[0])
for test in fold_idx:
train = np.setdiff1d(every, test)
dm = D[train].mean(0)
dc = D[train] - dm
dtd = dc.T @ dc
dty = dc.T @ X[train] # X is centred; dc sums to 0, so no need to centre X here
ym = X[train].mean(0)
yt = X[test]
sst = ((yt - yt.mean(0)) ** 2).sum(0)
dt = D[test] - dm
for m, cols in models.items():
coef = np.linalg.solve(dtd[np.ix_(cols, cols)] + alpha * np.eye(len(cols)), dty[cols])
pred = dt[:, cols] @ coef + ym
r2[m] += 1 - ((yt - pred) ** 2).sum(0) / sst
return {m: v / len(fold_idx) for m, v in r2.items()}
names = ["VF", "VL", "VR", "walk"]
t0 = time.time()
r2_models = demix(design("walk"), names)
unique = {n: r2_models["all"] - r2_models[f"without {n}"] for n in names}
print(f"Demixing of {X.shape[1]} blocks in {time.time() - t0:.0f} s")
UNIQUE_MIN = 0.01
rows = []
for n in names:
rows.append({"variable": n, "unique: blocks >= 0.01": int((unique[n] >= UNIQUE_MIN).sum()),
"unique: % of blocks": 100 * (unique[n] >= UNIQUE_MIN).mean(),
"unique: max": unique[n].max(), "unique: 99th pct": np.percentile(unique[n], 99),
"single: median R2": np.median(r2_models[f"only {n}"]), "single: max R2": r2_models[f"only {n}"].max()})
demix_table = pd.DataFrame(rows)
print(demix_table.round(4).to_string(index=False))
print(f"All four: median R2 {np.median(r2_models['all']):.4f}, max {r2_models['all'].max():.4f}; blocks with R2 <= 0: {(r2_models['all'] <= 0).sum()}")
print("Correlation over blocks of single-variable R2, VF against walking state:",
round(np.corrcoef(r2_models["only VF"], r2_models["only walk"])[0, 1], 3))
for thr in (0.005, 0.01, 0.02):
print(f" threshold {thr}: % of blocks encoding uniquely", {n: round(100 * (unique[n] >= thr).mean(), 2) for n in names})
fa = (unique["walk"] >= UNIQUE_MIN).mean(); ff = (unique["VF"] >= UNIQUE_MIN).mean()
print(f"\nCriteria: walking state <= 1%: {fa <= 0.01}; VF >= 5%: {ff >= 0.05}; VF >= 5 x walking: {ff >= 5 * fa}")
# Variant: the walking state as the code intended it (norm of both standardised components)
V["walk_intended"] = np.concatenate([(np.interp(times[f], bw.smoothed_velocity(by_fly[f])["t"],
bw.walking_state(bw.smoothed_velocity(by_fly[f])["Y"], bw.smoothed_velocity(by_fly[f])["Z"], as_intended=True)) >= 0.5).astype(float)
for f in flies])
r2_int = demix(design("walk_intended"), names)
u_int = {n: r2_int["all"] - r2_int[f"without {n}"] for n in names}
print("Variant, walking state as intended: % of blocks encoding uniquely",
{n: round(100 * (u_int[n] >= UNIQUE_MIN).mean(), 2) for n in names})
print("Walking fraction of time points: as run", round(V["walk"].mean(), 3), "; as intended", round(V["walk_intended"].mean(), 3))
# Per fly: the same model within each fly
rows = []
for j, f in enumerate(flies):
sl = slice(j * 3384, (j + 1) * 3384)
Df = design("walk")[sl]
Xf = X[sl] - X[sl].mean(0)
fifths = bw.contiguous_folds(3384, 5)
full = bw.multi_ridge_r2(Df, Xf, fifths, 1.0)
rows.append({"fly": f, **{f"{n} unique %": 100 * ((full - bw.multi_ridge_r2(Df[:, [i for i in range(4) if i != k]], Xf, fifths, 1.0)) >= UNIQUE_MIN).mean()
for k, n in enumerate(names)}})
print("\nPer fly (within-fly models):")
print(pd.DataFrame(rows).round(2).to_string(index=False))
Demixing of 5363 blocks in 10 s
variable unique: blocks >= 0.01 unique: % of blocks unique: max unique: 99th pct single: median R2 single: max R2
VF 40 0.7459 0.0315 0.0084 0.0010 0.0628
VL 4 0.0746 0.0117 0.0040 0.0003 0.0259
VR 8 0.1492 0.0159 0.0055 0.0001 0.0268
walk 5 0.0932 0.0119 0.0061 0.0015 0.0403
All four: median R2 0.0019, max 0.0651; blocks with R2 <= 0: 809
Correlation over blocks of single-variable R2, VF against walking state: 0.882
threshold 0.005: % of blocks encoding uniquely {'VF': np.float64(2.24), 'VL': np.float64(0.65), 'VR': np.float64(1.21), 'walk': np.float64(1.88)}
threshold 0.01: % of blocks encoding uniquely {'VF': np.float64(0.75), 'VL': np.float64(0.07), 'VR': np.float64(0.15), 'walk': np.float64(0.09)}
threshold 0.02: % of blocks encoding uniquely {'VF': np.float64(0.09), 'VL': np.float64(0.0), 'VR': np.float64(0.0), 'walk': np.float64(0.0)}
Criteria: walking state <= 1%: True; VF >= 5%: False; VF >= 5 x walking: True
Variant, walking state as intended: % of blocks encoding uniquely {'VF': np.float64(0.91), 'VL': np.float64(0.02), 'VR': np.float64(0.04), 'walk': np.float64(0.97)}
Walking fraction of time points: as run 0.472 ; as intended 0.573
Per fly (within-fly models):
fly VF unique % VL unique % VR unique % walk unique %
fly_087 4.77 1.03 3.11 0.00
fly_089 0.19 0.28 0.13 0.24
fly_094 2.74 0.88 0.00 8.35
fly_097 0.00 4.01 5.26 0.00
fly_098 0.02 1.04 0.07 0.52
fly_099 1.53 0.00 4.27 0.00
fly_100 0.58 0.37 0.82 0.06
fly_101 0.34 0.21 0.52 1.32
fly_105 0.02 0.06 0.06 0.07
Result: partially reproduced. At the 0.01 threshold, 5 blocks (0.09%) encode walking state uniquely and 40 (0.75%) encode VF uniquely, 8 times as many, but below the 5% that the criterion requires. Unique contributions are small everywhere (largest 0.032 for VF and 0.012 for walking state), and the full model has a median held-out R2 of 0.0019 per block, so the 0.01 threshold sits far above most blocks. The single-variable R2 maps of VF and walking state correlate at 0.882 over blocks, the similarity the paper describes. The VF advantage depends on the threshold and on the walking state: at a threshold of 0.005, walking state is unique in 1.88% of blocks and VF in 2.24%; with the walking state the code intended, walking state is unique in 0.97% and VF in 0.91%. Within single flies, walking state is unique in up to 8.35% of blocks (fly_094).
fig5c-global-below-threshold¶
The functional regions defined by the unique contributions of forward, left and right velocity are all correlated with velocity when the fly is below the movement threshold, a widespread increase in activity at the initiation of movement, while above the threshold each region correlates most with its own velocity component.
regions = {key: unique[n] >= UNIQUE_MIN for key, n in (("fwd", "VF"), ("left", "VL"), ("right", "VR"))}
print("Blocks per region:", {k: int(v.sum()) for k, v in regions.items()}, "; in two or more regions:",
int((sum(v.astype(int) for v in regions.values()) > 1).sum()))
reg_sig = {k: X[:, m].mean(1) for k, m in regions.items() if m.any()}
rows, checks5 = [], {}
for reading in ("a", "b"):
below_ok, above_ok = [], []
for rk, sig in reg_sig.items():
best = None
for key in ("fwd", "left", "right"):
sc = split_corr(sig, V[key], reading)
rows.append({"reading": reading, "region": rk, "component": key, "r below": sc["below"][0], "p below": sc["below"][1],
"r above": sc["above"][0], "p above": sc["above"][1]})
below_ok.append(sc["below"][0] > 0 and sc["below"][1] < 0.05)
if best is None or sc["above"][0] > best[1]:
best = (key, sc["above"][0])
above_ok.append(best[0] == rk)
checks5[reading] = {"all below-threshold r > 0, p < 0.05": all(below_ok) and len(below_ok) == 9,
"own component largest above": all(above_ok) and len(above_ok) == 3}
print(pd.DataFrame(rows).to_string(index=False, float_format=lambda v: f"{v:.3g}"))
print("\nCriteria by reading:", checks5)
rows = []
for j, f in enumerate(flies):
sl = slice(j * 3384, (j + 1) * 3384)
row = {"fly": f}
for rk, sig in reg_sig.items():
for key in ("fwd", "left", "right"):
comp = V[key][sl]
m = comp / V[key].std() <= 0.75
row[f"{rk}~{key} below"] = np.corrcoef(sig[sl][m], comp[m])[0, 1]
rows.append(row)
# Null for the below-threshold correlations: each fly's behaviour rolled circularly by 300 to 3084
# volumes (a new shift per fly and roll), 200 rolls; the threshold is applied to the rolled component.
rng_5c = np.random.default_rng(4)
shifts = rng_5c.integers(300, 3084, size=(200, len(flies)))
def rolled(v, sh):
return np.concatenate([np.roll(v[j * 3384:(j + 1) * 3384], s_) for j, s_ in enumerate(sh)])
rows_null = []
for reading in ("a", "b"):
for rk, sig in reg_sig.items():
for key in ("fwd", "left", "right"):
obs = split_corr(sig, V[key], reading)["below"][0]
nul = np.array([split_corr(sig, rolled(V[key], sh), reading)["below"][0] for sh in shifts])
rows_null.append({"reading": reading, "region": rk, "component": key, "r below": obs,
"null 95th pct": np.percentile(nul, 95), "null max": nul.max(),
"p (rolls >= r)": (1 + (nul >= obs).sum()) / (1 + len(nul))})
null5c = pd.DataFrame(rows_null)
print("\nBelow-threshold r against 200 rolls of the behaviour:")
print(null5c.to_string(index=False, float_format=lambda v: f"{v:.3g}"))
for reading in ("a", "b"):
sub = null5c[null5c.reading == reading]
print(f"Reading ({reading}): below-threshold r > 0 with roll p < 0.05 in {int(((sub['r below'] > 0) & (sub['p (rolls >= r)'] < 0.05)).sum())} of {len(sub)}")
print("\nPer fly, below-threshold r (reading a):")
print(pd.DataFrame(rows).round(2).to_string(index=False))
Blocks per region: {'fwd': 40, 'left': 4, 'right': 8} ; in two or more regions: 0
reading region component r below p below r above p above
a fwd fwd 0.157 9.7e-131 0.365 3.34e-215
a fwd left 0.0814 1.54e-39 0.0543 0.000312
a fwd right 0.122 6.9e-86 -0.0323 0.0281
a left fwd 0.0851 3.33e-39 -0.0026 0.83
a left left 0.0736 1.37e-32 0.196 1.56e-39
a left right -0.00798 0.2 -0.00171 0.907
a right fwd 0.134 1.27e-95 0.102 2.18e-17
a right left 0.0258 3.11e-05 0.0642 1.99e-05
a right right 0.125 5.97e-90 0.269 1.95e-77
b fwd fwd 0.221 1.88e-278 0.377 6.93e-167
b fwd left 0.117 2.84e-83 0.0466 0.0065
b fwd right 0.145 3.19e-126 -0.0249 0.134
b left fwd 0.118 5.42e-80 -0.0124 0.383
b left left 0.093 5.38e-53 0.17 2.04e-23
b left right 0.0131 0.0312 0.0101 0.542
b right fwd 0.193 2.02e-211 0.12 1.73e-17
b right left 0.0654 4.99e-27 0.0693 5.23e-05
b right right 0.16 8.13e-154 0.266 1.03e-59
Criteria by reading: {'a': {'all below-threshold r > 0, p < 0.05': False, 'own component largest above': True}, 'b': {'all below-threshold r > 0, p < 0.05': True, 'own component largest above': True}}
Below-threshold r against 200 rolls of the behaviour:
reading region component r below null 95th pct null max p (rolls >= r)
a fwd fwd 0.157 0.0123 0.0209 0.00498
a fwd left 0.0814 0.0107 0.0143 0.00498
a fwd right 0.122 0.0122 0.0205 0.00498
a left fwd 0.0851 0.0111 0.0213 0.00498
a left left 0.0736 0.0119 0.0184 0.00498
a left right -0.00798 0.0099 0.0225 0.886
a right fwd 0.134 0.00975 0.0172 0.00498
a right left 0.0258 0.00959 0.0154 0.00498
a right right 0.125 0.00988 0.015 0.00498
b fwd fwd 0.221 0.0116 0.0151 0.00498
b fwd left 0.117 0.0099 0.0152 0.00498
b fwd right 0.145 0.0108 0.0181 0.00498
b left fwd 0.118 0.0109 0.0175 0.00498
b left left 0.093 0.0109 0.0194 0.00498
b left right 0.0131 0.0103 0.0213 0.0249
b right fwd 0.193 0.0107 0.0253 0.00498
b right left 0.0654 0.0099 0.0144 0.00498
b right right 0.16 0.00952 0.0129 0.00498
Reading (a): below-threshold r > 0 with roll p < 0.05 in 8 of 9
Reading (b): below-threshold r > 0 with roll p < 0.05 in 9 of 9
Per fly, below-threshold r (reading a):
fly fwd~fwd below fwd~left below fwd~right below left~fwd below left~left below left~right below right~fwd below right~left below right~right below
fly_087 0.18 0.07 0.14 0.04 0.10 -0.12 0.11 -0.04 0.13
fly_089 0.19 0.11 0.10 0.15 0.09 0.02 0.22 0.09 0.13
fly_094 0.26 0.22 0.19 0.13 0.15 0.00 0.23 0.19 0.18
fly_097 0.11 0.10 0.17 0.07 0.04 0.05 0.12 0.03 0.16
fly_098 0.13 0.07 0.06 0.06 0.03 0.01 0.09 0.01 0.09
fly_099 0.07 0.03 0.14 0.02 0.04 0.01 0.05 -0.04 0.14
fly_100 0.18 0.07 0.06 0.05 0.09 -0.02 0.13 -0.01 0.09
fly_101 0.23 0.03 0.13 0.20 0.09 0.01 0.22 -0.02 0.14
fly_105 0.21 0.10 0.20 0.07 0.02 0.05 0.15 0.09 0.16
Result: reproduced. Under reading (b), all nine below-threshold correlations of the three regions (40, 4 and 8 blocks) with the three components are positive with p < 0.05 (0.013 to 0.221), and above the threshold each region correlates most with its own component (0.377, 0.170 and 0.266). The Pearson p-values treat time points as independent; against 200 rolls of each fly's behaviour in time, which keep the autocorrelation, the nine are also above chance under reading (b), the weakest (left region with VR, r = 0.013) with p = 0.025 against a null 95th percentile of 0.010. Under reading (a) one of the nine, the left region with VR, is -0.008 (Pearson p = 0.2, roll p = 0.89). The left region is only 4 blocks, and the regions are small because few blocks encode a component uniquely (fig5a-unique-variance).
fig5f-ipsi-contra¶
When the fly turns, neural activity increases in specific regions of the ipsilateral hemisphere and decreases in the same regions of the contralateral hemisphere.
fwd_mm = V["Y"] * bw.BALL_RADIUS_MM * bw.FICTRAC_HZ
rot_deg = V["Z"] * 180 / np.pi * bw.FICTRAC_HZ
slow = (fwd_mm >= 1) & (fwd_mm <= 3) & (np.abs(rot_deg) < 30)
conds = {"left": rot_deg > 100, "right": rot_deg < -100}
print("Time points: slow walking", int(slow.sum()), "; left turns", int(conds["left"].sum()), "; right turns", int(conds["right"].sum()))
print("By fly:", {f: (int(slow[j * 3384:(j + 1) * 3384].sum()), int(conds["left"][j * 3384:(j + 1) * 3384].sum()),
int(conds["right"][j * 3384:(j + 1) * 3384].sum())) for j, f in enumerate(flies)})
col_of = {(p, b): c for c, (p, b) in enumerate(zip(block_plane, block_index))}
partner = np.array([col_of.get((p, int(bw.mirror_block(np.array([b]))[0])), -1) for p, b in zip(block_plane, block_index)])
print("Blocks with a kept mirror partner:", int((partner >= 0).sum()), "of", len(partner))
base = X[slow].mean(0)
turn_res = {}
for turn, m in conds.items():
diff = X[m].mean(0) - base
top = np.flatnonzero(diff >= np.percentile(diff, 90))
top = top[partner[top] >= 0]
pchange = diff[partner[top]]
w = stats.wilcoxon(pchange)
in_a = (block_bx[top] < nbx // 2).mean()
turn_res[turn] = {"top blocks": len(top), "their change, median": float(np.median(diff[top])),
"mirror change, median": float(np.median(pchange)), "mirror change < 0 (%)": 100 * float((pchange < 0).mean()),
"Wilcoxon p": w.pvalue, "top blocks in half A (%)": 100 * in_a,
"all blocks: change range": (round(float(diff.min()), 3), round(float(diff.max()), 3))}
print(pd.DataFrame(turn_res).T.to_string())
opposite = (turn_res["left"]["top blocks in half A (%)"] > 50) != (turn_res["right"]["top blocks in half A (%)"] > 50)
ok = {t: r["mirror change, median"] < 0 and r["Wilcoxon p"] < 0.05 for t, r in turn_res.items()}
print("\nMirror decrease (median < 0, p < 0.05):", ok, "; increases in opposite halves:", opposite)
# Variant: other bounds for turning
for thr in (50, 150):
vals = {}
for turn, m in (("left", rot_deg > thr), ("right", rot_deg < -thr)):
diff = X[m].mean(0) - base
top = np.flatnonzero(diff >= np.percentile(diff, 90)); top = top[partner[top] >= 0]
vals[turn] = (round(float(np.median(diff[partner[top]])), 4), float(f"{stats.wilcoxon(diff[partner[top]]).pvalue:.3g}"))
print(f"Variant, turning threshold {thr} deg/s: mirror change median and p", vals)
# Variant: the contrast from the demixing model (ridge, alpha 1, of fig5a-unique-variance, fitted to all time points)
D = design("walk")
Dc = D - D.mean(0)
coef = np.linalg.solve(Dc.T @ Dc + np.eye(4), Dc.T @ X) # variable, block
for label, change_of in (("predicted turn - slow walking", lambda m: (D[m].mean(0) - D[slow].mean(0)) @ coef),
("VL or VR coefficient alone", lambda m: coef[1] if m is conds["left"] else coef[2])):
vals = {}
for turn, m in conds.items():
diff = change_of(m)
top = np.flatnonzero(diff >= np.percentile(diff, 90)); top = top[partner[top] >= 0]
pc = diff[partner[top]]
vals[turn] = {"top in half A (%)": round(100 * float((block_bx[top] < nbx // 2).mean()), 1),
"mirror median": float(f"{np.median(pc):.3g}"), "mirror < 0 (%)": round(100 * float((pc < 0).mean()), 1),
"Wilcoxon p": float(f"{stats.wilcoxon(pc).pvalue:.3g}")}
print(f"Variant, demixing model, {label}:", vals)
# Per fly
rows = []
for j, f in enumerate(flies):
sl = np.zeros(len(slow), bool); sl[j * 3384:(j + 1) * 3384] = True
row = {"fly": f}
for turn, m in conds.items():
if (m & sl).sum() < 5 or (slow & sl).sum() < 5:
row[turn] = np.nan
continue
diff = X[m & sl].mean(0) - X[slow & sl].mean(0)
top = np.flatnonzero(diff >= np.percentile(diff, 90)); top = top[partner[top] >= 0]
row[turn + " mirror median"] = float(np.median(diff[partner[top]]))
rows.append(row)
print("\nPer fly (flies with at least 5 points in each condition):")
print(pd.DataFrame(rows).round(3).to_string(index=False))
Time points: slow walking 1053 ; left turns 1774 ; right turns 2574
By fly: {'fly_087': (109, 273, 641), 'fly_089': (189, 238, 76), 'fly_094': (135, 143, 10), 'fly_097': (51, 78, 165), 'fly_098': (120, 233, 290), 'fly_099': (109, 25, 491), 'fly_100': (131, 307, 553), 'fly_101': (139, 380, 282), 'fly_105': (70, 97, 66)}
Blocks with a kept mirror partner: 5032 of 5363
top blocks their change, median mirror change, median mirror change < 0 (%) Wilcoxon p top blocks in half A (%) all blocks: change range
left 525 0.035227 0.00635 36.380952 0.0 81.714286 (-0.059, 0.126)
right 527 0.026062 0.002187 42.125237 0.000039 24.857685 (-0.075, 0.084)
Mirror decrease (median < 0, p < 0.05): {'left': False, 'right': False} ; increases in opposite halves: True
Variant, turning threshold 50 deg/s: mirror change median and p {'left': (0.0057, 3.67e-20), 'right': (0.0017, 0.000359)}
Variant, turning threshold 150 deg/s: mirror change median and p {'left': (0.0096, 3.15e-24), 'right': (0.0018, 0.000652)}
Variant, demixing model, predicted turn - slow walking: {'left': {'top in half A (%)': 77.8, 'mirror median': 0.0254, 'mirror < 0 (%)': 6.5, 'Wilcoxon p': 1.64e-83}, 'right': {'top in half A (%)': 17.0, 'mirror median': 0.0111, 'mirror < 0 (%)': 18.8, 'Wilcoxon p': 1.25e-51}}
Variant, demixing model, VL or VR coefficient alone: {'left': {'top in half A (%)': 76.1, 'mirror median': 0.00642, 'mirror < 0 (%)': 8.4, 'Wilcoxon p': 8.28e-80}, 'right': {'top in half A (%)': 18.7, 'mirror median': 0.00642, 'mirror < 0 (%)': 8.0, 'Wilcoxon p': 4.63e-78}}
Per fly (flies with at least 5 points in each condition):
fly left mirror median right mirror median
fly_087 0.015 0.013
fly_089 0.010 0.004
fly_094 0.010 0.047
fly_097 0.074 0.037
fly_098 0.024 0.007
fly_099 0.018 0.009
fly_100 0.018 0.009
fly_101 -0.005 0.003
fly_105 0.009 0.007
Result: not reproduced. Turning raises activity mostly on one side, opposite for the two directions (82% of the top 10% of blocks in half A for left turns, 25% for right turns), but their mirror blocks rise too, by a median 0.0064 for left turns and 0.0022 for right turns, with 36% and 42% of them falling, against slow walking. The same holds with turning thresholds of 50 and 150 deg/s, and in 8 of 9 flies for left turns and 9 of 9 for right turns. The paper derives Figure 5F from the demixing model; with the change predicted by the demixing model of fig5a-unique-variance, the mirror blocks rise as well (medians 0.0254 and 0.0111, 6.5% and 18.8% falling), and so do the VL and VR coefficients of the mirror blocks (median 0.0064, 8.4% and 8.0% negative).
Cross-correlation filters (Figure 6)¶
For each block, the behaviour at the neural time plus each shift (-5 to 5 s in 20 ms steps) is weighted by the block signal and summed over time points and flies (linear_filters.py), notch filtered at the volume rate, and deconvolved with the code's GCaMP6f kernel (deconvolution.py). The lead is the shift at the peak of the deconvolved filter; positive means activity precedes behaviour.
def filters_for(offset_s=0.0):
out = {key: np.zeros((len(bw.SHIFTS_S), X.shape[1])) for key in ("fwd", "left", "right")}
for j, f in enumerate(flies):
Xf = X[j * 3384:(j + 1) * 3384]
for key in out:
out[key] += bw.shifted_behaviour(by_fly[f], key, times[f] + offset_s) @ Xf
return out
t0 = time.time()
raw_filters = filters_for(0.0)
rate = float(np.mean([s["processed"]["rate"] for s in sessions]))
b_notch, a_notch = signal.iirnotch(rate, Q=2, fs=1 / 0.02)
toe = bw.gcamp_toeplitz(len(bw.SHIFTS_S))
def peaks(raw):
res = {}
for key, flt in raw.items():
notched = signal.filtfilt(b_notch, a_notch, flt, axis=0)
dec = bw.deconvolve(notched.T, toe) # block, shift
idx = dec.argmax(1)
top = dec.max(1) >= np.percentile(dec.max(1), 60) # top 40% by peak value
res[key] = {"dec": dec, "lead": bw.SHIFTS_S[idx], "peak": dec.max(1), "top": top,
"lead_raw": bw.SHIFTS_S[notched.T.argmax(1)]}
return res
pk = peaks(raw_filters)
print(f"Filters for {X.shape[1]} blocks in {time.time() - t0:.0f} s")
Filters for 5363 blocks in 11 s
fig6-peak-time-range¶
The time of peak cross-correlation between supervoxel activity and forward velocity ranges from about 300 ms before the change in behaviour, the earliest, to more than 1 s after it.
def lead_summary(res):
rows = []
for key, r in res.items():
L = r["lead"][r["top"]]
rows.append({"component": key, "blocks (top 40%)": int(r["top"].sum()),
"(a) earliest lead": L.max(), "(b) 99th pct lead": np.percentile(L, 99),
"median lead": np.median(L), "(b) 1st pct lead": np.percentile(L, 1), "(a) latest lead": L.min(),
"at window edge": int((np.abs(L) >= 4.9).sum()),
"median lead before deconvolution": np.median(r["lead_raw"][r["top"]])})
return pd.DataFrame(rows)
leads = lead_summary(pk)
print("Lead in s (positive: activity precedes behaviour). Expected for forward velocity: earliest 0.3 s (0.27 to 0.33), latest a lag over 1 s")
print(leads.round(3).to_string(index=False))
fw = leads.set_index("component").loc["fwd"]
for rd, e, l in (("a", "(a) earliest lead", "(a) latest lead"), ("b", "(b) 99th pct lead", "(b) 1st pct lead")):
print(f"Reading ({rd}): earliest within 0.27-0.33 s: {0.27 <= fw[e] <= 0.33}; latest lag > 1 s: {fw[l] < -1}")
hist = np.histogram(pk["fwd"]["lead"][pk["fwd"]["top"]], bins=np.arange(-5, 5.01, 0.25))
print("\nHistogram of VF leads (top 40%), 0.25 s bins, counts from -5 s:", hist[0].tolist())
# Variant: neural times at the start of each volume
pk_start = peaks(filters_for(-0.5 / rate))
print("\nVariant, neural times at the start of each volume:")
print(lead_summary(pk_start).round(3).to_string(index=False))
# Variant: top 20% and top 60%
for q in (80, 40):
L = pk["fwd"]["lead"][pk["fwd"]["peak"] >= np.percentile(pk["fwd"]["peak"], q)]
print(f"Variant, top {100 - q}% of VF filters: earliest {L.max():.2f}, 99th pct {np.percentile(L, 99):.2f}, "
f"1st pct {np.percentile(L, 1):.2f}, latest {L.min():.2f} s")
# Per plane, VF
print("\nVF lead by FDA plane (top 40%): median, 99th and 1st percentiles, range")
by_plane = pd.DataFrame({"plane": block_plane[pk["fwd"]["top"]], "lead": pk["fwd"]["lead"][pk["fwd"]["top"]]}).groupby("plane").lead
print(pd.concat([by_plane.agg(["size", "median"]), by_plane.quantile(0.99).rename("99th pct"), by_plane.quantile(0.01).rename("1st pct"),
by_plane.agg(["min", "max"])], axis=1).round(2).to_string())
L_top, P_top = pk["fwd"]["lead"][pk["fwd"]["top"]], block_plane[pk["fwd"]["top"]]
L_within = L_top - pd.Series(L_top).groupby(P_top).transform("median").values
print(f"Leads with each plane's median removed: 99th pct {np.percentile(L_within, 99):.2f} s, 1st pct {np.percentile(L_within, 1):.2f} s, "
f"spread (99th - 1st) {np.percentile(L_within, 99) - np.percentile(L_within, 1):.2f} s against {np.percentile(L_top, 99) - np.percentile(L_top, 1):.2f} s")
print(f"Volume period: {1 / rate:.3f} s")
rows = []
for j, f in enumerate(flies):
flt = bw.shifted_behaviour(by_fly[f], "fwd", times[f]) @ X[j * 3384:(j + 1) * 3384]
dec = bw.deconvolve(signal.filtfilt(b_notch, a_notch, flt, axis=0).T, toe)
t_ = dec.max(1) >= np.percentile(dec.max(1), 60)
L = bw.SHIFTS_S[dec.argmax(1)][t_]
rows.append({"fly": f, "99th pct lead": np.percentile(L, 99), "median lead": np.median(L), "1st pct lead": np.percentile(L, 1)})
print("\nPer fly, VF filters of that fly alone (top 40% within the fly):")
print(pd.DataFrame(rows).round(2).to_string(index=False))
fig, ax = plt.subplots(figsize=(6, 3))
ax.hist(pk["fwd"]["lead"][pk["fwd"]["top"]], bins=hist[1], color="0.4")
ax.set_xlabel("lead of activity over forward velocity (s)"); ax.set_ylabel("blocks")
plt.show()
Lead in s (positive: activity precedes behaviour). Expected for forward velocity: earliest 0.3 s (0.27 to 0.33), latest a lag over 1 s
component blocks (top 40%) (a) earliest lead (b) 99th pct lead median lead (b) 1st pct lead (a) latest lead at window edge median lead before deconvolution
fwd 2145 0.44 0.32 -0.22 -1.234 -1.44 0 -0.38
left 2145 0.74 0.30 -0.18 -1.240 -3.48 0 -0.46
right 2145 0.44 0.22 -0.16 -1.220 -1.94 0 -0.36
Reading (a): earliest within 0.27-0.33 s: False; latest lag > 1 s: True
Reading (b): earliest within 0.27-0.33 s: True; latest lag > 1 s: True
Histogram of VF leads (top 40%), 0.25 s bins, counts from -5 s: [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 22, 19, 225, 90, 598, 526, 509, 156, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
Variant, neural times at the start of each volume:
component blocks (top 40%) (a) earliest lead (b) 99th pct lead median lead (b) 1st pct lead (a) latest lead at window edge median lead before deconvolution
fwd 2145 0.7 0.58 0.04 -0.954 -1.16 0 -0.12
left 2145 1.0 0.56 0.08 -0.971 -3.20 0 -0.18
right 2145 0.7 0.48 0.10 -0.951 -1.68 0 -0.08
Variant, top 20% of VF filters: earliest 0.40, 99th pct 0.32, 1st pct -0.86, latest -1.20 s
Variant, top 60% of VF filters: earliest 0.44, 99th pct 0.32, 1st pct -1.38, latest -2.98 s
VF lead by FDA plane (top 40%): median, 99th and 1st percentiles, range
size median 99th pct 1st pct min max
plane
17 18 0.36 0.43 0.26 0.26 0.44
24 81 0.24 0.38 -0.15 -0.18 0.40
31 117 0.24 0.32 -0.24 -0.28 0.38
38 216 0.16 0.30 -0.84 -0.86 0.32
45 287 -0.12 0.28 -0.88 -1.18 0.40
52 345 -0.24 0.28 -0.88 -1.34 0.32
59 384 -0.26 0.26 -1.19 -1.38 0.32
66 377 -0.26 0.24 -1.36 -1.40 0.28
73 268 -0.28 0.22 -1.38 -1.44 0.26
80 52 -0.30 -0.10 -0.84 -0.88 -0.10
Leads with each plane's median removed: 99th pct 0.50 s, 1st pct -1.02 s, spread (99th - 1st) 1.52 s against 1.55 s
Volume period: 0.532 s
Per fly, VF filters of that fly alone (top 40% within the fly):
fly 99th pct lead median lead 1st pct lead
fly_087 0.36 -0.18 -1.45
fly_089 0.41 -0.68 -3.66
fly_094 0.42 -0.22 -1.35
fly_097 1.01 -0.50 -3.74
fly_098 0.70 -0.24 -2.05
fly_099 2.85 -0.02 -3.38
fly_100 0.36 -0.24 -2.41
fly_101 0.38 -0.22 -2.32
fly_105 4.55 -0.40 -4.12
Result: reproduced. Under reading (b), the 99th percentile of the VF lead is 0.32 s, within 0.27 to 0.33 s, and the 1st percentile is a lag of 1.23 s; under reading (a) the earliest lead is 0.44 s and the latest a lag of 1.44 s. For VL and VR the 99th percentiles are 0.30 and 0.22 s and the 1st percentiles lags of 1.24 and 1.22 s. The match within 0.03 s rests on the assumed neural times and is not robust to them: with the start of each volume as the neural time, the VF 99th percentile is 0.58 s. The median VF lead also changes steadily with FDA plane, from 0.36 s in plane 17 to -0.30 s in plane 80, a spread of 0.66 s against a volume period of 0.532 s; this is the pattern that the unknown acquisition time of each plane would produce, and the order in which planes were acquired is not in the dandiset. With each plane's median removed, the 99th percentile is 0.50 s and the 1st -1.02 s. In every plane from 38 to 80 the 1st percentile is a lag of at least 0.84 s, and in planes 59 to 73 a lag of more than 1 s. Fitted per fly, the 99th percentile is 0.36 to 4.55 s.
fig6h-peak-vs-width¶
Supervoxels whose peak correlation lags behaviour have broader temporal kernels, while those that anticipate behaviour have narrow kernels.
top = pk["fwd"]["top"]
width = bw.fwhm(pk["fwd"]["dec"][top])
lag = -pk["fwd"]["lead"][top]
r_w = stats.pearsonr(lag, width)
print(f"VF, top 40% ({top.sum()} blocks): lag {lag.min():.2f} to {lag.max():.2f} s, FWHM {width.min():.2f} to {width.max():.2f} s")
print(f"Pearson r(lag, FWHM) = {r_w.statistic:.3f}, p = {r_w.pvalue:.3g}; Spearman {stats.spearmanr(lag, width).statistic:.3f}")
inner = np.abs(lag) < 4.5
print(f"Without filters peaking within 0.5 s of the window edge ({(~inner).sum()} dropped): r = {stats.pearsonr(lag[inner], width[inner]).statistic:.3f}")
for key in ("left", "right"):
t_ = pk[key]["top"]
print(f"{key}: r(lag, FWHM) = {stats.pearsonr(-pk[key]['lead'][t_], bw.fwhm(pk[key]['dec'][t_])).statistic:.3f}")
tb = pk_start["fwd"]["top"]
print(f"Variant, volume start times: r = {stats.pearsonr(-pk_start['fwd']['lead'][tb], bw.fwhm(pk_start['fwd']['dec'][tb])).statistic:.3f}")
planes_top = block_plane[top]
rows = []
for z in Z:
m = planes_top == z
if m.sum() >= 10:
rr = stats.pearsonr(lag[m], width[m])
rows.append({"plane": z, "blocks": int(m.sum()), "r(lag, FWHM)": rr.statistic, "p": rr.pvalue})
within = pd.DataFrame(rows)
print("\nWithin each FDA plane (planes with at least 10 top blocks):")
print(within.to_string(index=False, float_format=lambda v: f"{v:.3g}"))
lag_res = lag - pd.Series(lag).groupby(planes_top).transform("mean").values
wid_res = width - pd.Series(width).groupby(planes_top).transform("mean").values
pr = stats.pearsonr(lag_res, wid_res)
print(f"Partial correlation controlling for plane (plane means removed): r = {pr.statistic:.3f}, p = {pr.pvalue:.3g}; "
f"plane-mean lag against plane-mean FWHM: r = {np.corrcoef(pd.Series(lag).groupby(planes_top).mean(), pd.Series(width).groupby(planes_top).mean())[0, 1]:.3f} over {len(set(planes_top))} planes")
print("Mean FWHM by lag bin (s):")
print(pd.Series(width).groupby(pd.cut(lag, [-5, -1, -0.5, 0, 0.5, 1, 5])).agg(["size", "mean"]).round(3).to_string())
VF, top 40% (2145 blocks): lag -0.44 to 1.44 s, FWHM 0.54 to 3.92 s Pearson r(lag, FWHM) = 0.672, p = 2.23e-281; Spearman 0.703 Without filters peaking within 0.5 s of the window edge (0 dropped): r = 0.672
left: r(lag, FWHM) = 0.442
right: r(lag, FWHM) = 0.449
Variant, volume start times: r = 0.673
Within each FDA plane (planes with at least 10 top blocks):
plane blocks r(lag, FWHM) p
17 18 0.259 0.3
24 81 0.259 0.0194
31 117 0.344 0.000148
38 216 0.734 9.55e-38
45 287 0.751 2.82e-53
52 345 0.681 2.27e-48
59 384 0.57 1.77e-34
66 377 0.607 2.13e-39
73 268 0.533 4.86e-21
80 52 0.222 0.114
Partial correlation controlling for plane (plane means removed): r = 0.623, p = 6.26e-231; plane-mean lag against plane-mean FWHM: r = 0.979 over 10 planes
Mean FWHM by lag bin (s):
size mean
(-0.5, 0.0] 665 1.384
(0.0, 0.5] 1124 1.892
(0.5, 1.0] 315 2.386
(1.0, 5.0] 41 2.617
Result: reproduced. Over the 2145 top VF blocks, lag and FWHM correlate at r = 0.672 (p = 2.2e-281; Spearman 0.703). Mean FWHM rises from 1.38 s for blocks leading behaviour by up to 0.5 s to 2.62 s for blocks lagging by more than 1 s. Plane accounts for part of it (the plane means of lag and FWHM correlate at 0.979 over 10 planes), but not for all: with plane means removed the correlation is 0.623 (p = 6.3e-231), and within single planes it is positive in all 10 planes with at least 10 blocks (0.222 to 0.751), with p < 0.05 in 8. VL and VR give r = 0.442 and 0.449, and volume-start times r = 0.673.
fig6fg-ipsi-leads-contra¶
Across nearly all brain regions, the peak correlation with angular velocity comes first in the hemisphere ipsilateral to the turn, then in the contralateral one.
res_ic = {}
for key, comp in (("left", "left"), ("right", "right")):
ipsi_is_a = ipsi_half[comp] == "A"
t_ = pk[key]["top"]
in_ipsi = (block_bx < nbx // 2) if ipsi_is_a else (block_bx >= nbx // 2)
pairs = [(c, partner[c]) for c in np.flatnonzero(in_ipsi & t_) if partner[c] >= 0 and t_[partner[c]]]
d = np.array([pk[key]["lead"][a] - pk[key]["lead"][b] for a, b in pairs])
w = stats.wilcoxon(d) if len(d) else None
res_ic[key] = {"ipsilateral half": ipsi_half[comp], "pairs": len(d), "median ipsi - contra lead (s)": float(np.median(d)),
"ipsi earlier (%)": 100 * float((d > 0).mean()), "ties (%)": 100 * float((d == 0).mean()),
"Wilcoxon p": w.pvalue}
print(pd.DataFrame(res_ic).T.to_string())
# By plane
for key in ("left", "right"):
ipsi_is_a = ipsi_half[key] == "A"
t_ = pk[key]["top"]
in_ipsi = (block_bx < nbx // 2) if ipsi_is_a else (block_bx >= nbx // 2)
rows = [(block_plane[c], pk[key]["lead"][c] - pk[key]["lead"][partner[c]])
for c in np.flatnonzero(in_ipsi & t_) if partner[c] >= 0 and t_[partner[c]]]
by = pd.DataFrame(rows, columns=["plane", "d"]).groupby("plane").d.agg(["size", "median"])
print(f"\n{key}: median ipsi - contra lead by plane")
print(by.round(3).T.to_string())
ipsilateral half pairs median ipsi - contra lead (s) ipsi earlier (%) ties (%) Wilcoxon p left A 868 0.16 67.281106 5.184332 0.0 right B 849 0.14 72.438163 6.007067 0.0 left: median ipsi - contra lead by plane plane 17 24 31 38 45 52 59 66 73 80 size 6.00 32.00 56.00 88.00 114.00 137.00 159.00 158.00 100.00 18.00 median 0.33 0.32 0.32 0.26 0.08 0.18 0.18 0.08 0.02 -0.01 right: median ipsi - contra lead by plane plane 17 24 31 38 45 52 59 66 73 80 size 4.00 28.00 51.00 92.0 118.00 129.00 152.00 155.00 102.00 18.00 median 0.29 0.26 0.26 0.2 0.22 0.22 0.19 0.02 0.03 -0.03
Result: reproduced. For VL, over 868 mirror pairs, the block in half A (the half with more VL correlation) leads its partner by a median 0.16 s and is earlier in 67% of pairs; for VR, over 849 pairs, the block in half B leads by 0.14 s and is earlier in 72%; both Wilcoxon p-values print as 0. The median difference is 0.26 to 0.33 s in planes 17 to 38 for VL and 0.19 to 0.29 s in planes 17 to 59 for VR, and near 0 in planes 73 and 80 (-0.03 to 0.03 s). Both blocks of a pair lie in the same FDA plane, so a timing offset between planes does not enter the difference; whether the two halves of a plane were acquired at the same time is not known.
figs2-tdtomato-control (blocked): a check on the raw channels¶
Correlations computed from the tdTomato structural channel show no behaviour-related signal, so motion artefacts are undetectably small.
The claim needs the motion-corrected tdTomato channel, which is not in the dandiset, so it stays blocked. The raw channels are there, without motion correction. This check, which does not decide the outcome, correlates raw tdTomato and raw GCaMP6f voxels with VF, VL and VR within each fly, on raw planes 18 and 30 of 49, after the code's high-pass filter. Voxels are those with mean tdTomato above the plane's median (a choice). Uncorrected tdTomato carries all the motion; little correlation there would bound the motion artefacts from above, while much correlation would not tell whether the paper's correction removed them.
thr_fly = bw.r_threshold(1e-4, 3384)
rows = []
for f in flies:
for z in bw.RAW_SLICES:
rz = raw_r[f][z]
row = {"fly": f, "raw plane": z, "voxels": rz["voxels"]}
for ch in ("green", "red"):
sig = np.abs(rz[ch]) > thr_fly
row[f"{ch}: correlated %"] = 100 * sig.any(0).mean()
row[f"{ch}: VF %"] = 100 * sig[0].mean()
row[f"{ch}: max |r|"] = np.abs(rz[ch]).max()
rows.append(row)
raw_tab = pd.DataFrame(rows)
print(f"Within-fly correlations, |r| threshold for p < 1e-4 with 3384 points: {thr_fly:.4f}")
print(raw_tab.round(2).to_string(index=False))
print("\nMedian over flies and planes:", raw_tab.drop(columns=["fly", "raw plane"]).median().round(2).to_dict())
print("tdTomato above GCaMP6f in correlated %:", int((raw_tab["red: correlated %"] > raw_tab["green: correlated %"]).sum()), "of", len(raw_tab))
Within-fly correlations, |r| threshold for p < 1e-4 with 3384 points: 0.0668
fly raw plane voxels green: correlated % green: VF % green: max |r| red: correlated % red: VF % red: max |r|
fly_087 18 16384 2.27 1.31 0.12 1.98 0.07 0.17
fly_087 30 16384 1.93 1.26 0.14 0.16 0.02 0.09
fly_089 18 16384 0.83 0.54 0.13 3.53 2.92 0.17
fly_089 30 16384 0.41 0.18 0.11 1.39 0.92 0.14
fly_094 18 16384 3.73 3.01 0.14 0.04 0.01 0.08
fly_094 30 16384 3.16 2.47 0.14 0.05 0.01 0.08
fly_097 18 16384 0.50 0.07 0.09 0.10 0.02 0.09
fly_097 30 16384 0.46 0.06 0.10 0.14 0.01 0.09
fly_098 18 16384 0.41 0.16 0.09 0.46 0.12 0.10
fly_098 30 16384 0.37 0.07 0.10 0.24 0.06 0.09
fly_099 18 16384 0.95 0.24 0.13 0.71 0.02 0.14
fly_099 30 16384 0.70 0.13 0.10 0.13 0.03 0.10
fly_100 18 16384 0.20 0.12 0.09 0.22 0.04 0.09
fly_100 30 16384 0.25 0.09 0.09 0.09 0.01 0.09
fly_101 18 16384 0.55 0.20 0.14 1.31 0.07 0.13
fly_101 30 16384 0.43 0.31 0.11 0.34 0.15 0.10
fly_105 18 16384 0.66 0.29 0.10 0.25 0.10 0.09
fly_105 30 16384 0.89 0.35 0.11 0.21 0.06 0.10
Median over flies and planes: {'voxels': 16384.0, 'green: correlated %': 0.6, 'green: VF %': 0.22, 'green: max |r|': 0.11, 'red: correlated %': 0.23, 'red: VF %': 0.05, 'red: max |r|': 0.09}
tdTomato above GCaMP6f in correlated %: 5 of 18
Check, not deciding the outcome. In the raw, uncorrected channels, the median share of voxels correlated with VF, VL or VR at p < 1e-4 within a fly is 0.23% for tdTomato and 0.60% for GCaMP6f over 9 flies and 2 planes, and tdTomato exceeds GCaMP6f in 5 of 18 fly-planes. The largest tdTomato shares are in fly_089 (3.53% and 1.39%, against 0.83% and 0.41% for GCaMP6f) and fly_101 plane 18 (1.31% against 0.55%). Even before motion correction, behaviour-related signal in tdTomato is small in most flies but not in all. Whether the paper's correction removed it in fly_089 cannot be tested without the corrected channel, so the claim stays blocked.
print(f"Notebook ran in {(time.time() - started) / 60:.1f} min")
Notebook ran in 13.0 min
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 | Outcome |
|---|---|
| dataset-structure | reproduced |
| fig2c-fraction-correlated | partially reproduced |
| fig2c-venn-categories | partially reproduced |
| fig2-symmetry | reproduced |
| fig2e-region-fractions | blocked: metadata insufficient |
| fig3-topographic-maps | blocked: metadata insufficient |
| fig4bcd-pc-structure | reproduced |
| fig4fg-ridge-r2 | partially reproduced |
| fig4h-r2-vs-movement | partially reproduced |
| fig4i-weights-stereotyped | reproduced |
| fig5a-unique-variance | partially reproduced |
| fig5c-global-below-threshold | reproduced |
| fig5f-ipsi-contra | not reproduced |
| fig6-peak-time-range | reproduced |
| fig6h-peak-vs-width | reproduced |
| fig6fg-ipsi-leads-contra | reproduced |
| fig6eg-region-sequence | blocked: metadata insufficient |
| fig7-within-region-gradients | blocked: metadata insufficient |
| figs2-tdtomato-control | blocked: data missing |
| figs4a-supervoxel-number | blocked: too large |
How the analysis differs from the paper's. It uses the authors' processed volumes in the FDA space, not the paper's own mean brain, and 12 of the 91 FDA planes. The paper's supervoxels (Ward clusters, 2000 per plane) are replaced by 6 x 6 voxel blocks. Behaviour is sampled at the midpoint of each volume, as the dandiset has no plane times; the paper sampled it at each plane's own time. Folds, the ridge penalty of the demixing models, the threshold for a unique contribution, the bounds for slow walking and turning, and the notch filter's width are choices made here, as neither the Methods nor the code fix them. The 1000 and 2000 PC models were not run. The Figure 5F test compares mean activity in turning and slow walking, where the paper uses the demixing model; the model's prediction is shown as a variant and gives the same result. The Figure 6 timing results depend on the assumed neural times and change with FDA plane; the lag-width relation holds within planes. The supervoxel-number sweep of Figure S4A was not attempted on the 12 planes: it needs every voxel's time series held for Ward clustering at seven sizes, while the streaming here keeps only block means and per-voxel sums.
What was found about the dandiset.
- It holds the raw and processed imaging and the ball tracking of all nine flies of the paper, and no recording twice. The processed volumes are in the Functional Drosophila Atlas, a space the preprint does not use, with voxel size given as 1, 1, 1 in an unknown unit.
- There are no neuropil labels for any of its spaces, and no mean brain, supervoxels or z-depth maps, so none of the paper's region-level results (Figures 2E to I, 3, 6E to G, 7) can be tested.
- Per-plane acquisition times are missing. The share of correlated voxels changes from 31% to 39% between the middle and the start of each volume, and the 99th percentile of forward-velocity lead from 0.32 to 0.58 s.
- Which half of the volume is the fly's left is not stated.
- The motion-corrected tdTomato channel, needed for the paper's motion control, is not included. In the raw, uncorrected tdTomato channel few voxels correlate with behaviour in most flies (median 0.23% of voxels on two raw planes), but up to 3.5% in fly_089.
- The nine behaviour videos (AVI) are deposited; the analysis does not use them.
- Every FicTrac series has unit "meters" and conversion 0.0045, though most hold angles in
radians;
SpatialSeriesRotationDeltaError, described in radians, has medians of 6067 to 10,435. - The paper's code computes the walking state with
np.sqrt(a, b), which uses only forward velocity; the line was written to combine forward and rotational velocity.