OSCR

The mPFC-reuniens-hippocampus pathway links brain circuitry and neural plasticity in antidepressant response.

Code ↔ Paper

4 matches between paragraphs of the paper and lines of its authors' code, computed by the harvester (lexical-v1). Click a colored paragraph or line to see its counterpart.

The 4 matches · 1 of them tie a paragraph to a whole file, not to given lines: a weak match, whose lines are not tinted
  1. [1] § Methods › In vivo Local Field Potential (LFP) recordings ↔ LFP_analysis.py, lines 126–153 · score 0.87 · 13–30 Hz, 4–12 Hz, high gamma, low gamma, downsampled, LFP
  2. [2] § Methods › Fiber photometry: data acquisition and analysis ↔ FiberPhotometryAnalysisTool/functions/analysis_data_loader.py, the whole file · a weak match · score 0.68 · find_peaks, peak amplitude, Episodes, GuPPy, Fiber photometry, duration
  3. [3] § Methods › Fiber photometry: data acquisition and analysis ↔ FiberPhotometryAnalysisTool/screens/plotting_screen.py, lines 16–106 · score 0.67 · find_peaks, peak amplitude, Episodes, GuPPy, Fiber photometry, duration
  4. [4] § Methods › Fiber photometry: data acquisition and analysis ↔ matlab_demodulation_script.m, lines 1–40 · score 0.57 · LabVIEW, fiber photometry, channel, signal

Paper

Loaded from Europe PMC by your browser, not stored by OSCR: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The authors' code

Python · 1,792 lines · 66 KB · MIT · 1 match

  2. """
  3. Source-data script for the quantitative panels of Supplementary Fig. 9.
  4. - Registration of channels and regions
  5. - Match behaviour and ephys files
  6. - NOE extraction and delta statistics
  7. - HAB running band-coherence extraction and statistics
  8. - HAB running power spectra
  10. """
  12. import os
  13. import re
  14. import glob
  15. import itertools
  16. from collections import defaultdict
  18. import matplotlib
  19. matplotlib.rcParams["pdf.fonttype"] = 42
  20. matplotlib.rcParams["svg.fonttype"] = "none"
  21. import matplotlib.pyplot as plt
  22. import numpy as np
  23. import pandas as pd
  24. from open_ephys.analysis import Session
  25. from scipy import signal
  26. from scipy.signal import hilbert, welch
  27. from scipy.signal.windows import dpss
  28. from scipy.stats import friedmanchisquare, wilcoxon
  31. np.random.seed(1)
  34. CHANNEL_TO_REGION = {
  35. "HCl": "HC",
  36. "HCr": "HC",
  37. "PFCl": "PFC",
  38. "PFCr": "PFC",
  39. "RE": "RE",
  40. }
  41. CHANNEL_TO_HEMI = {
  42. "HCl": "L",
  43. "HCr": "R",
  44. "PFCl": "L",
  45. "PFCr": "R",
  46. "RE": "M",
  47. }
  49. PAIR_POOL_MAP = {}
  50. for hemi_letter, suffix in [("L", "l"), ("R", "r")]:
  51. PAIR_POOL_MAP[tuple(sorted([f"HC{suffix}", f"PFC{suffix}"]))] = ("HC-PFC", hemi_letter)
  52. PAIR_POOL_MAP[tuple(sorted([f"HC{suffix}", "RE"]))] = ("HC-RE", hemi_letter)
  53. PAIR_POOL_MAP[tuple(sorted([f"PFC{suffix}", "RE"]))] = ("PFC-RE", hemi_letter)
  55. REGION_FIGURE_LABEL = {"HC": "vHIPP", "PFC": "mPFC", "RE": "RE"}
  56. PAIR_FIGURE_LABEL = {
  57. "HC-PFC": "mPFC-vHIPP",
  58. "HC-RE": "RE-vHIPP",
  59. "PFC-RE": "mPFC-RE",
  60. }
  61. PAC_DIRECTION_LABEL = {
  62. "HC-PFC": ("vHIPP", "mPFC"),
  63. "HC-RE": ("vHIPP", "RE"),
  64. "PFC-RE": ("mPFC", "RE"),
  65. }
  67. CONDITION_GROUPS = {
  68. "naive": ["naive1", "naive2"],
  69. "CDM": ["CDM_1", "CDM_2"],
  70. "C21": ["C21_1"],
  71. }
  72. CONDITION_LABEL = {"naive": "naive", "CDM": "post-CDM", "C21": "post-C21"}
  73. CONDITION_COLOR = {"naive": "#E69F00", "CDM": "#D55E00", "C21": "#56B4E9"}
  75. PSD_REGIONS = ["PFC", "RE", "HC"]
  76. COHERENCE_PAIRS = ["PFC-RE", "HC-RE", "HC-PFC"]
  78. BAR_PANELS_RUNNING = [
  79. {
  80. "panel": "m",
  81. "target": "HC-PFC",
  82. "metric": "coh",
  83. "band": "beta",
  84. "title": "Running - mPFC-vHIPP coherence β",
  85. },
  86. {
  87. "panel": "n",
  88. "target": "HC-RE",
  89. "metric": "coh",
  90. "band": "lowGamma",
  91. "title": "Running - RE-vHIPP coherence low-γ",
  92. },
  93. {
  94. "panel": "o",
  95. "target": "HC-PFC",
  96. "metric": "coh",
  97. "band": "lowGamma",
  98. "title": "Running - mPFC-vHIPP coherence low-γ",
  99. },
  100. ]
  102. BAR_PANELS_NOE = [
  103. {
  104. "panel": "p",
  105. "target": "PFC-RE",
  106. "metric": "pac_mi",
  107. "band": "theta_lowGamma",
  108. "title": "NOE - mPFC-RE θ-lowγ PAC-mi",
  109. },
  110. {
  111. "panel": "q",
  112. "target": "HC-RE",
  113. "metric": "coh",
  114. "band": "lowGamma",
  115. "title": "NOE - RE-vHIPP coherence lowγ",
  116. },
  117. {
  118. "panel": "r",
  119. "target": "HC-PFC",
  120. "metric": "pac_mi",
  121. "band": "theta_lowGamma",
  122. "title": "NOE - mPFC-vHIPP θ-lowγ PAC-mi",
  123. },
  124. ]
  126. def make_cfg():
  127. cfg = {
  128. "rootFoldEphys": r"D:\LTP_analysis\ProcessedData",
  129. "rootFoldVideo": r"D:\LTP_analysis\ProcessedData\VideoRec",
  130. "resultsDir": r"D:\LTP_analysis\results",
  131. "miceAvail": ["10084", "10085", "male_7pins", "male_9pins", "10418", "100134"],
  132. "batch": [1, 1, 1, 1, 2, 2],
  133. "conditions": ["naive1", "naive2", "CDM_1", "CDM_2", "C21_1", "C21_2"],
  134. "test": "NOE",
  135. "speedThresh_cm_s": 4.0,
  136. "px_to_cm": 1 / 8,
  137. "minEpochDur_s_NOE": 1.0,
  138. "minEpochDur_s_HAB": 0.5,
  139. "doDownsample": True,
  140. "fsTarget": 1000,
  141. "doNotch": False,
  142. "lineFreq": 50,
  143. "notchQ": 35,
  144. "freqRange": (1, 150),
  145. "bands": {
  146. "theta": (4, 12),
  147. "beta": (13, 30),
  148. "lowGamma": (30, 55),
  149. "highGamma": (65, 100),
  150. },
  151. }
  152. cfg["metaXlsx"] = os.path.join(cfg["rootFoldEphys"], "LFP_organization_bothbatches.xlsx")
  153. return cfg
  155. def pool_pair_name(a, b):
  156. return PAIR_POOL_MAP.get(tuple(sorted([a, b])), (None, None))
  159. def pool_channel_name(ch):
  160. return CHANNEL_TO_REGION.get(ch, ch), CHANNEL_TO_HEMI.get(ch, "?")
  163. def should_skip(mouse, cond, test):
  164. if mouse == "10418" and cond in ["C21_1", "C21_2"]:
  165. return True
  166. return False
  169. def _norm_str(s):
  170. return str(s).strip()
  173. def find_meta_row(meta, mouse, cond, test):
  174. m = meta.copy()
  175. m["mouse_id"] = m["mouse_id"].astype(str).str.strip()
  176. m["condition"] = m["condition"].astype(str).str.strip()
  177. m["test"] = m["test"].astype(str).str.strip()
  178. mask = (
  179. (m["mouse_id"] == _norm_str(mouse))
  180. & (m["condition"] == _norm_str(cond))
  181. & (m["test"] == _norm_str(test))
  182. )
  183. idx = np.where(mask.values)[0]
  184. if len(idx) == 0:
  185. return None, "no_metadata_match"
  186. return int(idx[0]), "ok"
  189. def parse_nose_xy(nose_series):
  190. x = np.zeros(len(nose_series), dtype=float)
  191. y = np.zeros(len(nose_series), dtype=float)
  192. for i, s in enumerate(nose_series):
  193. nums = [float(v) for v in re.findall(r"[-+]?\d*\.?\d+(?:[eE][-+]?\d+)?", str(s))]
  194. if len(nums) >= 2:
  195. x[i] = nums[0]
  196. y[i] = nums[1]
  197. return x, y
  200. def idx_to_epochs(time_vec, idx, min_dur):
  201. idx = np.asarray(idx, dtype=int)
  202. if idx.size == 0:
  203. return np.array([]), np.array([])
  204. d = np.diff(idx)
  205. split = np.where(d > 1)[0]
  206. starts = np.r_[0, split + 1]
  207. ends = np.r_[split, len(idx) - 1]
  208. ep_s = time_vec[idx[starts]]
  209. ep_e = time_vec[idx[ends]]
  210. keep = (ep_e - ep_s) >= min_dur
  211. return ep_s[keep], ep_e[keep]
  214. def resolve_video_dir(video_root, video_folder, batch_id):
  215. vf = str(video_folder)
  216. if os.path.isabs(vf) and os.path.exists(vf):
  217. return vf
  218. candidate = os.path.join(video_root, vf)
  219. if os.path.exists(candidate):
  220. return candidate
  221. base = os.path.basename(vf.rstrip("\\/"))
  222. candidate2 = os.path.join(video_root, f"Batch{batch_id}", base)
  223. if os.path.exists(candidate2):
  224. return candidate2
  225. return candidate
  228. def load_behavior_epochs(meta_row, test_name, cfg):
  229. video_dir = resolve_video_dir(cfg["rootFoldVideo"], meta_row["video_folder"], int(meta_row["batch_id"]))
  230. trial = str(meta_row["trial_number"])
  231. trial_suffix = trial[1:] if len(trial) > 1 else trial
  233. xls_name = f"processed_interpolated_T{trial_suffix}_V2.xlsx"
  234. xls_path = os.path.join(video_dir, xls_name)
  235. if not os.path.exists(xls_path):
  236. pattern = os.path.join(video_dir, f"processed_interpolated_T*{trial_suffix}*_V2.xlsx")
  237. cand = glob.glob(pattern)
  238. if cand:
  239. xls_path = cand[0]
  241. if not os.path.exists(xls_path):
  242. return None, None, None, f"behavior_file_missing:{xls_path}"
  244. b = pd.read_excel(xls_path)
  245. if "Timestamp_sec" not in b.columns:
  246. return None, None, None, "behavior_missing_Timestamp_sec"
  248. b.loc[b["Timestamp_sec"] <= 1, "Timestamp_sec"] = np.nan
  249. b = b.dropna(subset=["Timestamp_sec"]).copy()
  250. if len(b) == 0:
  251. return None, None, None, "behavior_empty_after_crop"
  253. b["Timestamp_sec"] = b["Timestamp_sec"] - b["Timestamp_sec"].iloc[0]
  255. if test_name == "NOE":
  256. if "region1_other" not in b.columns:
  257. return b, np.array([]), np.array([]), "missing_region1_other"
  258. idx = np.where(pd.to_numeric(b["region1_other"], errors="coerce").fillna(0).values > 0)[0]
  259. ep_s, ep_e = idx_to_epochs(b["Timestamp_sec"].values, idx, cfg["minEpochDur_s_NOE"])
  260. return b, ep_s, ep_e, "ok"
  262. stop_idx = np.where(b["Timestamp_sec"].values <= 600)[0]
  263. if len(stop_idx) == 0:
  264. return b, np.array([]), np.array([]), "hab_no_10min"
  265. stop10 = stop_idx[-1]
  267. if "nose" not in b.columns:
  268. return b, np.array([]), np.array([]), "missing_nose"
  270. noseX, noseY = parse_nose_xy(b["nose"].iloc[: stop10 + 1].values)
  271. dX = np.diff(noseX) * (40.0 / 680.0)
  272. dY = np.diff(noseY) * (30.0 / 481.0)
  273. dist = np.sqrt(dX ** 2 + dY ** 2)
  274. dt = np.diff(b["Timestamp_sec"].iloc[: stop10 + 1].values)
  275. dt[dt <= 0] = np.nan
  276. vel = dist / dt
  277. vel = np.nan_to_num(vel, nan=0.0, posinf=0.0, neginf=0.0)
  278. vel = np.r_[0, vel, 0]
  279. is_run = vel > cfg["speedThresh_cm_s"]
  281. idx = np.where(is_run)[0]
  282. ts = b["Timestamp_sec"].iloc[: len(vel)].values
  283. ep_s, ep_e = idx_to_epochs(ts, idx, cfg["minEpochDur_s_HAB"])
  284. return b, ep_s, ep_e, "ok"
  287. def preprocess_signal(x, fs, cfg):
  288. y = np.asarray(x, dtype=float)
  289. fs_out = float(fs)
  291. if cfg["doDownsample"] and fs_out > cfg["fsTarget"]:
  292. down = int(round(fs_out / cfg["fsTarget"]))
  293. y = signal.resample_poly(y, 1, down)
  294. fs_out = fs_out / down
  296. if cfg["doNotch"]:
  297. f0 = cfg["lineFreq"]
  298. q = cfg["notchQ"]
  299. if fs_out > (2 * (f0 + 2)):
  300. b, a = signal.iirnotch(w0=f0, Q=q, fs=fs_out)
  301. y = signal.filtfilt(b, a, y)
  303. y = signal.detrend(y, type="constant")
  304. return y, fs_out
  307. def reject_artifact_windows(x, centers, win_sec, fs, thresh_sd=4.0):
  308. if len(centers) == 0:
  309. return [], 0
  310. half = int(round(fs * win_sec / 2))
  311. n = len(x)
  312. amps = []
  313. valid_centers = []
  314. for c in centers:
  315. s, e = int(c - half), int(c + half)
  316. if s >= 0 and e <= n:
  317. amps.append(np.max(np.abs(x[s:e])))
  318. valid_centers.append(c)
  319. if len(amps) == 0:
  320. return [], 0
  321. amps = np.array(amps)
  322. med = np.median(amps)
  323. mad = np.median(np.abs(amps - med))
  324. robust_sd = 1.4826 * mad
  325. cutoff = med + thresh_sd * robust_sd
  326. keep = amps < cutoff
  327. clean = [c for c, k in zip(valid_centers, keep) if k]
  328. n_rejected = int(np.sum(~keep))
  329. return clean, n_rejected
  332. def compute_speed_on_lfp_time(bdf, t_lfp, px_to_cm=None):
  333. if "nose" not in bdf.columns:
  334. return np.zeros(len(t_lfp), dtype=float)
  336. px_to_cm_x = 40.0 / 680.0
  337. px_to_cm_y = 30.0 / 481.0
  339. bts = bdf["Timestamp_sec"].values.astype(float)
  340. nX, nY = parse_nose_xy(bdf["nose"].values)
  342. dX = np.diff(nX) * px_to_cm_x
  343. dY = np.diff(nY) * px_to_cm_y
  344. dt = np.diff(bts)
  345. dt[dt <= 0] = np.nan
  346. vel = np.sqrt(dX ** 2 + dY ** 2) / dt
  347. vel = np.nan_to_num(vel, nan=0.0, posinf=0.0, neginf=0.0)
  349. vel_ts = (bts[:-1] + bts[1:]) / 2.0
  350. speed_lfp = np.interp(t_lfp, vel_ts, vel, left=0.0, right=0.0)
  351. return speed_lfp
  354. def build_hab_speed_lfp_xy(bdf, t_lfp):
  355. px_to_cm_x = 40.0 / 680.0
  356. px_to_cm_y = 30.0 / 481.0
  358. bts = bdf["Timestamp_sec"].values.astype(float)
  359. nX, nY = parse_nose_xy(bdf["nose"].values)
  360. dX = np.diff(nX) * px_to_cm_x
  361. dY = np.diff(nY) * px_to_cm_y
  362. dt_b = np.diff(bts)
  363. dt_b[dt_b <= 0] = np.nan
  364. vel_b = np.sqrt(dX ** 2 + dY ** 2) / dt_b
  365. vel_b = np.nan_to_num(vel_b, nan=0.0, posinf=0.0, neginf=0.0)
  366. vel_ts = (bts[:-1] + bts[1:]) / 2.0
  367. return np.interp(t_lfp, vel_ts, vel_b, left=0.0, right=0.0)
  370. def tile_bout_to_centers(s_event, e_event, fs, n, win_sec, overlap=0.5):
  371. half = int(round(fs * win_sec / 2))
  372. if half < 2 or (2 * half + 1) >= n:
  373. return []
  375. bout_dur = float(e_event) - float(s_event)
  376. if bout_dur < win_sec:
  377. return []
  379. step = win_sec * (1.0 - overlap)
  380. centers = []
  381. t_center = float(s_event) + win_sec / 2.0
  383. while (t_center + win_sec / 2.0) <= (float(e_event) + 1e-9):
  384. t_lo = t_center - win_sec / 2.0
  385. t_hi = t_center + win_sec / 2.0
  387. if t_lo < float(s_event) - 1e-9 or t_hi > float(e_event) + 1e-9:
  388. t_center += step
  389. continue
  391. c = int(round(t_center * fs))
  392. if c >= half and c < (n - half - 1):
  393. centers.append(c)
  395. t_center += step
  397. return centers
  400. def _merge_close_epochs(ep_s, ep_e, max_gap_s, min_dur_s):
  401. if len(ep_s) == 0:
  402. return np.array([]), np.array([])
  403. order = np.argsort(ep_s)
  404. ep_s = np.asarray(ep_s[order], dtype=float)
  405. ep_e = np.asarray(ep_e[order], dtype=float)
  406. ms, me = [ep_s[0]], [ep_e[0]]
  407. for i in range(1, len(ep_s)):
  408. if ep_s[i] - me[-1] <= max_gap_s:
  409. me[-1] = max(me[-1], ep_e[i])
  410. else:
  411. ms.append(ep_s[i])
  412. me.append(ep_e[i])
  413. ms = np.array(ms)
  414. me = np.array(me)
  415. keep = (me - ms) >= min_dur_s
  416. return ms[keep], me[keep]
  419. def _interval_to_mask(time_vec, starts, ends):
  420. m = np.zeros(len(time_vec), dtype=bool)
  421. for s, e in zip(np.asarray(starts, dtype=float), np.asarray(ends, dtype=float)):
  422. m |= (time_vec >= s) & (time_vec <= e)
  423. return m
  426. def _object_mask_from_behavior(bdf):
  427. if "region1_other" in bdf.columns:
  428. return pd.to_numeric(bdf["region1_other"], errors="coerce").fillna(0).values > 0
  429. return np.zeros(len(bdf), dtype=bool)
  432. def build_speed_matched_baseline_mask(t_lfp, speed_lfp, event_mask, obj_mask_lfp, fs, win_sec, n_sd=1.0):
  433. n = len(t_lfp)
  434. half = int(round(fs * win_sec / 2))
  436. exp_speeds = speed_lfp[event_mask]
  437. if len(exp_speeds) < 5:
  438. return np.zeros(n, dtype=bool)
  440. mu = float(np.nanmean(exp_speeds))
  441. sd = float(np.nanstd(exp_speeds))
  442. if sd == 0:
  443. sd = max(0.5, 0.05 * mu)
  445. speed_lo = max(0.0, mu - n_sd * sd)
  446. speed_hi = mu + n_sd * sd
  448. speed_ok = (speed_lfp >= speed_lo) & (speed_lfp <= speed_hi)
  449. not_event = ~event_mask
  450. not_obj = ~obj_mask_lfp if obj_mask_lfp is not None else np.ones(n, dtype=bool)
  452. edge_ok = np.zeros(n, dtype=bool)
  453. edge_ok[half : n - half - 1] = True
  455. return speed_ok & not_event & not_obj & edge_ok
  458. def sample_baseline_centers(baseline_mask, n_needed, win_sec, fs, rng=None):
  459. if rng is None:
  460. rng = np.random.default_rng(seed=42)
  462. candidates = np.where(baseline_mask)[0]
  463. if len(candidates) == 0:
  464. return np.array([], dtype=int)
  466. min_sep = max(1, int(round(win_sec * fs)))
  467. order = rng.permutation(len(candidates))
  468. chosen = []
  469. for idx in order:
  470. c = int(candidates[idx])
  471. if all(abs(c - p) >= min_sep for p in chosen):
  472. chosen.append(c)
  473. if len(chosen) >= n_needed:
  474. break
  476. return np.array(sorted(chosen), dtype=int)
  479. def butter_bandpass(fs, band, order=3):
  480. lo = max(0.5, float(band[0]))
  481. hi = min(float(band[1]), fs / 2 - 1)
  482. if hi <= lo:
  483. return None, None
  484. b, a = signal.butter(order, [lo, hi], btype="bandpass", fs=fs)
  485. return b, a
  488. def window_slices(centers, half_win, n):
  489. out = []
  490. for c in centers:
  491. s = int(c - half_win)
  492. e = int(c + half_win + 1)
  493. if s >= 0 and e <= n:
  494. out.append((s, e))
  495. return out
  498. def coherence_band_mean(x, y, fs, centers, win_sec, band, NW=2):
  499. n = len(x)
  500. half = int(round(fs * win_sec / 2))
  501. vals = []
  502. for s, e in window_slices(centers, half, n):
  503. segx, segy = x[s:e], y[s:e]
  504. N = len(segx)
  505. if N < 32:
  506. continue
  507. K = max(1, 2 * NW - 1)
  508. tapers = dpss(N, NW, K)
  509. nfft = int(2 ** np.ceil(np.log2(N)))
  510. freqs = np.fft.rfftfreq(nfft, 1.0 / fs)
  511. Sxx = np.zeros(len(freqs))
  512. Syy = np.zeros(len(freqs))
  513. Sxy = np.zeros(len(freqs), dtype=complex)
  514. for taper in tapers:
  515. fx = np.fft.rfft(segx * taper, n=nfft)
  516. fy = np.fft.rfft(segy * taper, n=nfft)
  517. Sxx += np.abs(fx) ** 2
  518. Syy += np.abs(fy) ** 2
  519. Sxy += fx * np.conj(fy)
  520. Sxx /= K
  521. Syy /= K
  522. Sxy /= K
  523. denom = Sxx * Syy
  524. coh = np.where(denom > 0, np.abs(Sxy) ** 2 / denom, 0.0)
  525. keep = (freqs >= band[0]) & (freqs <= band[1])
  526. if np.any(keep):
  527. vals.append(float(np.nanmean(coh[keep])))
  528. return float(np.nanmean(vals)) if vals else np.nan
  531. def pac_mi(phase_sig, amp_sig, fs, centers, win_sec, phase_band, amp_band, n_bins=18):
  532. bp1 = butter_bandpass(fs, phase_band)
  533. bp2 = butter_bandpass(fs, amp_band)
  534. if bp1[0] is None or bp2[0] is None:
  535. return np.nan
  536. n = len(phase_sig)
  537. half = int(round(fs * win_sec / 2))
  538. bins = np.linspace(-np.pi, np.pi, n_bins + 1)
  539. vals = []
  540. for s, e in window_slices(centers, half, n):
  541. ph = np.angle(hilbert(signal.filtfilt(bp1[0], bp1[1], phase_sig[s:e])))
  542. amp = np.abs(hilbert(signal.filtfilt(bp2[0], bp2[1], amp_sig[s:e])))
  543. mean_amp = np.array(
  544. [
  545. np.mean(amp[(ph >= bins[bi]) & (ph < bins[bi + 1])])
  546. if np.any((ph >= bins[bi]) & (ph < bins[bi + 1]))
  547. else 0
  548. for bi in range(n_bins)
  549. ]
  550. )
  551. pa = mean_amp / (np.sum(mean_amp) + 1e-12)
  552. H = -np.sum(pa * np.log(pa + 1e-12))
  553. Hmax = np.log(n_bins)
  554. vals.append((Hmax - H) / (Hmax + 1e-12))
  555. return float(np.nanmean(vals)) if vals else np.nan
  558. def fdr_bh(pvals):
  559. p = np.array(pvals, dtype=float)
  560. n = len(p)
  561. if n == 0:
  562. return p
  563. order = np.argsort(p)
  564. ranked = p[order]
  565. q = np.empty(n, dtype=float)
  566. prev = 1.0
  567. for i in range(n - 1, -1, -1):
  568. prev = min(prev, ranked[i] * n / (i + 1))
  569. q[i] = prev
  570. q_out = np.empty(n, dtype=float)
  571. q_out[order] = np.minimum(q, 1.0)
  572. return q_out
  575. def load_channel_table(cfg):
  576. channel_data = pd.read_excel(cfg["metaXlsx"], sheet_name="electrodes")
  577. channel_data["animal"] = channel_data["animal"].astype(str).str.strip()
  578. return channel_data
  581. def load_meta_table(cfg, batch_id):
  582. return pd.read_excel(cfg["metaXlsx"], sheet_name=f"batch{batch_id}")
  585. def load_ephys_channels(cfg, meta_row, batch_id, channel_data, mouse):
  586. this_data_dir = os.path.join(
  587. cfg["rootFoldEphys"], "EphysRec", f"Batch{batch_id}", str(meta_row["ephy_folder"])
  588. )
  589. session_obj = Session(this_data_dir)
  590. rec = session_obj.recordnodes[0].recordings[0]
  591. stream_name = list(rec.continuous.keys())[0]
  592. cont = rec.continuous[stream_name]
  593. fs_raw = float(cont.metadata.sample_rate)
  594. t = np.asarray(cont.timestamps)
  596. eventline = rec.events.copy()
  597. st = pd.to_numeric(eventline["state"], errors="coerce").fillna(0).astype(int)
  598. smp = pd.to_numeric(eventline["sample_number"], errors="coerce")
  600. if np.sum(st == 1) < 1 or np.sum(st == 0) < 1:
  601. raise RuntimeError("bad_ttl_pattern")
  603. start_video = max(0, int(smp[st == 1].iloc[0] - np.min(t) * fs_raw + 1))
  604. stop_video = min(len(t) - 1, int(smp[st == 0].iloc[0] - np.min(t) * fs_raw + 1))
  605. if stop_video <= start_video:
  606. raise RuntimeError("empty_ttl_crop")
  608. ani_idx = np.where(channel_data["animal"].str.lower().values == str(mouse).strip().lower())[0]
  609. if len(ani_idx) == 0:
  610. raise RuntimeError("animal_not_in_electrode_sheet")
  611. ai = int(ani_idx[0])
  613. samples = np.asarray(cont.samples)
  614. time_first = samples.shape[0] == len(t)
  615. ref = (
  616. (samples[start_video:stop_video, 23] if time_first else samples[23, start_video:stop_video])
  617. if (batch_id == 1 and (samples.shape[1 if time_first else 0] > 23))
  618. else 0
  619. )
  621. channels = {}
  622. for name in ["HCl", "HCr", "RE", "PFCl", "PFCr"]:
  623. ch_id_raw = channel_data.loc[ai, name] if name in channel_data.columns else -1
  624. m_id = re.findall(r"[-+]?\d+", str(ch_id_raw))
  625. ch_id = int(m_id[0]) if m_id else -1
  626. if ch_id < 0:
  627. channels[name] = np.array([])
  628. continue
  629. ch_idx = ch_id - 1
  630. if time_first:
  631. channels[name] = (
  632. np.asarray(samples[start_video:stop_video, ch_idx], dtype=float) - ref
  633. if 0 <= ch_idx < samples.shape[1]
  634. else np.array([])
  635. )
  636. else:
  637. channels[name] = (
  638. np.asarray(samples[ch_idx, start_video:stop_video], dtype=float) - ref
  639. if 0 <= ch_idx < samples.shape[0]
  640. else np.array([])
  641. )
  642. return channels, fs_raw
  645. def preprocess_channels(channels, fs_raw, cfg, max_time_s=None):
  646. valid_names = [k for k, v in channels.items() if isinstance(v, np.ndarray) and v.size > 0]
  647. proc = {}
  648. fs = None
  649. for name in valid_names:
  650. proc[name], fs_this = preprocess_signal(channels[name], fs_raw, cfg)
  651. if fs is None:
  652. fs = fs_this
  653. if not proc:
  654. return {}, np.nan, np.array([])
  655. n = min(len(v) for v in proc.values())
  656. if max_time_s is not None:
  657. n = min(n, int(max_time_s * fs))
  658. for name in proc:
  659. proc[name] = proc[name][:n]
  660. t_lfp = np.arange(n) / fs
  661. return proc, fs, t_lfp
  664. def compute_group_summary(subject_df, value_col):
  665. rows = []
  666. group_cols = [c for c in subject_df.columns if c not in [value_col, "subject"]]
  667. for keys, g in subject_df.groupby(group_cols, sort=False):
  668. if not isinstance(keys, tuple):
  669. keys = (keys,)
  670. row = dict(zip(group_cols, keys))
  671. vals = g[value_col].values.astype(float)
  672. n = len(vals)
  673. row["n_subjects"] = n
  674. row["mean"] = float(np.nanmean(vals))
  675. row["sem"] = float(np.nanstd(vals, ddof=1) / np.sqrt(n)) if n > 1 else np.nan
  676. rows.append(row)
  677. return pd.DataFrame(rows)
  680. def save_dataframe(df, path):
  681. os.makedirs(os.path.dirname(path), exist_ok=True)
  682. df.to_csv(path, index=False)
  685. def compute_triplet_stats(
  686. results_long,
  687. out_dir,
  688. prefix,
  689. epoch_keep,
  690. domain_keep,
  691. metric_keep,
  692. bands_keep,
  693. win_sec,
  694. pool_hemispheres=False,
  695. require_both_hemis=False,
  696. c21_conds=None,
  697. normalize_method="zscore",
  698. apply_fdr=False,
  699. min_n_event_windows=0,
  700. ):
  701. if c21_conds is None:
  702. c21_conds = ["C21_1"]
  704. R0 = results_long.copy()
  706. if len(domain_keep) > 0:
  707. R0 = R0[R0["domain"].isin(domain_keep)].copy()
  708. if len(epoch_keep) > 0:
  709. R0 = R0[R0["epoch"].isin(epoch_keep)].copy()
  710. if len(metric_keep) > 0:
  711. R0 = R0[R0["metric"].isin(metric_keep)].copy()
  712. if bands_keep != "ALL":
  713. R0 = R0[R0["band"].isin(bands_keep)].copy()
  714. if win_sec > 0:
  715. R0 = R0[np.isclose(R0["win_sec"].astype(float), float(win_sec))].copy()
  716. if min_n_event_windows > 0 and "n_event_windows" in R0.columns:
  717. R0 = R0[
  718. pd.to_numeric(R0["n_event_windows"], errors="coerce").fillna(0).astype(int) >= int(min_n_event_windows)
  719. ].copy()
  721. if len(R0) == 0:
  722. raise RuntimeError(f"{prefix}: no rows left after filtering.")
  724. base_keys = ["win_sec", "epoch", "domain", "target", "metric", "band"]
  726. if pool_hemispheres:
  727. side = R0.groupby(["mouse", "condition"] + base_keys, as_index=False).agg(
  728. value=("value", "mean"),
  729. n_hemi=("hemisphere", "nunique"),
  730. )
  731. if require_both_hemis:
  732. is_lr_target = side["target"].isin(["HC", "PFC"])
  733. side = side[(~is_lr_target) | (side["n_hemi"] == 2)].copy()
  734. side["subject"] = side["mouse"].astype(str)
  735. else:
  736. side = R0.groupby(["mouse", "hemisphere", "condition"] + base_keys, as_index=False).agg(
  737. value=("value", "mean")
  738. )
  739. side["subject"] = side["mouse"].astype(str) + "|" + side["hemisphere"].astype(str)
  741. merge_keys = ["subject"] + base_keys
  743. def group_cond(df, conds, label):
  744. x = df[df["condition"].isin(conds)].copy()
  745. return x.groupby(merge_keys, as_index=False)["value"].mean().rename(columns={"value": label})
  747. naive = group_cond(side, CONDITION_GROUPS["naive"], "naive")
  748. cdm = group_cond(side, CONDITION_GROUPS["CDM"], "cdm")
  749. c21 = group_cond(side, c21_conds, "c21")
  751. merged_raw = naive.merge(cdm, on=merge_keys, how="inner").merge(c21, on=merge_keys, how="inner")
  752. if len(merged_raw) == 0:
  753. raise RuntimeError(f"{prefix}: no matched Naive/CDM/C21 rows after merging.")
  755. merged_norm = merged_raw.copy()
  756. if normalize_method == "zscore":
  757. for feat_keys, g in merged_norm.groupby(base_keys, sort=False):
  758. feat_vals = feat_keys if isinstance(feat_keys, tuple) else (feat_keys,)
  759. for subj in g["subject"].unique():
  760. mask = merged_norm["subject"] == subj
  761. for bk, bv in zip(base_keys, feat_vals):
  762. mask = mask & (merged_norm[bk] == bv)
  763. vals = merged_norm.loc[mask, ["naive", "cdm", "c21"]].values.flatten().astype(float)
  764. mu = np.nanmean(vals)
  765. sd = np.nanstd(vals)
  766. if sd > 0:
  767. merged_norm.loc[mask, "naive"] = (merged_norm.loc[mask, "naive"] - mu) / sd
  768. merged_norm.loc[mask, "cdm"] = (merged_norm.loc[mask, "cdm"] - mu) / sd
  769. merged_norm.loc[mask, "c21"] = (merged_norm.loc[mask, "c21"] - mu) / sd
  770. norm_label = "z-score"
  771. elif normalize_method == "none":
  772. norm_label = "raw"
  773. else:
  774. raise RuntimeError(f"Unsupported normalization: {normalize_method}")
  776. rows = []
  777. for keys, g in merged_norm.groupby(base_keys, sort=False):
  778. if not isinstance(keys, tuple):
  779. keys = (keys,)
  780. win_sec_i, epoch_i, domain_i, target_i, metric_i, band_i = keys
  781. vn = g["naive"].values.astype(float)
  782. vc = g["cdm"].values.astype(float)
  783. v2 = g["c21"].values.astype(float)
  784. n_subjects = len(g)
  785. if n_subjects < 3:
  786. continue
  788. try:
  789. chi2, p_friedman = friedmanchisquare(vn, vc, v2)
  790. except Exception:
  791. chi2, p_friedman = np.nan, np.nan
  793. try:
  794. _, p_naive_cdm = wilcoxon(vn, vc, alternative="two-sided")
  795. except Exception:
  796. p_naive_cdm = np.nan
  798. try:
  799. cdm_direction = np.nanmean(vc) - np.nanmean(vn)
  800. if cdm_direction > 0:
  801. _, p_cdm_c21 = wilcoxon(vc, v2, alternative="greater")
  802. else:
  803. _, p_cdm_c21 = wilcoxon(vc, v2, alternative="less")
  804. except Exception:
  805. p_cdm_c21 = np.nan
  807. try:
  808. _, p_naive_c21 = wilcoxon(vn, v2, alternative="two-sided")
  809. except Exception:
  810. p_naive_c21 = np.nan
  812. d1 = vc - vn
  813. d2 = v2 - vc
  814. rows.append(
  815. {
  816. "win_sec": float(win_sec_i),
  817. "epoch": epoch_i,
  818. "domain": domain_i,
  819. "target": target_i,
  820. "metric": metric_i,
  821. "band": band_i,
  822. "n_subjects": int(n_subjects),
  823. "mean_naive": float(np.nanmean(vn)),
  824. "mean_cdm": float(np.nanmean(vc)),
  825. "mean_c21": float(np.nanmean(v2)),
  826. "chi2_friedman": chi2,
  827. "p_friedman": p_friedman,
  828. "p_naive_cdm": p_naive_cdm,
  829. "p_cdm_c21": p_cdm_c21,
  830. "p_naive_c21": p_naive_c21,
  831. "restores_frac": float(np.mean(np.sign(d1) == -np.sign(d2))),
  832. "opposite_sign_group": bool(np.sign(np.nanmean(d1)) == -np.sign(np.nanmean(d2))),
  833. }
  834. )
  836. stats_all = pd.DataFrame(rows)
  837. if len(stats_all) == 0:
  838. raise RuntimeError(f"{prefix}: no analyzable groups.")
  840. if apply_fdr:
  841. for pcol in ["p_friedman", "p_naive_cdm", "p_cdm_c21", "p_naive_c21"]:
  842. valid = stats_all[pcol].notna()
  843. if valid.sum() > 0:
  844. stats_all.loc[valid, f"{pcol}_fdr"] = fdr_bh(stats_all.loc[valid, pcol].values)
  846. stats_all = stats_all.sort_values("p_friedman", na_position="last").reset_index(drop=True)
  848. suffix_hemi = "pooledHemi" if pool_hemispheres else "splitHemi"
  849. suffix_norm = normalize_method
  851. save_dataframe(stats_all, os.path.join(out_dir, f"{prefix}_stats_{suffix_hemi}_{suffix_norm}.csv"))
  852. save_dataframe(merged_raw, os.path.join(out_dir, f"{prefix}_plot_data_raw_{suffix_hemi}_{suffix_norm}.csv"))
  853. save_dataframe(merged_norm, os.path.join(out_dir, f"{prefix}_plot_data_norm_{suffix_hemi}_{suffix_norm}.csv"))
  855. return {
  856. "stats": stats_all,
  857. "plot_data_raw": merged_raw,
  858. "plot_data_norm": merged_norm,
  859. "norm_label": norm_label,
  860. "suffix_hemi": suffix_hemi,
  861. "suffix_norm": suffix_norm,
  862. }
  865. def run_noe_extraction(cfg, out_dir):
  866. channel_data = load_channel_table(cfg)
  867. win_sec = 0.5
  868. win_overlap = 0.5
  869. min_epoch_s = 0.5
  870. max_gap_s = 0.75
  871. min_exploration_s = 3.0
  872. min_bouts_per_session = 2
  873. baseline_n_sd = 1.0
  874. reject_baseline_if_object = True
  875. random_seed = 42
  876. artifact_thresh_sd = 4.0
  878. rng = np.random.default_rng(seed=random_seed)
  879. conn_bands = ["theta", "beta", "lowGamma"]
  880. pac_pairs = [("theta", "lowGamma")]
  882. all_rows = []
  883. logs = []
  884. bout_qc_rows = []
  886. orig_min_epoch = cfg["minEpochDur_s_NOE"]
  887. cfg["minEpochDur_s_NOE"] = 0.01
  889. for mi, mouse in enumerate(cfg["miceAvail"]):
  890. batch_id = int(cfg["batch"][mi])
  891. meta = load_meta_table(cfg, batch_id)
  893. for cond in cfg["conditions"]:
  894. test = "NOE"
  895. session_id = f"{mouse}|{cond}|{test}"
  897. if should_skip(mouse, cond, test):
  898. logs.append({"session_id": session_id, "mouse": mouse, "batch": batch_id, "condition": cond, "test": test, "status": "skipped_rule"})
  899. continue
  901. row_id, msg = find_meta_row(meta, mouse, cond, test)
  902. if row_id is None:
  903. logs.append({"session_id": session_id, "mouse": mouse, "batch": batch_id, "condition": cond, "test": test, "status": msg})
  904. continue
  906. md = meta.loc[row_id].copy()
  907. md["batch_id"] = batch_id
  909. try:
  910. channels, fs_raw = load_ephys_channels(cfg, md, batch_id, channel_data, mouse)
  911. except Exception as e:
  912. logs.append({"session_id": session_id, "mouse": mouse, "batch": batch_id, "condition": cond, "test": test, "status": f"load_ephys_failed:{e}"})
  913. continue
  915. bdf, ep_s, ep_e, bmsg = load_behavior_epochs(md, test, cfg)
  916. if bdf is None:
  917. logs.append({"session_id": session_id, "mouse": mouse, "batch": batch_id, "condition": cond, "test": test, "status": bmsg})
  918. continue
  920. n_raw = len(ep_s)
  921. ep_s, ep_e = _merge_close_epochs(ep_s, ep_e, max_gap_s, min_epoch_s)
  922. n_after = len(ep_s)
  923. total_exploration_s = float(np.sum(ep_e - ep_s)) if n_after > 0 else 0.0
  925. if n_after < min_bouts_per_session or total_exploration_s < min_exploration_s:
  926. reason = f"too_few_bouts:{n_after}" if n_after < min_bouts_per_session else f"too_little_exploration:{total_exploration_s:.1f}s"
  927. logs.append(
  928. {
  929. "session_id": session_id,
  930. "mouse": mouse,
  931. "batch": batch_id,
  932. "condition": cond,
  933. "test": test,
  934. "status": reason,
  935. "n_raw_bouts": n_raw,
  936. "n_merged_bouts": n_after,
  937. "total_exploration_s": total_exploration_s,
  938. }
  939. )
  940. continue
  942. proc, fs, t_lfp = preprocess_channels(channels, fs_raw, cfg)
  943. if not proc or len(t_lfp) < int(fs):
  944. logs.append({"session_id": session_id, "mouse": mouse, "batch": batch_id, "condition": cond, "test": test, "status": "too_short"})
  945. continue
  947. n = len(t_lfp)
  948. event_mask = _interval_to_mask(t_lfp, ep_s, ep_e)
  949. speed_lfp = compute_speed_on_lfp_time(bdf, t_lfp, cfg["px_to_cm"])
  951. obj_mask_lfp = None
  952. if reject_baseline_if_object:
  953. bts = bdf["Timestamp_sec"].values.astype(float)
  954. obj_beh = _object_mask_from_behavior(bdf).astype(float)
  955. obj_mask_lfp = (
  956. np.interp(t_lfp, bts, obj_beh, left=0.0, right=0.0) > 0.5
  957. if len(bts) > 1
  958. else np.zeros(n, dtype=bool)
  959. )
  961. event_centers_list = []
  962. per_bout_qc = []
  963. for s_evt, e_evt in zip(ep_s.astype(float), ep_e.astype(float)):
  964. centers_this_bout = tile_bout_to_centers(s_evt, e_evt, fs, n, win_sec=win_sec, overlap=win_overlap)
  965. event_centers_list.extend(centers_this_bout)
  966. per_bout_qc.append(
  967. {
  968. "session_id": session_id,
  969. "mouse": mouse,
  970. "condition": cond,
  971. "bout_start_s": float(s_evt),
  972. "bout_end_s": float(e_evt),
  973. "bout_duration_s": float(e_evt - s_evt),
  974. "n_windows_tiled": len(centers_this_bout),
  975. }
  976. )
  978. event_centers = np.array(sorted(event_centers_list), dtype=int)
  979. if len(event_centers) == 0:
  980. logs.append({"session_id": session_id, "mouse": mouse, "batch": batch_id, "condition": cond, "test": test, "status": "no_event_windows"})
  981. continue
  983. artifact_keep = np.ones(len(event_centers), dtype=bool)
  984. for _, x_ch in proc.items():
  985. clean_centers_ch, _ = reject_artifact_windows(x_ch, list(event_centers), win_sec, fs, thresh_sd=artifact_thresh_sd)
  986. clean_set = set(clean_centers_ch)
  987. for ci, c in enumerate(event_centers):
  988. if c not in clean_set:
  989. artifact_keep[ci] = False
  990. event_centers = event_centers[artifact_keep]
  991. n_artifact_rejected = int(np.sum(~artifact_keep))
  992. n_event_windows = len(event_centers)
  993. if n_event_windows == 0:
  994. logs.append({"session_id": session_id, "mouse": mouse, "batch": batch_id, "condition": cond, "test": test, "status": "all_event_windows_rejected"})
  995. continue
  997. baseline_mask = build_speed_matched_baseline_mask(
  998. t_lfp,
  999. speed_lfp,
  1000. event_mask,
  1001. obj_mask_lfp,
  1002. fs=fs,
  1003. win_sec=win_sec,
  1004. n_sd=baseline_n_sd,
  1005. )
  1006. n_pool = int(np.sum(baseline_mask))
  1007. baseline_centers = sample_baseline_centers(baseline_mask, n_needed=n_event_windows, win_sec=win_sec, fs=fs, rng=rng)
  1009. baseline_keep = np.ones(len(baseline_centers), dtype=bool)
  1010. for _, x_ch in proc.items():
  1011. clean_centers_ch, _ = reject_artifact_windows(x_ch, list(baseline_centers), win_sec, fs, thresh_sd=artifact_thresh_sd)
  1012. clean_set = set(clean_centers_ch)
  1013. for ci, c in enumerate(baseline_centers):
  1014. if c not in clean_set:
  1015. baseline_keep[ci] = False
  1016. baseline_centers = baseline_centers[baseline_keep]
  1017. n_baseline_windows = len(baseline_centers)
  1018. n_valid_pairs = min(n_event_windows, n_baseline_windows)
  1020. bout_qc_rows.extend(per_bout_qc)
  1021. epoch_centers = {"event": event_centers, "baseline": baseline_centers}
  1022. epoch_metrics = {}
  1024. for epoch_name, centers in epoch_centers.items():
  1025. m = {}
  1026. if len(centers) == 0:
  1027. epoch_metrics[epoch_name] = m
  1028. continue
  1030. for a, b in itertools.combinations(sorted(proc.keys()), 2):
  1031. canonical, pair_hemi = pool_pair_name(a, b)
  1032. if canonical is None:
  1033. continue
  1034. xa, xb = proc[a], proc[b]
  1036. for bname in conn_bands:
  1037. val = coherence_band_mean(xa, xb, fs, centers, win_sec, cfg["bands"][bname])
  1038. m[("pair", canonical, pair_hemi, "coh", bname)] = val
  1040. for phase_band_name, amp_band_name in pac_pairs:
  1041. mi_val = pac_mi(
  1042. xa,
  1043. xb,
  1044. fs,
  1045. centers,
  1046. win_sec,
  1047. cfg["bands"][phase_band_name],
  1048. cfg["bands"][amp_band_name],
  1049. )
  1050. m[("pair", canonical, pair_hemi, "pac_mi", f"{phase_band_name}_{amp_band_name}")] = mi_val
  1052. epoch_metrics[epoch_name] = m
  1054. for key, val in m.items():
  1055. domain, target, hemi, metric, band = key
  1056. phase_region, amp_region = PAC_DIRECTION_LABEL.get(target, ("", ""))
  1057. all_rows.append(
  1058. {
  1059. "session_id": session_id,
  1060. "mouse": mouse,
  1061. "batch": batch_id,
  1062. "condition": cond,
  1063. "test": test,
  1064. "epoch": epoch_name,
  1065. "win_sec": win_sec,
  1066. "domain": domain,
  1067. "target": target,
  1068. "target_figure_label": PAIR_FIGURE_LABEL.get(target, target),
  1069. "hemisphere": hemi,
  1070. "metric": metric,
  1071. "band": band,
  1072. "phase_region": phase_region if metric == "pac_mi" else "",
  1073. "amplitude_region": amp_region if metric == "pac_mi" else "",
  1074. "value": float(val) if np.isfinite(val) else np.nan,
  1075. "notes": "ok",
  1076. "n_raw_bouts": int(n_raw),
  1077. "n_merged_bouts": int(n_after),
  1078. "total_exploration_s": float(total_exploration_s),
  1079. "n_event_windows": int(n_event_windows),
  1080. "n_baseline_windows": int(n_baseline_windows),
  1081. "n_valid_pairs": int(n_valid_pairs),
  1082. "baseline_pool_size": int(n_pool),
  1083. "n_artifact_rejected": int(n_artifact_rejected),
  1084. }
  1085. )
  1087. me = epoch_metrics.get("event", {})
  1088. mb = epoch_metrics.get("baseline", {})
  1089. for kk in set(me.keys()) & set(mb.keys()):
  1090. ve, vb = me[kk], mb[kk]
  1091. vd = ve - vb if np.isfinite(ve) and np.isfinite(vb) else np.nan
  1092. domain, target, hemi, metric, band = kk
  1093. phase_region, amp_region = PAC_DIRECTION_LABEL.get(target, ("", ""))
  1094. all_rows.append(
  1095. {
  1096. "session_id": session_id,
  1097. "mouse": mouse,
  1098. "batch": batch_id,
  1099. "condition": cond,
  1100. "test": test,
  1101. "epoch": "delta",
  1102. "win_sec": win_sec,
  1103. "domain": domain,
  1104. "target": target,
  1105. "target_figure_label": PAIR_FIGURE_LABEL.get(target, target),
  1106. "hemisphere": hemi,
  1107. "metric": metric,
  1108. "band": band,
  1109. "phase_region": phase_region if metric == "pac_mi" else "",
  1110. "amplitude_region": amp_region if metric == "pac_mi" else "",
  1111. "value": float(vd) if np.isfinite(vd) else np.nan,
  1112. "notes": "event_minus_speedMatchedBaseline",
  1113. "n_raw_bouts": int(n_raw),
  1114. "n_merged_bouts": int(n_after),
  1115. "total_exploration_s": float(total_exploration_s),
  1116. "n_event_windows": int(n_event_windows),
  1117. "n_baseline_windows": int(n_baseline_windows),
  1118. "n_valid_pairs": int(n_valid_pairs),
  1119. "baseline_pool_size": int(n_pool),
  1120. "n_artifact_rejected": int(n_artifact_rejected),
  1121. }
  1122. )
  1124. logs.append(
  1125. {
  1126. "session_id": session_id,
  1127. "mouse": mouse,
  1128. "batch": batch_id,
  1129. "condition": cond,
  1130. "test": test,
  1131. "status": "ok",
  1132. "n_raw_bouts": int(n_raw),
  1133. "n_merged_bouts": int(n_after),
  1134. "total_exploration_s": float(total_exploration_s),
  1135. "n_event_windows": int(n_event_windows),
  1136. "n_baseline_windows": int(n_baseline_windows),
  1137. "n_valid_pairs": int(n_valid_pairs),
  1138. "baseline_pool_size": int(n_pool),
  1139. "n_artifact_rejected": int(n_artifact_rejected),
  1140. }
  1141. )
  1143. cfg["minEpochDur_s_NOE"] = orig_min_epoch
  1145. results_long = pd.DataFrame(all_rows)
  1146. session_log = pd.DataFrame(logs)
  1147. bout_qc = pd.DataFrame(bout_qc_rows)
  1149. save_dataframe(results_long, os.path.join(out_dir, "NOE_results_long.csv"))
  1150. save_dataframe(session_log, os.path.join(out_dir, "NOE_session_log.csv"))
  1151. save_dataframe(bout_qc, os.path.join(out_dir, "NOE_bout_qc.csv"))
  1153. return results_long, session_log, bout_qc
  1156. def run_hab_running_band_coherence(cfg, out_dir):
  1157. channel_data = load_channel_table(cfg)
  1159. win_sec = 0.5
  1160. win_overlap = 0.5
  1161. max_time_s = 600
  1162. speed_run_thresh = 5.0
  1163. min_windows = 3
  1164. artifact_thresh = 4.0
  1166. long_rows = []
  1167. session_log = []
  1169. for mi, mouse in enumerate(cfg["miceAvail"]):
  1170. batch_id = int(cfg["batch"][mi])
  1171. meta = load_meta_table(cfg, batch_id)
  1173. for cond in cfg["conditions"]:
  1174. test = "HAB"
  1175. session_id = f"{mouse}|{cond}|{test}"
  1177. if should_skip(mouse, cond, test):
  1178. session_log.append({"session_id": session_id, "status": "skipped_rule"})
  1179. continue
  1181. row_id, msg = find_meta_row(meta, mouse, cond, test)
  1182. if row_id is None:
  1183. session_log.append({"session_id": session_id, "status": msg})
  1184. continue
  1186. md = meta.loc[row_id].copy()
  1187. md["batch_id"] = batch_id
  1189. try:
  1190. channels, fs_raw = load_ephys_channels(cfg, md, batch_id, channel_data, mouse)
  1191. except Exception as e:
  1192. session_log.append({"session_id": session_id, "status": f"ephys_failed:{e}"})
  1193. continue
  1195. video_dir = resolve_video_dir(cfg["rootFoldVideo"], md["video_folder"], batch_id)
  1196. trial = str(md["trial_number"])
  1197. trial_suffix = trial[1:] if len(trial) > 1 else trial
  1198. xls_name = f"processed_interpolated_T{trial_suffix}_V2.xlsx"
  1199. xls_path = os.path.join(video_dir, xls_name)
  1200. if not os.path.exists(xls_path):
  1201. pattern = os.path.join(video_dir, f"processed_interpolated_T*{trial_suffix}*_V2.xlsx")
  1202. cand = glob.glob(pattern)
  1203. if cand:
  1204. xls_path = cand[0]
  1206. if not os.path.exists(xls_path):
  1207. session_log.append({"session_id": session_id, "status": "behavior_missing"})
  1208. continue
  1210. bdf = pd.read_excel(xls_path)
  1211. if "Timestamp_sec" not in bdf.columns or "nose" not in bdf.columns:
  1212. session_log.append({"session_id": session_id, "status": "missing_columns"})
  1213. continue
  1215. bdf.loc[bdf["Timestamp_sec"] <= 1, "Timestamp_sec"] = np.nan
  1216. bdf = bdf.dropna(subset=["Timestamp_sec"]).copy()
  1217. if len(bdf) == 0:
  1218. session_log.append({"session_id": session_id, "status": "empty_behavior"})
  1219. continue
  1220. bdf["Timestamp_sec"] = bdf["Timestamp_sec"] - bdf["Timestamp_sec"].iloc[0]
  1221. bdf = bdf[bdf["Timestamp_sec"] <= max_time_s].copy()
  1222. if len(bdf) < 100:
  1223. session_log.append({"session_id": session_id, "status": "too_short"})
  1224. continue
  1226. proc, fs, t_lfp = preprocess_channels(channels, fs_raw, cfg, max_time_s=max_time_s)
  1227. if not proc or len(t_lfp) < int(fs * 2):
  1228. session_log.append({"session_id": session_id, "status": "too_short_lfp"})
  1229. continue
  1231. speed_lfp = build_hab_speed_lfp_xy(bdf, t_lfp)
  1232. n = len(t_lfp)
  1233. half = int(round(fs * win_sec / 2))
  1234. step = int(win_sec * (1 - win_overlap) * fs)
  1235. all_centers = np.arange(half, n - half, step, dtype=int)
  1236. if len(all_centers) < 10:
  1237. session_log.append({"session_id": session_id, "status": "too_few_windows"})
  1238. continue
  1240. win_speeds = np.array([np.mean(speed_lfp[max(0, c - half) : min(n, c + half)]) for c in all_centers])
  1242. clean_mask = np.ones(len(all_centers), dtype=bool)
  1243. for _, x_ch in proc.items():
  1244. cl, _ = reject_artifact_windows(x_ch, list(all_centers), win_sec, fs, artifact_thresh)
  1245. cl_set = set(cl)
  1246. for j, c in enumerate(all_centers):
  1247. if c not in cl_set:
  1248. clean_mask[j] = False
  1250. clean_centers = all_centers[clean_mask]
  1251. clean_speeds = win_speeds[clean_mask]
  1252. run_centers = clean_centers[clean_speeds > speed_run_thresh]
  1254. if len(run_centers) < min_windows:
  1255. session_log.append({"session_id": session_id, "status": "too_few_running_windows"})
  1256. continue
  1258. valid_names = [k for k, v in proc.items() if isinstance(v, np.ndarray) and v.size > 0]
  1259. pair_list = list(itertools.combinations(valid_names, 2))
  1261. for chA, chB in pair_list:
  1262. pair_name, pair_hemi = pool_pair_name(chA, chB)
  1263. if pair_name is None:
  1264. continue
  1265. xA, xB = proc[chA], proc[chB]
  1266. for bname in ["theta", "beta", "lowGamma"]:
  1267. val = coherence_band_mean(xA, xB, fs, run_centers, win_sec, cfg["bands"][bname])
  1268. if np.isfinite(val):
  1269. long_rows.append(
  1270. {
  1271. "session_id": session_id,
  1272. "mouse": mouse,
  1273. "condition": cond,
  1274. "hemisphere": pair_hemi,
  1275. "domain": "pair",
  1276. "target": pair_name,
  1277. "target_figure_label": PAIR_FIGURE_LABEL.get(pair_name, pair_name),
  1278. "epoch": "running",
  1279. "win_sec": win_sec,
  1280. "n_windows": len(run_centers),
  1281. "metric": "coh",
  1282. "band": bname,
  1283. "value": float(val),
  1284. }
  1285. )
  1287. session_log.append(
  1288. {
  1289. "session_id": session_id,
  1290. "mouse": mouse,
  1291. "condition": cond,
  1292. "status": "ok",
  1293. "n_run": len(run_centers),
  1294. "mean_speed": float(np.mean(clean_speeds)) if len(clean_speeds) else np.nan,
  1295. }
  1296. )
  1298. results_long = pd.DataFrame(long_rows)
  1299. session_log = pd.DataFrame(session_log)
  1301. save_dataframe(results_long, os.path.join(out_dir, "HAB_running_band_coherence_results_long.csv"))
  1302. save_dataframe(session_log, os.path.join(out_dir, "HAB_running_band_coherence_session_log.csv"))
  1304. return results_long, session_log
  1307. def run_hab_psd_and_coherence_spectra(cfg, out_dir):
  1308. channel_data = load_channel_table(cfg)
  1310. win_sec_welch = 1.0
  1311. overlap_welch = 0.5
  1312. freq_range = (1, 80)
  1313. speed_thresh = 5.0
  1315. psd_store = defaultdict(lambda: defaultdict(lambda: defaultdict(list)))
  1316. coh_store = defaultdict(lambda: defaultdict(lambda: defaultdict(list)))
  1317. freqs_ref_psd = None
  1318. freqs_ref_coh = None
  1320. cond_map = {c: label for label, conds in CONDITION_GROUPS.items() for c in conds}
  1322. for mi, mouse in enumerate(cfg["miceAvail"]):
  1323. batch_id = int(cfg["batch"][mi])
  1324. meta = load_meta_table(cfg, batch_id)
  1326. for cond in cfg["conditions"]:
  1327. if cond not in cond_map:
  1328. continue
  1330. cond_label = cond_map[cond]
  1331. test = "HAB"
  1332. session_id = f"{mouse}|{cond}|{test}"
  1334. if should_skip(mouse, cond, test):
  1335. continue
  1337. row_id, _ = find_meta_row(meta, mouse, cond, test)
  1338. if row_id is None:
  1339. continue
  1341. md = meta.loc[row_id].copy()
  1342. md["batch_id"] = batch_id
  1344. try:
  1345. channels, fs_raw = load_ephys_channels(cfg, md, batch_id, channel_data, mouse)
  1346. except Exception:
  1347. continue
  1349. bdf, _, _, bmsg = load_behavior_epochs(md, test, cfg)
  1350. if bdf is None:
  1351. continue
  1353. proc, fs, t_lfp = preprocess_channels(channels, fs_raw, cfg)
  1354. if not proc:
  1355. continue
  1357. speed_lfp = compute_speed_on_lfp_time(bdf, t_lfp, cfg["px_to_cm"])
  1358. running_mask = speed_lfp >= speed_thresh
  1359. if running_mask.sum() < int(win_sec_welch * fs):
  1360. continue
  1362. nperseg = int(win_sec_welch * fs)
  1363. noverlap = int(nperseg * overlap_welch)
  1365. for ch_name, x in proc.items():
  1366. if ch_name not in CHANNEL_TO_REGION:
  1367. continue
  1368. region, hemi = pool_channel_name(ch_name)
  1369. if region not in PSD_REGIONS:
  1370. continue
  1372. store_key = f"{mouse}|{hemi}"
  1373. x_run = x[running_mask]
  1374. if len(x_run) < nperseg:
  1375. continue
  1377. freqs, psd = welch(
  1378. x_run,
  1379. fs=fs,
  1380. nperseg=nperseg,
  1381. noverlap=noverlap,
  1382. window="hann",
  1383. scaling="density",
  1384. )
  1385. if freqs_ref_psd is None:
  1386. freqs_ref_psd = freqs
  1387. psd_store[cond_label][region][store_key].append(psd)
  1389. for a, b in itertools.combinations(sorted(proc.keys()), 2):
  1390. pair_name, pair_hemi = pool_pair_name(a, b)
  1391. if pair_name is None:
  1392. continue
  1393. xa = proc[a][running_mask]
  1394. xb = proc[b][running_mask]
  1395. if len(xa) < nperseg or len(xb) < nperseg:
  1396. continue
  1397. freqs_c, coh = signal.coherence(
  1398. xa,
  1399. xb,
  1400. fs=fs,
  1401. nperseg=nperseg,
  1402. noverlap=noverlap,
  1403. window="hann",
  1404. )
  1405. if freqs_ref_coh is None:
  1406. freqs_ref_coh = freqs_c
  1407. store_key = f"{mouse}|{pair_hemi}"
  1408. coh_store[cond_label][pair_name][store_key].append(coh)
  1410. if freqs_ref_psd is None:
  1411. raise RuntimeError("No HAB running PSD could be computed.")
  1412. if freqs_ref_coh is None:
  1413. raise RuntimeError("No HAB running coherence spectrum could be computed.")
  1415. freq_mask_psd = (freqs_ref_psd >= freq_range[0]) & (freqs_ref_psd <= freq_range[1])
  1416. freq_mask_coh = (freqs_ref_coh >= freq_range[0]) & (freqs_ref_coh <= freq_range[1])
  1418. psd_subject_rows = []
  1419. for cond_label, reg_dict in psd_store.items():
  1420. for region, subj_dict in reg_dict.items():
  1421. for subj_key, psd_list in subj_dict.items():
  1422. subj_mean = np.nanmean(np.stack(psd_list, axis=0), axis=0)
  1423. for f, v in zip(freqs_ref_psd[freq_mask_psd], subj_mean[freq_mask_psd]):
  1424. psd_subject_rows.append(
  1425. {
  1426. "analysis": "running_psd",
  1427. "condition_group": cond_label,
  1428. "condition_label": CONDITION_LABEL[cond_label],
  1429. "region": region,
  1430. "region_figure_label": REGION_FIGURE_LABEL[region],
  1431. "subject": subj_key,
  1432. "frequency_hz": float(f),
  1433. "value": float(v),
  1434. }
  1435. )
  1437. psd_subject = pd.DataFrame(psd_subject_rows)
  1438. psd_summary = compute_group_summary(psd_subject, "value")
  1440. coh_subject_rows = []
  1441. for cond_label, pair_dict in coh_store.items():
  1442. for pair_name, subj_dict in pair_dict.items():
  1443. for subj_key, coh_list in subj_dict.items():
  1444. subj_mean = np.nanmean(np.stack(coh_list, axis=0), axis=0)
  1445. for f, v in zip(freqs_ref_coh[freq_mask_coh], subj_mean[freq_mask_coh]):
  1446. coh_subject_rows.append(
  1447. {
  1448. "analysis": "running_coherence_spectrum",
  1449. "condition_group": cond_label,
  1450. "condition_label": CONDITION_LABEL[cond_label],
  1451. "target": pair_name,
  1452. "target_figure_label": PAIR_FIGURE_LABEL[pair_name],
  1453. "subject": subj_key,
  1454. "frequency_hz": float(f),
  1455. "value": float(v),
  1456. }
  1457. )
  1459. coh_subject = pd.DataFrame(coh_subject_rows)
  1460. coh_summary = compute_group_summary(coh_subject, "value")
  1462. save_dataframe(psd_subject, os.path.join(out_dir, "running_psd_subject_curves.csv"))
  1463. save_dataframe(psd_summary, os.path.join(out_dir, "running_psd_group_summary.csv"))
  1464. save_dataframe(coh_subject, os.path.join(out_dir, "running_coherence_spectrum_subject_curves.csv"))
  1465. save_dataframe(coh_summary, os.path.join(out_dir, "running_coherence_spectrum_group_summary.csv"))
  1467. return psd_subject, psd_summary, coh_subject, coh_summary
  1470. def _star(p):
  1471. if not np.isfinite(p):
  1472. return "ns"
  1473. if p < 0.001:
  1474. return "***"
  1475. if p < 0.01:
  1476. return "**"
  1477. if p < 0.05:
  1478. return "*"
  1479. return "ns"
  1482. def export_selected_bar_panel_csvs(panel_specs, stats_dict, out_dir, prefix):
  1483. stats_all = stats_dict["stats"]
  1484. plot_data_norm = stats_dict["plot_data_norm"]
  1485. plot_data_raw = stats_dict["plot_data_raw"]
  1487. combined_rows = []
  1489. for spec in panel_specs:
  1490. mask = (
  1491. (stats_all["domain"] == "pair")
  1492. & (stats_all["target"] == spec["target"])
  1493. & (stats_all["metric"] == spec["metric"])
  1494. & (stats_all["band"] == spec["band"])
  1495. )
  1496. stat_row = stats_all.loc[mask].copy()
  1497. if len(stat_row) == 0:
  1498. continue
  1499. stat_row = stat_row.iloc[[0]].copy()
  1500. stat_row["panel"] = spec["panel"]
  1501. stat_row["figure_title"] = spec["title"]
  1503. plot_mask = (
  1504. np.isclose(plot_data_norm["win_sec"].astype(float), float(stat_row["win_sec"].iloc[0]))
  1505. & (plot_data_norm["epoch"] == stat_row["epoch"].iloc[0])
  1506. & (plot_data_norm["domain"] == stat_row["domain"].iloc[0])
  1507. & (plot_data_norm["target"] == stat_row["target"].iloc[0])
  1508. & (plot_data_norm["metric"] == stat_row["metric"].iloc[0])
  1509. & (plot_data_norm["band"] == stat_row["band"].iloc[0])
  1510. )
  1511. d_norm = plot_data_norm.loc[plot_mask].copy()
  1512. d_raw = plot_data_raw.loc[plot_mask].copy()
  1514. panel_rows = []
  1515. for _, row in d_norm.iterrows():
  1516. raw_row = d_raw.loc[d_raw["subject"] == row["subject"]]
  1517. raw_vals = raw_row.iloc[0] if len(raw_row) else None
  1518. for cond_key in ["naive", "cdm", "c21"]:
  1519. out_row = {
  1520. "panel": spec["panel"],
  1521. "figure_title": spec["title"],
  1522. "subject": row["subject"],
  1523. "target": row["target"],
  1524. "target_figure_label": PAIR_FIGURE_LABEL.get(row["target"], row["target"]),
  1525. "metric": row["metric"],
  1526. "band": row["band"],
  1527. "condition_group": cond_key.upper() if cond_key != "naive" else "naive",
  1528. "condition_label": CONDITION_LABEL["naive"] if cond_key == "naive" else CONDITION_LABEL[cond_key.upper()],
  1529. "value_norm": float(row[cond_key]),
  1530. "value_raw": float(raw_vals[cond_key]) if raw_vals is not None else np.nan,
  1531. }
  1532. if row["metric"] == "pac_mi":
  1533. phase_region, amp_region = PAC_DIRECTION_LABEL.get(row["target"], ("", ""))
  1534. out_row["phase_region"] = phase_region
  1535. out_row["amplitude_region"] = amp_region
  1536. panel_rows.append(out_row)
  1538. panel_df = pd.DataFrame(panel_rows)
  1539. save_dataframe(panel_df, os.path.join(out_dir, f"{prefix}_panel_{spec['panel']}.csv"))
  1540. save_dataframe(stat_row, os.path.join(out_dir, f"{prefix}_panel_{spec['panel']}_stats.csv"))
  1541. combined_rows.append(panel_df)
  1543. if combined_rows:
  1544. save_dataframe(pd.concat(combined_rows, ignore_index=True), os.path.join(out_dir, f"{prefix}_all_panels_long.csv"))
  1547. def export_curve_panel_csvs(panel_key, subject_df, summary_df, out_dir):
  1548. save_dataframe(subject_df, os.path.join(out_dir, f"{panel_key}_subject_curves.csv"))
  1549. save_dataframe(summary_df, os.path.join(out_dir, f"{panel_key}_group_summary.csv"))
  1552. def add_pairwise_brackets(ax, stat_row, y_top, y_step):
  1553. p_nc = float(stat_row.get("p_naive_cdm", np.nan))
  1554. p_cc = float(stat_row.get("p_cdm_c21", np.nan))
  1555. p_n2 = float(stat_row.get("p_naive_c21", np.nan))
  1557. level = y_top
  1558. if np.isfinite(p_nc) and p_nc < 0.05:
  1559. ax.plot([0, 0, 1, 1], [level, level + 0.04 * y_step, level + 0.04 * y_step, level], "k-", lw=0.8)
  1560. ax.text(0.5, level + 0.05 * y_step, _star(p_nc), ha="center", va="bottom", fontsize=9)
  1561. level += 0.12 * y_step
  1563. if np.isfinite(p_cc) and p_cc < 0.05:
  1564. ax.plot([1, 1, 2, 2], [level, level + 0.04 * y_step, level + 0.04 * y_step, level], "k-", lw=0.8)
  1565. ax.text(1.5, level + 0.05 * y_step, _star(p_cc), ha="center", va="bottom", fontsize=9)
  1566. level += 0.12 * y_step
  1568. if np.isfinite(p_n2) and p_n2 < 0.05:
  1569. ax.plot([0, 0, 2, 2], [level, level + 0.04 * y_step, level + 0.04 * y_step, level], "k-", lw=0.8)
  1570. ax.text(1.0, level + 0.05 * y_step, _star(p_n2), ha="center", va="bottom", fontsize=9)
  1573. def plot_selected_bar_panel(ax, spec, stats_dict):
  1574. stats_all = stats_dict["stats"]
  1575. plot_data_norm = stats_dict["plot_data_norm"]
  1577. mask = (
  1578. (stats_all["domain"] == "pair")
  1579. & (stats_all["target"] == spec["target"])
  1580. & (stats_all["metric"] == spec["metric"])
  1581. & (stats_all["band"] == spec["band"])
  1582. )
  1583. if mask.sum() == 0:
  1584. ax.set_axis_off()
  1585. ax.text(0.5, 0.5, f"Missing panel {spec['panel']}", ha="center", va="center")
  1586. return
  1588. stat_row = stats_all.loc[mask].iloc[0]
  1589. plot_mask = (
  1590. np.isclose(plot_data_norm["win_sec"].astype(float), float(stat_row["win_sec"]))
  1591. & (plot_data_norm["epoch"] == stat_row["epoch"])
  1592. & (plot_data_norm["domain"] == stat_row["domain"])
  1593. & (plot_data_norm["target"] == stat_row["target"])
  1594. & (plot_data_norm["metric"] == stat_row["metric"])
  1595. & (plot_data_norm["band"] == stat_row["band"])
  1596. )
  1597. d = plot_data_norm.loc[plot_mask].copy()
  1598. if len(d) == 0:
  1599. ax.set_axis_off()
  1600. ax.text(0.5, 0.5, f"No data for panel {spec['panel']}", ha="center", va="center")
  1601. return
  1603. x = np.array([0, 1, 2], dtype=float)
  1604. values = np.c_[d["naive"].values.astype(float), d["cdm"].values.astype(float), d["c21"].values.astype(float)]
  1606. means = np.nanmean(values, axis=0)
  1607. ax.bar(x, means, width=0.55, color=[CONDITION_COLOR["naive"], CONDITION_COLOR["CDM"], CONDITION_COLOR["C21"]], edgecolor="black", linewidth=0.8, zorder=1)
  1609. for row in values:
  1610. ax.plot(x, row, "-o", color="0.35", lw=0.8, ms=3.5, mfc="white", zorder=2)
  1611. ax.scatter([0], [row[0]], s=12, color=CONDITION_COLOR["naive"], zorder=3, edgecolor="black", linewidth=0.2)
  1612. ax.scatter([1], [row[1]], s=12, color=CONDITION_COLOR["CDM"], zorder=3, edgecolor="black", linewidth=0.2)
  1613. ax.scatter([2], [row[2]], s=12, color=CONDITION_COLOR["C21"], zorder=3, edgecolor="black", linewidth=0.2)
  1615. all_vals = values.flatten()
  1616. ymax = float(np.nanmax(all_vals))
  1617. ymin = float(np.nanmin(all_vals))
  1618. dy = (ymax - ymin) if ymax > ymin else (abs(ymax) * 0.2 + 1e-3)
  1619. add_pairwise_brackets(ax, stat_row, ymax + 0.18 * dy, dy)
  1621. ax.axhline(0, color="black", lw=0.9)
  1622. ax.set_xlim(-0.6, 2.6)
  1623. ax.set_xticks([0, 1, 2])
  1624. ax.set_xticklabels([CONDITION_LABEL["naive"], CONDITION_LABEL["CDM"], CONDITION_LABEL["C21"]], rotation=0, fontsize=8)
  1625. ax.set_ylabel("z-score", fontsize=8)
  1626. ax.set_title(spec["title"], fontsize=9, fontweight="bold")
  1627. ax.spines["top"].set_visible(False)
  1628. ax.spines["right"].set_visible(False)
  1629. ax.tick_params(labelsize=8)
  1632. def add_band_guides(ax, freq_max):
  1633. guides = [(4, "θ"), (12, "β"), (30, "low-γ"), (60, "high-γ")]
  1634. for x, label in guides:
  1635. if x <= freq_max:
  1636. ax.axvline(x, ls=":", lw=0.8, color="0.55", zorder=0)
  1637. labels = [(2.5, "δ"), (8, "θ"), (21, "β"), (45, "low-γ"), (70, "high-γ")]
  1638. for x, label in labels:
  1639. if x <= freq_max:
  1640. ax.text(x, 1.01, label, transform=ax.get_xaxis_transform(), ha="center", va="bottom", fontsize=8)
  1643. def plot_curve_panel(ax, summary_df, value_col, x_col, title, ylabel, ylim=None, log_y=False):
  1644. for cond_key in ["naive", "CDM", "C21"]:
  1645. cond_data = summary_df[summary_df["condition_group"] == cond_key].copy()
  1646. if len(cond_data) == 0:
  1647. continue
  1648. x = cond_data[x_col].values.astype(float)
  1649. y = cond_data["mean"].values.astype(float)
  1650. sem = cond_data["sem"].values.astype(float)
  1651. order = np.argsort(x)
  1652. x = x[order]
  1653. y = y[order]
  1654. sem = sem[order]
  1655. if log_y:
  1656. ax.semilogy(x, y, color=CONDITION_COLOR[cond_key], lw=1.8, label=CONDITION_LABEL[cond_key])
  1657. ax.fill_between(x, np.clip(y - sem, 1e-20, None), y + sem, color=CONDITION_COLOR[cond_key], alpha=0.22)
  1658. else:
  1659. ax.plot(x, y, color=CONDITION_COLOR[cond_key], lw=1.8, label=CONDITION_LABEL[cond_key])
  1660. ax.fill_between(x, y - sem, y + sem, color=CONDITION_COLOR[cond_key], alpha=0.22)
  1662. add_band_guides(ax, float(np.nanmax(summary_df[x_col].values.astype(float))))
  1663. ax.set_xlabel("Frequency (Hz)", fontsize=8)
  1664. ax.set_ylabel(ylabel, fontsize=8)
  1665. ax.set_title(title, fontsize=9, fontweight="bold")
  1666. ax.spines["top"].set_visible(False)
  1667. ax.spines["right"].set_visible(False)
  1668. ax.tick_params(labelsize=8)
  1669. if ylim is not None:
  1670. ax.set_ylim(ylim)
  1673. def make_quantitative_figure(psd_summary, coh_summary, hab_stats, noe_stats, out_dir):
  1674. fig, axes = plt.subplots(4, 3, figsize=(12.5, 14.5), constrained_layout=True)
  1676. for j, region in enumerate(["PFC", "RE", "HC"]):
  1677. d = psd_summary[psd_summary["region"] == region].copy()
  1678. plot_curve_panel(
  1679. axes[0, j],
  1680. d,
  1681. value_col="mean",
  1682. x_col="frequency_hz",
  1683. title=f"Power spectrum - {REGION_FIGURE_LABEL[region]}",
  1684. ylabel="PSD (V²/Hz)",
  1685. log_y=True,
  1686. )
  1687. axes[0, j].set_xlim(1, 80)
  1688. axes[0, j].legend(frameon=False, fontsize=7)
  1690. for j, pair_name in enumerate(["PFC-RE", "HC-RE", "HC-PFC"]):
  1691. d = coh_summary[coh_summary["target"] == pair_name].copy()
  1692. plot_curve_panel(
  1693. axes[1, j],
  1694. d,
  1695. value_col="mean",
  1696. x_col="frequency_hz",
  1697. title=f"Coherence - {PAIR_FIGURE_LABEL[pair_name]}",
  1698. ylabel="Coherence",
  1699. log_y=False,
  1700. )
  1701. axes[1, j].set_xlim(1, 80)
  1702. axes[1, j].set_ylim(0, 1)
  1703. axes[1, j].legend(frameon=False, fontsize=7)
  1705. for j, spec in enumerate(BAR_PANELS_RUNNING):
  1706. plot_selected_bar_panel(axes[2, j], spec, hab_stats)
  1708. for j, spec in enumerate(BAR_PANELS_NOE):
  1709. plot_selected_bar_panel(axes[3, j], spec, noe_stats)
  1711. panel_letters = [
  1712. ["g", "h", "i"],
  1713. ["j", "k", "l"],
  1714. ["m", "n", "o"],
  1715. ["p", "q", "r"],
  1716. ]
  1717. for i in range(4):
  1718. for j in range(3):
  1719. axes[i, j].text(-0.15, 1.06, panel_letters[i][j], transform=axes[i, j].transAxes, fontsize=12, fontweight="bold")
  1721. for ext in ["png", "pdf", "svg"]:
  1722. path = os.path.join(out_dir, f"Figure5_quantitative_panels.{ext}")
  1723. fig.savefig(path, dpi=300 if ext == "png" else None, bbox_inches="tight")
  1725. plt.close(fig)
  1728. def main():
  1729. cfg = make_cfg()
  1730. output_root = os.path.join(cfg["resultsDir"], "figure5_source_data")
  1731. noe_dir = os.path.join(output_root, "NOE")
  1732. hab_dir = os.path.join(output_root, "HAB_running")
  1733. panel_dir = os.path.join(output_root, "panels")
  1734. os.makedirs(noe_dir, exist_ok=True)
  1735. os.makedirs(hab_dir, exist_ok=True)
  1736. os.makedirs(panel_dir, exist_ok=True)
  1738. noe_results_long, noe_session_log, noe_bout_qc = run_noe_extraction(cfg, noe_dir)
  1739. noe_stats = compute_triplet_stats(
  1740. results_long=noe_results_long,
  1741. out_dir=noe_dir,
  1742. prefix="NOE_delta_pair_coh_pac",
  1743. epoch_keep=["delta"],
  1744. domain_keep=["pair"],
  1745. metric_keep=["coh", "pac_mi"],
  1746. bands_keep=["theta", "beta", "lowGamma", "theta_lowGamma"],
  1747. win_sec=0.5,
  1748. pool_hemispheres=False,
  1749. require_both_hemis=False,
  1750. c21_conds=["C21_1"],
  1751. normalize_method="zscore",
  1752. apply_fdr=False,
  1753. min_n_event_windows=2,
  1754. )
  1756. hab_running_results_long, hab_running_session_log = run_hab_running_band_coherence(cfg, hab_dir)
  1757. hab_stats = compute_triplet_stats(
  1758. results_long=hab_running_results_long,
  1759. out_dir=hab_dir,
  1760. prefix="HAB_running_pair_coherence",
  1761. epoch_keep=["running"],
  1762. domain_keep=["pair"],
  1763. metric_keep=["coh"],
  1764. bands_keep=["theta", "beta", "lowGamma"],
  1765. win_sec=0.5,
  1766. pool_hemispheres=False,
  1767. require_both_hemis=False,
  1768. c21_conds=["C21_1"],
  1769. normalize_method="zscore",
  1770. apply_fdr=False,
  1771. min_n_event_windows=0,
  1772. )
  1774. psd_subject, psd_summary, coh_subject, coh_summary = run_hab_psd_and_coherence_spectra(cfg, hab_dir)
  1776. export_selected_bar_panel_csvs(BAR_PANELS_RUNNING, hab_stats, panel_dir, "running")
  1777. export_selected_bar_panel_csvs(BAR_PANELS_NOE, noe_stats, panel_dir, "noe")
  1779. export_curve_panel_csvs("power_spectrum", psd_subject, psd_summary, panel_dir)
  1780. export_curve_panel_csvs("coherence_spectrum", coh_subject, coh_summary, panel_dir)
  1782. make_quantitative_figure(psd_summary, coh_summary, hab_stats, noe_stats, panel_dir)
  1784. print(f"NOE rows: {len(noe_results_long)}")
  1785. print(f"HAB running rows: {len(hab_running_results_long)}")
  1786. print(f"PSD subject rows: {len(psd_subject)}")
  1787. print(f"Coherence spectrum subject rows: {len(coh_subject)}")
  1788. print(f"Output root: {output_root}")
  1791. if __name__ == "__main__":
  1792. main()

LFP_analysis.py at commit 8ce9fb2, under MIT · at the source

Overview

Authors: Maxime Veleanu1, Louise Schuberth1,2,3, Antje Kilias4, Jakob Weber1, Jan Warneke1, Lovis Würz1, Tim Schwär1, Rebecca Heck1, Joelle Müller1, David H. Sarrazin5, Marguerite Anselin1, Lukas Rutke1, Guillermo Suarez Marchi1, Stella Zimmermann1, Zoe Borgeest1, Martin Balzinger5,6, Alina Blendinger5, Anna Catarata1, Samira Assaad Dib1, Yaroslav Sych5
and 6 other authorsThibault Cholvin4, Marlene Bartos4, Katharina Domschke1,7, Claus Normann1, Stefan Vestring1, Tsvetan Serchov1,5,6
ORCID iDs: Maxime Veleanu, Louise Schuberth, Antje Kilias, David H. Sarrazin, Guillermo Suarez Marchi, Stella Zimmermann, Yaroslav Sych, Thibault Cholvin, Marlene Bartos, Katharina Domschke, Claus Normann, Stefan Vestring, Tsvetan Serchov
  1. Department of Psychiatry and Psychotherapy, Medical Center—University of Freiburg, Faculty of Medicine, University of Freiburg,Freiburg, Germany
  2. Spemann Graduate School of Biology and Medicine (SGBM), Freiburg, Germany
  3. Faculty of Biology, University of Freiburg,Freiburg, Germany
  4. Institute for Physiology I, Medical Faculty, University of Freiburg,Freiburg, Germany
  5. Centre National de la Recherche Scientifique (CNRS), Institute of Cellular and Integrative Neurosciences (INCI) UPR, University of Strasbourg,Strasbourg, 3212 France
  6. University of Strasbourg Institute for Advanced Study (USIAS),Strasbourg, France
  7. German Center for Mental Health (DZPG), Partner Site Berlin/Potsdam,Berlin, Germany
Institutions: University of Freiburg (Germany); Université de Strasbourg (France); Deutsches Zentrum für Psychische Gesundheit (Germany)
Journal: Nature communications, volume 17, issue 1, article 9361
Dates: received 31 March 2025; accepted 17 August 2026; published online 1 September 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.1038/s41467-026-77248-y · PMID 42680757 · PMCID PMC13534506 · OpenAlex W7204912409
Open access: gold, a free copy (OpenAlex)
Status: code verified
Categories: mouse (organism), depression (population)
Methods: Spectral & time-frequency, Connectivity, Statistics, Preprocessing, Evoked potentials, fMRI & imaging, Smoothing, state filtering, decompositions
Keywords: Neural circuits, Depression, Synaptic plasticity, Stress and resilience
MeSH: Antidepressive Agents*, Depression*, Hippocampus*, Midline Thalamic Nuclei*, Neuronal Plasticity*, Prefrontal Cortex*, Animals, Disease Models, Animal, Ketamine, Long-Term Potentiation, Male, Mice, Mice, Inbred C57BL, Neural Pathways (* major topic)
Topic: Neurotransmitter Receptor Influence on Behavior (Cellular and Molecular Neuroscience, Neuroscience), according to OpenAlex
Funding: Berta-Ottenstein Program for Clinician Scientists, Faculty of Medicine, University of Freiburg, Germany; Hans A. Krebs Medical Scientist Fellowship Program, University of Freiburg, Germany; Deutsche Forschungsgemeinschaft (DFG) (DFG NO370/7-1, SE 2666/2-1, SE 2666/2-3); Centre National de la Recherche Scientifique (CNRS) (CNRS UPR3212); Université de Strasbourg (Unistra) (2021-2028)
Citations: not cited yet (Europe PMC); 103 references in the paper
Research resources: Wild-type C57BL/6 J mice RRID:IMSR_JAX:000664

Abstract

The pathophysiology of depression involves multiple biological processes, including circuit dysfunction and impaired neuroplasticity, yet an integrative view linking these processes remains elusive. Here, we identify a convergent circuit for antidepressant response and plasticity modulation. We demonstrate that chemogenetic activation of the infralimbic cortex (IL) exerts rapid antidepressant-like effects across multiple behavioral domains in a mouse model of stress-induced depression. IL stimulation exerts top-down control over the hippocampus, enhancing structural plasticity, restoring long-term potentiation deficits and improving state-dependent network dynamics in the ventral hippocampus (vHIPP). We identify the thalamic nucleus reuniens (RE) as a necessary mediator of these effects. Notably, direct inhibition of RE, its inputs from IL or projections to vHIPP, blocks both IL stimulation-induced antidepressant response and the therapeutic and neuroplastic effects of ketamine. Our findings demonstrate that the functional IL → RE→vHIPP circuit plays a central role in the antidepressant response, linking circuit activity, hippocampal plasticity, and depressive-like behaviors.

Reproduced under the paper's license (CC BY), from the paper cited above.

Repositories

Its files are read in the Code ↔ Paper reader above, with 4 matches between paragraphs and lines of code.

mveleanu/veleanu_et_al

License: MIT
State: the link answers, verified on 26 September 2026
Evidence: files inventoried
Commit: 8ce9fb28c1b42d135b796a2780b46aee923ac44b, 30 March 2026
Languages: Python (7), MATLAB (2)
Size: 98 files, 9 scripts
Software Heritage: not archived
Found in: “Code availability”
Holds: README, license file
Not found: CITATION.cff, environment file, tests, continuous integration, documentation
Tools: NumPy (6 files), pandas (6 files), Matplotlib (3 files), SciPy (3 files), h5py (1 file), Signal Processing Toolbox (1 file), Open Ephys analysis tools (1 file)
Availability: 1 check, the latest on 26 September 2026: the link answers
  • 26 September 2026: the link answers
11 files
At the source: github.com/mveleanu/veleanu_et_al

Zenodo 21111701

License: MIT
State: the link answers, verified on 26 September 2026
Evidence: files inventoried
Size: 1 file
Software Heritage: not checked
Found in: “Code availability”
Not found: README, license file, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: NumPy (6 files), pandas (6 files), Matplotlib (3 files), SciPy (3 files), h5py (1 file), Signal Processing Toolbox (1 file), Open Ephys analysis tools (1 file)
Availability: 1 check, the latest on 26 September 2026: the link answers (HTTP 200)
  • 26 September 2026: the link answers (HTTP 200)
11 files
At the source:

Code availability

Custom Python scripts used for analyzing fiber photometry and LFP data are available in our GitHub repository (https://github.com/mveleanu/veleanu_et_al) and archived on Zenodo (https://doi.org/10.5281/zenodo.21111701)103.

Reproduced under the paper's license (CC BY), from the paper cited above.

Tracing map

Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.

What the map holds:

  • 2 repositories of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
  • 18 scripts, each with its path and the digest of its content;
  • 4 matches between paragraphs of the paper and lines of the code (method lexical-v1);
  • neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

No dataset and no data link were found in the paper.

Data availability

The data generated in this study are provided in the Supplementary Information and the Source Data file. Source data are provided with this paper.

Reproduced under the paper's license (CC BY), from the paper cited above.

Versions

The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.

Version 1, 27 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 26 authors, 4 keywords, 14 MeSH terms, 5 funders, 102 references, 1 RRID.

Cite

This paper

Veleanu, M., Schuberth, L., Kilias, A., Weber, J., Warneke, J., Würz, L., Schwär, T., Heck, R., Müller, J., Sarrazin, D. H., Anselin, M., Rutke, L., Suarez Marchi, G., Zimmermann, S., Borgeest, Z., Balzinger, M., Blendinger, A., Catarata, A., Dib, S. A., . . . Serchov, T. (2026). The mPFC-reuniens-hippocampus pathway links brain circuitry and neural plasticity in antidepressant response. Nature communications, 17(1), 9361. https://doi.org/10.1038/s41467-026-77248-y

BibTeX

@article{veleanu2026mpfc,
author = {Veleanu, Maxime and Schuberth, Louise and Kilias, Antje and Weber, Jakob and Warneke, Jan and Würz, Lovis and Schwär, Tim and Heck, Rebecca and Müller, Joelle and Sarrazin, David H. and Anselin, Marguerite and Rutke, Lukas and Suarez Marchi, Guillermo and Zimmermann, Stella and Borgeest, Zoe and Balzinger, Martin and Blendinger, Alina and Catarata, Anna and Dib, Samira Assaad and Sych, Yaroslav and Cholvin, Thibault and Bartos, Marlene and Domschke, Katharina and Normann, Claus and Vestring, Stefan and Serchov, Tsvetan},
title = {{The mPFC-reuniens-hippocampus pathway links brain circuitry and neural plasticity in antidepressant response}},
journal = {Nature communications},
year = {2026},
month = sep,
volume = {17},
number = {1},
pages = {9361},
publisher = {Nature Publishing Group},
issn = {2041-1723},
doi = {10.1038/s41467-026-77248-y},
url = {https://doi.org/10.1038/s41467-026-77248-y},
pmid = {42680757},
pmcid = {PMC13534506}
}

RIS

TY - JOUR
AU - Veleanu, Maxime
AU - Schuberth, Louise
AU - Kilias, Antje
AU - Weber, Jakob
AU - Warneke, Jan
AU - Würz, Lovis
AU - Schwär, Tim
AU - Heck, Rebecca
AU - Müller, Joelle
AU - Sarrazin, David H.
AU - Anselin, Marguerite
AU - Rutke, Lukas
AU - Suarez Marchi, Guillermo
AU - Zimmermann, Stella
AU - Borgeest, Zoe
AU - Balzinger, Martin
AU - Blendinger, Alina
AU - Catarata, Anna
AU - Dib, Samira Assaad
AU - Sych, Yaroslav
AU - Cholvin, Thibault
AU - Bartos, Marlene
AU - Domschke, Katharina
AU - Normann, Claus
AU - Vestring, Stefan
AU - Serchov, Tsvetan
TI - The mPFC-reuniens-hippocampus pathway links brain circuitry and neural plasticity in antidepressant response
T2 - Nature communications
J2 - Nat Commun
PY - 2026
DA - 2026/09/01
VL - 17
IS - 1
SP - 9361
SN - 2041-1723
PB - Nature Publishing Group
DO - 10.1038/s41467-026-77248-y
UR - https://doi.org/10.1038/s41467-026-77248-y
LA - en
ER -

CSL-JSON

{
"id": "10.1038/s41467-026-77248-y",
"type": "article-journal",
"title": "The mPFC-reuniens-hippocampus pathway links brain circuitry and neural plasticity in antidepressant response",
"container-title": "Nature communications",
"author": [
{
"family": "Veleanu",
"given": "Maxime"
},
{
"family": "Schuberth",
"given": "Louise"
},
{
"family": "Kilias",
"given": "Antje"
},
{
"family": "Weber",
"given": "Jakob"
},
{
"family": "Warneke",
"given": "Jan"
},
{
"family": "Würz",
"given": "Lovis"
},
{
"family": "Schwär",
"given": "Tim"
},
{
"family": "Heck",
"given": "Rebecca"
},
{
"family": "Müller",
"given": "Joelle"
},
{
"family": "Sarrazin",
"given": "David H."
},
{
"family": "Anselin",
"given": "Marguerite"
},
{
"family": "Rutke",
"given": "Lukas"
},
{
"family": "Suarez Marchi",
"given": "Guillermo"
},
{
"family": "Zimmermann",
"given": "Stella"
},
{
"family": "Borgeest",
"given": "Zoe"
},
{
"family": "Balzinger",
"given": "Martin"
},
{
"family": "Blendinger",
"given": "Alina"
},
{
"family": "Catarata",
"given": "Anna"
},
{
"family": "Dib",
"given": "Samira Assaad"
},
{
"family": "Sych",
"given": "Yaroslav"
},
{
"family": "Cholvin",
"given": "Thibault"
},
{
"family": "Bartos",
"given": "Marlene"
},
{
"family": "Domschke",
"given": "Katharina"
},
{
"family": "Normann",
"given": "Claus"
},
{
"family": "Vestring",
"given": "Stefan"
},
{
"family": "Serchov",
"given": "Tsvetan"
}
],
"container-title-short": "Nat Commun",
"volume": "17",
"issue": "1",
"page": "9361",
"DOI": "10.1038/s41467-026-77248-y",
"PMID": "42680757",
"PMCID": "PMC13534506",
"ISSN": "2041-1723",
"publisher": "Nature Publishing Group",
"URL": "https://doi.org/10.1038/s41467-026-77248-y",
"language": "en",
"issued": {
"date-parts": [
[
2026,
9,
1
]
]
}
}

The tracing map gets a citation of its own once an author has validated it and it has a DOI.

Similar papers

The papers with a page that share the most with this one: the tools found in their code, their categories, datasets, cited references and authors, the rarest counting most.

[1] doi:10.1038/s41467-026-76574-5 [code]
Prefrontal-raphe-habenular circuit drives fear memory recall and updating.
Journal: Nature communications
In common: Signal Processing Toolbox, mouse, 5 references
[2] doi:10.1038/s41593-026-02262-8 [code]
Cheese3D enables sensitive detection and analysis of whole-face movement in mice.
Journal: Nature neuroscience
In common: Open Ephys analysis tools, h5py, pandas, 3 other tools, mouse, 1 reference
[3] doi:10.1016/j.isci.2026.116825 [code]
Social hierarchy shapes behavioral and transcriptional responses to chronic stress and ketamine in male mice.
Journal: iScience
In common: h5py, pandas, SciPy, 2 other tools, depression, mouse, 2 references
[4] doi:10.1038/s41398-026-04128-w
A nucleus accumbens-projecting prefrontal cortex circuit underlies chronic social stress-induced depression-like behaviors.
Journal: Translational psychiatry
In common: depression, mouse, 4 references
[5] doi:10.1016/j.celrep.2026.117590 [code]
Impaired behavioral inhibition in Fmr1 KO mice is linked to disrupted visual cortex theta oscillations.
Journal: Cell reports
In common: Open Ephys analysis tools, pandas, SciPy, 2 other tools, mouse, 1 reference
[6] doi:10.1186/s12974-026-03890-4 [code]
Neuronal toll-like receptor-4 regulation of matrix metalloproteinase-9 activity mediates dentate circuit dysfunction after traumatic brain injury.
Journal: Journal of neuroinflammation
In common: Open Ephys analysis tools, h5py, pandas, 3 other tools, mouse
[7] doi:10.1038/s41593-026-02255-7 [code]
Neural circuits encode prior knowledge of temporal statistics.
Journal: Nature neuroscience
In common: Open Ephys analysis tools, Signal Processing Toolbox, pandas, 3 other tools, mouse
[8] doi:10.1038/s41398-026-04025-2 [code]
Brain energetic landscapes shape state dysregulation in major depressive disorder: a morphological network controllability perspective.
Journal: Translational psychiatry
In common: h5py, pandas, SciPy, 2 other tools, depression, 2 references
[9] doi:10.1038/s44220-026-00680-y [code]
The neuroimaging correlates of depression established across six large-scale population datasets.
Journal: Nature. Mental health
In common: pandas, SciPy, Matplotlib, 1 other tool, depression, 3 references
[10] doi:10.1038/s41467-026-76581-6 [code]
Thalamocortical bursts encode reward contingencies and drive associative learning.
Journal: Nature communications
In common: h5py, Signal Processing Toolbox, pandas, 3 other tools, mouse, 1 reference

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: questions and error reports about this paper and its code, from signed-in readers and its authors. It opens with sign-in.

Reproductions: reports from readers who ran the authors' code: what they reproduced, with which environment, commit and data. It opens with sign-in.

Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.