Reproducing Lohani, Moberly et al. (2022) from dandiset 001172¶
Paper. Lohani S, Moberly AH, Benisty H, Landa B, Jing M, Li Y, Higley MJ, Cardin JA (2022). Spatiotemporally heterogeneous coordination of cholinergic and neocortical activity. Nature Neuroscience 25, 1706-1713. doi:10.1038/s41593-022-01202-6. This attempt used the full text of the accepted author manuscript in PubMed Central (PMC10661869); the version of record is not open access.
Dandiset. DANDI:001172, version 0.260129.0829. The paper cites no DANDI version (its data statement offers the data on request; the dandiset was published later). This is the only published version. The draft holds the same 32 assets (same blobs and sizes) and differs only in metadata (anatomy terms, author roles), none of which the analysis uses.
What the dandiset holds. Mesoscopic imaging of ACh3.0 (blue excitation) and jRCaMP1b (green excitation) with a 395 nm control channel, raw and as atlas-registered ΔF/F, face video with FaceMap output, pupil area, wheel speed and synchronisation pulses, for 4 of the paper's 6 mice (grabAM05 to grabAM08): 4 spontaneous, 6 air-puff and 6 visual stimulus sessions (data check). It holds no ECoG, no drug sessions, no ChAT-GCaMP6s mice and no control (GFP, mCherry) mice.
Outside inputs. The Allen parcel masks (parcells_updated121519.mat) and the brain mask
(brainMask.mat) of the paper's code (GitHub cardin-higley-lab/Lohani_Moberly_et_al_2022, commit
c337686), which name the pixels of the dandiset's ΔF/F. They bring no measurements; the claims that
use them are with_labels and depend on allen-parcel-link. The procedures follow that code where
the Methods leave a choice open (see meso.py). The same repository's
Figure2Regression/Fig2RawData.mat, the paper's per-mouse R2 values, is used only for a comparison
in the Figure 2 section; it decides no outcome.
Timing. The files' ΔF/F series start 49.9 s too early and their FaceMap and pupil timestamps are wrong (data check). The analysis times ΔF/F frame i at blue LED pulse i + 499 and FaceMap frame i at camera pulse i, as the paper's code does.
| Claim | Tested | |
|---|---|---|
| allen-parcel-link | yes | link of the code's parcels to the data |
| fig2a-r2-differs-across-areas | yes | |
| fig2b-ap-gradient-inverted | yes | |
| fig2c-face-best-predictor-ach | yes | |
| fig2c-face-best-predictor-calcium | yes | |
| s7-face-onset-transients | yes | |
| s7-locomotion-onset-transients | yes | |
| s7-ach-transients-same-across-states | yes | |
| s7-calcium-transients-largest-at-locomotion | yes | |
| s7-ap-gradient-transients | yes | |
| fig3b-ecog-desynchronization | blocked | no ECoG in the dandiset |
| fig3c-ach-sustained-face | yes | 3 mice |
| fig3c-ach-locomotion-no-further-increase | yes | 3 mice |
| fig3c-calcium-sustained | yes | 3 mice |
| s7-ap-gradient-sustained | yes | 3 mice |
| fig3d-ach-frontal-above-visual | yes | 3 mice |
| fig3e-scopolamine-motor-not-visual | blocked | no drug sessions |
| fig4a-ach-correlations-up-with-face | yes | 3 mice |
| fig4a-ach-correlations-down-with-locomotion | yes | 3 mice |
| fig4b-calcium-correlations-up-with-face | yes | 3 mice |
| fig4b-calcium-correlations-down-posterior | yes | 3 mice |
| fig5-coupling-up-with-face | yes | 3 mice |
| fig5-coupling-down-frontal-with-locomotion | yes | 3 mice |
| s3-hemodynamic-correction | blocked | no control mice |
| s2-ach-sensor-validation | blocked | no slice or stimulation data |
| s8-cholinergic-axons | blocked | no ChAT-GCaMP6s mice |
The notebook takes about 17 minutes, most of it reading the ΔF/F movies.
import itertools
import time
import warnings
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from scipy import stats
from scipy.ndimage import uniform_filter1d
from statsmodels.stats.anova import AnovaRM
from statsmodels.stats.multitest import multipletests
import meso
pd.set_option("display.width", 220)
pd.set_option("display.max_columns", 60)
pd.set_option("display.max_rows", 300)
warnings.filterwarnings("ignore", category=RuntimeWarning)
started = time.time()
parcels = meso.load_parcels()
W, frac_in_mask = meso.parcel_weights(parcels)
sel = pd.DataFrame(parcels["sel"])
LEFT = sel.index[sel.hemi == "L"].to_numpy()
left_areas = sel.area[LEFT].to_numpy()
ap_left = sel.ap_rank[LEFT].to_numpy()
sessions = []
for mouse, names in meso.SESSIONS.items():
for name in names:
if (mouse, name) == meso.SESSION_WITHOUT_PUFFS:
continue
sessions.append(meso.load_session(mouse, name, W))
print(f"loaded {len(sessions)} sessions from dandiset {meso.DANDISET_ID} version {meso.DANDISET_VERSION} "
f"in {time.time() - started:.0f} s")
MICE = list(meso.SESSIONS)
loaded 9 sessions from dandiset 001172 version 0.260129.0829 in 607 s
Data check¶
Each check of the metadata against the values is stated under its output. Values of ΔF/F are shown in percent, as in the paper.
rows = []
for s in sessions:
rows.append({"mouse": s["mouse"], "session": s["session"], "raw_frames": s["n_raw"],
"dff_frames": s["ach"].shape[0], "blue_pulses": s["ttl_counts"]["TTLSignalBlueLED"],
"green_pulses": s["ttl_counts"]["TTLSignalGreenLED"],
"camera_pulses": s["ttl_counts"]["TTLSignalPupilCamera"], "face_frames": len(s["face_norm"]),
"minutes": round((s["t"][-1] - s["t"][0]) / 60, 1),
"frame_interval_ms": round(1000 * np.median(np.diff(s["t"])), 2),
"airpuffs": None if s["airpuff"] is None else len(s["airpuff"]),
"wheel_activation": len(s["wheel_activation"]),
"parcellated": str(s["parcellated_shape"])})
info = pd.DataFrame(rows)
print(info.to_string(index=False))
print("\nMice:", info.mouse.nunique(), "(the paper and its code use 6, grabAM05 to grabAM10); sessions used:", len(info))
print("raw frames - ΔF/F frames per session:", sorted(set(info.raw_frames - info.dff_frames)))
url, f = meso.open_nwb(meso.asset_path(*meso.SESSION_WITHOUT_PUFFS))
print("Left out,", meso.SESSION_WITHOUT_PUFFS, "intervals in file:", list(f["intervals"].keys()))
f.close()
mouse session raw_frames dff_frames blue_pulses green_pulses camera_pulses face_frames minutes frame_interval_ms airpuffs wheel_activation parcellated grabAM05 11222019-spont 18098 17099 18098 18097 18105 18078 28.5 100.0 NaN 8 (16799, 93) grabAM05 11232019-airpuffs 26203 25204 26203 26202 26205 26178 42.0 100.0 7.0 6 (24904, 96) grabAM05 11262019-airpuffs 25264 24265 25264 25264 25290 25262 40.4 100.0 7.0 3 (23965, 90) grabAM06 11222019-spont 15930 14931 15930 15929 15965 15938 24.9 100.0 NaN 30 (14631, 107) grabAM06 11242019-airpuffs 24676 23677 24676 24676 24677 24650 39.5 100.0 7.0 13 (23377, 90) grabAM07 11232019-spont 15341 14342 15341 15340 15328 15301 23.9 100.0 NaN 14 (14042, 98) grabAM07 11242019-airpuffs 24685 23686 24685 24685 24685 24657 39.5 100.0 2.0 26 (23386, 89) grabAM08 11232019-spont 15958 14959 15958 15957 15972 15944 24.9 100.0 NaN 18 (14659, 100) grabAM08 11242019-airpuffs 22250 21251 22250 22249 22270 22243 35.4 100.0 6.0 10 (20951, 97) Mice: 4 (the paper and its code use 6, grabAM05 to grabAM10); sessions used: 9 raw frames - ΔF/F frames per session: [999]
Left out, ('grabAM07', '11262019-airpuffs') intervals in file: ['WheelActivation']
# Sessions in the dandiset by mouse and type, from the archive's asset list.
import re
from dandi.dandiapi import DandiAPIClient
with DandiAPIClient() as client:
paths = [a.path for a in client.get_dandiset(meso.DANDISET_ID, meso.DANDISET_VERSION).get_assets()
if a.path.endswith(".nwb")]
ses = pd.DataFrame([re.match(r"sub-(?P<mouse>\w+)/.*_ses-\d+-(?P<type>[a-z-]+)_", p).groupdict() for p in paths])
print(f"{len(paths)} NWB files")
print(ses.groupby(["mouse", "type"]).size().unstack(fill_value=0).to_string())
print("total per type:", ses.type.value_counts().to_dict())
16 NWB files
type airpuffs spont vis-stim
mouse
grabAM05 2 1 2
grabAM06 1 1 2
grabAM07 2 1 1
grabAM08 1 1 1
total per type: {'airpuffs': 6, 'vis-stim': 6, 'spont': 4}
The dandiset holds 16 NWB files from 4 of the 6 mice of the paper's code (grabAM05 to grabAM08): per mouse one spontaneous session and one or two air-puff sessions (the code's state analyses use these), and one or two visual stimulus sessions, 6 in all, of which one per mouse is opened here. Each ΔF/F movie has 999 frames fewer than its raw movie, as the code's detrending removes them. The frame interval is 100 ms (10 Hz per colour, Methods). The air-puff session grabAM07 11262019 has no Airpuff table, so it is left out.
# The ΔF/F pixels against the paper code's brain mask, under the 8 rotations and reflections.
def orientations(img):
for k in range(4):
for flip in (False, True):
m = np.rot90(img, k)
yield (k * 90, flip), (m.T if flip else m)
agree = []
for s in sessions:
for colour in ("ach", "ca"):
nz = s[colour + "_nonzero"]
agree.append({"session": f"{s['mouse']} {s['session']}", "colour": colour, "nonzero_pixels": int(nz.sum()),
**{f"rot{k}{'_T' if fl else ''}": round((nz == m).mean(), 4) for (k, fl), m in orientations(parcels["mask"])}})
agree = pd.DataFrame(agree)
print("Reading 1, each movie separately (pixels non-zero in any frame of the movie):")
print(agree.to_string(index=False))
union = []
for s in sessions:
nz = s["ach_nonzero"] | s["ca_nonzero"]
union.append({"session": f"{s['mouse']} {s['session']}", "nonzero_pixels": int(nz.sum()),
**{f"rot{k}{'_T' if fl else ''}": round((nz == m).mean(), 4) for (k, fl), m in orientations(parcels["mask"])}})
union = pd.DataFrame(union)
print("\nReading 2, each session (pixels non-zero in the blue or the green movie):")
print(union.to_string(index=False))
odd = [(s["mouse"], s["session"], c) for s in sessions for c in ("ach", "ca") if s[c + "_nonzero"].sum() != parcels["mask"].sum()]
for mouse, name, c in odd:
s = next(x for x in sessions if x["mouse"] == mouse and x["session"] == name)
nz = s[c + "_nonzero"]
m = meso.orient(parcels["mask"])
print(f"{mouse} {name} {c}: zero in every frame inside the mask {int((m & ~nz).sum())} pixels; non-zero outside the mask {int((~m & nz).sum())}")
hit = {p["name"]: int((meso.orient(parcels["indicators"][:, :, p["index"]]) & m & ~nz).sum()) for p in parcels["sel"]}
in_any = np.any([meso.orient(parcels["indicators"][:, :, p["index"]]) for p in parcels["sel"]], axis=0)
print(" of these, in each selected parcel (parcels overlap at borders):", {k: v for k, v in hit.items() if v},
"; in no selected parcel:", int((m & ~nz & ~in_any).sum()))
print("code brainMask pixels:", int(parcels["mask"].sum()))
print("parcels with fewer than 100% of pixels in the mask:",
{sel.name[i]: round(frac_in_mask[i], 2) for i in range(len(sel)) if frac_in_mask[i] < 1})
Reading 1, each movie separately (pixels non-zero in any frame of the movie):
session colour nonzero_pixels rot0 rot0_T rot90 rot90_T rot180 rot180_T rot270 rot270_T
grabAM05 11222019-spont ach 25414 0.8317 0.7681 0.7708 0.8407 1.000 0.7739 0.7708 0.9818
grabAM05 11222019-spont ca 25414 0.8317 0.7681 0.7708 0.8407 1.000 0.7739 0.7708 0.9818
grabAM05 11232019-airpuffs ach 25414 0.8317 0.7681 0.7708 0.8407 1.000 0.7739 0.7708 0.9818
grabAM05 11232019-airpuffs ca 25414 0.8317 0.7681 0.7708 0.8407 1.000 0.7739 0.7708 0.9818
grabAM05 11262019-airpuffs ach 25414 0.8317 0.7681 0.7708 0.8407 1.000 0.7739 0.7708 0.9818
grabAM05 11262019-airpuffs ca 25414 0.8317 0.7681 0.7708 0.8407 1.000 0.7739 0.7708 0.9818
grabAM06 11222019-spont ach 25414 0.8317 0.7681 0.7708 0.8407 1.000 0.7739 0.7708 0.9818
grabAM06 11222019-spont ca 25414 0.8317 0.7681 0.7708 0.8407 1.000 0.7739 0.7708 0.9818
grabAM06 11242019-airpuffs ach 25414 0.8317 0.7681 0.7708 0.8407 1.000 0.7739 0.7708 0.9818
grabAM06 11242019-airpuffs ca 25414 0.8317 0.7681 0.7708 0.8407 1.000 0.7739 0.7708 0.9818
grabAM07 11232019-spont ach 25414 0.8317 0.7681 0.7708 0.8407 1.000 0.7739 0.7708 0.9818
grabAM07 11232019-spont ca 25414 0.8317 0.7681 0.7708 0.8407 1.000 0.7739 0.7708 0.9818
grabAM07 11242019-airpuffs ach 25414 0.8317 0.7681 0.7708 0.8407 1.000 0.7739 0.7708 0.9818
grabAM07 11242019-airpuffs ca 25414 0.8317 0.7681 0.7708 0.8407 1.000 0.7739 0.7708 0.9818
grabAM08 11232019-spont ach 25414 0.8317 0.7681 0.7708 0.8407 1.000 0.7739 0.7708 0.9818
grabAM08 11232019-spont ca 25151 0.8357 0.7721 0.7748 0.8447 0.996 0.7779 0.7748 0.9792
grabAM08 11242019-airpuffs ach 25414 0.8317 0.7681 0.7708 0.8407 1.000 0.7739 0.7708 0.9818
grabAM08 11242019-airpuffs ca 25414 0.8317 0.7681 0.7708 0.8407 1.000 0.7739 0.7708 0.9818
Reading 2, each session (pixels non-zero in the blue or the green movie):
session nonzero_pixels rot0 rot0_T rot90 rot90_T rot180 rot180_T rot270 rot270_T
grabAM05 11222019-spont 25414 0.8317 0.7681 0.7708 0.8407 1.0 0.7739 0.7708 0.9818
grabAM05 11232019-airpuffs 25414 0.8317 0.7681 0.7708 0.8407 1.0 0.7739 0.7708 0.9818
grabAM05 11262019-airpuffs 25414 0.8317 0.7681 0.7708 0.8407 1.0 0.7739 0.7708 0.9818
grabAM06 11222019-spont 25414 0.8317 0.7681 0.7708 0.8407 1.0 0.7739 0.7708 0.9818
grabAM06 11242019-airpuffs 25414 0.8317 0.7681 0.7708 0.8407 1.0 0.7739 0.7708 0.9818
grabAM07 11232019-spont 25414 0.8317 0.7681 0.7708 0.8407 1.0 0.7739 0.7708 0.9818
grabAM07 11242019-airpuffs 25414 0.8317 0.7681 0.7708 0.8407 1.0 0.7739 0.7708 0.9818
grabAM08 11232019-spont 25414 0.8317 0.7681 0.7708 0.8407 1.0 0.7739 0.7708 0.9818
grabAM08 11242019-airpuffs 25414 0.8317 0.7681 0.7708 0.8407 1.0 0.7739 0.7708 0.9818
grabAM08 11232019-spont ca: zero in every frame inside the mask 263 pixels; non-zero outside the mask 0
of these, in each selected parcel (parcels overlap at borders): {'R-TEa': 8, 'R-AUD': 263, 'R-SSs': 9} ; in no selected parcel: 0
code brainMask pixels: 25414
parcels with fewer than 100% of pixels in the mask: {'R-RSPagl': np.float64(0.99), 'L-RSPagl': np.float64(1.0), 'R-RSPd': np.float64(0.55), 'L-RSPd': np.float64(0.64), 'R-MOs': np.float64(0.71), 'L-MOs': np.float64(0.73)}
In every session the pixels that are non-zero in the ΔF/F are the 25414 pixels of the code's
brainMask.mat rotated by 180 degrees, and no other orientation comes close (at most 0.98). In
one movie (grabAM08 11232019 spont, calcium) 263 pixels inside the mask are zero in every frame,
all of them in the right AUD parcel (8 and 9 of them also in the overlapping right TEa and SSs
borders); no pixel outside the mask is non-zero. The left-hemisphere summaries do not use these
parcels. So the dandiset's ΔF/F is the
code's atlas-registered final_dFoF.mat, and the code's parcel masks apply after the same
rotation. The NWB description of these series mentions only detrending and ΔF/F; the hemodynamic
correction and atlas registration that the code applies before saving them are not stated (the
next checks bear on the correction).
# Frame offset between the raw blue movie and the blue ΔF/F, from the frame-to-frame changes of
# the mean over the frame (raw) and over the brain mask (ΔF/F), for every session used here.
def raw_means(path, n=1100, series="OnePhotonSeriesBlueExcitationGreenChannel"):
url, f = meso.open_nwb(path)
d = f["acquisition/" + series + "/data"]
parts = meso.reduce_many([(url, d, lambda b, k0: b.reshape(b.shape[0], -1).mean(1), range(n))])[0]
f.close()
return np.concatenate(parts)[:n]
def best_offset(raw, dff, n_dff=400):
r, d = np.diff(raw), np.diff(dff[:n_dff])
cs = np.array([np.corrcoef(r[o:o + n_dff - 1], d)[0, 1] for o in range(len(r) - n_dff + 2)])
return int(cs.argmax()), cs.max(), cs[meso.DETREND_DROP]
offsets = []
for s in sessions:
o, c, c2 = best_offset(raw_means(meso.asset_path(s["mouse"], s["session"])), s["ach"][:, -1])
offsets.append({"session": f"{s['mouse']} {s['session']}", "offset_frames": o, "r": round(c, 3), "r_at_499": round(c2, 3),
"file_start_minus_pulse0_s": round(s["ach_file_start"] - s["blue_ttl"][0], 3),
"pulse_499_minus_pulse0_s": round(s["blue_ttl"][meso.DETREND_DROP] - s["blue_ttl"][0], 2)})
offsets = pd.DataFrame(offsets)
print(offsets.to_string(index=False))
session offset_frames r r_at_499 file_start_minus_pulse0_s pulse_499_minus_pulse0_s grabAM05 11222019-spont 499 0.944 0.944 0.0 49.9 grabAM05 11232019-airpuffs 499 0.941 0.941 0.0 49.9 grabAM05 11262019-airpuffs 499 0.923 0.923 0.0 49.9 grabAM06 11222019-spont 499 0.938 0.938 0.0 49.9 grabAM06 11242019-airpuffs 499 0.916 0.916 0.0 49.9 grabAM07 11232019-spont 499 0.912 0.912 0.0 49.9 grabAM07 11242019-airpuffs 499 0.920 0.920 0.0 49.9 grabAM08 11232019-spont 499 0.805 0.805 0.0 49.9 grabAM08 11242019-airpuffs 499 0.846 0.846 0.0 49.9
In every session ΔF/F frame 0 is raw frame 499 (the code's detrending drops the first 500 samples). The files give the ΔF/F series the starting time of the first blue LED pulse, the time of raw frame 0, so their ΔF/F times are about 49.9 s early. The analysis times ΔF/F frame i at blue pulse i + 499. The next check shows the same offset from behaviour.
# Is the archived blue ΔF/F corrected for hemodynamics? In one session, the blue (470 nm) and violet
# (395 nm, the hemodynamic reference) signals are compared before correction (raw movies, native
# space, 8 x 8 pixel blocks) and in the archived ΔF/F (atlas space, pixels in the brain mask), over
# the first 2000 ΔF/F frames (raw frames 499 to 2498), after removing a 30 s moving mean.
hs = sessions[0]
path = meso.asset_path(hs["mouse"], hs["session"])
n_h = 2000
url, f = meso.open_nwb(path)
def blocks(b, k0):
n, r, c = b.shape
return b.reshape(n, r // 8, 8, c // 8, 8).astype(float).mean((2, 4)).reshape(n, -1)
mpix = meso.orient(parcels["mask"]).ravel()
raw_jobs = [(url, f["acquisition/" + nm + "/data"], blocks, range(meso.DETREND_DROP + n_h))
for nm in ("OnePhotonSeriesBlueExcitationGreenChannel", "OnePhotonSeriesVioletExcitationGreenChannel")]
dff_jobs = [(url, f["processing/ophys/" + nm + "/data"], lambda b, k0: b.reshape(b.shape[0], -1)[:, mpix], range(n_h))
for nm in ("DFFOnePhotonSeriesBlue", "DFFOnePhotonSeriesViolet")]
print("raw shapes:", [j[1].shape for j in raw_jobs], "; ΔF/F shapes:", [j[1].shape for j in dff_jobs])
res = meso.reduce_many(raw_jobs + dff_jobs)
f.close()
rb, rv = (np.concatenate(r)[meso.DETREND_DROP:meso.DETREND_DROP + n_h] for r in res[:2])
db_, dv_ = (np.concatenate(r)[:n_h] for r in res[2:])
def colcorr(a, b):
a = a - uniform_filter1d(a, 300, axis=0, mode="nearest")
b = b - uniform_filter1d(b, 300, axis=0, mode="nearest")
a, b = a - a.mean(0), b - b.mean(0)
with np.errstate(invalid="ignore", divide="ignore"):
return (a * b).sum(0) / np.sqrt((a * a).sum(0) * (b * b).sum(0))
bright = rb.mean(0) > np.median(rb.mean(0))
r_raw = colcorr(rb[:, bright], rv[:, bright])
r_dff = colcorr(db_, dv_)
print(f"{hs['mouse']} {hs['session']}, 30 s moving mean removed, blue against violet:")
print(f" raw, mean over the frame: r = {colcorr(rb.mean(1, keepdims=True), rv.mean(1, keepdims=True))[0]:.3f}; "
f"per 8 x 8 block (the brighter half, {bright.sum()} blocks): median r = {np.nanmedian(r_raw):.3f} (IQR {np.nanpercentile(r_raw, 25):.3f} to {np.nanpercentile(r_raw, 75):.3f})")
print(f" archived ΔF/F, mean over the brain mask: r = {colcorr(db_.mean(1, keepdims=True), dv_.mean(1, keepdims=True))[0]:.3f}; "
f"per pixel ({np.isfinite(r_dff).sum()} pixels): median r = {np.nanmedian(r_dff):.3f} (IQR {np.nanpercentile(r_dff, 25):.3f} to {np.nanpercentile(r_dff, 75):.3f})")
mean = lambda x: x.mean(1, keepdims=True)
print(f" archived blue ΔF/F mean against raw blue mean r = {colcorr(mean(db_), mean(rb))[0]:.3f}, against raw violet mean "
f"r = {colcorr(mean(db_), mean(rv))[0]:.3f}; archived violet ΔF/F mean against raw violet mean r = {colcorr(mean(dv_), mean(rv))[0]:.3f}")
print(f" violet ΔF/F non-zero pixels: {int((dv_ != 0).any(0).sum())} of {mpix.sum()} in the mask")
del db_, dv_
raw shapes: [(18098, 512, 512), (18098, 512, 512)] ; ΔF/F shapes: [(17099, 256, 256), (17099, 256, 256)]
grabAM05 11222019-spont, 30 s moving mean removed, blue against violet: raw, mean over the frame: r = 0.576; per 8 x 8 block (the brighter half, 2048 blocks): median r = 0.298 (IQR 0.167 to 0.434) archived ΔF/F, mean over the brain mask: r = 0.070; per pixel (25414 pixels): median r = 0.007 (IQR -0.013 to 0.028) archived blue ΔF/F mean against raw blue mean r = 0.825, against raw violet mean r = 0.032; archived violet ΔF/F mean against raw violet mean r = 0.969 violet ΔF/F non-zero pixels: 25414 of 25414 in the mask
Before correction, with slow drift removed, the raw blue and violet signals share much of their fluctuation (r = 0.576 for the frame means, median 0.298 over 8 x 8 blocks). In the archived ΔF/F they share almost none (r = 0.070 for the brain means, median 0.007 per pixel), while the archived blue still follows the raw blue (r = 0.825) and not the raw violet (r = 0.032), and the archived violet follows the raw violet (r = 0.969). This is consistent with the archived blue ΔF/F having the violet-shared (hemodynamic) component regressed out, as the code does; it does not show which regression was used, since the registration transforms needed to repeat it in native space are not in the files. One session, the first 2000 frames.
def best_lag(a, b, maxlag):
a = (a - np.nanmean(a)) / np.nanstd(a)
b = (b - np.nanmean(b)) / np.nanstd(b)
n, out = len(a), []
for L in range(-maxlag, maxlag + 1):
x, y = a[max(0, L): n + min(0, L)], b[max(0, -L): n - max(0, L)]
out.append(np.mean(x * y))
out = np.array(out)
return int(out.argmax()) - maxlag, out.max(), out[maxlag]
clock = []
for s in sessions:
brain_ca = s["ca"][:, -1]
row = {"session": f"{s['mouse']} {s['session']}"}
for label, t in (("file", s["ca_file_start"] + np.arange(len(brain_ca)) / 10), ("pulse", s["t"])):
sp = np.interp(t, s["speed_100hz_t"], np.abs(s["speed_100hz"]))
L, r, r0 = best_lag(brain_ca, sp, 800)
row.update({f"ca_wheel_lag_{label}_frames": L, f"r_{label}": round(r, 3), f"r_at_0_{label}": round(r0, 3)})
for label, t in (("file", s["face_file_t"]), ("pulse", s["face_t"])):
ok = np.r_[np.diff(t) > 0, False]
tt = np.arange(max(t[0], s["t"][0]), min(t[ok][-1], s["t"][-1]), 0.1)
fa = np.interp(tt, t[ok], np.abs(s["face_norm"][ok]))
sp = np.interp(tt, s["speed_100hz_t"], np.abs(s["speed_100hz"]))
L, r, r0 = best_lag(fa, sp, 600)
row.update({f"face_wheel_lag_{label}_s": L / 10, f"rf_{label}": round(r, 3)})
ft = s["face_file_t"]
row.update({"face_file_dt_ms": round(1000 * np.median(np.diff(ft)), 1), "face_file_nonmonotonic": int((np.diff(ft) <= 0).sum()),
"face_file_span_s": round(ft[ft > ft[0]].max() - ft[0], 1),
"last_two_equal_pulse0": bool(np.allclose(ft[-2:], s["face_t"][0])), "camera_span_s": round(s["face_t"][-1] - s["face_t"][0], 1),
"camera_dt_ms": round(1000 * np.median(np.diff(s["face_t"])), 1)})
clock.append(row)
clock = pd.DataFrame(clock)
print(clock.to_string(index=False))
session ca_wheel_lag_file_frames r_file r_at_0_file ca_wheel_lag_pulse_frames r_pulse r_at_0_pulse face_wheel_lag_file_s rf_file face_wheel_lag_pulse_s rf_pulse face_file_dt_ms face_file_nonmonotonic face_file_span_s last_two_equal_pulse0 camera_span_s camera_dt_ms grabAM05 11222019-spont -499 0.319 -0.015 0 0.312 0.312 60.0 0.058 0.0 0.281 111.1 2 2008.3 True 1807.9 100.0 grabAM05 11232019-airpuffs -498 0.333 -0.001 0 0.346 0.346 60.0 0.044 0.0 0.191 111.1 2 2908.3 True 2618.0 100.0 grabAM05 11262019-airpuffs -498 0.292 -0.022 1 0.340 0.339 44.5 0.142 1.3 0.350 111.1 2 2806.6 True 2526.4 100.0 grabAM06 11222019-spont -500 0.493 0.017 -2 0.493 0.489 -42.1 0.098 -0.3 0.515 111.1 2 1770.6 True 1593.9 100.0 grabAM06 11242019-airpuffs -499 0.378 -0.003 0 0.436 0.436 53.0 0.055 0.4 0.240 111.1 2 2738.6 True 2465.2 100.0 grabAM07 11232019-spont -500 0.352 -0.024 -1 0.370 0.365 31.9 0.158 0.8 0.358 111.1 2 1699.8 True 1530.2 100.0 grabAM07 11242019-airpuffs -500 0.352 -0.034 -2 0.359 0.345 58.2 0.135 1.4 0.329 111.1 2 2739.3 True 2465.9 100.0 grabAM08 11232019-spont -512 0.501 -0.002 -13 0.514 0.462 23.1 0.127 0.0 0.745 111.1 2 1771.2 True 1594.5 100.0 grabAM08 11242019-airpuffs -504 0.405 -0.012 -4 0.424 0.407 -54.3 0.161 0.5 0.590 111.1 2 2471.1 True 2224.5 100.0
Calcium (brain mean) and wheel speed are best aligned at lag -498 to -512 frames under the files' ΔF/F times and at -13 to 1 frames with ΔF/F frame i at blue pulse i + 499, which confirms the offset from the data. The FaceMap and pupil timestamps in the files are spaced 111 ms (9 Hz), span longer than the session and are not monotonic at the end; the camera pulses are 100 ms apart (10 Hz, Methods). With frame i at camera pulse i, as the code does, face movement and wheel speed align within 1.4 s; with the files' timestamps the best lag is 23 to 60 s either way and the correlation weaker. The analysis uses the pulses. In grabAM08 11232019 spont the best calcium-wheel lag under pulse timing is -13 frames (r = 0.514 against 0.462 at lag 0), calcium leading wheel speed by 1.3 s; in the same session the ΔF/F frame offset from the raw movie is 499 (previous check) and face movement and wheel speed align at lag 0, so its ΔF/F timing does not rest on the wheel. The lead is not explained here. Wheel speed is near 0 at rest and high during running, so its cross-correlation with calcium contrasts the two states; a peak within 13 frames of lag 0 under pulse timing, in every session, is what shows that the ΔF/F and the wheel share a clock.
# The locomotion change points re-derived here against the file's WheelActivation table (the
# code's onsets after its own selection); and the ranges and units of the variables used.
rng_rows = []
for s in sessions:
wa = s["wheel_activation"]
d_on = np.array([np.min(np.abs(s["wheel_on_all"] - x)) for x in wa[:, 0]])
d_off = np.array([np.min(np.abs(s["wheel_off_all"] - x)) for x in wa[:, 1]])
rng_rows.append({"session": f"{s['mouse']} {s['session']}", "WheelActivation": len(wa), "rederived_epochs": len(s["wheel_on_all"]),
"max_on_diff_s": d_on.max(), "max_off_diff_s": d_off.max(),
"running_fraction": round(np.sum(s["wheel_off_all"] - s["wheel_on_all"]) / (s["t"][-1] - s["t"][0]), 2),
"speed_min": round(s["speed_100hz"].min(), 3), "speed_max": round(s["speed_100hz"].max(), 3),
"speed_unit": s["wheel_unit"], "pupil_range": f"{s['pupil'].min():.0f} to {s['pupil'].max():.0f}",
"pupil_unit": s["pupil_unit"], "faceNorm_range": f"{s['face_norm'].min():.1f} to {s['face_norm'].max():.1f}",
"ach_pct_range": f"{100 * np.nanmin(s['ach'][:, LEFT]):.1f} to {100 * np.nanmax(s['ach'][:, LEFT]):.1f}",
"ca_pct_range": f"{100 * np.nanmin(s['ca'][:, LEFT]):.1f} to {100 * np.nanmax(s['ca'][:, LEFT]):.1f}"})
print(pd.DataFrame(rng_rows).to_string(index=False))
session WheelActivation rederived_epochs max_on_diff_s max_off_diff_s running_fraction speed_min speed_max speed_unit pupil_range pupil_unit faceNorm_range ach_pct_range ca_pct_range grabAM05 11222019-spont 8 12 2.273737e-13 2.273737e-13 0.08 -0.140 0.186 unknown 834 to 7375 unknown -3.3 to 62.9 -1.0 to 1.9 -3.3 to 7.2 grabAM05 11232019-airpuffs 6 16 2.273737e-13 5.684342e-14 0.05 -0.088 0.266 unknown 23 to 4257 unknown -2.6 to 203.6 -1.2 to 2.1 -3.5 to 5.5 grabAM05 11262019-airpuffs 3 18 0.000000e+00 0.000000e+00 0.05 -0.128 0.235 unknown 12 to 2408 unknown -2.8 to 634.0 -1.3 to 2.3 -3.5 to 5.8 grabAM06 11222019-spont 30 61 2.273737e-13 2.273737e-13 0.39 -0.107 0.354 unknown 679 to 6290 unknown -3.0 to 158.8 -1.4 to 2.8 -3.3 to 4.8 grabAM06 11242019-airpuffs 13 23 2.273737e-13 4.547474e-13 0.08 -0.086 0.253 unknown 175 to 7997 unknown -2.7 to 496.3 -1.8 to 3.1 -4.0 to 5.6 grabAM07 11232019-spont 14 39 1.136868e-13 2.273737e-13 0.31 -0.102 0.331 unknown 85 to 2226 unknown -2.8 to 165.5 -1.3 to 2.3 -3.0 to 5.9 grabAM07 11242019-airpuffs 26 91 4.547474e-13 2.273737e-13 0.50 -0.113 0.375 unknown 175 to 3896 unknown -2.3 to 158.2 -1.3 to 2.0 -3.5 to 4.8 grabAM08 11232019-spont 18 22 2.273737e-13 2.273737e-13 0.16 -0.077 0.233 unknown 77 to 3606 unknown -3.2 to 612.5 -2.8 to 4.5 -3.7 to 8.4 grabAM08 11242019-airpuffs 10 18 2.273737e-13 2.273737e-13 0.10 -0.070 0.239 unknown 39 to 2539 unknown -4.1 to 617.9 -2.7 to 5.0 -4.6 to 9.0
Every WheelActivation onset and offset in the files is one of the epochs re-derived here, to the sample, so the port of the code's change-point method matches the code's output. The wheel speed is signed (its minimum is -0.07 to -0.14, backward turns) and stated in "unknown" units; the code treats it as m/s (thresholds of 0.02 and 0.03 m/s) and its range is consistent with that, not proven. Pupil area has no unit (pixels by FaceMap).
# The parcellated series in the files: which ΔF/F it averages, with which masks, and its frame offset.
s0 = sessions[0]
url, f = meso.open_nwb(meso.asset_path(s0["mouse"], s0["session"]))
masks = f["processing/ophys/ImageSegmentation/ParcellatedPlaneSegmentation/image_mask"][:] > 0
parc = f["processing/ophys/Fluorescence/RoiResponseSeriesParcellated/data"][:1000, :10]
linked = f["processing/ophys/ImageSegmentation/ParcellatedPlaneSegmentation/imaging_plane"].name
start_parc = float(f["processing/ophys/Fluorescence/RoiResponseSeriesParcellated/starting_time"][()])
res = []
for colour, name in (("ach", "DFFOnePhotonSeriesBlue"), ("ca", "DFFOnePhotonSeriesGreen")):
d = f["processing/ophys/" + name + "/data"]
mov = np.concatenate(meso.reduce_many([(url, d, lambda b, k0: b.reshape(b.shape[0], -1), range(1300))])[0])[:1300]
for label, m in (("as stored", masks[:10]), ("transposed", np.transpose(masks[:10], (0, 2, 1)))):
tr = mov @ (m.reshape(10, -1) / m.reshape(10, -1).sum(1, keepdims=True)).T
cs = [np.mean([np.corrcoef(tr[o:o + 1000, k], parc[:, k])[0, 1] for k in range(10)]) for o in range(0, 301)]
o = int(np.argmax(cs))
res.append({"dff": colour, "masks": label, "best_offset_frames": o, "mean_r": round(cs[o], 4),
"sd_ratio": round(float(np.mean(tr[o:o + 1000].std(0) / parc.std(0))), 3)})
f.close()
print(pd.DataFrame(res).to_string(index=False))
print("imaging plane linked to the parcellation:", linked, "; starting time", round(start_parc, 3),
"s; green pulse 0 at", round(s0["green_ttl"][0], 3), "s; ROI counts per session:", sorted({str(s["parcellated_shape"][1]) for s in sessions}))
dff masks best_offset_frames mean_r sd_ratio ach as stored 150 0.6800 0.268 ach transposed 150 0.7153 0.406 ca as stored 150 0.8538 0.792 ca transposed 150 0.9988 0.993 imaging plane linked to the parcellation: /processing/ophys/ImageSegmentation/ParcellatedPlaneSegmentation/imaging_plane ; starting time 10.115 s; green pulse 0 at 10.115 s; ROI counts per session: ['100', '107', '89', '90', '93', '96', '97', '98']
The RoiResponseSeriesParcellated series is the calcium (green-excitation) ΔF/F averaged over
its ROI masks transposed, 150 frames after the ΔF/F's first frame (r near 1, SD ratio near 1). The
file links its segmentation to the ACh imaging plane, stores the masks transposed relative to the
ΔF/F, and gives it the starting time of the first green pulse, 64.9 s before its first frame. The
ROIs are a functional parcellation whose number differs between sessions, not the paper's Allen
parcels; it is not used here.
# Visual stimulus table of one session and the Methods' stimulus; duplicate check across sessions.
url, f = meso.open_nwb(meso.asset_path("grabAM05", meso.VIS_SESSIONS["grabAM05"]))
g = f["intervals/VisualStimulus"]
vis_tab = pd.DataFrame({k: list(g[k][:5]) for k in ("contrast", "orientation", "stimulus_frequency", "spatial_frequency", "stimulus_size")})
vis_tab["screen_coordinates"] = [list(np.round(r, 2)) for r in g["screen_coordinates"][:5]]
vis_tab["duration_s"] = np.round(g["stop_time"][:5] - g["start_time"][:5], 2)
print(vis_tab.to_string(index=False))
print("contrasts:", sorted(set(g["contrast"][:])))
f.close()
for mouse, ses in meso.VIS_SESSIONS.items():
url, f = meso.open_nwb(meso.asset_path(mouse, ses))
g = f["intervals/VisualStimulus"]
sc = g["screen_coordinates"]
print(f"{mouse} {ses}: columns {list(g.keys())}; screen_coordinates description {sc.attrs.get('description')!r}; "
f"distinct rows {np.unique(np.round(sc[:], 2), axis=0).tolist()}")
f.close()
brain = {f"{s['mouse']} {s['session']}": s["ca"][:2000, -1] for s in sessions}
keys = list(brain)
pair_r = max(abs(np.corrcoef(brain[a], brain[b])[0, 1]) for i, a in enumerate(keys) for b in keys[i + 1:])
print("largest |r| between the first 2000 frames of the calcium brain mean of two sessions:", round(pair_r, 3))
print("distinct ΔF/F frame counts:", info.dff_frames.nunique(), "of", len(info))
puff_run = []
for s in sessions:
if s["airpuff"] is None:
continue
for p in s["airpuff"][:, 0]:
puff_run.append(np.any((s["wheel_on_all"] > p - 2) & (s["wheel_on_all"] < p + 5)))
print("air puffs in the tables:", len(puff_run), "; with a locomotion onset from -2 to +5 s:", int(np.sum(puff_run)))
contrast orientation stimulus_frequency spatial_frequency stimulus_size screen_coordinates duration_s
5 0 2 0.04 2 [20.0, 667.97, 317.97, 932.03] 2.01
10 0 2 0.04 2 [20.0, 667.97, 317.97, 932.03] 2.00
50 0 2 0.04 2 [20.0, 667.97, 317.97, 932.03] 2.00
10 0 2 0.04 2 [20.0, 667.97, 317.97, 932.03] 2.00
50 0 2 0.04 2 [20.0, 667.97, 317.97, 932.03] 2.00
contrasts: [np.int64(2), np.int64(5), np.int64(10), np.int64(20), np.int64(50), np.int64(100)]
grabAM05 11222019-vis-stim: columns ['contrast', 'id', 'orientation', 'screen_coordinates', 'spatial_frequency', 'start_time', 'stimulus_frequency', 'stimulus_size', 'stop_time']; screen_coordinates description 'Screen coordinates that define the position of the stimulus on the monitor, i.e. [x1,y1,x2,y2].Left (x1): The x-coordinate of the left edge of the rectangle.Top (y1): The y-coordinate of the top edge of the rectangle.Right (x2): The x-coordinate of the right edge of the rectangle.Bottom (y2): The y-coordinate of the bottom edge of the rectangle.'; distinct rows [[20.0, 667.97, 317.97, 932.03]]
grabAM06 11222019-vis-stim: columns ['contrast', 'id', 'orientation', 'screen_coordinates', 'spatial_frequency', 'start_time', 'stimulus_frequency', 'stimulus_size', 'stop_time']; screen_coordinates description 'Screen coordinates that define the position of the stimulus on the monitor, i.e. [x1,y1,x2,y2].Left (x1): The x-coordinate of the left edge of the rectangle.Top (y1): The y-coordinate of the top edge of the rectangle.Right (x2): The x-coordinate of the right edge of the rectangle.Bottom (y2): The y-coordinate of the bottom edge of the rectangle.'; distinct rows [[20.0, 667.97, 317.97, 932.03]]
grabAM07 11232019-vis-stim: columns ['contrast', 'id', 'orientation', 'screen_coordinates', 'spatial_frequency', 'start_time', 'stimulus_frequency', 'stimulus_size', 'stop_time']; screen_coordinates description 'Screen coordinates that define the position of the stimulus on the monitor, i.e. [x1,y1,x2,y2].Left (x1): The x-coordinate of the left edge of the rectangle.Top (y1): The y-coordinate of the top edge of the rectangle.Right (x2): The x-coordinate of the right edge of the rectangle.Bottom (y2): The y-coordinate of the bottom edge of the rectangle.'; distinct rows [[20.0, 667.97, 317.97, 932.03]]
grabAM08 11232019-vis-stim: columns ['contrast', 'id', 'orientation', 'screen_coordinates', 'spatial_frequency', 'start_time', 'stimulus_frequency', 'stimulus_size', 'stop_time']; screen_coordinates description 'Screen coordinates that define the position of the stimulus on the monitor, i.e. [x1,y1,x2,y2].Left (x1): The x-coordinate of the left edge of the rectangle.Top (y1): The y-coordinate of the top edge of the rectangle.Right (x2): The x-coordinate of the right edge of the rectangle.Bottom (y2): The y-coordinate of the bottom edge of the rectangle.'; distinct rows [[20.0, 667.97, 317.97, 932.03]] largest |r| between the first 2000 frames of the calcium brain mean of two sessions: 0.211 distinct ΔF/F frame counts: 9 of 9 air puffs in the tables: 29 ; with a locomotion onset from -2 to +5 s: 0
In the VisualStimulus table the stimulus size is 2 and the first screen coordinate is 20: the
Methods give 20 degree gratings shown for 1 to 2 s, so the columns from stimulus size on are
shifted by one (size holds the duration, the first coordinate holds the size). No two sessions
share a recording. No air puff in the tables is followed by a locomotion onset within 5 s, as the
code's selection of puffs would give (process_spike2_two_cams.m drops those); the tables then
hold only the puffs the code kept, which the files do not say (consistent with the data, not
shown).
# States and events per session, with the code's parameters.
state_rows, matched, dur_rows = [], {}, []
for s in sessions:
st = meso.sustained_states(s)
m = meso.match_epochs(st)
matched[(s["mouse"], s["session"])] = m
s["sit"] = st["sit"]
s["face_onsets"] = meso.face_onsets(s)
s["face_onsets_all"] = meso.face_onsets(s, include_last=True)
for k in ("loc", "face_high", "face_low"):
dur_rows.append({"mouse": s["mouse"], "state": k, "epochs": len(st[k]), "seconds": float(np.sum(st[k][:, 1] - st[k][:, 0]))})
s["loc_onsets"] = meso.locomotion_onsets(s)
state_rows.append({"session": f"{s['mouse']} {s['session']}", "loc": len(st["loc"]), "face_high": len(st["face_high"]),
"face_low": len(st["face_low"]), "matched_per_state": 0 if m is None else len(m["loc"]),
"matched_seconds": 0 if m is None else round(float(np.sum(m["loc"][:, 1] - m["loc"][:, 0])), 1),
"face_onsets": len(s["face_onsets"]), "face_onsets_last_kept": len(s["face_onsets_all"]),
"locomotion_onsets": len(s["loc_onsets"])})
state_rows = pd.DataFrame(state_rows)
print(state_rows.to_string(index=False))
print("\nper mouse:")
state_rows["mouse"] = state_rows.session.str.split().str[0]
print(state_rows.groupby("mouse")[["matched_per_state", "matched_seconds", "face_onsets", "face_onsets_last_kept", "locomotion_onsets"]].sum().to_string())
print("\nmatched epochs have the same duration in each state, so matched_seconds holds for every state.")
print("all sustained epochs per mouse and state, before matching (count and seconds):")
print(pd.DataFrame(dur_rows).groupby(["mouse", "state"]).sum().unstack("state").round(1).to_string())
session loc face_high face_low matched_per_state matched_seconds face_onsets face_onsets_last_kept locomotion_onsets
grabAM05 11222019-spont 4 4 6 4 29.2 8 9 8
grabAM05 11232019-airpuffs 3 17 30 3 20.7 12 15 6
grabAM05 11262019-airpuffs 3 28 35 3 21.8 11 23 3
grabAM06 11222019-spont 21 0 0 0 0.0 0 0 28
grabAM06 11242019-airpuffs 2 11 21 2 14.6 25 25 11
grabAM07 11232019-spont 17 0 12 0 0.0 11 11 13
grabAM07 11242019-airpuffs 29 0 6 0 0.0 1 1 25
grabAM08 11232019-spont 5 5 23 5 39.2 15 14 17
grabAM08 11242019-airpuffs 8 10 27 8 54.5 15 15 8
per mouse:
matched_per_state matched_seconds face_onsets face_onsets_last_kept locomotion_onsets
mouse
grabAM05 10 71.7 31 47 17
grabAM06 2 14.6 25 25 39
grabAM07 0 0.0 12 12 38
grabAM08 13 93.7 30 29 25
matched epochs have the same duration in each state, so matched_seconds holds for every state.
all sustained epochs per mouse and state, before matching (count and seconds):
epochs seconds
state face_high face_low loc face_high face_low loc
mouse
grabAM05 49 71 10 515.4 589.8 91.4
grabAM06 11 21 23 73.8 136.4 191.6
grabAM07 0 18 46 0.0 136.2 751.0
grabAM08 15 50 13 119.1 410.3 106.7
grabAM07 runs during 31% and 50% of its two sessions and has no sustained high face epoch, so
the matched-state analyses (Figures 3c, 3d, 4 and 5) rest on 3 mice, and grabAM06 contributes 2
epochs per state (14.6 s per state, against 71.7 s for grabAM05 and 93.7 s for grabAM08), with the
same weight as the other mice in every test across mice. A section at the end repeats those tests
without grabAM06. Face onsets drop the quiescent period after the last locomotion epoch, as the
code's extractFacePupilTStamps.m does (sitOn = [0; wheelOff(1:end-1)]), while its
getSustainedStates.m keeps it; keeping it adds 16 onsets in grabAM05, none in grabAM06 and grabAM07, and in
grabAM08 gives one fewer (29 against 30; 14 against 15 in its spontaneous session). An added period
can remove an onset because the thresholds are the 60th and 40th percentiles of face movement over
all quiescent frames (meso.face_onsets, as in the code), so the added frames move them. The
variant is under s7-face-onset-transients.
allen-parcel-link¶
The paper's Allen parcel masks, from its code, can be laid on the dandiset's ΔF/F pixels, so that the areas the paper names can be found in the data.
vis_rows, vis_maps, vis_agree = [], {}, []
for mouse in MICE:
path = meso.asset_path(mouse, meso.VIS_SESSIONS[mouse])
url, f = meso.open_nwb(path)
blue = f["acquisition/TTLsTable/timestamp"][:][f["acquisition/TTLsTable/ttl_type"][:] == 0]
d = f["processing/ophys/DFFOnePhotonSeriesGreen/data"]
n = d.shape[0]
t = blue[meso.DETREND_DROP:meso.DETREND_DROP + n]
g = f["intervals/VisualStimulus"]
on = g["start_time"][:][g["contrast"][:] == 100]
raw = raw_means(path, n=1300, series="OnePhotonSeriesGreenExcitationRedChannel")
w = np.nan_to_num(W)
db = f["processing/ophys/DFFOnePhotonSeriesBlue/data"]
def red_green(b, k0):
flat = b.reshape(b.shape[0], -1)
return k0, flat @ w, (flat != 0).any(0)
parts, parts_blue = meso.reduce_many([(url, d, red_green, None), (url, db, lambda b, k0: (b.reshape(b.shape[0], -1) != 0).any(0), None)])
f.close()
x = np.full((n, W.shape[1]), np.nan)
for k0, v, _nz in parts:
x[k0:k0 + len(v)] = v
for colour, nzs in (("ach", parts_blue), ("ca", [p[2] for p in parts])):
nz = np.any(nzs, axis=0).reshape(256, 256)
vis_agree.append({"session": f"{mouse} {meso.VIS_SESSIONS[mouse]}", "colour": colour, "nonzero_pixels": int(nz.sum()),
**{f"rot{k}{'_T' if fl else ''}": round((nz == m).mean(), 4) for (k, fl), m in orientations(parcels["mask"])}})
o, c, c499 = best_offset(raw, x[:, -1], n_dff=700)
resp = []
for e in on:
pre = x[(t >= e - 2) & (t < e)]
post = x[(t >= e) & (t < e + 1)]
if len(pre) == 20 and len(post) == 10:
resp.append(post.mean(0) - pre.mean(0))
n_used = len(resp)
resp = 100 * np.nanmean(resp, axis=0)[:-1]
resp[np.isnan(W[0, :-1])] = np.nan
vis_maps[mouse] = resp
order = np.argsort(-np.nan_to_num(resp, nan=-np.inf))
lv = sel.index[(sel.hemi == "L") & (sel.area == "VISp")][0]
rv = sel.index[(sel.hemi == "R") & (sel.area == "VISp")][0]
vis_rows.append({"mouse": mouse, "session": meso.VIS_SESSIONS[mouse], "trials_100pct": len(on), "trials_used": n_used, "dff_offset_frames": o, "offset_r": round(c, 3), "r_at_499": round(c499, 3),
"top3": ", ".join(sel.name[order[:3]]), "L_VISp": round(resp[lv], 3), "R_VISp": round(resp[rv], 3),
"rank_L_VISp": int(np.flatnonzero(order == lv)[0]) + 1,
"top_is_left_visual": sel.hemi[order[0]] == "L" and sel.area[order[0]] in meso.GROUPS["visual"]})
vis_rows = pd.DataFrame(vis_rows)
vis_agree = pd.DataFrame(vis_agree)
print("Visual stimulus sessions, each movie separately (pixels non-zero in any frame):")
print(vis_agree.to_string(index=False))
all_agree = pd.concat([agree, vis_agree])
others = all_agree.drop(columns=["session", "colour", "nonzero_pixels", "rot180"]).max().max()
print(f"brain mask agreement, rotation 180, per movie over all 13 sessions: min {all_agree.rot180.min()}, movies at 1.0: "
f"{(all_agree.rot180 == 1).sum()} of {len(all_agree)}; per session (union, the 9 state sessions): min {union.rot180.min()}; "
f"best other orientation {others}")
print(vis_rows.to_string(index=False))
print("dff_offset_frames: frame offset of the green ΔF/F against the raw green movie, from frame-to-frame changes of their means")
print("response range over parcels and mice (percent ΔF/F):", round(np.nanmin(list(vis_maps.values())), 3), "to", round(np.nanmax(list(vis_maps.values())), 3))
Visual stimulus sessions, each movie separately (pixels non-zero in any frame):
session colour nonzero_pixels rot0 rot0_T rot90 rot90_T rot180 rot180_T rot270 rot270_T
grabAM05 11222019-vis-stim ach 25414 0.8317 0.7681 0.7708 0.8407 1.0 0.7739 0.7708 0.9818
grabAM05 11222019-vis-stim ca 25414 0.8317 0.7681 0.7708 0.8407 1.0 0.7739 0.7708 0.9818
grabAM06 11222019-vis-stim ach 25414 0.8317 0.7681 0.7708 0.8407 1.0 0.7739 0.7708 0.9818
grabAM06 11222019-vis-stim ca 25414 0.8317 0.7681 0.7708 0.8407 1.0 0.7739 0.7708 0.9818
grabAM07 11232019-vis-stim ach 25414 0.8317 0.7681 0.7708 0.8407 1.0 0.7739 0.7708 0.9818
grabAM07 11232019-vis-stim ca 25414 0.8317 0.7681 0.7708 0.8407 1.0 0.7739 0.7708 0.9818
grabAM08 11232019-vis-stim ach 25414 0.8317 0.7681 0.7708 0.8407 1.0 0.7739 0.7708 0.9818
grabAM08 11232019-vis-stim ca 25414 0.8317 0.7681 0.7708 0.8407 1.0 0.7739 0.7708 0.9818
brain mask agreement, rotation 180, per movie over all 13 sessions: min 0.996, movies at 1.0: 25 of 26; per session (union, the 9 state sessions): min 1.0; best other orientation 0.9818
mouse session trials_100pct trials_used dff_offset_frames offset_r r_at_499 top3 L_VISp R_VISp rank_L_VISp top_is_left_visual
grabAM05 11222019-vis-stim 50 47 499 0.981 0.981 L-VlSp, L-VlSam, L-VlSpm 0.211 0.039 1 True
grabAM06 11222019-vis-stim 50 49 499 0.972 0.972 L-VlSp, L-Tea, L-VlSli 0.108 0.065 1 True
grabAM07 11232019-vis-stim 50 49 499 0.978 0.978 L-VlSpm, L-VlSam, L-VlSp 0.216 0.054 3 True
grabAM08 11232019-vis-stim 50 49 499 0.946 0.946 L-VlSal, L-VlSp, L-Tea 0.215 -0.024 2 True
dff_offset_frames: frame offset of the green ΔF/F against the raw green movie, from frame-to-frame changes of their means
response range over parcels and mice (percent ΔF/F): -0.11 to 0.291
fig, axes = plt.subplots(1, 4, figsize=(14, 3.5))
for ax, mouse in zip(axes, MICE):
img = np.full((256, 256), np.nan)
for i, p in enumerate(parcels["sel"]):
img[meso.orient(parcels["indicators"][:, :, p["index"]])] = vis_maps[mouse][i]
im = ax.imshow(img, cmap="viridis")
ax.set_title(f"{mouse}: calcium, 100% contrast")
ax.axis("off")
plt.colorbar(im, ax=ax, fraction=0.046, label="% ΔF/F")
plt.tight_layout()
plt.show()
Result: partially reproduced. Read per session (pixels non-zero in either movie), the non-zero pixels equal the code's brain mask rotated by 180 degrees in all 9 state sessions (agreement 1.0, best other orientation 0.98). Read per movie, 25 of the 26 movies of the 13 sessions (including both movies of the 4 visual stimulus sessions) agree at 1.0; one calcium movie agrees at 0.996, below the criterion's 99.9%, because 263 pixels inside the mask, all in the right AUD parcel, are zero in every frame. This attempt takes the per-session reading (a deviation); under the per-movie reading the mask part fails on that one movie. With 100% contrast gratings the largest calcium response is in a left visual parcel in all 4 mice and left VISp exceeds right VISp in all 4, but left VISp itself ranks first in 2 of 4 mice (ranks 1, 1, 3 and 2), short of the 3 the criterion asks for. The parcels line up with the data; their exact borders are not tested beyond this. The gratings are at the same screen coordinates in every session; the files do not say on which side of the mouse the monitor stood, so the expected hemisphere is not checked from the data.
Regression of ACh and calcium on behaviour (Figure 2)¶
fig2a: the variance explained by FaceMap PCs 1 to 25, locomotion and pupil differs across cortical areas. fig2b: the anterior-posterior gradients are inverted between ACh and calcium. fig2c: face PC1 is the best single predictor of both signals.
Per session, the predictors are taken at the ΔF/F frame times (FaceMap and pupil interpolated from
the camera pulses, wheel speed averaged to 100 Hz and interpolated); cv_r2 scales them with the
training folds' mean and SD. The outputs are the 23 left parcels. Shuffled variables are permuted in
time with a fixed seed. The wheel predictor is the stored WheelSpeed, which is signed (backward turns
are negative, data check). The Methods compute locomotion speed "as a differential of distance", and
the code's wheel_changepoints.m keeps its sign (db1SpeedMpS), taking the absolute value only for
the threshold on mean epoch speed; the change points re-derived from WheelSpeed with that method
match the file's WheelActivation table (data check). The regression code itself is not public.
|speed| is shown as a variant after the fig2c claims, followed by a comparison with the paper's own
per-mouse values.
rng_reg = np.random.default_rng(2023)
rng_abs = np.random.default_rng(2022)
reg = []
for s in sessions:
keep = (s["t"] >= s["face_t"][0]) & (s["t"] <= s["face_t"][-1])
t = s["t"][keep]
face = np.column_stack([np.interp(t, s["face_t"], s["face_pcs"][:, k]) for k in range(25)])
pupil = np.interp(t, s["face_t"], s["pupil"])
wheel = s["wheel_frame"][keep]
X = np.column_stack([face, pupil, wheel])
Xv = np.column_stack([face, pupil, np.abs(wheel)])
cols = {"face_pc1": [0], "pupil": [25], "wheel": [26]}
for colour in ("ach", "ca"):
Y = s[colour][keep][:, LEFT]
row = {"mouse": s["mouse"], "session": s["session"], "colour": colour, "frames": int(keep.sum())}
row["full"] = meso.cv_r2(X, Y)
for name, c in cols.items():
row["single_" + name] = meso.cv_r2(X[:, c], Y)
Xs = X.copy()
Xs[:, c] = Xs[rng_reg.permutation(len(Xs))][:, c]
row["delta_" + name] = row["full"] - meso.cv_r2(Xs, Y)
# variant: |speed| as the wheel predictor
row["full_abs"] = meso.cv_r2(Xv, Y)
row["single_wheel_abs"] = meso.cv_r2(Xv[:, [26]], Y)
for name, c in cols.items():
Xs = Xv.copy()
Xs[:, c] = Xs[rng_abs.permutation(len(Xs))][:, c]
row[f"delta_{name}_abs"] = row["full_abs"] - meso.cv_r2(Xs, Y)
reg.append(row)
print("sessions fitted:", len(reg) // 2, "; frames per session:", sorted({r["frames"] for r in reg}))
def per_mouse(colour, key, clip=False):
out = {}
for mouse in MICE:
v = np.mean([r[key] for r in reg if r["mouse"] == mouse and r["colour"] == colour], axis=0)
out[mouse] = np.clip(v, 0, None) if clip else v
return pd.DataFrame(out, index=left_areas)
full = {c: per_mouse(c, "full") for c in ("ach", "ca")}
for c in ("ach", "ca"):
print(f"\nfull-model cross-validated R2, {c}, per left parcel and mouse (sessions averaged):")
print(full[c].round(3).assign(ap_rank=ap_left).sort_values("ap_rank").to_string())
print("range", round(full[c].values.min(), 3), "to", round(full[c].values.max(), 3))
sessions fitted: 9 ; frames per session: [14342, 14931, 14959, 17099, 21251, 23677, 23686, 24265, 25204]
full-model cross-validated R2, ach, per left parcel and mouse (sessions averaged):
grabAM05 grabAM06 grabAM07 grabAM08 ap_rank
MOs 0.323 0.397 0.286 0.215 1
MOp 0.303 0.372 0.241 0.227 2
SSp-m 0.119 0.292 0.116 0.236 3
SSp-ul 0.227 0.294 0.151 0.231 4
SSp-n 0.116 0.285 0.076 0.234 5
SSp-ll 0.249 0.289 0.158 0.176 6
SSp-un 0.208 0.294 0.131 0.213 7
SSp-tr 0.250 0.268 0.168 0.139 8
SSp-bfd 0.204 0.280 0.135 0.172 9
SSs 0.111 0.185 0.053 0.143 10
PTLp 0.243 0.268 0.155 0.120 11
VISam 0.243 0.255 0.141 0.100 12
AUD 0.101 0.155 0.031 0.082 13
RSPagl 0.236 0.249 0.157 0.121 14
RSPd 0.245 0.268 0.137 0.123 15
VISal 0.155 0.184 0.112 0.088 16
VISpm 0.230 0.255 0.148 0.130 17
TEa 0.092 0.138 0.030 0.066 18
VISli 0.098 0.165 0.052 0.065 19
VISl 0.127 0.164 0.114 0.085 20
VISp 0.187 0.272 0.168 0.138 21
VISpor 0.064 0.118 0.056 0.039 22
VISpl 0.045 0.083 0.083 0.077 23
range 0.03 to 0.397
full-model cross-validated R2, ca, per left parcel and mouse (sessions averaged):
grabAM05 grabAM06 grabAM07 grabAM08 ap_rank
MOs 0.226 0.304 0.181 0.175 1
MOp 0.203 0.253 0.188 0.187 2
SSp-m 0.111 0.153 0.114 0.093 3
SSp-ul 0.196 0.246 0.187 0.193 4
SSp-n 0.113 0.199 0.128 0.084 5
SSp-ll 0.254 0.306 0.231 0.291 6
SSp-un 0.188 0.275 0.171 0.144 7
SSp-tr 0.269 0.319 0.246 0.285 8
SSp-bfd 0.192 0.282 0.214 0.137 9
SSs 0.186 0.285 0.187 0.172 10
PTLp 0.261 0.320 0.253 0.254 11
VISam 0.240 0.292 0.251 0.273 12
AUD 0.215 0.361 0.184 0.228 13
RSPagl 0.246 0.374 0.268 0.310 14
RSPd 0.217 0.360 0.227 0.308 15
VISal 0.221 0.276 0.210 0.217 16
VISpm 0.226 0.332 0.263 0.295 17
TEa 0.197 0.348 0.204 0.246 18
VISli 0.195 0.343 0.214 0.252 19
VISl 0.203 0.298 0.244 0.246 20
VISp 0.229 0.387 0.284 0.304 21
VISpor 0.181 0.363 0.243 0.257 22
VISpl 0.185 0.332 0.269 0.275 23
range 0.084 to 0.387
for c in ("ach", "ca"):
fr = stats.friedmanchisquare(*full[c].values)
print(f"{c}: Friedman across 23 parcels, mice as blocks (n = {full[c].shape[1]}): chi2 = {fr.statistic:.2f}, p = {fr.pvalue:.2g}; "
f"mean R2 over parcels per mouse {full[c].mean().round(3).to_dict()}")
ach: Friedman across 23 parcels, mice as blocks (n = 4): chi2 = 67.98, p = 1.4e-06; mean R2 over parcels per mouse {'grabAM05': 0.182, 'grabAM06': 0.24, 'grabAM07': 0.126, 'grabAM08': 0.14}
ca: Friedman across 23 parcels, mice as blocks (n = 4): chi2 = 63.51, p = 6.7e-06; mean R2 over parcels per mouse {'grabAM05': 0.207, 'grabAM06': 0.305, 'grabAM07': 0.216, 'grabAM08': 0.227}
Result: reproduced. Full-model R2 differs across the 23 left parcels for ACh (Friedman p = 1.4e-06) and calcium (p = 6.7e-06) with 4 mice.
for c in ("ach", "ca"):
rho, p = stats.spearmanr(ap_left, full[c].mean(axis=1))
per = {m: round(stats.spearmanr(ap_left, full[c][m])[0], 2) for m in MICE}
print(f"{c}: Spearman of mean R2 with anterior-posterior rank (1 = MOs): rho = {rho:.3f}, p = {p:.2g}; per mouse rho {per}")
fig, axes = plt.subplots(1, 2, figsize=(10, 3.5), sharey=True)
for ax, c in zip(axes, ("ach", "ca")):
o = np.argsort(ap_left)
for m in MICE:
ax.plot(ap_left[o], full[c][m].values[o], ".", alpha=0.5, label=m)
ax.plot(ap_left[o], full[c].mean(axis=1).values[o], "k-", label="mean")
ax.set_xlabel("anterior-posterior rank (1 = MOs)")
ax.set_title(c)
axes[0].set_ylabel("cross-validated R2")
axes[1].legend(fontsize=7)
plt.tight_layout()
plt.show()
ach: Spearman of mean R2 with anterior-posterior rank (1 = MOs): rho = -0.761, p = 2.5e-05; per mouse rho {'grabAM05': np.float64(-0.58), 'grabAM06': np.float64(-0.85), 'grabAM07': np.float64(-0.44), 'grabAM08': np.float64(-0.87)}
ca: Spearman of mean R2 with anterior-posterior rank (1 = MOs): rho = 0.555, p = 0.0059; per mouse rho {'grabAM05': np.float64(-0.01), 'grabAM06': np.float64(0.7), 'grabAM07': np.float64(0.65), 'grabAM08': np.float64(0.57)}
# Variant: spontaneous sessions only (one per mouse).
for c in ("ach", "ca"):
sp = pd.DataFrame({m: [r["full"] for r in reg if r["mouse"] == m and r["colour"] == c and "spont" in r["session"]][0] for m in MICE},
index=left_areas)
fr = stats.friedmanchisquare(*sp.values)
rho, p = stats.spearmanr(ap_left, sp.mean(axis=1))
print(f"spontaneous only, {c}: Friedman chi2 = {fr.statistic:.2f}, p = {fr.pvalue:.2g}; Spearman with AP rank rho = {rho:.3f}, p = {p:.2g}")
spontaneous only, ach: Friedman chi2 = 70.27, p = 6e-07; Spearman with AP rank rho = -0.897, p = 6.7e-09 spontaneous only, ca: Friedman chi2 = 42.77, p = 0.005; Spearman with AP rank rho = 0.655, p = 0.00069
Result: reproduced. Mean R2 falls from anterior to posterior for ACh (rho = -0.761, p = 2.5e-05) and rises for calcium (rho = 0.555, p = 0.0059). Spontaneous sessions alone give the same signs (rho -0.897 and 0.655). For calcium one mouse (grabAM05) has no gradient (rho -0.01).
def variable_table(colour, prefix, abs_wheel=False):
rows = {}
for name in ("face_pc1", "pupil", "wheel"):
key = f"{prefix}_{name}"
if abs_wheel and (prefix == "delta" or name == "wheel"):
key += "_abs"
rows[name] = per_mouse(colour, key, clip=(prefix == "single")).mean()
return pd.DataFrame(rows)
def friedman_exact_p(values):
"""Exact p of Friedman's test without ties, enumerating every assignment of ranks within blocks."""
n, k = values.shape
def stat(R):
return 12 / (n * k * (k + 1)) * (R.sum(0) ** 2).sum() - 3 * n * (k + 1)
obs = stat(stats.rankdata(values, axis=1))
perms = np.array(list(itertools.permutations(range(1, k + 1))))
hits = [stat(perms[list(c)]) >= obs - 1e-9 for c in itertools.product(range(len(perms)), repeat=n)]
return float(np.mean(hits))
def judge_variables(colour, abs_wheel=False):
for prefix in ("single", "delta"):
tab = variable_table(colour, prefix, abs_wheel)
fr = stats.friedmanchisquare(*tab.values.T)
best = tab.idxmax(axis=1)
print(f"{colour} {prefix} R2 (mean over 23 left parcels{', negatives set to 0' if prefix == 'single' else ''}"
f"{'; wheel is |speed|' if abs_wheel else ''}):")
print(tab.round(4).to_string())
print(f" mean over mice: {tab.mean().round(4).to_dict()}; best variable per mouse: {best.to_dict()}; "
f"Friedman chi2 = {fr.statistic:.2f}, p = {fr.pvalue:.3g} (exact p = {friedman_exact_p(tab.values):.4f})")
for other in ("pupil", "wheel"):
w = stats.wilcoxon(tab.face_pc1, tab[other])
print(f" face_pc1 above {other} in {(tab.face_pc1 > tab[other]).sum()} of {len(tab)} mice "
f"(margin per mouse {(tab.face_pc1 - tab[other]).round(3).to_dict()}); signed-rank p = {w.pvalue:.3f}")
judge_variables("ach")
ach single R2 (mean over 23 left parcels, negatives set to 0):
face_pc1 pupil wheel
grabAM05 0.0866 0.0001 0.0017
grabAM06 0.0684 0.0079 0.0296
grabAM07 0.0475 0.0000 0.0115
grabAM08 0.1213 0.0287 0.0601
mean over mice: {'face_pc1': 0.0809, 'pupil': 0.0092, 'wheel': 0.0257}; best variable per mouse: {'grabAM05': 'face_pc1', 'grabAM06': 'face_pc1', 'grabAM07': 'face_pc1', 'grabAM08': 'face_pc1'}; Friedman chi2 = 8.00, p = 0.0183 (exact p = 0.0046)
face_pc1 above pupil in 4 of 4 mice (margin per mouse {'grabAM05': 0.087, 'grabAM06': 0.06, 'grabAM07': 0.047, 'grabAM08': 0.093}); signed-rank p = 0.125
face_pc1 above wheel in 4 of 4 mice (margin per mouse {'grabAM05': 0.085, 'grabAM06': 0.039, 'grabAM07': 0.036, 'grabAM08': 0.061}); signed-rank p = 0.125
ach delta R2 (mean over 23 left parcels):
face_pc1 pupil wheel
grabAM05 0.0728 0.0157 0.0024
grabAM06 0.1067 0.0303 0.0052
grabAM07 0.0803 0.0157 0.0047
grabAM08 0.0710 0.0015 0.0025
mean over mice: {'face_pc1': 0.0827, 'pupil': 0.0158, 'wheel': 0.0037}; best variable per mouse: {'grabAM05': 'face_pc1', 'grabAM06': 'face_pc1', 'grabAM07': 'face_pc1', 'grabAM08': 'face_pc1'}; Friedman chi2 = 6.50, p = 0.0388 (exact p = 0.0417)
face_pc1 above pupil in 4 of 4 mice (margin per mouse {'grabAM05': 0.057, 'grabAM06': 0.076, 'grabAM07': 0.065, 'grabAM08': 0.069}); signed-rank p = 0.125
face_pc1 above wheel in 4 of 4 mice (margin per mouse {'grabAM05': 0.07, 'grabAM06': 0.102, 'grabAM07': 0.076, 'grabAM08': 0.068}); signed-rank p = 0.125
Result: reproduced. Face PC1 has the highest single-variable R2 in all 4 mice (mean 0.081, against 0.026 for wheel speed and 0.009 for pupil), Friedman p = 0.018 (exact p = 0.0046). Its unique contribution is also the largest in all 4 mice (Friedman p = 0.039, exact p = 0.042). The signed-rank tests give p = 0.125, the smallest possible with 4 mice.
judge_variables("ca")
ca single R2 (mean over 23 left parcels, negatives set to 0):
face_pc1 pupil wheel
grabAM05 0.0780 0.0008 0.0253
grabAM06 0.1419 0.0258 0.1263
grabAM07 0.0895 0.0000 0.0638
grabAM08 0.1386 0.0140 0.1211
mean over mice: {'face_pc1': 0.112, 'pupil': 0.0102, 'wheel': 0.0841}; best variable per mouse: {'grabAM05': 'face_pc1', 'grabAM06': 'face_pc1', 'grabAM07': 'face_pc1', 'grabAM08': 'face_pc1'}; Friedman chi2 = 8.00, p = 0.0183 (exact p = 0.0046)
face_pc1 above pupil in 4 of 4 mice (margin per mouse {'grabAM05': 0.077, 'grabAM06': 0.116, 'grabAM07': 0.089, 'grabAM08': 0.125}); signed-rank p = 0.125
face_pc1 above wheel in 4 of 4 mice (margin per mouse {'grabAM05': 0.053, 'grabAM06': 0.016, 'grabAM07': 0.026, 'grabAM08': 0.017}); signed-rank p = 0.125
ca delta R2 (mean over 23 left parcels):
face_pc1 pupil wheel
grabAM05 0.0300 0.0281 0.0123
grabAM06 0.0527 0.0062 0.0206
grabAM07 0.0717 0.0231 0.0154
grabAM08 0.0763 0.0441 0.0076
mean over mice: {'face_pc1': 0.0577, 'pupil': 0.0254, 'wheel': 0.014}; best variable per mouse: {'grabAM05': 'face_pc1', 'grabAM06': 'face_pc1', 'grabAM07': 'face_pc1', 'grabAM08': 'face_pc1'}; Friedman chi2 = 6.50, p = 0.0388 (exact p = 0.0417)
face_pc1 above pupil in 4 of 4 mice (margin per mouse {'grabAM05': 0.002, 'grabAM06': 0.047, 'grabAM07': 0.049, 'grabAM08': 0.032}); signed-rank p = 0.125
face_pc1 above wheel in 4 of 4 mice (margin per mouse {'grabAM05': 0.018, 'grabAM06': 0.032, 'grabAM07': 0.056, 'grabAM08': 0.069}); signed-rank p = 0.125
Variant: the wheel predictor as |speed| (backward turns counted as running forward). The single-variable and unique-contribution tables are repeated with it in place of the signed speed (the face and pupil single-variable models do not change).
for c in ("ach", "ca"):
judge_variables(c, abs_wheel=True)
fa = per_mouse(c, "full_abs")
fr = stats.friedmanchisquare(*fa.values)
rho, p = stats.spearmanr(ap_left, fa.mean(axis=1))
print(f"{c} full model with |speed|: Friedman across parcels p = {fr.pvalue:.2g}; Spearman with AP rank rho = {rho:.3f}, p = {p:.2g}; "
f"mean R2 change from signed speed {(fa - full[c]).values.mean():.4f}\n")
ach single R2 (mean over 23 left parcels, negatives set to 0; wheel is |speed|):
face_pc1 pupil wheel
grabAM05 0.0866 0.0001 0.0031
grabAM06 0.0684 0.0079 0.0349
grabAM07 0.0475 0.0000 0.0119
grabAM08 0.1213 0.0287 0.0628
mean over mice: {'face_pc1': 0.0809, 'pupil': 0.0092, 'wheel': 0.0282}; best variable per mouse: {'grabAM05': 'face_pc1', 'grabAM06': 'face_pc1', 'grabAM07': 'face_pc1', 'grabAM08': 'face_pc1'}; Friedman chi2 = 8.00, p = 0.0183 (exact p = 0.0046)
face_pc1 above pupil in 4 of 4 mice (margin per mouse {'grabAM05': 0.087, 'grabAM06': 0.06, 'grabAM07': 0.047, 'grabAM08': 0.093}); signed-rank p = 0.125
face_pc1 above wheel in 4 of 4 mice (margin per mouse {'grabAM05': 0.083, 'grabAM06': 0.034, 'grabAM07': 0.036, 'grabAM08': 0.058}); signed-rank p = 0.125
ach delta R2 (mean over 23 left parcels; wheel is |speed|):
face_pc1 pupil wheel
grabAM05 0.0698 0.0160 0.0029
grabAM06 0.1062 0.0316 0.0064
grabAM07 0.0791 0.0159 0.0050
grabAM08 0.0692 0.0015 0.0023
mean over mice: {'face_pc1': 0.0811, 'pupil': 0.0163, 'wheel': 0.0042}; best variable per mouse: {'grabAM05': 'face_pc1', 'grabAM06': 'face_pc1', 'grabAM07': 'face_pc1', 'grabAM08': 'face_pc1'}; Friedman chi2 = 6.50, p = 0.0388 (exact p = 0.0417)
face_pc1 above pupil in 4 of 4 mice (margin per mouse {'grabAM05': 0.054, 'grabAM06': 0.075, 'grabAM07': 0.063, 'grabAM08': 0.068}); signed-rank p = 0.125
face_pc1 above wheel in 4 of 4 mice (margin per mouse {'grabAM05': 0.067, 'grabAM06': 0.1, 'grabAM07': 0.074, 'grabAM08': 0.067}); signed-rank p = 0.125
ach full model with |speed|: Friedman across parcels p = 1.3e-06; Spearman with AP rank rho = -0.758, p = 2.8e-05; mean R2 change from signed speed 0.0004
ca single R2 (mean over 23 left parcels, negatives set to 0; wheel is |speed|):
face_pc1 pupil wheel
grabAM05 0.0780 0.0008 0.0483
grabAM06 0.1419 0.0258 0.1474
grabAM07 0.0895 0.0000 0.0687
grabAM08 0.1386 0.0140 0.1255
mean over mice: {'face_pc1': 0.112, 'pupil': 0.0102, 'wheel': 0.0974}; best variable per mouse: {'grabAM05': 'face_pc1', 'grabAM06': 'wheel', 'grabAM07': 'face_pc1', 'grabAM08': 'face_pc1'}; Friedman chi2 = 6.50, p = 0.0388 (exact p = 0.0417)
face_pc1 above pupil in 4 of 4 mice (margin per mouse {'grabAM05': 0.077, 'grabAM06': 0.116, 'grabAM07': 0.089, 'grabAM08': 0.125}); signed-rank p = 0.125
face_pc1 above wheel in 3 of 4 mice (margin per mouse {'grabAM05': 0.03, 'grabAM06': -0.005, 'grabAM07': 0.021, 'grabAM08': 0.013}); signed-rank p = 0.250
ca delta R2 (mean over 23 left parcels; wheel is |speed|):
face_pc1 pupil wheel
grabAM05 0.0285 0.0335 0.0210
grabAM06 0.0511 0.0084 0.0268
grabAM07 0.0692 0.0236 0.0161
grabAM08 0.0723 0.0445 0.0080
mean over mice: {'face_pc1': 0.0553, 'pupil': 0.0275, 'wheel': 0.018}; best variable per mouse: {'grabAM05': 'pupil', 'grabAM06': 'face_pc1', 'grabAM07': 'face_pc1', 'grabAM08': 'face_pc1'}; Friedman chi2 = 4.50, p = 0.105 (exact p = 0.1250)
face_pc1 above pupil in 3 of 4 mice (margin per mouse {'grabAM05': -0.005, 'grabAM06': 0.043, 'grabAM07': 0.046, 'grabAM08': 0.028}); signed-rank p = 0.250
face_pc1 above wheel in 4 of 4 mice (margin per mouse {'grabAM05': 0.008, 'grabAM06': 0.024, 'grabAM07': 0.053, 'grabAM08': 0.064}); signed-rank p = 0.125
ca full model with |speed|: Friedman across parcels p = 7.7e-06; Spearman with AP rank rho = 0.558, p = 0.0056; mean R2 change from signed speed 0.0040
Result (variant, does not change an outcome): With |speed| as the wheel predictor, face PC1 stays the best single variable and has the largest unique contribution in all 4 mice for ACh. For calcium, wheel speed is higher than face PC1 in grabAM06 (0.147 against 0.142), so face is best in 3 of 4 mice (Friedman p = 0.039), and its unique contribution is the largest in 3 of 4 (in grabAM05 pupil is higher, 0.0335 against 0.0285; Friedman p = 0.105); under this reading the calcium claim would be partially reproduced. Full-model R2 changes by 0.0004 (ACh) and 0.0040 (calcium) on average, and the fig2a and fig2b results keep their signs and significance (Friedman p = 1.3e-06 and 7.7e-06; rho = -0.758 and 0.558).
Comparison with the paper's own per-mouse values. The code repository holds
Figure2Regression/Fig2RawData.mat, which MakeFig2.m plots: cross-validated R2 for 6 mice and
every parcel, for the full model, each single variable and the model without each variable. The
file does not name the mice of its rows, so the rows are matched to the mice here by the full-model
R2 across the 23 left parcels. The paper's unique contribution is the mean full-model R2 minus the
mean R2 without the variable (negatives set to 0 before averaging, as MakeFig2.m does); here it is
the drop when the variable is shuffled, averaged over parcels.
paper = meso.paper_r2()
rows = [f"row{i + 1}" for i in range(6)]
P = {c: {k: pd.DataFrame(v.T, index=meso.AP_ORDER, columns=rows).loc[left_areas] for k, v in paper[c].items()}
for c in paper}
print("Spearman correlation across the 23 left parcels of full-model R2, this attempt (rows) against the file's rows:")
for c in ("ach", "ca"):
corr = pd.DataFrame({r: [stats.spearmanr(full[c][m], P[c]["full"][r])[0] for m in MICE] for r in rows}, index=MICE)
print(f"{c}:\n{corr.round(2).to_string()}")
print(f" best row per mouse {corr.idxmax(axis=1).to_dict()}; mean full R2 here {full[c].mean().round(3).to_dict()}; "
f"in the file {P[c]['full'].mean().round(3).to_dict()}")
def paper_variable_table(c, kind):
pos = {k: v.clip(lower=0).mean() for k, v in P[c].items()}
if kind == "single":
return pd.DataFrame({"face_pc1": pos["face"], "pupil": pos["pupil"], "wheel": pos["wheel"]})
return pd.DataFrame({"face_pc1": pos["full"] - pos["noFace"], "pupil": pos["full"] - pos["noPupil"],
"wheel": pos["full"] - pos["noWheel"]})
for c in ("ach", "ca"):
for kind in ("single", "delta"):
tab = paper_variable_table(c, kind)
fr = stats.friedmanchisquare(*tab.values.T)
print(f"\nthe file's {c} {kind} R2 (mean over 23 left parcels), 6 mice:\n{tab.round(4).to_string()}")
print(f" best variable per row {tab.idxmax(axis=1).to_dict()}; face_pc1 best in {(tab.idxmax(axis=1) == 'face_pc1').sum()} of 6; "
f"Friedman p = {fr.pvalue:.3g}")
print(f"\nwheel single-variable R2 per mouse here, {c}: signed {variable_table(c, 'single').wheel.round(4).to_dict()}, "
f"|speed| {variable_table(c, 'single', abs_wheel=True).wheel.round(4).to_dict()}")
Spearman correlation across the 23 left parcels of full-model R2, this attempt (rows) against the file's rows:
ach:
row1 row2 row3 row4 row5 row6
grabAM05 0.98 0.88 0.90 0.47 0.33 0.67
grabAM06 0.63 0.89 0.64 0.87 0.69 0.93
grabAM07 0.85 0.82 0.98 0.49 0.21 0.73
grabAM08 0.42 0.72 0.46 0.95 0.75 0.91
best row per mouse {'grabAM05': 'row1', 'grabAM06': 'row6', 'grabAM07': 'row3', 'grabAM08': 'row4'}; mean full R2 here {'grabAM05': 0.182, 'grabAM06': 0.24, 'grabAM07': 0.126, 'grabAM08': 0.14}; in the file {'row1': 0.166, 'row2': 0.278, 'row3': 0.126, 'row4': 0.142, 'row5': 0.066, 'row6': 0.113}
ca:
row1 row2 row3 row4 row5 row6
grabAM05 0.94 0.83 0.52 0.43 0.45 0.60
grabAM06 0.52 0.72 0.72 0.79 0.68 0.38
grabAM07 0.68 0.57 0.98 0.85 0.47 0.72
grabAM08 0.74 0.75 0.87 0.95 0.62 0.55
best row per mouse {'grabAM05': 'row1', 'grabAM06': 'row4', 'grabAM07': 'row3', 'grabAM08': 'row4'}; mean full R2 here {'grabAM05': 0.207, 'grabAM06': 0.305, 'grabAM07': 0.216, 'grabAM08': 0.227}; in the file {'row1': 0.193, 'row2': 0.284, 'row3': 0.229, 'row4': 0.254, 'row5': 0.033, 'row6': 0.127}
the file's ach single R2 (mean over 23 left parcels), 6 mice:
face_pc1 pupil wheel
row1 0.0690 0.0000 0.0014
row2 0.1034 0.0003 0.0041
row3 0.0460 0.0000 0.0129
row4 0.1260 0.0097 0.0581
row5 0.1216 0.0003 0.0000
row6 0.0350 0.0000 0.0009
best variable per row {'row1': 'face_pc1', 'row2': 'face_pc1', 'row3': 'face_pc1', 'row4': 'face_pc1', 'row5': 'face_pc1', 'row6': 'face_pc1'}; face_pc1 best in 6 of 6; Friedman p = 0.0057
the file's ach delta R2 (mean over 23 left parcels), 6 mice:
face_pc1 pupil wheel
row1 0.0634 0.0094 0.0029
row2 0.1875 0.0394 0.0078
row3 0.0769 0.0169 0.0054
row4 0.0699 0.0051 0.0038
row5 0.0418 -0.0028 -0.0092
row6 0.0853 0.0305 0.0036
best variable per row {'row1': 'face_pc1', 'row2': 'face_pc1', 'row3': 'face_pc1', 'row4': 'face_pc1', 'row5': 'face_pc1', 'row6': 'face_pc1'}; face_pc1 best in 6 of 6; Friedman p = 0.00248
wheel single-variable R2 per mouse here, ach: signed {'grabAM05': 0.0017, 'grabAM06': 0.0296, 'grabAM07': 0.0115, 'grabAM08': 0.0601}, |speed| {'grabAM05': 0.0031, 'grabAM06': 0.0349, 'grabAM07': 0.0119, 'grabAM08': 0.0628}
the file's ca single R2 (mean over 23 left parcels), 6 mice:
face_pc1 pupil wheel
row1 0.0758 0.0008 0.0246
row2 0.1253 0.0264 0.0770
row3 0.0783 0.0000 0.0781
row4 0.1794 0.0109 0.1306
row5 0.0750 0.0239 0.0124
row6 0.0697 0.0188 0.0347
best variable per row {'row1': 'face_pc1', 'row2': 'face_pc1', 'row3': 'face_pc1', 'row4': 'face_pc1', 'row5': 'face_pc1', 'row6': 'face_pc1'}; face_pc1 best in 6 of 6; Friedman p = 0.0057
the file's ca delta R2 (mean over 23 left parcels), 6 mice:
face_pc1 pupil wheel
row1 0.0230 0.0260 0.0136
row2 0.0842 0.0055 0.0241
row3 0.0644 0.0236 0.0256
row4 0.0912 0.0426 0.0051
row5 0.0038 -0.0002 -0.0025
row6 0.0487 0.0147 -0.0004
best variable per row {'row1': 'pupil', 'row2': 'face_pc1', 'row3': 'face_pc1', 'row4': 'face_pc1', 'row5': 'face_pc1', 'row6': 'face_pc1'}; face_pc1 best in 5 of 6; Friedman p = 0.0302
wheel single-variable R2 per mouse here, ca: signed {'grabAM05': 0.0253, 'grabAM06': 0.1263, 'grabAM07': 0.0638, 'grabAM08': 0.1211}, |speed| {'grabAM05': 0.0483, 'grabAM06': 0.1474, 'grabAM07': 0.0687, 'grabAM08': 0.1255}
Result (comparison, decides no outcome): Three of the file's rows match a mouse here closely in full-model R2 across parcels, for both signals: row 1 grabAM05 (Spearman 0.98 for ACh, 0.94 for calcium), row 3 grabAM07 (0.98 and 0.98) and row 4 grabAM08 (0.95 and 0.95). grabAM06 has no clear match: its best rows are 6 for ACh (0.93) and 4 for calcium (0.79, which grabAM08 matches better), and row 2 gives 0.89 and 0.72. So rows 1, 3 and 4 are very likely grabAM05, grabAM07 and grabAM08, which fits rows 1 to 4 being grabAM05 to grabAM08, but grabAM06's row is not established. For those three mice the mean full-model R2 here is close to the file's (ACh 0.182, 0.126 and 0.140 against 0.166, 0.126 and 0.142; calcium 0.207, 0.216 and 0.227 against 0.193, 0.229 and 0.254), a first outside check of the ridge regression used here. On the wheel predictor grabAM05 is the mouse where the two readings differ most: its calcium wheel R2 is 0.0253 with signed speed and 0.0483 with |speed|, against 0.0246 in row 1, and its ACh wheel R2 0.0017 and 0.0031 against 0.0014. In grabAM07 and grabAM08 the two readings are closer to each other than to the file (calcium 0.0638 and 0.0687 against 0.0781; 0.1211 and 0.1255 against 0.1306), so they do not choose. In the file face PC1 is the best single variable in all 6 mice for both signals, and its calcium unique contribution is the largest in 5 of 6 (in row 1 pupil is higher, 0.0260 against 0.0230, as in grabAM05 here with |speed|; with signed speed face is higher there by 0.002).
Result: reproduced. The first revision, which used |speed|, called this partially reproduced; the third verification round agreed with the change to reproduced. With the signed wheel speed, which the Methods and the code define, face PC1 has the highest single-variable R2 in all 4 mice (mean 0.112, against 0.084 for wheel speed and 0.010 for pupil), Friedman p = 0.018 (exact p = 0.0046); its margin over wheel speed is 0.016 to 0.053. Its unique contribution is also the largest in all 4 mice (Friedman p = 0.039, exact p = 0.042), in grabAM05 by 0.002 over pupil. With |speed| face is best in 3 of 4 mice (variant above). The paper's own values side with the signed reading in grabAM05, the one mouse where the two readings differ clearly.
Transients at face movement and locomotion onsets (Figure S7)¶
Face movement onset: ACh 0.64 ± 0.09, calcium 0.69 ± 0.15 % ΔF/F. Locomotion onset: ACh 0.68 ± 0.17, calcium 1.50 ± 0.11 % ΔF/F. ACh transients do not differ between onset types; calcium transients are largest at locomotion onset. ACh, but not calcium, transients have an anterior-posterior gradient.
WINDOWS = {"face": (4, 4, (-4, -2)), "loc": (5, 5, (-5, -3))} # pre, post, baseline (s)
def peaks_from(kind, events):
"""events[(mouse, session)]: event times. Per colour and mouse, the peak in 0 to 1 s of the
baseline-subtracted trial mean (percent ΔF/F, 46 parcels), and the trials used per mouse."""
pre, post, base = WINDOWS[kind]
out, counts = {}, {}
for colour in ("ach", "ca"):
out[colour] = {}
for mouse in MICE:
trials = []
for s in sessions:
if s["mouse"] != mouse:
continue
tr = meso.event_traces(s[colour], s["t"], events[(s["mouse"], s["session"])], pre, post)
if len(tr):
trials.append(tr)
trials = np.concatenate(trials)
i0, i1 = int((pre + base[0]) * 10), int((pre + base[1]) * 10)
norm = 100 * (trials - trials[:, i0:i1].mean(1, keepdims=True))
mean = norm.mean(0)
out[colour][mouse] = mean[pre * 10: pre * 10 + 10].max(0)[:-1]
counts[mouse] = len(trials)
return out, counts
def transient_peaks(kind):
key = "face_onsets" if kind == "face" else "loc_onsets"
return peaks_from(kind, {(s["mouse"], s["session"]): s[key] for s in sessions})
peaks = {}
for kind in ("face", "loc"):
peaks[kind], counts = transient_peaks(kind)
print(kind, "onsets used per mouse:", counts)
# Null: the same peak at random times in sustained quiescence (the code's quiescent epochs, 10 s
# from their ends), with the same number of events per session as onsets, 200 draws.
rng_null = np.random.default_rng(7)
N_NULL = 200
def random_quiescent_times(s, n, rng):
ep = s["sit"]
if n == 0 or len(ep) == 0:
return np.array([])
length = ep[:, 1] - ep[:, 0]
k = rng.choice(len(ep), size=n, p=length / length.sum())
return ep[k, 0] + rng.random(n) * length[k]
null_wc = {}
for kind in ("face", "loc"):
key = "face_onsets" if kind == "face" else "loc_onsets"
draws = {"ach": [], "ca": []}
for _ in range(N_NULL):
ev = {(s["mouse"], s["session"]): random_quiescent_times(s, len(s[key]), rng_null) for s in sessions}
pk, _c = peaks_from(kind, ev)
for colour in draws:
draws[colour].append({m: np.nanmean(pk[colour][m][LEFT]) for m in MICE})
null_wc[kind] = {colour: pd.DataFrame(d) for colour, d in draws.items()}
def whole_cortex(kind, colour):
return pd.Series({m: np.nanmean(peaks[kind][colour][m][LEFT]) for m in MICE})
def report_transient(kind, expected):
for colour in ("ach", "ca"):
v = whole_cortex(kind, colour)
tt = stats.ttest_1samp(v, 0)
exp, tol = expected[colour]
print(f"{kind} onset, {colour}: per mouse {v.round(3).to_dict()}; mean {v.mean():.3f} ± {v.sem():.3f} (SEM) % ΔF/F; "
f"paper {exp} ± {tol} ({exp - tol:.2f} to {exp + tol:.2f}), within tolerance: {abs(v.mean() - exp) <= tol + 1e-12}; positive in {(v > 0).sum()} of {len(v)} mice; one-sample t = {tt.statistic:.2f}, p = {tt.pvalue:.3g}")
nm = null_wc[kind][colour]
tn = stats.ttest_rel(v, nm.mean())
print(f" null peak at random quiescent times ({N_NULL} draws): per mouse mean {nm.mean().round(3).to_dict()}, 95th percentile "
f"{nm.quantile(0.95).round(3).to_dict()}; observed above the 95th percentile in {(v > nm.quantile(0.95)).sum()} of {len(v)} mice; "
f"observed minus null mean {(v - nm.mean()).mean():.3f} % ΔF/F, paired t = {tn.statistic:.2f}, p = {tn.pvalue:.3g}")
print("parcel peak range over mice (percent):", {c: (round(np.nanmin([peaks[kind][c][m][LEFT] for m in MICE]), 2),
round(np.nanmax([peaks[kind][c][m][LEFT] for m in MICE]), 2)) for c in ("ach", "ca")})
face onsets used per mouse: {'grabAM05': 31, 'grabAM06': 25, 'grabAM07': 12, 'grabAM08': 30}
loc onsets used per mouse: {'grabAM05': 17, 'grabAM06': 39, 'grabAM07': 38, 'grabAM08': 25}
report_transient("face", {"ach": (0.64, 0.09), "ca": (0.69, 0.15)})
# Variant: face onsets with the quiescent period after the last locomotion epoch kept (the code's
# extractFacePupilTStamps.m leaves it out; getSustainedStates.m keeps it).
alt, alt_counts = peaks_from("face", {(s["mouse"], s["session"]): s["face_onsets_all"] for s in sessions})
print("variant, last quiescent period kept: onsets used per mouse", alt_counts)
for colour in ("ach", "ca"):
v = pd.Series({m: np.nanmean(alt[colour][m][LEFT]) for m in MICE})
print(f" {colour}: per mouse {v.round(3).to_dict()}; mean {v.mean():.3f} ± {v.sem():.3f} % ΔF/F; t-test p = {stats.ttest_1samp(v, 0).pvalue:.3g}")
face onset, ach: per mouse {'grabAM05': 0.506, 'grabAM06': 0.939, 'grabAM07': 0.536, 'grabAM08': 0.929}; mean 0.728 ± 0.119 (SEM) % ΔF/F; paper 0.64 ± 0.09 (0.55 to 0.73), within tolerance: True; positive in 4 of 4 mice; one-sample t = 6.09, p = 0.00888
null peak at random quiescent times (200 draws): per mouse mean {'grabAM05': 0.049, 'grabAM06': 0.088, 'grabAM07': 0.076, 'grabAM08': 0.07}, 95th percentile {'grabAM05': 0.111, 'grabAM06': 0.209, 'grabAM07': 0.178, 'grabAM08': 0.227}; observed above the 95th percentile in 4 of 4 mice; observed minus null mean 0.657 % ΔF/F, paired t = 5.74, p = 0.0105
face onset, ca: per mouse {'grabAM05': 0.888, 'grabAM06': 1.029, 'grabAM07': 0.572, 'grabAM08': 1.19}; mean 0.920 ± 0.131 (SEM) % ΔF/F; paper 0.69 ± 0.15 (0.54 to 0.84), within tolerance: False; positive in 4 of 4 mice; one-sample t = 7.00, p = 0.00598
null peak at random quiescent times (200 draws): per mouse mean {'grabAM05': 0.113, 'grabAM06': 0.157, 'grabAM07': 0.168, 'grabAM08': 0.177}, 95th percentile {'grabAM05': 0.261, 'grabAM06': 0.316, 'grabAM07': 0.34, 'grabAM08': 0.377}; observed above the 95th percentile in 4 of 4 mice; observed minus null mean 0.766 % ΔF/F, paired t = 5.88, p = 0.00983
parcel peak range over mice (percent): {'ach': (np.float64(0.1), np.float64(1.8)), 'ca': (np.float64(0.14), np.float64(2.03))}
variant, last quiescent period kept: onsets used per mouse {'grabAM05': 47, 'grabAM06': 25, 'grabAM07': 12, 'grabAM08': 29}
ach: per mouse {'grabAM05': 0.486, 'grabAM06': 0.939, 'grabAM07': 0.536, 'grabAM08': 0.952}; mean 0.728 ± 0.126 % ΔF/F; t-test p = 0.0103
ca: per mouse {'grabAM05': 0.815, 'grabAM06': 1.029, 'grabAM07': 0.572, 'grabAM08': 1.067}; mean 0.871 ± 0.114 % ΔF/F; t-test p = 0.00467
Result: partially reproduced. Both transients are positive in all 4 mice (ACh p = 0.0089, calcium p = 0.0060). The ACh mean, 0.728 % ΔF/F, is inside the paper's 0.55 to 0.73, by 0.002; the calcium mean, 0.920, is above the paper's 0.54 to 0.84. The same peak taken at random times in sustained quiescence is 0.049 to 0.088 (ACh) and 0.113 to 0.177 (calcium) per mouse, and every mouse's transient is above the 95th percentile of its null. Keeping the last quiescent period gives ACh 0.728 and calcium 0.871, still outside tolerance.
report_transient("loc", {"ach": (0.68, 0.17), "ca": (1.50, 0.11)})
loc onset, ach: per mouse {'grabAM05': 0.345, 'grabAM06': 0.654, 'grabAM07': 0.491, 'grabAM08': 1.555}; mean 0.761 ± 0.272 (SEM) % ΔF/F; paper 0.68 ± 0.17 (0.51 to 0.85), within tolerance: True; positive in 4 of 4 mice; one-sample t = 2.80, p = 0.0679
null peak at random quiescent times (200 draws): per mouse mean {'grabAM05': 0.062, 'grabAM06': 0.062, 'grabAM07': 0.044, 'grabAM08': 0.082}, 95th percentile {'grabAM05': 0.16, 'grabAM06': 0.135, 'grabAM07': 0.088, 'grabAM08': 0.243}; observed above the 95th percentile in 4 of 4 mice; observed minus null mean 0.699 % ΔF/F, paired t = 2.63, p = 0.0783
loc onset, ca: per mouse {'grabAM05': 1.336, 'grabAM06': 1.696, 'grabAM07': 1.472, 'grabAM08': 3.436}; mean 1.985 ± 0.489 (SEM) % ΔF/F; paper 1.5 ± 0.11 (1.39 to 1.61), within tolerance: False; positive in 4 of 4 mice; one-sample t = 4.06, p = 0.027
null peak at random quiescent times (200 draws): per mouse mean {'grabAM05': 0.141, 'grabAM06': 0.163, 'grabAM07': 0.111, 'grabAM08': 0.212}, 95th percentile {'grabAM05': 0.357, 'grabAM06': 0.27, 'grabAM07': 0.191, 'grabAM08': 0.413}; observed above the 95th percentile in 4 of 4 mice; observed minus null mean 1.828 % ΔF/F, paired t = 3.89, p = 0.0302
parcel peak range over mice (percent): {'ach': (np.float64(0.09), np.float64(3.03)), 'ca': (np.float64(0.7), np.float64(5.38))}
Result: partially reproduced. Both transients are positive in all 4 mice. The ACh mean, 0.761 % ΔF/F, is inside the paper's 0.51 to 0.85 but its t-test gives p = 0.068, driven by one mouse (grabAM08, 1.555); the calcium mean, 1.985, is above the paper's 1.39 to 1.61 (p = 0.027). Against the null peak at random quiescent times (per mouse means 0.044 to 0.082 for ACh, 0.111 to 0.212 for calcium), every mouse's transient is above the null's 95th percentile.
def rm_anova(values, factor_name, levels):
"""values[level][mouse] = array over the 23 left parcels. Two-way repeated measures ANOVA."""
rows = []
for lev in levels:
for m in MICE:
v = values[lev].get(m)
if v is None:
continue
for a, x in zip(left_areas, v):
rows.append({"mouse": m, "area": a, factor_name: lev, "y": x})
df = pd.DataFrame(rows)
complete = df.groupby("mouse").y.apply(lambda y: y.notna().all() and len(y) == 23 * len(levels))
df = df[df.mouse.isin(complete.index[complete])]
return AnovaRM(df, "y", "mouse", within=["area", factor_name]).fit().anova_table, df.mouse.nunique()
def transient_anova(colour):
vals = {kind: {m: peaks[kind][colour][m][LEFT] for m in MICE} for kind in ("face", "loc")}
tab, n = rm_anova(vals, "onset", ["face", "loc"])
print(f"{colour}: two-way repeated measures ANOVA, n = {n} mice")
print(tab.round(4).to_string())
diff = whole_cortex("loc", colour) - whole_cortex("face", colour)
print(f"whole-cortex locomotion minus face per mouse {diff.round(3).to_dict()}, mean {diff.mean():.3f} % ΔF/F")
return tab
tab_ach_tr = transient_anova("ach")
ach: two-way repeated measures ANOVA, n = 4 mice
F Value Num DF Den DF Pr > F
area 9.3928 22.0 66.0 0.0000
onset 0.0276 1.0 3.0 0.8787
area:onset 0.8544 22.0 66.0 0.6498
whole-cortex locomotion minus face per mouse {'grabAM05': -0.161, 'grabAM06': -0.285, 'grabAM07': -0.045, 'grabAM08': 0.625}, mean 0.034 % ΔF/F
Result: reproduced. The onset-type effect on ACh transients is not significant (F = 0.03, p = 0.88; mean locomotion minus face 0.034 % ΔF/F, per mouse -0.285 to 0.625, three mice one way and one the other), while the area effect is (p < 0.0001). This is an absence result with 4 mice and a test with 3 denominator degrees of freedom, so it carries little weight.
tab_ca_tr = transient_anova("ca")
ca: two-way repeated measures ANOVA, n = 4 mice
F Value Num DF Den DF Pr > F
area 31.5903 22.0 66.0 0.0000
onset 6.9360 1.0 3.0 0.0781
area:onset 2.7762 22.0 66.0 0.0007
whole-cortex locomotion minus face per mouse {'grabAM05': 0.447, 'grabAM06': 0.667, 'grabAM07': 0.9, 'grabAM08': 2.245}, mean 1.065 % ΔF/F
Result: partially reproduced. The calcium transient is larger at locomotion onset than at face onset in all 4 mice (by 0.447 to 2.245 % ΔF/F), but the onset-type effect has p = 0.078; the area effect has p < 0.0001.
for colour in ("ach", "ca"):
for kind in ("face", "loc"):
mean_peak = np.nanmean([peaks[kind][colour][m][LEFT] for m in MICE], axis=0)
rho, p = stats.spearmanr(ap_left, mean_peak)
print(f"{colour} {kind} onset: Spearman of mean peak with anterior-posterior rank (1 = MOs): rho = {rho:.3f}, p = {p:.3g}")
ach face onset: Spearman of mean peak with anterior-posterior rank (1 = MOs): rho = -0.605, p = 0.00224 ach loc onset: Spearman of mean peak with anterior-posterior rank (1 = MOs): rho = -0.752, p = 3.51e-05 ca face onset: Spearman of mean peak with anterior-posterior rank (1 = MOs): rho = 0.014, p = 0.95 ca loc onset: Spearman of mean peak with anterior-posterior rank (1 = MOs): rho = 0.106, p = 0.631
Result: reproduced. ACh transients are larger in anterior areas at face onset (rho = -0.605, p = 0.0022) and locomotion onset (rho = -0.752, p = 3.5e-05); calcium transients have no gradient (p = 0.95 and 0.63).
Sustained states (Figure 3c, 3d)¶
ACh: low face -0.13 ± 0.04, high face 0.18 ± 0.08, locomotion 0.16 ± 0.06 % ΔF/F. Calcium: low face -0.02 ± 0.02, high face 0.17 ± 0.04, locomotion 0.63 ± 0.19 % ΔF/F. ACh is higher in frontal-motor than visual areas during high face movement.
STATES = ("face_low", "face_high", "loc")
amp = {c: {st: {} for st in STATES} for c in ("ach", "ca")}
for mouse in MICE:
for c in ("ach", "ca"):
sess_means = {st: [] for st in STATES}
for s in sessions:
m = matched[(s["mouse"], s["session"])]
if s["mouse"] != mouse or m is None:
continue
for st in STATES:
ep = [s[c][(s["t"] > a) & (s["t"] < b)].mean(0) for a, b in m[st]]
sess_means[st].append(np.nanmean(ep, axis=0))
for st in STATES:
if sess_means[st]:
amp[c][st][mouse] = 100 * np.nanmean(sess_means[st], axis=0)[LEFT]
def lsd(c, a, b, mice=None):
"""Fisher's LSD on whole-cortex state means with the state-by-mouse error of the RM ANOVA."""
mice = mice or [m for m in MICE if all(m in amp[c][st] for st in STATES)]
Y = np.array([[np.nanmean(amp[c][st][m]) for st in STATES] for m in mice])
n, k = Y.shape
resid = Y - Y.mean(1, keepdims=True) - Y.mean(0, keepdims=True) + Y.mean()
mse = (resid ** 2).sum() / ((n - 1) * (k - 1))
ia, ib = STATES.index(a), STATES.index(b)
tval = (Y[:, ia].mean() - Y[:, ib].mean()) / np.sqrt(2 * mse / n)
return tval, 2 * stats.t.sf(abs(tval), (n - 1) * (k - 1)), Y[:, ia] - Y[:, ib]
def sustained_report(c, expected):
tab, n = rm_anova(amp[c], "state", list(STATES))
print(f"{c}: two-way repeated measures ANOVA (area by state), n = {n} mice")
print(tab.round(4).to_string())
for st in STATES:
v = pd.Series({m: np.nanmean(x) for m, x in amp[c][st].items()})
exp, tol = expected[st]
print(f" {st}: per mouse {v.round(3).to_dict()}; mean {v.mean():.3f} ± {v.sem():.3f} % ΔF/F; paper {exp} ± {tol} ({exp - tol:.2f} to {exp + tol:.2f}); "
f"within tolerance: {abs(v.mean() - exp) <= tol + 1e-12}")
for a, b in (("face_high", "face_low"), ("loc", "face_high")):
tval, p, d = lsd(c, a, b)
print(f" LSD {a} - {b}: per mouse {np.round(d, 3)}, mean {d.mean():.3f}, t = {tval:.2f}, p = {p:.3f}")
sustained_report("ach", {"face_low": (-0.13, 0.04), "face_high": (0.18, 0.08), "loc": (0.16, 0.06)})
ach: two-way repeated measures ANOVA (area by state), n = 3 mice
F Value Num DF Den DF Pr > F
area 8.1731 22.0 44.0 0.0000
state 11.0809 2.0 4.0 0.0234
area:state 8.4597 44.0 88.0 0.0000
face_low: per mouse {'grabAM05': -0.069, 'grabAM06': -0.167, 'grabAM08': -0.152}; mean -0.129 ± 0.031 % ΔF/F; paper -0.13 ± 0.04 (-0.17 to -0.09); within tolerance: True
face_high: per mouse {'grabAM05': 0.07, 'grabAM06': 0.173, 'grabAM08': 0.231}; mean 0.158 ± 0.047 % ΔF/F; paper 0.18 ± 0.08 (0.10 to 0.26); within tolerance: True
loc: per mouse {'grabAM05': 0.097, 'grabAM06': 0.242, 'grabAM08': 0.411}; mean 0.250 ± 0.091 % ΔF/F; paper 0.16 ± 0.06 (0.10 to 0.22); within tolerance: False
LSD face_high - face_low: per mouse [0.138 0.34 0.384], mean 0.287, t = 3.42, p = 0.027
LSD loc - face_high: per mouse [0.027 0.069 0.179], mean 0.092, t = 1.09, p = 0.336
Result: reproduced. With the 3 mice that have matched epochs, whole-cortex ACh is -0.129 % ΔF/F in low and 0.158 in high face states, both inside the paper's ranges; high exceeds low in all 3 mice (LSD p = 0.027) and the state effect has p = 0.023. Without grabAM06 (2 epochs per state) high still exceeds low in both other mice, but the LSD test has p = 0.207 (robustness section).
tval, p, d = lsd("ach", "loc", "face_high")
v = pd.Series({m: np.nanmean(x) for m, x in amp["ach"]["loc"].items()})
print(f"ACh locomotion mean {v.mean():.3f} % ΔF/F (paper 0.16 ± 0.06); locomotion minus high face per mouse {np.round(d, 3)}, "
f"mean {d.mean():.3f}; LSD t = {tval:.2f}, p = {p:.3f}")
ACh locomotion mean 0.250 % ΔF/F (paper 0.16 ± 0.06); locomotion minus high face per mouse [0.027 0.069 0.179], mean 0.092; LSD t = 1.09, p = 0.336
Result: reproduced. Locomotion does not differ significantly from high face (LSD p = 0.34), as the criterion asks. The data point the other way from the claim: the locomotion mean, 0.250 % ΔF/F, is above the paper's 0.10 to 0.22, and it is higher than high face in all 3 mice (by 0.027 to 0.179). With 3 mice the test has little power, so this absence result carries little weight.
sustained_report("ca", {"face_low": (-0.02, 0.02), "face_high": (0.17, 0.04), "loc": (0.63, 0.19)})
ca: two-way repeated measures ANOVA (area by state), n = 3 mice
F Value Num DF Den DF Pr > F
area 15.9902 22.0 44.0 0.0000
state 48.3029 2.0 4.0 0.0016
area:state 14.9958 44.0 88.0 0.0000
face_low: per mouse {'grabAM05': -0.052, 'grabAM06': 0.013, 'grabAM08': 0.032}; mean -0.002 ± 0.025 % ΔF/F; paper -0.02 ± 0.02 (-0.04 to 0.00); within tolerance: True
face_high: per mouse {'grabAM05': 0.17, 'grabAM06': 0.166, 'grabAM08': 0.156}; mean 0.164 ± 0.004 % ΔF/F; paper 0.17 ± 0.04 (0.13 to 0.21); within tolerance: True
loc: per mouse {'grabAM05': 0.827, 'grabAM06': 1.365, 'grabAM08': 1.126}; mean 1.106 ± 0.156 % ΔF/F; paper 0.63 ± 0.19 (0.44 to 0.82); within tolerance: False
LSD face_high - face_low: per mouse [0.222 0.152 0.124], mean 0.166, t = 1.37, p = 0.244
LSD loc - face_high: per mouse [0.656 1.2 0.969], mean 0.942, t = 7.75, p = 0.001
Result: partially reproduced. Low face (-0.002 % ΔF/F) and high face (0.164) are inside the paper's ranges; locomotion (1.106) is above its 0.44 to 0.82. Both orderings hold in all 3 mice; locomotion above high face has LSD p = 0.001, but high above low has p = 0.24.
for c, (a, b) in (("ach", ("face_high", "face_low")), ("ca", ("loc", "face_high"))):
mice = [m for m in MICE if m in amp[c][a] and m in amp[c][b]]
diff = np.mean([amp[c][a][m] - amp[c][b][m] for m in mice], axis=0)
rho, p = stats.spearmanr(ap_left, diff)
print(f"{c} {a} minus {b} (n = {len(mice)} mice): Spearman with anterior-posterior rank (1 = MOs) rho = {rho:.3f}, p = {p:.3g}; "
f"difference range {diff.min():.3f} to {diff.max():.3f} % ΔF/F")
ach face_high minus face_low (n = 3 mice): Spearman with anterior-posterior rank (1 = MOs) rho = -0.801, p = 4.31e-06; difference range 0.037 to 0.546 % ΔF/F ca loc minus face_high (n = 3 mice): Spearman with anterior-posterior rank (1 = MOs) rho = 0.329, p = 0.125; difference range 0.470 to 1.643 % ΔF/F
Result: partially reproduced. The ACh increase from low to high face is larger in anterior areas (rho = -0.801, p = 4.3e-06); the calcium increase from high face to locomotion has no significant gradient (rho = 0.329, p = 0.13).
fm = np.isin(left_areas, ["MOp", "MOs"])
vis = np.isin(left_areas, meso.GROUPS["visual"])
mice = list(amp["ach"]["face_high"])
fr_ = np.array([np.nanmean(amp["ach"]["face_high"][m][fm]) for m in mice])
vi_ = np.array([np.nanmean(amp["ach"]["face_high"][m][vis]) for m in mice])
tt = stats.ttest_rel(fr_, vi_)
print(f"high face ACh, mice {mice}: frontal-motor {np.round(fr_, 3)}, visual {np.round(vi_, 3)} % ΔF/F; "
f"difference {np.round(fr_ - vi_, 3)}; paired t = {tt.statistic:.2f}, p = {tt.pvalue:.3f}")
high face ACh, mice ['grabAM05', 'grabAM06', 'grabAM08']: frontal-motor [0.166 0.341 0.427], visual [0.032 0.092 0.108] % ΔF/F; difference [0.134 0.249 0.319]; paired t = 4.34, p = 0.049
Result: reproduced. In high face states ACh is higher in frontal-motor than visual areas in all 3 mice (by 0.134 to 0.319 % ΔF/F, paired p = 0.049). The p-value is just below 0.05 and rests on grabAM06: without it the difference holds in both other mice but p = 0.247 with 1 degree of freedom.
Correlations between areas and between signals (Figures 4 and 5)¶
Facial movement increases ACh correlations among anterior somatosensory and motor areas and calcium correlations broadly; locomotion decreases ACh correlations across the cortex and calcium correlations mainly in posterior areas. Facial movement increases ACh-calcium coupling broadly; locomotion decreases it in frontal cortex.
rng_perm = np.random.default_rng(5)
iu = np.tril_indices(23, -1)
z = {key: {st: {} for st in STATES} for key in ("ach", "ca", "bg")}
for mouse in MICE:
for st in STATES:
acc = {key: [] for key in z}
for s in sessions:
m = matched[(s["mouse"], s["session"])]
if s["mouse"] != mouse or m is None:
continue
for a, b in m[st]:
A = meso.epoch_data(s["ach"][:, LEFT], s["t"], a, b)
C = meso.epoch_data(s["ca"][:, LEFT], s["t"], a, b)
acc["ach"].append(np.arctanh(meso.max_crosscorr(A, A)[iu]))
acc["ca"].append(np.arctanh(meso.max_crosscorr(C, C)[iu]))
acc["bg"].append(np.arctanh(np.diag(meso.max_crosscorr(A, C))))
for key in z:
if acc[key]:
z[key][st][mouse] = np.array(acc[key])
print("epochs per mouse and state:", {m: len(z["ach"]["loc"].get(m, [])) for m in MICE})
pairs = pd.DataFrame({"a": left_areas[iu[0]], "b": left_areas[iu[1]], "rank_a": ap_left[iu[0]], "rank_b": ap_left[iu[1]]})
group_of = {a: g for g, areas in meso.GROUPS.items() for a in areas}
def compare(key, s1, s2, n_perm=10000):
mice = [m for m in MICE if m in z[key][s1] and m in z[key][s2]]
obs, p = meso.permutation_test([z[key][s1][m] for m in mice], [z[key][s2][m] for m in mice], n_perm, rng_perm)
sig = multipletests(p, alpha=0.01, method="fdr_by")[0]
up, down = sig & (obs > 0), sig & (obs < 0)
mean_r = {st: np.tanh(np.mean([z[key][st][m].mean(0) for m in mice], axis=0)) for st in (s1, s2)}
print(f"{key} {s1} - {s2}: mice {mice}; values {len(obs)}; significant (BY q = 0.01) increases {up.sum()}, decreases {down.sum()}; "
f"mean difference in z {obs.mean():.3f} (range {obs.min():.3f} to {obs.max():.3f}); uncorrected p < 0.05: {(p < 0.05).sum()}; "
f"smallest p {p.min():.4f}")
print(f" mean r over values: {s1} {mean_r[s1].mean():.3f}, {s2} {mean_r[s2].mean():.3f}")
per = {m: round(float((z[key][s1][m].mean(0) - z[key][s2][m].mean(0)).mean()), 3) for m in mice}
print(f" mean difference in z per mouse (epochs {[len(z[key][s1][m]) for m in mice]}): {per}")
sig05 = multipletests(p, alpha=0.05, method="fdr_by")[0]
print(f" variant, BY q = 0.05: increases {(sig05 & (obs > 0)).sum()}, decreases {(sig05 & (obs < 0)).sum()}")
return obs, p, up, down
def group_cover(mask):
hit = set(pairs.a[mask]) | set(pairs.b[mask])
return {g: any(a in hit for a in areas) for g, areas in meso.GROUPS.items()}
epochs per mouse and state: {'grabAM05': 10, 'grabAM06': 2, 'grabAM07': 0, 'grabAM08': 13}
obs, p, up, down = compare("ach", "face_high", "face_low")
ms = pairs.a.str.match(r"MO|SS") & pairs.b.str.match(r"MO|SS")
print(f"fraction with a significant increase: motor-somatosensory pairs {up[ms].mean():.3f} ({up[ms].sum()} of {ms.sum()}), "
f"other pairs {up[~ms].mean():.3f} ({up[~ms].sum()} of {(~ms).sum()})")
print(f"mean difference in z: motor-somatosensory {obs[ms].mean():.3f}, other {obs[~ms].mean():.3f}")
ach face_high - face_low: mice ['grabAM05', 'grabAM06', 'grabAM08']; values 253; significant (BY q = 0.01) increases 64, decreases 0; mean difference in z 0.257 (range -0.027 to 0.549); uncorrected p < 0.05: 159; smallest p 0.0001
mean r over values: face_high 0.791, face_low 0.691
mean difference in z per mouse (epochs [10, 2, 13]): {'grabAM05': 0.095, 'grabAM06': 0.314, 'grabAM08': 0.362}
variant, BY q = 0.05: increases 98, decreases 0
fraction with a significant increase: motor-somatosensory pairs 0.867 (39 of 45), other pairs 0.120 (25 of 208)
mean difference in z: motor-somatosensory 0.419, other 0.222
Result: reproduced. 64 of 253 pairs increase significantly and none decreases; 87% of motor-somatosensory pairs increase significantly against 12% of the other pairs. The mean change is positive in each of the 3 mice.
obs, p, up, down = compare("ach", "loc", "face_high")
print(f"fraction of pairs with a significant decrease {down.mean():.3f}; area groups in a decreased pair {group_cover(down)}")
ach loc - face_high: mice ['grabAM05', 'grabAM06', 'grabAM08']; values 253; significant (BY q = 0.01) increases 0, decreases 243; mean difference in z -0.531 (range -0.954 to -0.189); uncorrected p < 0.05: 253; smallest p 0.0001
mean r over values: loc 0.539, face_high 0.791
mean difference in z per mouse (epochs [10, 2, 13]): {'grabAM05': -0.337, 'grabAM06': -0.648, 'grabAM08': -0.608}
variant, BY q = 0.05: increases 0, decreases 250
fraction of pairs with a significant decrease 0.960; area groups in a decreased pair {'visual': True, 'retrosplenial': True, 'auditory': True, 'parietal': True, 'somatosensory': True, 'frontal-motor': True}
Result: reproduced. 243 of 253 pairs (96%) decrease significantly, none increases, and every area group is involved.
obs, p, up, down = compare("ca", "face_high", "face_low")
print(f"fraction of pairs with a significant increase {up.mean():.3f}; area groups in an increased pair {group_cover(up)}")
ca face_high - face_low: mice ['grabAM05', 'grabAM06', 'grabAM08']; values 253; significant (BY q = 0.01) increases 161, decreases 0; mean difference in z 0.397 (range 0.107 to 0.755); uncorrected p < 0.05: 224; smallest p 0.0001
mean r over values: face_high 0.843, face_low 0.672
mean difference in z per mouse (epochs [10, 2, 13]): {'grabAM05': 0.395, 'grabAM06': 0.567, 'grabAM08': 0.227}
variant, BY q = 0.05: increases 189, decreases 0
fraction of pairs with a significant increase 0.636; area groups in an increased pair {'visual': True, 'retrosplenial': True, 'auditory': True, 'parietal': True, 'somatosensory': True, 'frontal-motor': True}
Result: reproduced. 161 of 253 pairs (64%) increase significantly, none decreases, and every area group is involved.
obs, p, up, down = compare("ca", "loc", "face_high")
post = (pairs.rank_a >= 12) & (pairs.rank_b >= 12)
print(f"fraction with a significant decrease: posterior pairs {down[post].mean():.3f} ({down[post].sum()} of {post.sum()}), "
f"other pairs {down[~post].mean():.3f} ({down[~post].sum()} of {(~post).sum()})")
print(f"mean difference in z: posterior {obs[post].mean():.3f}, other {obs[~post].mean():.3f}")
ca loc - face_high: mice ['grabAM05', 'grabAM06', 'grabAM08']; values 253; significant (BY q = 0.01) increases 0, decreases 100; mean difference in z -0.302 (range -0.641 to 0.194); uncorrected p < 0.05: 149; smallest p 0.0001
mean r over values: loc 0.744, face_high 0.843
mean difference in z per mouse (epochs [10, 2, 13]): {'grabAM05': -0.183, 'grabAM06': -0.513, 'grabAM08': -0.21}
variant, BY q = 0.05: increases 0, decreases 117
fraction with a significant decrease: posterior pairs 0.833 (55 of 66), other pairs 0.241 (45 of 187)
mean difference in z: posterior -0.456, other -0.248
Result: reproduced. 100 pairs decrease significantly and none increases; 83% of posterior pairs decrease significantly against 24% of the other pairs.
obs, p, up, down = compare("bg", "face_high", "face_low")
print(pd.DataFrame({"area": left_areas, "diff_z": obs.round(3), "p": p.round(4), "up": up, "down": down}).to_string(index=False))
bg face_high - face_low: mice ['grabAM05', 'grabAM06', 'grabAM08']; values 23; significant (BY q = 0.01) increases 10, decreases 0; mean difference in z 0.471 (range 0.173 to 0.737); uncorrected p < 0.05: 22; smallest p 0.0001
mean r over values: face_high 0.573, face_low 0.185
mean difference in z per mouse (epochs [10, 2, 13]): {'grabAM05': 0.182, 'grabAM06': 0.694, 'grabAM08': 0.536}
variant, BY q = 0.05: increases 20, decreases 0
area diff_z p up down
VISp 0.484 0.0013 False False
VISpm 0.677 0.0001 True False
VISam 0.593 0.0001 True False
VISl 0.366 0.0015 False False
VISli 0.407 0.0127 False False
VISpl 0.255 0.0042 False False
VISpor 0.273 0.0175 False False
VISal 0.583 0.0004 True False
RSPagl 0.737 0.0001 True False
RSPd 0.632 0.0001 True False
TEa 0.445 0.0031 False False
AUD 0.362 0.0007 True False
PTLp 0.597 0.0001 True False
SSp-bfd 0.480 0.0088 False False
SSs 0.327 0.0043 False False
SSp-tr 0.511 0.0001 True False
SSp-n 0.173 0.2225 False False
SSp-un 0.340 0.0074 False False
SSp-ul 0.366 0.0035 False False
SSp-ll 0.539 0.0001 True False
SSp-m 0.328 0.0110 False False
MOp 0.697 0.0006 True False
MOs 0.655 0.0034 False False
Result: partially reproduced. ACh-calcium coupling increases in all 23 parcels and significantly in 10 (none decreases), short of the half the criterion asks for; 22 of 23 have uncorrected p < 0.05.
obs, p, up, down = compare("bg", "loc", "face_high")
frontal = np.isin(left_areas, ["MOs", "MOp"])
print(pd.DataFrame({"area": left_areas, "diff_z": obs.round(3), "p": p.round(4), "up": up, "down": down}).to_string(index=False))
print(f"significant decreases: frontal {down[frontal].sum()} of {frontal.sum()}, other {down[~frontal].sum()} of {(~frontal).sum()}")
bg loc - face_high: mice ['grabAM05', 'grabAM06', 'grabAM08']; values 23; significant (BY q = 0.01) increases 0, decreases 0; mean difference in z -0.188 (range -0.370 to -0.037); uncorrected p < 0.05: 8; smallest p 0.0006
mean r over values: loc 0.436, face_high 0.573
mean difference in z per mouse (epochs [10, 2, 13]): {'grabAM05': -0.138, 'grabAM06': -0.269, 'grabAM08': -0.158}
variant, BY q = 0.05: increases 0, decreases 2
area diff_z p up down
VISp -0.169 0.1043 False False
VISpm -0.247 0.0052 False False
VISam -0.233 0.0126 False False
VISl -0.132 0.1572 False False
VISli -0.169 0.0558 False False
VISpl -0.111 0.2001 False False
VISpor -0.194 0.0916 False False
VISal -0.262 0.0088 False False
RSPagl -0.200 0.0485 False False
RSPd -0.037 0.7998 False False
TEa -0.183 0.0377 False False
AUD -0.115 0.1833 False False
PTLp -0.370 0.0007 False False
SSp-bfd -0.084 0.4769 False False
SSs -0.085 0.4328 False False
SSp-tr -0.344 0.0006 False False
SSp-n -0.232 0.4000 False False
SSp-un -0.157 0.2607 False False
SSp-ul -0.127 0.2270 False False
SSp-ll -0.216 0.0083 False False
SSp-m -0.193 0.2051 False False
MOp -0.204 0.0761 False False
MOs -0.268 0.0720 False False
significant decreases: frontal 0 of 2, other 0 of 21
Result: not reproduced. No parcel decreases significantly after FDR correction. Coupling is lower with locomotion in all 23 parcels (by 0.037 to 0.370 in z) and in each mouse; the largest decreases are in PTLp and SSp-tr, not in MOs and MOp (uncorrected p 0.072 and 0.076).
Robustness of the matched-state results to grabAM06¶
grabAM06 contributes 2 matched epochs per state (data check) but has the same weight as the other two mice in the ANOVA, the LSD tests, the paired t-test and the permutation tests. The cells below repeat the tests of Figures 3c, 3d, 4 and 5 with grabAM05 and grabAM08 only. The outcomes are judged on the 3 mice, as the criteria fix; these are variants. With 2 mice the LSD tests have 2 degrees of freedom and the paired t-test 1.
rng_rob = np.random.default_rng(11)
two = [m for m in MICE if m in amp["ach"]["face_high"] and m != "grabAM06"]
print("mice:", two)
for c in ("ach", "ca"):
for st in STATES:
v = pd.Series({m: np.nanmean(amp[c][st][m]) for m in two})
print(f"{c} {st}: per mouse {v.round(3).to_dict()}, mean {v.mean():.3f} % ΔF/F")
for a, b in (("face_high", "face_low"), ("loc", "face_high")):
tval, p, d = lsd(c, a, b, mice=two)
print(f" {c} LSD {a} - {b}: per mouse {np.round(d, 3)}, t = {tval:.2f}, p = {p:.3f}")
fr2 = np.array([np.nanmean(amp["ach"]["face_high"][m][fm]) for m in two])
vi2 = np.array([np.nanmean(amp["ach"]["face_high"][m][vis]) for m in two])
tt2 = stats.ttest_rel(fr2, vi2)
print(f"fig3d, high face ACh frontal-motor minus visual: {np.round(fr2 - vi2, 3)}; paired t = {tt2.statistic:.2f}, p = {tt2.pvalue:.3f}")
for c, (a, b) in (("ach", ("face_high", "face_low")), ("ca", ("loc", "face_high"))):
rho, p = stats.spearmanr(ap_left, np.mean([amp[c][a][m] - amp[c][b][m] for m in two], axis=0))
print(f"s7 gradient, {c} {a} minus {b}: rho = {rho:.3f}, p = {p:.3g}")
ms = pairs.a.str.match(r"MO|SS") & pairs.b.str.match(r"MO|SS")
post = (pairs.rank_a >= 12) & (pairs.rank_b >= 12)
frontal = np.isin(left_areas, ["MOs", "MOp"])
for key, s1, s2 in (("ach", "face_high", "face_low"), ("ach", "loc", "face_high"), ("ca", "face_high", "face_low"),
("ca", "loc", "face_high"), ("bg", "face_high", "face_low"), ("bg", "loc", "face_high")):
obs, p = meso.permutation_test([z[key][s1][m] for m in two], [z[key][s2][m] for m in two], 10000, rng_rob)
sig = multipletests(p, alpha=0.01, method="fdr_by")[0]
up, down = sig & (obs > 0), sig & (obs < 0)
line = f"{key} {s1} - {s2}: significant increases {up.sum()}, decreases {down.sum()} of {len(obs)}"
if key == "ach" and s1 == "face_high":
line += f"; increase in motor-somatosensory pairs {up[ms].mean():.3f}, other pairs {up[~ms].mean():.3f}"
if key == "ca" and s1 == "loc":
line += f"; decrease in posterior pairs {down[post].mean():.3f}, other pairs {down[~post].mean():.3f}"
if key == "bg" and s1 == "loc":
line += f"; frontal decreases {down[frontal].sum()} of 2, MOs and MOp p {p[frontal].round(4)}"
print(line)
mice: ['grabAM05', 'grabAM08']
ach face_low: per mouse {'grabAM05': -0.069, 'grabAM08': -0.152}, mean -0.111 % ΔF/F
ach face_high: per mouse {'grabAM05': 0.07, 'grabAM08': 0.231}, mean 0.151 % ΔF/F
ach loc: per mouse {'grabAM05': 0.097, 'grabAM08': 0.411}, mean 0.254 % ΔF/F
ach LSD face_high - face_low: per mouse [0.138 0.384], t = 1.84, p = 0.207
ach LSD loc - face_high: per mouse [0.027 0.179], t = 0.73, p = 0.543
ca face_low: per mouse {'grabAM05': -0.052, 'grabAM08': 0.032}, mean -0.010 % ΔF/F
ca face_high: per mouse {'grabAM05': 0.17, 'grabAM08': 0.156}, mean 0.163 % ΔF/F
ca loc: per mouse {'grabAM05': 0.827, 'grabAM08': 1.126}, mean 0.976 % ΔF/F
ca LSD face_high - face_low: per mouse [0.222 0.124], t = 1.53, p = 0.266
ca LSD loc - face_high: per mouse [0.656 0.969], t = 7.18, p = 0.019
fig3d, high face ACh frontal-motor minus visual: [0.134 0.319]; paired t = 2.45, p = 0.247
s7 gradient, ach face_high minus face_low: rho = -0.874, p = 5.28e-08
s7 gradient, ca loc minus face_high: rho = 0.437, p = 0.0372
ach face_high - face_low: significant increases 57, decreases 0 of 253; increase in motor-somatosensory pairs 0.778, other pairs 0.106
ach loc - face_high: significant increases 0, decreases 234 of 253
ca face_high - face_low: significant increases 128, decreases 0 of 253
ca loc - face_high: significant increases 0, decreases 85 of 253; decrease in posterior pairs 0.727, other pairs 0.198 bg face_high - face_low: significant increases 7, decreases 0 of 23
bg loc - face_high: significant increases 0, decreases 0 of 23; frontal decreases 0 of 2, MOs and MOp p [0.246 0.424]
Result (variant, does not change an outcome): Without grabAM06 every direction holds in both remaining mice. Most tests across mice lose significance: ACh high minus low face LSD p = 0.207, locomotion minus high face p = 0.543, calcium high minus low p = 0.266, frontal-motor minus visual paired p = 0.247. Calcium locomotion minus high face stays significant (LSD p = 0.019). The ACh sustained-state gradient stays (rho = -0.874, p = 5.3e-08), and the calcium one (locomotion minus high face) becomes significant (rho = 0.437, p = 0.037, against p = 0.125 with 3 mice), which would meet the calcium half of the s7-ap-gradient-sustained criterion. The permutation tests give similar patterns with fewer pairs: ACh correlations rise in 57 pairs (0.778 of motor-somatosensory pairs, 0.106 of others) and fall in 234 with locomotion; calcium correlations rise in 128 and fall in 85 (0.727 of posterior pairs, 0.198 of others); ACh-calcium coupling rises in 7 of 23 parcels and falls in none with locomotion. So the significance of fig3c-ach-sustained-face and fig3d-ach-frontal-above-visual depends on grabAM06, the calcium gradient of s7-ap-gradient-sustained is significant only without it, and the counts of Figures 4 and 5 keep their pattern.
print(f"notebook ran in {(time.time() - started) / 60:.1f} min")
notebook ran in 16.6 min
Summary¶
| Claim | Outcome |
|---|---|
| allen-parcel-link | partially reproduced |
| fig2a-r2-differs-across-areas | reproduced |
| fig2b-ap-gradient-inverted | reproduced |
| fig2c-face-best-predictor-ach | reproduced |
| fig2c-face-best-predictor-calcium | reproduced |
| s7-face-onset-transients | partially reproduced |
| s7-locomotion-onset-transients | partially reproduced |
| s7-ach-transients-same-across-states | reproduced |
| s7-calcium-transients-largest-at-locomotion | partially reproduced |
| s7-ap-gradient-transients | reproduced |
| fig3b-ecog-desynchronization | blocked: data missing |
| fig3c-ach-sustained-face | reproduced |
| fig3c-ach-locomotion-no-further-increase | reproduced |
| fig3c-calcium-sustained | partially reproduced |
| s7-ap-gradient-sustained | partially reproduced |
| fig3d-ach-frontal-above-visual | reproduced |
| fig3e-scopolamine-motor-not-visual | blocked: data missing |
| fig4a-ach-correlations-up-with-face | reproduced |
| fig4a-ach-correlations-down-with-locomotion | reproduced |
| fig4b-calcium-correlations-up-with-face | reproduced |
| fig4b-calcium-correlations-down-posterior | reproduced |
| fig5-coupling-up-with-face | partially reproduced |
| fig5-coupling-down-frontal-with-locomotion | not reproduced |
| s3-hemodynamic-correction | blocked: data missing |
| s2-ach-sensor-validation | blocked: data missing |
| s8-cholinergic-axons | blocked: data missing |
A failure to reproduce here may be a failure of this attempt, not of the paper: see the reasoning and the differences from the paper's method.
How the analysis differs from the paper's. 4 of the paper's 6 mice are in the dandiset, and only 3 have matched sustained states (grabAM07 has no sustained high face epoch), so Figures 3c, 3d, 4 and 5 rest on 3 mice. One air-puff session without an air-puff table is left out. The areas come from the code's parcel masks. The ridge regression code is not public; its penalty is chosen here by marginal likelihood over a grid. The wheel predictor is the signed speed, as the Methods and code define it (with |speed| the calcium fig2c claim would be partial); for the three mice that can be matched to rows of the paper's own per-mouse R2 file, the full-model R2 here is close to the paper's. With 4 mice the paper's signed-rank tests cannot reach p < 0.05, so the direction in every mouse is used. Calcium amplitudes at locomotion are larger here than in the paper. grabAM06 contributes 2 matched epochs per state; without it the directions of Figures 3c and 3d hold and most of their tests across mice are not significant, but calcium locomotion above high face stays significant (p = 0.019) and the calcium sustained-state gradient becomes significant (p = 0.037).
What was found about the dandiset. The ΔF/F series start 49.9 s too early (frame 0 is raw frame 499); the FaceMap and pupil timestamps are spaced 111 ms where the camera ran at 10 Hz and are not monotonic; the parcellated series is calcium but is linked to the ACh imaging plane, its masks are transposed and its start is 64.9 s early; the ΔF/F descriptions omit the hemodynamic correction and atlas registration (the data are consistent with a correction having been applied); one calcium movie has 263 pixels in the right AUD parcel that are zero in every frame; WheelSpeed is signed; the Allen parcels used by the paper are not in the files; the visual stimulus table's columns from stimulus size on are shifted by one; one air-puff session has no air-puff table; the air-puff tables appear to hold only the puffs without running; wheel speed and pupil area have no units; two of the paper's six mice, the ECoG, the drug sessions and the control and ChAT mice are not in the dandiset.