Pipeline

  1"""S00 — Data provenance, quality control and session mapping."""
  2
  3from __future__ import annotations
  4
  5import math
  6
  7import matplotlib.pyplot as plt
  8import numpy as np
  9import pandas as pd
 10from scipy import stats
 11
 12from ..behavior import clock_shift_s
 13from ..data import CONDITION_COLORS
 14from ..stats import fmt_p, perm_p
 15from .base import Study
 16
 17
 18def _eta2(x: np.ndarray, lab: np.ndarray) -> float:
 19    keep = lab != "none"
 20    x, lab = x[keep], lab[keep]
 21    if len(lab) < 30 or len(np.unique(lab)) < 2:
 22        return np.nan
 23    mu = x.mean(0)
 24    tot = ((x - mu) ** 2).sum(0)
 25    btw = np.zeros(x.shape[1])
 26    for u in np.unique(lab):
 27        m = lab == u
 28        btw += m.sum() * (x[m].mean(0) - mu) ** 2
 29    return float(np.mean(btw / np.maximum(tot, 1e-12)))
 30
 31
 32class DataQC(Study):
 33    id = "s00"
 34    number = 0
 35    title = "Data provenance, quality control and session mapping"
 36    short = "Recovering session boundaries from concatenated recordings and validating them against behavior."
 37    question = ("How are the calcium recordings organised in time, and can each imaging block be assigned to its "
 38                "behavioral session with evidence rather than assumption?")
 39    background = """
 40**Acquisition and preprocessing (as reported by the lab).** Six male C57BL/6N mice (6–7 weeks at arrival) were
 41imaged in prelimbic mPFC (AP +1.80 mm, ML +0.35 mm, DV −2.30 mm, no angle) with a Doric 2-color fluorescence
 42miniscope (0.5 mm lens SICL_D_500_80) and AAV2/9-Syn-GCaMP6f. Frame rate is 6.666 frames/s. Imaging days are
 43CSDS days 1, 2, 3, 5, 7 and 9, plus day 11. Each defeat day is a 15-min recording split into baseline (BL),
 44social defeat (SD) and sensory stress (SS). Epoch length is nominally 2,000 frames and varies from 1,998 to
 452,001. Raw frames (630 × 630 px) had black frames replaced by the mean of neighbouring frames and were
 46spatially downsampled by a factor of 3 (210 × 210 px). Slow baseline drift in mean frame intensity was removed
 47with a rubberband (convex-hull) baseline fitted by piecewise-linear and smoothing splines. Minian could not
 48reliably detect neurons in the defeat data, so ROIs were drawn manually on a maximum projection after epochs
 49were merged within each day and fields of view were checked across days, then aligned by hand across days:
 50848 neurons in total. ΔF/F was computed on the full concatenated recording (F₀ = 10th percentile in each 12-s
 51window; the memoir approximated the frame rate as 7 Hz, about 84 frames). CASCADE converted each ΔF/F sample
 52to a spike probability (Rupprecht et al. 2021; model and noise level not recorded). ANY-maze triggered the
 53miniscope, so behavior and calcium start together. Annotation files of about 305 s (155 s for each SIT half)
 54use a clock that is 5 s ahead of calcium time.
 55
 56**Why a QC study.** The exported `TS_{id}.xlsx` workbooks hold one long neuron × frame matrix per animal with no
 57session markers. Every downstream contrast depends on assigning each ~2,000-frame block to the right session.
 58The lab's description fixes the order within a day, and this study verifies it from the data. Day 11
 59(BL → SIT → POST) is in the workbooks; the SIT block is two halves of about 1,000 frames (empty enclosure,
 60then CD-1) and was not used in the previous article.
 61
 62**Day 7 and incomplete epochs.** Day 7 was not recorded for mouse 6597, and sensory stress on day 9 was not
 63included for that mouse. Mouse 6669 has no SIT imaging block. Day-7 data exist for the other five mice.
 64Primary analyses still omit day 7 so that every animal is compared on the same days.
 65"""
 66    methods = """
 67**Data format.** Each workbook holds two signal families, `dff_mc` (motion-corrected ΔF/F) and
 68`predictions_dff_mc` (inferred spike probability). Rows are neurons (row 0 is the frame index) and columns are
 69frames; time is split into `part1/2/3` sheets only because Excel limits a sheet to 16,384 columns. Parts were
 70concatenated in order. The first and last 33 frames of each recording are constant. CASCADE was used only to
 71write the spike-probability columns, so those samples are excluded as unusable frames rather than as
 72documented CASCADE padding.
 73
 74**Block structure.** Behavioral sessions last ≈300 s. Each SIT half lasts ≈150 s of calcium time inside one
 75block (annotation files read ≈155 s because of the 5 s clock). We tested whether the
 76recording is composed of fixed 2,000-frame blocks (300 s at 6.666 Hz). The lab reports that epoch length can
 77be 1,998–2,001 frames, so a boundary can sit a few frames off this grid. For every candidate boundary
 78*t* = 2,000·*k*, a boundary score was computed as the Euclidean distance between the population mean ΔF/F
 79vectors in the 30 frames before and after *t*. This was compared with the same score at 400 random time
 80points per animal, and each boundary is reported as a percentile of that null distribution. True session
 81boundaries should sit in the upper tail, because focal plane, baseline fluorescence and behavioral context
 82change between sessions. This test uses the nominal grid, so it is not biased by the correction below.
 83
 84**Epoch-length correction.** Around each nominal cut *k* (frame 2,000·*k*), the step score (20 frames each
 85side) was computed at every frame within ±15 frames. A step was *reliable* if it exceeded the 95th percentile
 86of step scores at 400 random times and did not sit at the edge of the search window. Because offsets
 87accumulate from frame 0, offset_k = k·δ, where δ is the mean deviation of epoch length from 2,000 frames. δ was
 88estimated per animal as the median of offset_k / k over reliable cuts (≥ 4 required; otherwise δ = 0), limited
 89to the lab's reported 1,998–2,001-frame range. All analyses cut epoch *k* at round(k·(2,000 + δ)). One
 90parameter per animal is estimated, so individual noisy detections do not move individual boundaries.
 91
 92**Session order hypothesis.** The protocol implies, for each of days 1–9, BL → SD → SS. On day 11 it implies
 93BL → SIT (empty enclosure, then CD1 target, ≈1,000 frames each) → POST. The lab identified the missing epoch
 94for the two short recordings (6597: SS on day 9; 6669: SIT). Every single-slot deletion is still scored below
 95as a check on that assignment.
 96
 97**Behavior-encoding alignment test.** If block *b* truly corresponds to session *s*, neural activity in *b*
 98should be better explained by the ethogram of *s* than by the ethogram of any other session. For every
 99(block, session) pair we:
1001. Binned inferred spikes in 1-s bins and z-scored each neuron.
1012. Assigned each bin the annotated behavior of the candidate session.
1023. Computed η², the fraction of each neuron's variance explained by behavior labels, averaged over neurons.
1034. Expressed η² as a z-score against 40 circular shifts of the label sequence, which preserves bout structure
104   and destroys alignment.
105
106The mean z on the hypothesised diagonal was compared with 5,000 random block→session permutations, per animal
107and pooled across animals with complete mappings. For animals with one block too few, the mean diagonal z was
108computed for each possible missing slot, giving a likelihood ranking of candidate mappings.
109
110**Frame-to-video timing.** Calcium and behavior start together. Frame *i* of a block maps to annotation time
111*i*/6.666 s, plus 5 s when that file's duration is ≈305 s or ≈155 s.
112"""
113    limitations = """
114* The drift model assumes a constant mean epoch length per animal. Individual epochs vary by 1–2 frames
115  (≈ 0.2 s) around it, which the model does not recover.
116* ΔF/F baselines were estimated on the concatenated recording, so the 12-s F₀ window can cross an epoch boundary.
117* The 33 constant frames at each end of a workbook are excluded. Their origin is not documented.
118* CASCADE model and noise level were not recorded. Spike probability is paired one-to-one with ΔF/F.
119* Mouse 6597 has no day-9 sensory-stress block, and mouse 6669 has no SIT imaging block.
120"""
121    references = [
122        "Rupprecht P. et al. (2021) A database and deep learning toolbox for noise-optimized, generalized spike inference from calcium imaging. Nat Neurosci 24:1324–1337.",
123        "Dong Z. et al. (2022) Minian, an open-source miniscope analysis pipeline. eLife 11:e70661.",
124        "Optimal choice of baseline correction for multivariate calibration of spectra (rubberband baseline method; cited by the lab).",
125        "Golden S.A. et al. (2011) A standardized protocol for repeated social defeat stress in mice. Nat Protoc 6:1183–1191.",
126    ]
127
128    def run(self) -> dict:
129        ds = self.ds
130        rng = np.random.default_rng(0)
131        fs, block = ds.fs, ds.block
132
133        summary, boundary_rows = [], []
134        for sid in ds.subjects:
135            d = ds.neural(sid, "dff")
136            total = d.shape[1]
137            n_blocks = math.ceil(total / block)
138            sess = ds.all_sessions[ds.all_sessions.subject == sid]
139            summary.append({"subject": sid, "neurons": d.shape[0], "frames": total,
140                            "minutes": round(total / fs / 60, 1), "blocks": n_blocks,
141                            "protocol_slots": len(ds._slots(sid)), "sessions_mapped": len(sess),
142                            "mapping": sess.mapping.iloc[0]})
143            score = lambda t, d=d: np.linalg.norm(d[:, t - 30:t].mean(1) - d[:, t:t + 30].mean(1))
144            rand = np.array([score(t) for t in rng.integers(100, total - 100, 400)])
145            for t in range(block, total - 50, block):
146                boundary_rows.append({"subject": sid, "frame": t, "percentile": float((rand < score(t)).mean())})
147        summary = pd.DataFrame(summary)
148        bnd = pd.DataFrame(boundary_rows)
149        frac90 = float((bnd.percentile > 0.9).mean())
150        binom_p = float(stats.binomtest(int((bnd.percentile > 0.9).sum()), len(bnd), 0.1, alternative="greater").pvalue)
151
152        # Behavior-encoding alignment matrices
153        align, drop_rows, perm_rows = {}, [], []
154        for sid in ds.subjects:
155            p = ds.neural(sid, "spikes")
156            total = p.shape[1]
157            n_blocks = math.ceil(total / block)
158            bounds = ds.epoch_boundaries(sid)
159            slots = ds._slots(sid)
160            labels = {}
161            for cond, day in slots:
162                labels[(cond, day)] = self._slot_labels(sid, cond, day, n_bins=300)
163            m = np.full((n_blocks, len(slots)), np.nan)
164            for b in range(n_blocks):
165                seg = p[:, bounds[b]:bounds[b + 1]]
166                n_bins = int(seg.shape[1] / fs)
167                edges = (np.arange(n_bins + 1) * fs).astype(int)
168                x = np.stack([seg[:, edges[i]:edges[i + 1]].mean(1) for i in range(n_bins)])
169                x = (x - x.mean(0)) / (x.std(0) + 1e-9)
170                for j, sl in enumerate(slots):
171                    lab = labels[sl][: len(x)]
172                    e = _eta2(x, lab)
173                    if np.isnan(e):
174                        continue
175                    null = [_eta2(x, np.roll(lab, rng.integers(20, len(lab) - 20))) for _ in range(40)]
176                    m[b, j] = (e - np.nanmean(null)) / (np.nanstd(null) + 1e-12)
177            align[sid] = (m, slots)
178            if n_blocks == len(slots):
179                diag = np.nanmean(np.diag(m))
180                null = [np.nanmean(m[np.arange(n_blocks), rng.permutation(n_blocks)]) for _ in range(5000)]
181                perm_rows.append({"subject": sid, "diagonal_z": diag, "shuffled_z_mean": float(np.mean(null)),
182                                  "p_perm": perm_p(diag, np.array(null))})
183            else:
184                for j, sl in enumerate(slots):
185                    keep = [k for k in range(len(slots)) if k != j]
186                    drop_rows.append({"subject": sid, "missing_slot": f"{sl[0]}{sl[1]}",
187                                      "diagonal_z": float(np.nanmean([m[b, keep[b]] for b in range(n_blocks)]))})
188        perm_df = pd.DataFrame(perm_rows)
189        drop_df = pd.DataFrame(drop_rows)
190        complete = [sid for sid in ds.subjects if align[sid][0].shape[0] == len(align[sid][1])]
191        obs = np.mean([np.nanmean(np.diag(align[s][0])) for s in complete])
192        null = [np.mean([np.nanmean(align[s][0][np.arange(len(align[s][1])), rng.permutation(len(align[s][1]))])
193                         for s in complete]) for _ in range(5000)]
194        pooled_p = perm_p(obs, np.array(null))
195
196        bnd_rows, cut_rows = [], []
197        for sid in ds.subjects:
198            info = ds.boundary_info[sid]
199            rel = [o for o in info["offsets"] if o["reliable"]]
200            n_cuts = info["n_cuts"]
201            bnd_rows.append({"subject": sid, "cuts": n_cuts, "reliable_cuts": info["n_reliable"],
202                             "epoch_frames": round(info["epoch_frames"], 3),
203                             "max_shift_vs_fixed_grid": int(round(n_cuts * info["delta_frames"])),
204                             "median_abs_residual": float(np.median([abs(o["offset"] - o["cut"] * info["delta_frames"])
205                                                                     for o in rel])) if rel else np.nan,
206                             "last_epoch_frames": int(ds.neural(sid).shape[1] - round(n_cuts * info["epoch_frames"]))})
207            cut_rows += [{"subject": sid, **o, "predicted_offset": round(o["cut"] * info["delta_frames"], 2)}
208                         for o in info["offsets"]]
209        bnd_df, cut_df = pd.DataFrame(bnd_rows), pd.DataFrame(cut_rows)
210
211        # Figures
212        self._fig_population(bnd)
213        self._fig_cuts(cut_df, bnd_df)
214        self._fig_alignment(align, complete)
215        if not drop_df.empty:
216            self._fig_drop(drop_df)
217        self._fig_raster()
218
219        sessions = ds.all_sessions[["subject", "condition", "day", "block", "start", "stop", "n_frames", "obs_id", "mapping", "included"]]
220        self.table(summary, "dataset_summary", "Recording summary per animal",
221                   "Neurons, frames and the number of 2,000-frame blocks versus protocol slots.")
222        self.table(bnd_df, "epoch_lengths", "Estimated epoch length per animal",
223                   "Per-epoch length from the drift of reliable boundary steps away from the 2,000-frame grid; "
224                   "max shift = correction applied at the last cut.")
225        self.table(cut_df, "boundary_offsets", "Detected boundary offsets",
226                   "Offset (frames) of the largest population ΔF/F step within ±15 frames of each nominal cut.", show=20)
227        self.table(perm_df, "alignment_permutation", "Behavior-encoding alignment (complete mappings)",
228                   "Mean z of the hypothesised block→session diagonal vs 5,000 random assignments.")
229        if not drop_df.empty:
230            self.table(drop_df.sort_values(["subject", "diagonal_z"], ascending=[True, False]), "missing_slot_ranking",
231                       "Candidate missing sessions (animals with one block too few)",
232                       "Higher diagonal z = better behavioral fit of the remaining mapping.", show=12)
233        self.table(sessions, "session_map", "Session map used by all studies",
234                   "Frame ranges (in the concatenated recording) for every session.", show=30)
235
236        n_neurons = int(summary.neurons.sum())
237        ranked = perm_df.sort_values("diagonal_z", ascending=False)
238        strong = ", ".join(str(s) for s in ranked[ranked.p_perm < 0.05].subject) or "none individually"
239        weak = ", ".join(str(s) for s in ranked[ranked.p_perm >= 0.05].subject) or "none"
240        late = []
241        for sid, d in drop_df.groupby("subject") if not drop_df.empty else []:
242            top = d.sort_values("diagonal_z", ascending=False).head(4)
243            late.append(f"{sid}: " + ", ".join(f"{r.missing_slot} ({r.diagonal_z:.2f})" for r in top.itertuples()))
244        self.highlight("Cohort", f"{len(ds.subjects)} mice · {n_neurons} neurons",
245                       f"{summary.neurons.min()}–{summary.neurons.max()} neurons per animal tracked across "
246                       f"{int(ds.all_sessions.shape[0])} mapped sessions ({int(ds.sessions.shape[0])} analysed after excluding day 7).", "info")
247        self.highlight("Epoch length", f"{bnd_df.epoch_frames.min():.2f}–{bnd_df.epoch_frames.max():.2f} frames",
248                       f"Estimated from {int(bnd_df.reliable_cuts.sum())} reliable boundary steps. Correcting the fixed "
249                       f"2,000-frame grid shifts late cuts by up to {int(bnd_df.max_shift_vs_fixed_grid.max())} frames "
250                       f"({bnd_df.max_shift_vs_fixed_grid.max() / fs:.1f} s).", "info")
251        self.highlight("Session boundaries", f"{frac90:.0%} of boundaries > 90th pct",
252                       f"Population activity changes abruptly every 2,000 frames (binomial vs 10% chance, {fmt_p(binom_p)}).",
253                       "significant" if binom_p < 0.05 else "null")
254        self.highlight("Session order validated" if pooled_p < 0.05 else "Session order not confirmed", f"pooled {fmt_p(pooled_p)}",
255                       f"Neural blocks are better explained by their own session's ethogram than by shuffled assignments "
256                       f"(mean diagonal z = {obs:.2f}, {len(complete)} animals).",
257                       "significant" if pooled_p < 0.05 else "null")
258        missing = {int(k): v for k, v in (ds.cfg.get("session_mapping", {}).get("missing_slots") or {}).items()}
259        if not drop_df.empty:
260            named = ", ".join(f"{sid} {slot}" for sid, slot in missing.items())
261            self.highlight("Missing epochs", named,
262                           "Identified by the lab and omitted from the session map. Figure 4 ranks every other "
263                           "single-slot deletion as a check.", "info")
264
265        results = f"""
266**Recording structure.** The {len(ds.subjects)} animals contributed {n_neurons} neurons in total
267({summary.neurons.min()}–{summary.neurons.max()} per animal), each tracked continuously for
268{summary.minutes.min()}–{summary.minutes.max()} min of imaging. Recordings split into ≈2,000-frame blocks.
269{frac90:.0%} of nominal block edges exceeded the 90th percentile of random-time boundary scores (expected 10%;
270binomial {fmt_p(binom_p)}). Figure 1 shows the abrupt population-level steps at block edges.
271
272**Epoch length.** Reliable boundary steps ({int(bnd_df.reliable_cuts.sum())} of {int(bnd_df.cuts.sum())} cuts) sat
273progressively later than the 2,000-frame grid, which means epochs are slightly longer than 2,000 frames. The
274estimated length was {bnd_df.epoch_frames.min():.2f}–{bnd_df.epoch_frames.max():.2f} frames per animal, inside the lab's
275reported 1,998–2,001 range (Table 2, Figure 2). The median absolute residual of detected cuts around the drift
276model was {np.nanmedian(bnd_df.median_abs_residual):.1f} frames. Using the corrected cuts moves the last session
277boundary by up to {int(bnd_df.max_shift_vs_fixed_grid.max())} frames ({bnd_df.max_shift_vs_fixed_grid.max() / fs:.1f} s), which would
278otherwise misalign late sessions with their behavior annotations. The final epoch of every recording is short
279({bnd_df.last_epoch_frames.min()}–{bnd_df.last_epoch_frames.max()} frames), consistent with a truncated last epoch.
280
281**Session order.** Under the protocol order (BL → SD → SS on days 1–9; BL → SIT → POST on day 11), the
282behavior-encoding alignment z on the hypothesised diagonal averaged {obs:.2f} across the {len(complete)} animals
283with complete mappings, versus {np.mean(null):.2f} for 5,000 random block assignments (pooled permutation
284{fmt_p(pooled_p)}; per-animal values in Table 4, Figure 3). Individually significant animals: {strong}. Not
285individually significant: {weak}.
286
287**Incomplete recordings.** The lab identified the missing epochs: mouse 6597 has no day-7 recording and no
288sensory-stress block on day 9; mouse 6669 has no SIT imaging block. Alignment rankings of single-slot
289deletions were {'; '.join(late) if late else 'n/a'}. Deleting an early slot fits clearly worse (Figure 4).
290The map uses the lab's assignment, not the top-scoring candidate.
291"""
292        return {"results": results, "parameters": {"block_frames": block, "fs_hz": fs, "boundary_window": 30,
293                                                   "alignment_bin_s": 1.0, "alignment_null_shifts": 40,
294                                                   "permutations": 5000}}
295
296    def _slot_labels(self, sid: int, cond: str, day: int, n_bins: int) -> np.ndarray:
297        ds = self.ds
298        lab = np.full(n_bins, "none", dtype=object)
299        def paint(obs, offset):
300            ev = ds.behavior[ds.behavior["Observation id"] == obs]
301            if ev.empty:
302                return
303            shift = clock_shift_s(float(ev["Total length"].iloc[0]), float(ev["Start (s)"].min()))
304            for a, b, beh in zip(ev["Start (s)"], ev["Stop (s)"], ev["Behavior"]):
305                i0, i1 = int(np.floor(a - shift + offset)), int(np.ceil(b - shift + offset))
306                lab[max(i0, 0):min(i1, n_bins)] = beh
307        if cond == "SIT":
308            no = ds._obs_id(sid, "SIT", day, "NoCD1")
309            paint(no, 0.0)
310            ev_no = ds.behavior[ds.behavior["Observation id"] == no]
311            length = float(ev_no["Total length"].iloc[0]) - clock_shift_s(
312                float(ev_no["Total length"].iloc[0]), float(ev_no["Start (s)"].min()))
313            cd1 = np.full(n_bins, "none", dtype=object)
314            saved = lab.copy()
315            lab[:] = "none"
316            paint(ds._obs_id(sid, "SIT", day, "CD1"), length)
317            cd1[:] = lab
318            lab[:] = saved
319            lab[int(length):] = cd1[int(length):]
320        else:
321            obs = ds._obs_id(sid, cond, day)
322            if obs:
323                paint(obs, 0.0)
324        return lab
325
326    def _fig_population(self, bnd: pd.DataFrame) -> None:
327        ds = self.ds
328        fig, axes = plt.subplots(len(ds.subjects), 1, figsize=(11, 1.6 * len(ds.subjects)), sharex=False)
329        for ax, sid in zip(axes, ds.subjects):
330            d = ds.neural(sid, "dff")
331            m = np.convolve(d.mean(0), np.ones(25) / 25, "same")
332            t = np.arange(len(m)) / ds.fs / 60
333            for r in ds.all_sessions[ds.all_sessions.subject == sid].itertuples():
334                ax.axvspan(r.start / ds.fs / 60, r.stop / ds.fs / 60, color=CONDITION_COLORS[r.condition], alpha=0.18, lw=0)
335            ax.plot(t, m, color="k", lw=0.4)
336            ax.set_ylabel(f"{sid}\nΔF/F", fontsize=8)
337            ax.set_xlim(0, t[-1])
338        axes[-1].set_xlabel("Imaging time (min, concatenated)")
339        handles = [plt.Rectangle((0, 0), 1, 1, color=CONDITION_COLORS[c], alpha=0.4) for c in CONDITION_COLORS]
340        axes[0].legend(handles, list(CONDITION_COLORS), ncol=6, fontsize=7, loc="upper center", bbox_to_anchor=(0.5, 1.6))
341        self.figure(fig, "population_blocks", "Population activity and recovered sessions",
342                    "Population-mean ΔF/F (25-frame smoothing) for each animal across the full concatenated recording. "
343                    "Shading marks the recovered sessions, coloured by condition. Steps in baseline fluorescence coincide "
344                    "with 2,000-frame block edges.")
345
346    def _fig_cuts(self, cuts: pd.DataFrame, bnd: pd.DataFrame) -> None:
347        subs = list(bnd.subject)
348        fig, axes = plt.subplots(1, len(subs), figsize=(2.4 * len(subs), 2.8), sharey=True)
349        for ax, sid in zip(np.atleast_1d(axes), subs):
350            d = cuts[cuts.subject == sid]
351            rel = d[d.reliable]
352            ax.scatter(d[~d.reliable].cut, d[~d.reliable].offset, s=10, facecolor="none", edgecolor="0.6")
353            ax.scatter(rel.cut, rel.offset, s=14, color="k")
354            ax.plot(d.cut, d.predicted_offset, color="#b5523b", lw=1.2)
355            ax.axhline(0, color="0.5", ls=":", lw=0.7)
356            ep = float(bnd.loc[bnd.subject == sid, "epoch_frames"].iloc[0])
357            ax.set_title(f"{sid}: {ep:.2f} fr/epoch", fontsize=8)
358            ax.set_xlabel("cut k")
359        np.atleast_1d(axes)[0].set_ylabel("offset vs 2,000·k (frames)")
360        self.figure(fig, "boundary_drift", "Epoch boundaries drift from the 2,000-frame grid",
361                    "Offset of the largest population ΔF/F step within ±15 frames of each nominal cut. Filled points are "
362                    "reliable steps (> 95th percentile of random-time steps, not at the window edge); open points are not. "
363                    "Red line: drift model k·δ used to place the cuts.")
364
365    def _fig_alignment(self, align: dict, complete: list[int]) -> None:
366        fig, axes = plt.subplots(1, len(complete), figsize=(3.3 * len(complete), 3.3))
367        axes = np.atleast_1d(axes)
368        for ax, sid in zip(axes, complete):
369            m, slots = align[sid]
370            im = ax.imshow(m, cmap="RdBu_r", vmin=-4, vmax=4)
371            ax.plot(np.arange(len(slots)), np.arange(len(slots)), "k.", ms=2)
372            ax.set_title(str(sid))
373            ax.set_xticks(range(len(slots)), [f"{c}{d}" for c, d in slots], rotation=90, fontsize=5)
374            ax.set_yticks(range(m.shape[0]), [str(i) for i in range(m.shape[0])], fontsize=5)
375            ax.set_xlabel("Behavioral session")
376        axes[0].set_ylabel("Imaging block")
377        fig.colorbar(im, ax=axes, shrink=0.7, label="alignment z")
378        self.figure(fig, "alignment_matrices", "Behavior-encoding alignment matrices",
379                    "z-scored variance in neural activity explained by each session's ethogram (columns) for each imaging "
380                    "block (rows). Black dots mark the hypothesised protocol order; positive values along the diagonal "
381                    "support the mapping.")
382
383    def _fig_drop(self, drop_df: pd.DataFrame) -> None:
384        subs = drop_df.subject.unique()
385        fig, axes = plt.subplots(1, len(subs), figsize=(5.5 * len(subs), 3))
386        for ax, sid in zip(np.atleast_1d(axes), subs):
387            d = drop_df[drop_df.subject == sid]
388            missing = {int(k): v for k, v in (self.ds.cfg.get("session_mapping", {}).get("missing_slots") or {}).items()}
389            ax.bar(range(len(d)), d.diagonal_z, color=["#b5523b" if s == missing.get(int(sid)) else "0.6" for s in d.missing_slot])
390            ax.set_xticks(range(len(d)), d.missing_slot, rotation=90, fontsize=7)
391            ax.set_title(f"{sid}: fit if this slot is missing")
392            ax.set_ylabel("mean diagonal z")
393        self.figure(fig, "missing_slot_candidates", "Which session is missing?",
394                    "Mean alignment z of the remaining mapping when each protocol slot is assumed missing. Early deletions "
395                    "fit worse. Red marks the epoch the lab identified as absent (6597: SS day 9; 6669: SIT).")
396
397    def _fig_raster(self) -> None:
398        ds = self.ds
399        sid = 6675
400        sess = list(ds.iter_sessions(sid, ["BL", "SD", "SS"]))
401        sess = [s for s in sess if s.day == 5]
402        fig, axes = plt.subplots(2, len(sess), figsize=(13, 5), gridspec_kw={"height_ratios": [5, 1]})
403        behs = sorted({b for s in sess for b in ds.labels(s)} - {"none"})
404        cmap = plt.get_cmap("tab20")
405        colors = {b: cmap(i % 20) for i, b in enumerate(behs)}
406        for j, s in enumerate(sess):
407            x = ds.activity(s)
408            order = np.argsort(-x.mean(1))
409            axes[0, j].imshow(x[order], aspect="auto", cmap="magma", vmax=np.percentile(x, 99),
410                              extent=[0, s.n_frames / ds.fs, x.shape[0], 0])
411            axes[0, j].set_title(f"{sid} · {s.condition} · day {s.day}")
412            lab = ds.labels(s)
413            rgb = np.array([colors.get(b, (1, 1, 1, 1)) for b in lab])[None]
414            axes[1, j].imshow(rgb, aspect="auto", extent=[0, s.n_frames / ds.fs, 0, 1])
415            axes[1, j].set_yticks([])
416            axes[1, j].set_xlabel("time (s)")
417        axes[0, 0].set_ylabel("neuron (sorted by rate)")
418        handles = [plt.Rectangle((0, 0), 1, 1, color=colors[b]) for b in behs]
419        fig.legend(handles, behs, ncol=min(7, len(behs)), fontsize=7, loc="lower center", bbox_to_anchor=(0.5, -0.06))
420        self.figure(fig, "raster_example", "Example day: inferred spikes with the ethogram",
421                    "Inferred spike probability for every neuron of animal 6675 on CSDS day 5 (baseline, defeat, post-defeat "
422                    "sensory contact), with the aligned behavioral annotation below each raster.")

  1"""S01 — Behavioral phenotyping across CSDS."""
  2
  3from __future__ import annotations
  4
  5import matplotlib.pyplot as plt
  6import numpy as np
  7import pandas as pd
  8
  9from ..data import CONDITION_COLORS
 10from ..stats import fdr_bh, fmt_p
 11from ..stats.compare import day_trend
 12from .base import Study
 13
 14DEFEAT = ["Physical interaction", "Attack", "Submissive", "Freezing", "Escape", "Danger"]
 15SOCIAL_SIT = ["Interaction", "Approach"]
 16
 17
 18class Behavior(Study):
 19    id = "s01"
 20    number = 1
 21    title = "Behavioral phenotyping across chronic social defeat"
 22    short = "Time budgets, the evolution of defeat responses across days, and a social-interaction phenotype."
 23    question = ("How does behavior evolve across repeated defeat, and how does social-test behavior relate to the "
 24                "lab's susceptible/resilient classification?")
 25    background = """
 26In the CSDS paradigm, C57BL/6 mice are defeated daily by a novel CD-1 aggressor and then housed in sensory
 27contact across a perforated divider (Golden et al. 2011). Susceptibility is classically defined in the social
 28interaction test (SIT) by a social-interaction ratio below 1: time near a CD-1 target divided by time near an
 29empty enclosure (Krishnan et al. 2007). Defeat episodes are characterised by attacks, submissive postures,
 30freezing and escape. How these responses change across days is itself informative about coping strategy
 31(Hollis & Kabbaj 2014). Here the ethograms are hand-scored state annotations (BORIS-style) for every session.
 32"""
 33    methods = """
 34**Annotations.** State events (start/stop, seconds) were loaded for every session. Labels were canonicalised:
 35`standing`→Standing, `exploring`→Exploring, `grooming`→Grooming, `escape/huir`→Escape. Observation IDs encode
 36animal, condition (BL, SD, SS, SIT_NoCD1, SIT_CD1, POST_SIT), CSDS day and date.
 37
 38**Time budgets.** For each session, the time in each behavior was divided by total annotated time. Budgets
 39were averaged per animal within condition, then across animals.
 40
 41**Defeat dynamics.** Within SD sessions (days 1–9 without day 7), the fraction of time spent in each defeat-related state
 42(physical interaction, attack, submissive, freezing, escape, danger) was regressed on CSDS day separately for
 43each animal. The animal is the unit of inference: the six per-animal slopes were tested against 0 with an exact
 44two-sided Wilcoxon signed-rank test (summary-statistic random-effects approach; Holmes & Friston 1998). The
 45reported slope is the pooled within-animal estimate from `fraction ~ day + animal` (OLS with animal fixed
 46effects). The same analysis was applied to *non-physical interaction* in BL and SS sessions, i.e. time spent
 47investigating the aggressor across the divider. A decline would indicate emerging social avoidance. P-values
 48were FDR-corrected (Benjamini–Hochberg) across all trend tests. With six animals, the smallest attainable
 49two-sided p is 0.031.
 50
 51**SIT phenotype.** The susceptible/resilient label of each animal was provided by the lab. To check how much of
 52it is visible in the annotations, social engagement was quantified as the seconds annotated as *Interaction* or
 53*Approach* in each SIT phase (≈150 s each). The social-interaction ratio proxy is SI = engagement(CD1) / engagement(empty), and
 54animals with SI < 1 are labelled *susceptible-like*. This is a proxy: the canonical ratio uses tracked time in an
 55interaction zone, which is not available in these files.
 56"""
 57    limitations = """
 58* The SI ratio is derived from ethogram states, not zone tracking. The lab pointed to the interaction- and
 59  corner-zone times (Shared_CSDS); those scores are not in this repository. n = 6, with no non-defeated control
 60  group. Phenotype labels are descriptive and are not used as a grouping factor for inferential tests.
 61* "Attack" refers to the annotated attack state in the focal animal's ethogram. Whether it codes being attacked
 62  or retaliation should be confirmed with the scorer.
 63* Observer reliability (inter-rater agreement) is not available in these files.
 64"""
 65    references = [
 66        "Golden S.A. et al. (2011) A standardized protocol for repeated social defeat stress in mice. Nat Protoc 6:1183–1191.",
 67        "Krishnan V. et al. (2007) Molecular adaptations underlying susceptibility and resistance to social defeat in brain reward regions. Cell 131:391–404.",
 68        "Hollis F. & Kabbaj M. (2014) Social defeat as an animal model for depression. ILAR J 55:221–232.",
 69        "Friard O. & Gamba M. (2016) BORIS: a free, versatile open-source event-logging software for video/audio coding and live observations. Methods Ecol Evol 7:1325–1330.",
 70        "Holmes A.P. & Friston K.J. (1998) Generalisability, random effects and population inference. NeuroImage 7:S754.",
 71    ]
 72
 73    def run(self) -> dict:
 74        b = self.ds.behavior[~self.ds.behavior.day.isin(self.ds.exclude_days)].copy()
 75        b["cond"] = np.where(b.condition == "SIT", "SIT_" + b.sit_target.fillna(""),
 76                             np.where(b.condition == "POST", "POST", b.condition))
 77        tot = b.groupby("Observation id")["Duration (s)"].sum().rename("total")
 78        frac = (b.groupby(["subject", "cond", "day", "Observation id", "Behavior"])["Duration (s)"].sum()
 79                .reset_index().join(tot, on="Observation id"))
 80        frac["fraction"] = frac["Duration (s)"] / frac["total"]
 81        wide = frac.pivot_table(index=["subject", "cond", "day", "Observation id"], columns="Behavior",
 82                                values="fraction", fill_value=0).reset_index()
 83        wide = wide.rename(columns={"cond": "condition"})
 84        self.ctx.cache["behavior_fractions"] = wide
 85
 86        # Time budgets
 87        budget = wide.drop(columns=["day", "Observation id"]).groupby(["subject", "condition"]).mean()
 88        budget = budget.groupby("condition").mean()
 89        order = [c for c in ["BL", "SD", "SS", "SIT_NoCD1", "SIT_CD1", "POST"] if c in budget.index]
 90        budget = budget.loc[order]
 91        self._fig_budget(budget)
 92
 93        # Defeat dynamics
 94        trend_rows = []
 95        sd = wide[(wide.condition == "SD") & (wide.day <= 9)]
 96        for beh in DEFEAT:
 97            if beh in sd:
 98                trend_rows.append(day_trend(sd.rename(columns={beh: "y"}), "y", "SD") | {"behavior": beh})
 99        for cond in ["BL", "SS"]:
100            d = wide[(wide.condition == cond) & (wide.day <= 9)].rename(columns={"Non-physical interaction": "y"})
101            trend_rows.append(day_trend(d, "y", cond) | {"behavior": "Non-physical interaction"})
102        trends = pd.DataFrame(trend_rows)[["behavior", "condition", "slope_per_day", "median_animal_slope", "n_positive",
103                                           "n_animals", "p", "n_sessions"]]
104        trends["q_fdr"] = fdr_bh(trends.p)
105        self._fig_defeat(wide)
106
107        # SIT
108        sit_ev = b[b.cond.str.startswith("SIT")]
109        eng = (sit_ev[sit_ev.Behavior.isin(SOCIAL_SIT)].groupby(["subject", "cond"])["Duration (s)"].sum()
110               .reindex(pd.MultiIndex.from_product([sorted(sit_ev.subject.unique()), ["SIT_CD1", "SIT_NoCD1"]],
111                                                   names=["subject", "cond"]), fill_value=0.0))
112        si = eng.unstack("cond").rename(columns={"SIT_CD1": "CD1_engagement_s", "SIT_NoCD1": "empty_engagement_s"})
113        si["SI_ratio_proxy"] = si["CD1_engagement_s"] / si["empty_engagement_s"].replace(0, np.nan)
114        si["phenotype_proxy"] = np.where(si.SI_ratio_proxy < 1, "susceptible-like", "resilient-like")
115        si = si.reset_index()
116        si["phenotype_lab"] = si.subject.map(self.ds.phenotype)
117        si["proxy_agrees"] = (si.phenotype_proxy.str.replace("-like", "") == si.phenotype_lab)
118        defeat_idx = sd.groupby("subject")[["Submissive", "Freezing", "Escape"]].mean().sum(1).rename("defeat_response_index")
119        si = si.join(defeat_idx, on="subject")
120        self.ctx.cache["phenotype"] = si
121        self._fig_sit(si)
122
123        self.table(trends, "defeat_trends", "Change in behavior across CSDS days",
124                   "Pooled within-animal slope (fraction of session time per day); p = exact Wilcoxon of per-animal slopes "
125                   "(animal = unit of inference); q = Benjamini–Hochberg.")
126        self.table(si, "sit_phenotype", "SIT engagement and social-interaction ratio proxy",
127                   "Engagement = annotated Interaction + Approach seconds in each SIT phase. Defeat-response index = mean fraction of SD "
128                   "time in submissive + freezing + escape.")
129        self.table(budget.reset_index().round(4), "time_budgets", "Mean time budget per condition",
130                   "Fraction of annotated time, averaged within animal then across animals.")
131
132        n_sus = int((si.phenotype_proxy == "susceptible-like").sum())
133        sig = trends[trends.q_fdr < 0.05]
134        for _, r in trends.sort_values("p").head(3).iterrows():
135            kind = "significant" if r.q_fdr < 0.05 else "trend" if r.p < 0.05 else "null"
136            self.highlight(f"{r.behavior} ({r.condition})", f"{r.slope_per_day * 100:+.2f} %/day",
137                           f"Change in time fraction per CSDS day; {r.n_positive}/{r.n_animals} animals positive, "
138                           f"exact Wilcoxon {fmt_p(r.p)} (q = {r.q_fdr:.3f}).", kind)
139        n_agree = int(si.proxy_agrees.sum())
140        self.highlight("Lab phenotype", f"{int((si.phenotype_lab == 'susceptible').sum())} susceptible · "
141                       f"{int((si.phenotype_lab == 'resilient').sum())} resilient",
142                       "Classification provided by the lab; used for all phenotype analyses (S10).", "info")
143        self.highlight("Ethogram SI proxy vs lab phenotype", f"{n_agree}/{len(si)} agree",
144                       f"{n_sus}/{len(si)} animals have an annotation-based SI ratio < 1. The proxy does not reproduce the "
145                       "lab classification well and is not used to label animals.", "caution")
146
147        sd_budget = budget.loc["SD"].sort_values(ascending=False)
148        top = lambda c: ", ".join(f"{k.lower()} {v:.0%}" for k, v in budget.loc[c].sort_values(ascending=False).head(3).items())
149        results = f"""
150**Time budgets.** The three dominant states were {top('BL')} in baseline sessions and {top('SS')} in
151post-defeat sensory-contact sessions. Defeat sessions carried a large defeat repertoire: physical interaction
152{sd_budget.get('Physical interaction', 0):.0%}, submissive {sd_budget.get('Submissive', 0):.0%}, freezing
153{sd_budget.get('Freezing', 0):.0%}, escape {sd_budget.get('Escape', 0):.0%} and attack
154{sd_budget.get('Attack', 0):.0%} of annotated time (Figure 1).
155
156**Defeat dynamics.** {len(sig)} of {len(trends)} day-trend tests survived FDR correction. The strongest
157effects are listed in Table 1 and plotted in Figure 2. Positive slopes mean the behavior occupies more of each
158session as defeats accumulate.
159
160**SIT phenotype.** The lab classified 4703, 6597 and 6669 as resilient and 6445, 6672 and 6675 as susceptible.
161An annotation-based proxy (Interaction + Approach time with the CD-1 vs the empty enclosure) gave SI ratios of
162{si.SI_ratio_proxy.min():.2f}–{si.SI_ratio_proxy.max():.2f}, with {n_sus} of {len(si)} animals below 1. The proxy matched the
163lab classification for {n_agree} of {len(si)} animals (Table 2, Figure 3). Ethogram states are not a substitute for
164zone-based tracking, so the lab classification is used for all phenotype analyses (S10).
165"""
166        return {"results": results, "parameters": {"sit_engagement_states": SOCIAL_SIT, "defeat_states": DEFEAT}}
167
168    def _fig_budget(self, budget: pd.DataFrame) -> None:
169        cols = budget.sum().sort_values(ascending=False).index
170        fig, ax = plt.subplots(figsize=(8, 3.6))
171        bottom = np.zeros(len(budget))
172        cmap = plt.get_cmap("tab20")
173        for i, c in enumerate(cols):
174            ax.bar(budget.index, budget[c], bottom=bottom, color=cmap(i % 20), label=c, width=0.7)
175            bottom += budget[c].to_numpy()
176        ax.set_ylabel("fraction of annotated time")
177        ax.legend(ncol=1, fontsize=7, bbox_to_anchor=(1.01, 1), loc="upper left")
178        self.figure(fig, "time_budgets", "Time budget by condition",
179                    "Mean fraction of annotated time spent in each behavior, per condition (animal-averaged).")
180
181    def _fig_defeat(self, wide: pd.DataFrame) -> None:
182        sd = wide[(wide.condition == "SD") & (wide.day <= 9)]
183        fig, axes = plt.subplots(1, 2, figsize=(10, 3.4))
184        cmap = plt.get_cmap("tab10")
185        for i, beh in enumerate([b for b in DEFEAT if b in sd]):
186            g = sd.groupby("day")[beh].agg(["mean", "sem"])
187            axes[0].errorbar(g.index, g["mean"], yerr=g["sem"], marker="o", ms=3, capsize=2, label=beh, color=cmap(i))
188        axes[0].set_xlabel("CSDS day")
189        axes[0].set_ylabel("fraction of SD session")
190        axes[0].set_title("Defeat repertoire across days")
191        axes[0].legend(fontsize=7)
192        for cond in ["BL", "SS"]:
193            d = wide[(wide.condition == cond) & (wide.day <= 9)].groupby("day")["Non-physical interaction"].agg(["mean", "sem"])
194            axes[1].errorbar(d.index, d["mean"], yerr=d["sem"], marker="o", ms=3, capsize=2, label=cond,
195                             color=CONDITION_COLORS[cond])
196        axes[1].set_xlabel("CSDS day")
197        axes[1].set_ylabel("fraction of session")
198        axes[1].set_title("Investigation across the divider")
199        axes[1].legend(fontsize=7)
200        self.figure(fig, "defeat_dynamics", "Evolution of defeat and social behavior",
201                    "Left: fraction of each defeat session in defeat-related states across CSDS days (mean ± SEM, n = 6). "
202                    "Right: non-physical interaction in baseline and post-defeat sessions.")
203
204    def _fig_sit(self, si: pd.DataFrame) -> None:
205        fig, axes = plt.subplots(1, 2, figsize=(8.5, 3.3))
206        for _, r in si.iterrows():
207            axes[0].plot([0, 1], [r.empty_engagement_s, r.CD1_engagement_s], marker="o",
208                         color="#b5523b" if r.phenotype_lab == "susceptible" else "#3b7d6e")
209            axes[0].text(1.05, r.CD1_engagement_s, str(int(r.subject)), fontsize=7, va="center")
210        axes[0].set_xticks([0, 1], ["empty", "CD1"])
211        axes[0].set_ylabel("engagement (s)")
212        axes[0].set_title("SIT engagement")
213        s = si.sort_values("SI_ratio_proxy")
214        axes[1].bar(s.subject.astype(str), s.SI_ratio_proxy,
215                    color=["#b5523b" if p == "susceptible" else "#3b7d6e" for p in s.phenotype_lab])
216        axes[1].axhline(1, color="k", ls="--", lw=0.8)
217        axes[1].set_ylabel("SI ratio proxy (CD1 / empty)")
218        axes[1].set_title("Social-interaction ratio proxy")
219        self.figure(fig, "sit_phenotype", "Social interaction test phenotype",
220                    "Left: time engaged with the enclosure when empty versus containing a CD-1 (per animal). Right: ratio "
221                    "proxy (dashed line = 1). Colour = lab phenotype: red = susceptible, green = resilient.")

  1"""S02 — Single-neuron activity across CSDS phases."""
  2
  3from __future__ import annotations
  4
  5import matplotlib.pyplot as plt
  6import numpy as np
  7import pandas as pd
  8
  9from ..data import CORE_CONDITIONS
 10from ..stats import animal_means, fdr_bh, fe_model, fmt_p, friedman, paired
 11from ..stats.compare import day_trend
 12from ..viz import day_lines, paired_conditions, star
 13from .base import Study
 14
 15
 16class Activity(Study):
 17    id = "s02"
 18    number = 2
 19    title = "Neural activity across CSDS phases"
 20    short = "Inferred firing, active-neuron fraction and behavior-locked population responses."
 21    question = ("Does the overall level of neural activity differ between baseline, defeat and post-defeat sensory "
 22                "contact, does it change across defeat days, and which behaviors drive population activity?")
 23    background = """
 24Chronic social defeat changes the excitability of stress-related circuits. CSDS increases activity of
 25prelimbic glutamatergic neurons and reduces PV interneuron activity (Nat Commun 2023). Resilient animals
 26show elevated dCA1 inferred spike rates (Biol Psychiatry Glob Open Sci 2025), and the pre-stress activity of NAc D1-MSNs predicts
 27later resilience (Muir et al. 2018). Changes in mean activity also matter methodologically: firing rate
 28influences pairwise correlations and therefore network measures (de la Rocha et al. 2007; Cohen & Kohn
 292011). Establishing activity effects is a prerequisite for interpreting the network studies that follow.
 30"""
 31    methods = """
 32**Signal.** CASCADE-inferred spike probability (expected spikes per frame) multiplied by the frame rate gives
 33an estimated event rate in Hz (Rupprecht et al. 2021). ΔF/F was analysed in parallel as a sensitivity check.
 34
 35**Session measures.** For each neuron and session: mean event rate; *active* if mean rate exceeded 0.1 Hz. Session
 36values are the mean over neurons and the fraction of active neurons.
 37
 38**Condition contrasts.** Primary contrasts use days 1–9, day 7 excluded, where BL, SD and SS share a day, an
 39animal and a recording length (≈2,000 frames). Only days on which an animal has all three conditions enter its
 40means (this drops BL and SD of day 9 for 6597, which has no SS that day), so every condition mean averages the
 41same days. Per-animal condition means were compared with a Friedman test followed by exact Wilcoxon
 42signed-rank tests (SD vs BL, SS vs BL, SS vs SD). A session-level OLS with animal × day block fixed effects
 43(`measure ~ condition + block`) gives within-day effect estimates; its session-level p-values are descriptive
 44because sessions of one animal are not independent. The day trend within each condition was estimated per
 45animal and the per-animal slopes tested against 0 (exact Wilcoxon). P-values were FDR-corrected across
 46measures × contrasts.
 47
 48**Behavior-locked population responses.** Within each session, the population rate (mean over neurons) was
 49z-scored, and its mean computed within each annotated behavior. Per-animal means (behaviors with ≥ 100 frames)
 50were tested against 0 with exact Wilcoxon tests, separately for defeat sessions and for baseline/post-defeat
 51sessions.
 52
 53**SIT.** Population rate during the CD-1 phase was compared with the empty-enclosure phase (paired Wilcoxon).
 54"""
 55    limitations = """
 56* Inferred spike rates depend on CASCADE model choice and noise level. Comparisons are within animal and
 57  within the same neurons, which limits but does not remove this dependence.
 58* With n = 6 the smallest attainable exact two-sided Wilcoxon p is 0.031. A null result is weak evidence of
 59  no effect.
 60* BL, SD and SS are always recorded in that order within a day. Condition contrasts are therefore confounded
 61  with position in the recording (photobleaching, within-day drift, time since handling).
 62* Slow fluorescence drift across the recording can mimic day trends. ΔF/F was computed on the concatenated
 63  recording (F₀ = 10th percentile in a 12-s window), so that window can cross an epoch boundary.
 64* The CASCADE model and noise level were not recorded. Spike probability is paired one-to-one with ΔF/F.
 65"""
 66    references = [
 67        "Rupprecht P. et al. (2021) CASCADE spike inference. Nat Neurosci 24:1324–1337.",
 68        "CSDS decreases claustrum glutamatergic and prelimbic PV interneuron activity and increases prelimbic glutamatergic activity (2023) Nat Commun. doi:10.1038/s41467-023-43636-x",
 69        "Diminished social memory and hippocampal correlates of social interactions in chronic social defeat stress susceptibility (2025) Biol Psychiatry Glob Open Sci. https://www.bpsgos.org/article/S2667-1743(25)00009-6/fulltext",
 70        "Muir J. et al. (2018) In vivo fiber photometry reveals signature of future stress susceptibility in nucleus accumbens. Neuropsychopharmacology 43:255–263.",
 71        "de la Rocha J. et al. (2007) Correlation between neural spike trains increases with firing rate. Nature 448:802–806.",
 72        "Cohen M.R. & Kohn A. (2011) Measuring and interpreting neuronal correlations. Nat Neurosci 14:811–819.",
 73    ]
 74
 75    def run(self) -> dict:
 76        ds = self.ds
 77        rows, beh_rows = [], []
 78        for s in ds.iter_sessions():
 79            x = ds.activity(s, "spikes")
 80            rate = x.mean(1) * ds.fs
 81            dff = ds.activity(s, "dff").mean(1)
 82            rows.append({"key": s.key, "subject": s.subject, "condition": s.condition, "day": s.day,
 83                         "mapping": s.mapping, "rate_hz": rate.mean(), "active_fraction": (rate > 0.1).mean(),
 84                         "dff": dff.mean()})
 85            pop = x.mean(0)
 86            z = (pop - pop.mean()) / (pop.std() + 1e-12)
 87            lab = ds.labels(s)
 88            for b in np.unique(lab):
 89                m = lab == b
 90                if b != "none" and m.sum() >= 1:
 91                    beh_rows.append({"subject": s.subject, "condition": s.condition, "behavior": b,
 92                                     "frames": int(m.sum()), "z_sum": float(z[m].sum())})
 93        df = pd.DataFrame(rows)
 94        self.ctx.cache["activity"] = df
 95
 96        core = df[df.condition.isin(CORE_CONDITIONS) & (df.day <= 9)]
 97        stat_rows, fe_rows = [], []
 98        pivots = {}
 99        for metric in ["rate_hz", "active_fraction", "dff"]:
100            pv = animal_means(core, metric, CORE_CONDITIONS)
101            pivots[metric] = pv
102            fr = friedman(pv)
103            for a, b in [("BL", "SD"), ("BL", "SS"), ("SD", "SS")]:
104                stat_rows.append({"metric": metric, **paired(pv, a, b), **fr})
105            fe_rows.append(fe_model(core, metric, CORE_CONDITIONS))
106        stats_df = pd.DataFrame(stat_rows)
107        stats_df["q_fdr"] = fdr_bh(stats_df.p)
108        fe_df = pd.concat(fe_rows)
109        trend = pd.DataFrame([day_trend(core, m, c) for m in ["rate_hz", "active_fraction"] for c in CORE_CONDITIONS])
110        trend["q_fdr"] = fdr_bh(trend.p)
111
112        # Behavior-locked
113        bl = pd.DataFrame(beh_rows)
114        bl["group"] = np.where(bl.condition == "SD", "SD", np.where(bl.condition.isin(["BL", "SS"]), "BL+SS", "other"))
115        agg = bl[bl.group != "other"].groupby(["group", "subject", "behavior"]).agg(frames=("frames", "sum"), z_sum=("z_sum", "sum"))
116        agg["z"] = agg.z_sum / agg.frames
117        agg = agg[agg.frames >= 100].reset_index()
118        from scipy import stats as st
119        brow = []
120        for (g, b), d in agg.groupby(["group", "behavior"]):
121            if len(d) >= 4:
122                p = st.wilcoxon(d.z, method="exact").pvalue if np.any(d.z != 0) else np.nan
123                brow.append({"group": g, "behavior": b, "n_animals": len(d), "mean_z": d.z.mean(), "p": p})
124        beh_stats = pd.DataFrame(brow)
125        beh_stats["q_fdr"] = fdr_bh(beh_stats.p)
126
127        sit = animal_means(df, "rate_hz", ["SIT_NoCD1", "SIT_CD1"])
128        sit_test = paired(sit, "SIT_NoCD1", "SIT_CD1") if len(sit) >= 3 else {}
129
130        self._fig_conditions(pivots, stats_df)
131        self._fig_days(core)
132        self._fig_behavior(agg, beh_stats)
133
134        self.table(stats_df[["metric", "contrast", "n", "median_diff", "dz", "p", "q_fdr", "friedman_chi2", "friedman_p"]],
135                   "condition_tests", "Condition contrasts (per-animal means, days 1–9, day 7 excluded)",
136                   "Exact Wilcoxon signed-rank on per-animal means; Friedman across BL/SD/SS; q = BH-FDR.")
137        self.table(fe_df, "fixed_effects", "Session-level OLS with animal × day fixed effects",
138                   "measure ~ condition + block (animal × day). p_session is descriptive (sessions are not independent).")
139        self.table(trend, "day_trends", "Day trends within condition",
140                   "slope_per_day = pooled within-animal slope; p = exact Wilcoxon of per-animal slopes.")
141        self.table(beh_stats.sort_values("p"), "behavior_locked", "Behavior-locked population rate",
142                   "Mean z-scored population rate during each behavior (per-animal mean tested against 0).")
143        self.table(df, "session_activity", "Session-level activity", show=15)
144
145        r = stats_df[(stats_df.metric == "rate_hz")].set_index("contrast")
146        sdbl = r.loc["SD vs BL"]
147        self.highlight("Firing during defeat", f"{sdbl.median_diff:+.3f} Hz (SD − BL)",
148                       f"Median paired change in inferred event rate, {sdbl.n_increase}/{sdbl.n} animals increased, "
149                       f"{fmt_p(sdbl.p)}; Friedman {fmt_p(sdbl.friedman_p)}.",
150                       "significant" if sdbl.q_fdr < 0.05 else "trend" if sdbl.p < 0.05 else "null")
151        top = beh_stats.sort_values("p").head(2)
152        for _, b in top.iterrows():
153            self.highlight(f"Population response: {b.behavior}", f"z = {b.mean_z:+.2f}",
154                           f"{b.group} sessions, {b.n_animals} animals, {fmt_p(b.p)} (q = {b.q_fdr:.3f}).",
155                           "significant" if b.q_fdr < 0.05 else "trend" if b.p < 0.05 else "null")
156        if sit_test:
157            self.highlight("SIT CD-1 vs empty", f"{sit_test['median_diff']:+.3f} Hz",
158                           f"Population rate change when the CD-1 is present, {fmt_p(sit_test['p'])}.",
159                           "trend" if sit_test["p"] < 0.05 else "null")
160
161        means = pivots["rate_hz"].mean()
162        results = f"""
163**Activity levels.** Mean inferred event rates on days 1–9, day 7 excluded were BL {means['BL']:.3f} Hz, SD {means['SD']:.3f} Hz and
164SS {means['SS']:.3f} Hz (animal means; Friedman {fmt_p(r.loc['SD vs BL'].friedman_p)}). For SD vs BL, the median
165paired difference was {sdbl.median_diff:+.3f} Hz ({sdbl.n_increase}/{sdbl.n} animals higher during defeat,
166{fmt_p(sdbl.p)}, q = {sdbl.q_fdr:.3f}). All contrasts for rate, active fraction and ΔF/F are in Table 1;
167fixed-effect estimates are in Table 2.
168
169**Across days.** Day trends are reported in Table 3 (Figure 2). {int((trend.q_fdr < 0.05).sum())} of {len(trend)}
170trends survived FDR correction.
171
172**Behavior-locked responses.** Figure 3 and Table 4 show the population z-score during each behavior.
173{int((beh_stats.q_fdr < 0.05).sum())} of {len(beh_stats)} behavior × context combinations differed from the session
174mean after FDR correction{'' if (beh_stats.q_fdr < 0.05).any() else '. At the population-average level, behavior did not consistently shift activity across animals. Single-neuron selectivity is examined in S03 and decoding in S08'}.
175"""
176        return {"results": results, "parameters": {"active_threshold_hz": 0.1, "min_behavior_frames": 100,
177                                                   "fs_hz": ds.fs}}
178
179    def _fig_conditions(self, pivots: dict, stats_df: pd.DataFrame) -> None:
180        fig, axes = plt.subplots(1, 3, figsize=(10, 3.3))
181        labels = {"rate_hz": "event rate (Hz)", "active_fraction": "fraction active (>0.1 Hz)", "dff": "mean ΔF/F"}
182        for ax, (m, pv) in zip(axes, pivots.items()):
183            st = stats_df[stats_df.metric == m].set_index("contrast")
184            stars = {("BL", "SD"): star(st.loc["SD vs BL", "p"]), ("BL", "SS"): star(st.loc["SS vs BL", "p"])}
185            paired_conditions(ax, pv, labels[m], labels[m].split(" (")[0], stars)
186        self.figure(fig, "activity_by_condition", "Activity by CSDS phase",
187                    "Per-animal means (grey lines) and group mean ± SEM for days 1–9, day 7 excluded. Stars: exact Wilcoxon vs BL "
188                    "(* p<0.05; n.s. not significant; with n = 6 the minimum p is 0.031).")
189
190    def _fig_days(self, core: pd.DataFrame) -> None:
191        fig, axes = plt.subplots(1, 2, figsize=(9, 3.2))
192        day_lines(axes[0], core, "rate_hz", CORE_CONDITIONS, "event rate (Hz)", "Event rate across days")
193        day_lines(axes[1], core, "active_fraction", CORE_CONDITIONS, "fraction active", "Active neurons across days")
194        self.figure(fig, "activity_by_day", "Activity across CSDS days",
195                    "Mean ± SEM across animals for each condition and day.")
196
197    def _fig_behavior(self, agg: pd.DataFrame, beh_stats: pd.DataFrame) -> None:
198        groups = ["SD", "BL+SS"]
199        fig, axes = plt.subplots(1, 2, figsize=(11, 3.4), gridspec_kw={"width_ratios": [1.4, 1]})
200        for ax, g in zip(axes, groups):
201            d = agg[agg.group == g]
202            s = beh_stats[beh_stats.group == g].set_index("behavior")
203            order = d.groupby("behavior").z.mean().sort_values().index
204            order = [b for b in order if b in s.index]
205            for i, b in enumerate(order):
206                v = d[d.behavior == b].z
207                ax.bar(i, v.mean(), color="#b5523b" if v.mean() > 0 else "#3b7d6e", alpha=0.75)
208                ax.scatter(np.full(len(v), i) + np.random.default_rng(i).uniform(-0.15, 0.15, len(v)), v, s=8, c="k")
209                ax.text(i, max(v.max(), 0) + 0.05, star(s.loc[b, "p"]), ha="center", fontsize=7)
210            ax.axhline(0, color="k", lw=0.6)
211            ax.set_xticks(range(len(order)), order, rotation=45, ha="right", fontsize=7)
212            ax.set_ylabel("population rate (z)")
213            ax.set_title(f"{g} sessions")
214        self.figure(fig, "behavior_locked_population", "Behavior-locked population activity",
215                    "Mean z-scored population rate during each annotated behavior; dots are animals. Stars: exact Wilcoxon vs 0.")