Procognitive restoration of PV neuron plasticity in neurodevelopmental disorders.

Code ↔ Paper

3 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 3 matches

[1] § Methods › Ripple analysis ↔ Ripples_Z-scored_Single Trial Script.py, lines 877–938 · score 0.91 · detected ripple event, peak LFP amplitude, bandpass filtered, ripple band, peak power, CA1 LFP
[2] § Methods › Extracting NREM periods ↔ Ripples_Z-scored_Single Trial Script.py, lines 99–102 · score 0.76 · 2–16 Hz, 0–300 Hz, theta dominance, 5–10 Hz, sleep, power
[3] § Extended Data ↔ Ripples_Z-scored_Single Trial Script.py, lines 110–114 · score 0.68 · 100–250 Hz, ripple duration, ripple frequency, envelope, peak, filtered

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 paper is loaded by your browser, with JavaScript: read it at doi.org or on Europe PMC.

The authors' code

Python · 938 lines · 52 KB · CC-BY-4.0 · 3 matches

# -*- coding: utf-8 -*-
"""
Created on Tue May 20 16:30:53 2025
@author: HT_bo
"""
import numpy as np
import scipy.io
import scipy.signal
from scipy.stats import zscore
import pandas as pd
import matplotlib.pyplot as plt
# from open_ephys.analysis import Session # Bypassing for manual load
import os
from pathlib import Path
import sys
import json # For reading structure.oebin
# --- Add Tkinter for file dialogs ---
import tkinter as tk
from tkinter import filedialog
# ------------------------------------------------------------------------------
# Configuration
# ------------------------------------------------------------------------------
# --- Setup Tkinter for file dialogs (main window won't be shown) ---
root_tk = tk.Tk()
root_tk.withdraw() # Hide the main tkinter window
# --- Prompt user for input paths ---
print("Please select the main recording session directory (e.g., '.../experiment1/recording1')")
session_dir_str = filedialog.askdirectory(title="Select Recording Session Directory (e.g., .../experiment1/recording1)")
if not session_dir_str:
print("No session directory selected. Exiting.")
sys.exit()
session_dir = Path(session_dir_str)
print(f"Selected session directory: {session_dir}")
# --- Extract path part for filenames ---
try:
# Assuming session_dir is like .../Record Node 101/experiment1/recording1
# We want the name of the folder 3 levels up (e.g., "WT,tdTomato,264,pre")
path_part_for_filename = session_dir.parent.parent.parent.name
if not path_part_for_filename:
path_part_for_filename = session_dir.name # Fallback 1
print(f"Using path part for filenames: {path_part_for_filename}")
except AttributeError:
print("Warning: Could not derive desired path part from 3 levels up. Using recording name as fallback.")
path_part_for_filename = session_dir.name # Fallback 2
except Exception as e:
print(f"Warning: Error deriving path part: {e}. Using recording name as fallback.")
path_part_for_filename = session_dir.name
try:
RECORDING_NAME = session_dir.name # e.g., "recording1"
EXPERIMENT_NAME = session_dir.parent.name # e.g., "experiment1"
except Exception: # Fallback if path is too short
RECORDING_NAME = path_part_for_filename # If session_dir itself was the base element for path_part
EXPERIMENT_NAME = "experiment"
print("Please select the MAT file containing electrode mapping (e.g., ..._Behavior_and_Optogenetics_TimeStamps.mat)")
mat_file_path_str = filedialog.askopenfilename(
title="Select Electrode Mapping MAT File",
filetypes=(("MAT files", "*.mat"), ("All files", "*.*"))
)
if not mat_file_path_str:
print("No MAT file selected. Exiting.")
sys.exit()
MAT_FILE_PATH = Path(mat_file_path_str)
print(f"Selected MAT file: {MAT_FILE_PATH}")
OUTPUT_DIR_BASE = session_dir.resolve() # Output will be inside the selected session_dir
os.makedirs(OUTPUT_DIR_BASE, exist_ok=True)
print(f"Output base directory: {OUTPUT_DIR_BASE}")
# --- Channel & Region Mapping ---
CA1_REGION_IDS = [281, 282, 283, 284] # Example: Your Region IDs for CA1 areas from MAT file
DG_REGION_ID = 261 # Example: Your Region ID for the DG noise channel from MAT file. Set to None if no DG/noise channel.
# !!! IMPORTANT: OFFSET FOR LOGICAL TO PHYSICAL CHANNEL MAPPING !!!
# Based on: "lfp channels were recorded from 9 to 24 which would be considered Channel 1 -16 in timestamps.mat file"
# This means logical channel 1 (from MAT file) is physical channel 9 (CH9 in structure.oebin).
# So, physical_channel = logical_channel + LOGICAL_TO_PHYSICAL_OFFSET
LOGICAL_TO_PHYSICAL_OFFSET = 8
# --- LFP Processing Parameters ---
TARGET_FS_DOWNSAMPLED = 1000.0 # Hz - Desired sampling rate after decimation
# --- Sleep State Detection Parameters ---
SPECTROGRAM_WINDOW_SEC = 10.0 # s
SPECTROGRAM_STEP_SEC = 1.0 # s
SCORING_BUFFER_SEC = 5 # s, min duration to confirm state change
# --- Frequency Bands for Sleep Scoring Features ---
PCA_SPECTROGRAM_FREQ_RANGE = (0, 300) # Hz, for PCA input
THETA_DOMINANCE_THETA_BAND = (5, 10) # Hz
THETA_DOMINANCE_TOTAL_BAND = (2, 16) # Hz (for normalization of theta power)
# --- NREM Estimation Options ---
USE_DT_RATIO_FOR_NREM = True # <--- Set to True to use Delta/Theta ratio for NREM detection
DT_RATIO_THRESHOLD = 1.5 # Threshold for Delta/Theta ratio
DT_DELTA_BAND = (0.5, 4) # Hz, Delta band for ratio calculation
DT_THETA_BAND = (5, 10) # Hz, Theta band for ratio calculation
# --- Ripple Detection Parameters ---
RIPPLE_THRESHOLDS = (2, 5) # (low_z_power, high_z_power_peak)
RIPPLE_DURATIONS = (20, 20, 150) # (min_isi_ms, min_dur_ms, max_dur_ms)
RIPPLE_FREQ_RANGE = (100, 250) # For ripple *detection*
RIPPLE_ENVELOPE_FILTER_HZ = (1, 20)
# --- Plotting Parameters ---
RIPPLE_PLOT_WINDOW_MS = 50 # Window around ripple peak for spectrogram/PSD (+/- ms)
SPECTROGRAM_CMAP = 'viridis'
# ------------------------------------------------------------------------------
# Helper Functions
# ------------------------------------------------------------------------------
def load_electrode_mapping(mat_file_path_local):
"""Loads ElectrodeVsRegisteredAreasNum from a .mat file."""
try:
mat_data = scipy.io.loadmat(mat_file_path_local)
if 'ElectrodeVsRegisteredAreasNum' in mat_data:
mapping = mat_data['ElectrodeVsRegisteredAreasNum']
return pd.DataFrame(mapping, columns=['Channel', 'RegionID'])
else:
print(f"Error: 'ElectrodeVsRegisteredAreasNum' not found in {mat_file_path_local}")
print(f"Available keys: {list(mat_data.keys())}")
return None
except FileNotFoundError:
print(f"Error: MAT file not found at {mat_file_path_local}")
return None
except Exception as e:
print(f"Error loading MAT file: {e}")
return None
def butter_bandpass_filter(data, lowcut, highcut, fs, order=3, axis=0):
nyq = 0.5 * fs
low = lowcut / nyq
high = highcut / nyq
if low <= 0: low = 1e-6
if high >= 1: high = 1 - 1e-6
if low >= high:
if lowcut == 0 and highcut > 0 and highcut < nyq :
b, a = scipy.signal.butter(order, high, btype='lowpass')
elif lowcut > 0 and lowcut < nyq and highcut >= nyq * 0.999 :
b, a = scipy.signal.butter(order, low, btype='highpass')
else:
print(f" Problematic band definition for bandpass: low {lowcut}, high {highcut}. Returning data copy.")
return data.copy()
else:
b, a = scipy.signal.butter(order, [low, high], btype='band')
try:
y = scipy.signal.filtfilt(b, a, data.astype(np.float64), axis=axis)
return y
except ValueError as e:
print(f"Filter error in butter_bandpass_filter: {e} with low={low}, high={high}. Returning unfiltered data.")
return data.copy()
def compute_spectrogram_custom(data, fs, window_size_sec, step_size_sec):
print(f"Computing spectrogram: fs={fs:.2f}, window={window_size_sec}s, step={step_size_sec}s")
if data is None or data.ndim != 1 or len(data) == 0:
print("Error: Spectrogram input data must be 1D and not empty.")
return None, None, None
window_samples = int(round(window_size_sec * fs))
step_samples = int(round(step_size_sec * fs))
if window_samples <= 0 or step_samples <= 0:
print("Error: Window or step size results in non-positive samples.")
return None, None, None
if window_samples > len(data):
# print(f"Warning: Window size ({window_samples}) > data length ({len(data)}). Adjusting window to data length.")
window_samples = len(data) # Adjust window if too long
noverlap = window_samples - step_samples
noverlap = max(0, noverlap)
if noverlap >= window_samples and window_samples > 0 :
noverlap = window_samples - 1
elif noverlap >= window_samples and window_samples == 0:
print("Error: window_samples is 0 in spectrogram.")
return None, None, None
frequencies, times, Sxx = scipy.signal.spectrogram(
data.astype(np.float64), fs=fs, window='hann',
nperseg=window_samples, noverlap=noverlap,
scaling='density', mode='psd'
)
# print(f"Spectrogram computed. Shape: {Sxx.shape}")
return Sxx, frequencies, times
def compute_pca_custom(spectrogram_data_for_pca):
# print(f"Computing PCA on spectrogram of shape: {spectrogram_data_for_pca.shape}")
if spectrogram_data_for_pca is None or spectrogram_data_for_pca.shape[0] < 2 or spectrogram_data_for_pca.shape[1] < 2:
print("Error: Invalid data for PCA (None, empty, or too small).")
return None
data_for_pca = spectrogram_data_for_pca.T
zscored_data = zscore(data_for_pca, axis=0, nan_policy='omit')
zscored_data = np.nan_to_num(zscored_data, nan=0.0, posinf=np.nanmax(zscored_data[np.isfinite(zscored_data)]) if np.any(np.isfinite(zscored_data)) else 0.0, neginf=np.nanmin(zscored_data[np.isfinite(zscored_data)]) if np.any(np.isfinite(zscored_data)) else 0.0)
from sklearn.decomposition import PCA
pca = PCA(n_components=1)
try:
pc1 = pca.fit_transform(zscored_data)
except ValueError as e:
print(f"Error during PCA fit_transform: {e}")
# print(" Data for PCA (first 5 rows after nan_to_num): \n", zscored_data[:5,:])
return None
# print(f"PCA computed. PC1 shape: {pc1.shape}, Explained var: {pca.explained_variance_ratio_[0]:.4f}")
return pc1.squeeze()
def compute_theta_dominance_custom(spectrogram_data, frequencies, theta_band_config, total_band_config):
# print(f"Computing Theta Dominance. Spec shape: {spectrogram_data.shape}")
if spectrogram_data is None or frequencies is None or spectrogram_data.size == 0 or frequencies.size == 0:
print("Error: Invalid input for theta dominance calculation.")
return None
theta_band_indices = np.where((frequencies >= theta_band_config[0]) & (frequencies <= theta_band_config[1]))[0]
total_band_indices = np.where((frequencies >= total_band_config[0]) & (frequencies <= total_band_config[1]))[0]
if len(theta_band_indices) == 0 or len(total_band_indices) == 0:
print(f"Warning: Theta band ({theta_band_config} Hz) or total power band ({total_band_config} Hz) not found in frequencies for theta dominance.")
return np.full(spectrogram_data.shape[1], np.nan)
epsilon = 1e-12
theta_power = np.nanmean(spectrogram_data[theta_band_indices, :], axis=0)
total_power = np.nanmean(spectrogram_data[total_band_indices, :], axis=0)
theta_dominance = np.full_like(theta_power, np.nan)
valid_mask = (total_power > epsilon) & np.isfinite(theta_power) & np.isfinite(total_power)
theta_dominance[valid_mask] = theta_power[valid_mask] / total_power[valid_mask]
return theta_dominance
def score_sleep_states_custom(pc1, theta_dominance, times, step_size_sec, buffer_sec):
print("\n--- Scoring Sleep States (LFP-based only) ---")
if pc1 is None or theta_dominance is None or times is None:
print("Error: Missing PC1, Theta Dominance, or Times for scoring.")
return None, {}
num_time_points = len(pc1)
if not (len(theta_dominance) == num_time_points and len(times) == num_time_points):
print(f"Error: Length mismatch in scoring inputs. PC1:{len(pc1)}, Theta:{len(theta_dominance)}, Times:{len(times)}")
return None, {}
pc1_finite = pc1[np.isfinite(pc1)]
theta_finite = theta_dominance[np.isfinite(theta_dominance)]
nrem_threshold_pc1 = np.nanpercentile(pc1_finite, 75) if len(pc1_finite) > 0 else np.nan
rem_threshold_theta = np.nanpercentile(theta_finite, 75) if len(theta_finite) > 0 else np.nan
if np.isnan(nrem_threshold_pc1) or np.isnan(rem_threshold_theta):
print("Error: Could not calculate NREM/REM thresholds (NaN).")
print(f" PC1 finite len: {len(pc1_finite)}, Theta finite len: {len(theta_finite)}")
return None, {}
print(f"NREM Threshold (PC1 >): {nrem_threshold_pc1:.3f}")
print(f"REM Threshold (Theta Dominance >): {rem_threshold_theta:.3f}")
sleep_states = np.zeros(num_time_points, dtype=int)
current_state = 0
buffer_samples = int(round(buffer_sec / step_size_sec))
if buffer_samples < 1: buffer_samples = 1
nrem_counter, rem_counter, awake_counter = 0, 0, 0
for i in range(num_time_points):
pc1_val, theta_val = pc1[i], theta_dominance[i]
is_nrem_like = np.isfinite(pc1_val) and pc1_val > nrem_threshold_pc1
is_rem_like = np.isfinite(theta_val) and theta_val > rem_threshold_theta
if current_state == 0: # Awake
if is_nrem_like: nrem_counter += 1; rem_counter = 0; awake_counter=0
elif is_rem_like: rem_counter += 1; nrem_counter = 0; awake_counter=0
else: nrem_counter = 0; rem_counter = 0;
if nrem_counter >= buffer_samples: current_state = 1; nrem_counter=0; rem_counter=0; awake_counter=0
elif rem_counter >= buffer_samples: current_state = 2; nrem_counter=0; rem_counter=0; awake_counter=0
elif current_state == 1: # NREM
if is_rem_like and (not is_nrem_like): rem_counter +=1; nrem_counter=0; awake_counter = 0
elif not is_nrem_like : awake_counter += 1; rem_counter = 0; nrem_counter=0;
else: nrem_counter +=1; awake_counter = 0; rem_counter = 0
if rem_counter >= buffer_samples : current_state = 2; rem_counter=0; nrem_counter=0; awake_counter=0
elif awake_counter >= buffer_samples : current_state = 0; awake_counter=0; nrem_counter=0; rem_counter=0
elif current_state == 2: # REM
if not is_rem_like: awake_counter += 1; rem_counter=0; nrem_counter=0;
else: rem_counter +=1; awake_counter = 0; nrem_counter=0;
if awake_counter >= buffer_samples: current_state = 0; awake_counter=0; rem_counter=0; nrem_counter=0
sleep_states[i] = current_state
print("Sleep scoring complete.")
thresholds_used = {'nrem_pc1': nrem_threshold_pc1, 'rem_theta': rem_threshold_theta}
return sleep_states, thresholds_used
def find_swr_custom(lfp, timestamps, fs, thresholds, durations, freq_range, envelope_filter_hz, noise_lfp=None):
if lfp is None or len(lfp) == 0:
# print("Warning: Empty LFP segment passed to find_swr_custom.")
return pd.DataFrame()
low_thresh, high_thresh = thresholds
min_isi_ms, min_dur_ms, max_dur_ms = durations
lfp_float64 = lfp.astype(np.float64)
noise_lfp_float64 = noise_lfp.astype(np.float64) if noise_lfp is not None else None
filtered_lfp = butter_bandpass_filter(lfp_float64, freq_range[0], freq_range[1], fs, order=3)
rectified = filtered_lfp ** 2
envelope_raw = butter_bandpass_filter(rectified, envelope_filter_hz[0], envelope_filter_hz[1], fs, order=3)
if len(envelope_raw) > 0 and not np.all(np.isnan(envelope_raw)):
median_env = np.nanmedian(envelope_raw)
mad_env = np.nanmedian(np.abs(envelope_raw - median_env))
if mad_env == 0 or np.isnan(mad_env): mad_env = 1e-9
envelope_z = 0.6745 * (envelope_raw - median_env) / mad_env
else:
return pd.DataFrame()
above_thresh = envelope_z > low_thresh
rising = np.where(np.diff(above_thresh.astype(int)) == 1)[0] + 1
falling = np.where(np.diff(above_thresh.astype(int)) == -1)[0] + 1
if len(rising) == 0 or len(falling) == 0: return pd.DataFrame()
if falling[0] < rising[0]:
first_valid_falling_idx = np.searchsorted(falling, rising[0])
if first_valid_falling_idx == len(falling): return pd.DataFrame()
falling = falling[first_valid_falling_idx:]
if len(rising) == 0 or len(falling) == 0: return pd.DataFrame()
if rising[-1] > falling[-1]:
last_valid_rising_idx = np.searchsorted(rising, falling[-1], side='right')
if last_valid_rising_idx == 0 : return pd.DataFrame()
rising = rising[:last_valid_rising_idx]
min_len_edges = min(len(rising), len(falling))
rising, falling = rising[:min_len_edges], falling[:min_len_edges]
if len(rising) == 0: return pd.DataFrame()
events = np.column_stack((rising, falling))
valid_event_mask = events[:,1] > events[:,0]
events = events[valid_event_mask]
if events.shape[0] == 0: return pd.DataFrame()
merged_events = []
if len(events) > 0:
current_event = list(events[0])
min_isi_samples = int(min_isi_ms / 1000 * fs)
for next_start, next_end in events[1:]:
if current_event[1] >= next_start :
current_event[1] = max(current_event[1], next_end)
elif (next_start - current_event[1]) < min_isi_samples:
current_event[1] = next_end
else:
if current_event[1] > current_event[0]: merged_events.append(current_event)
current_event = [next_start, next_end]
if current_event[1] > current_event[0]: merged_events.append(current_event)
if not merged_events: return pd.DataFrame()
merged_events = np.array(merged_events)
final_ripples = []
z_filtered_lfp_for_amplitude = np.array([])
if len(filtered_lfp) > 0 and not np.all(np.isnan(filtered_lfp)):
try:
z_filtered_lfp_for_amplitude = zscore(filtered_lfp, nan_policy='omit')
z_filtered_lfp_for_amplitude = np.nan_to_num(z_filtered_lfp_for_amplitude)
except Exception as e_zscore_filt:
print(f"Warning: zscore failed for filtered_lfp: {e_zscore_filt}")
z_filtered_lfp_for_amplitude = filtered_lfp # Fallback to non-zscored if error
noise_envelope_z = None
if noise_lfp_float64 is not None and len(noise_lfp_float64) == len(lfp_float64):
noise_filtered = butter_bandpass_filter(noise_lfp_float64, freq_range[0], freq_range[1], fs, order=3)
noise_rectified = noise_filtered**2
noise_envelope_raw = butter_bandpass_filter(noise_rectified, envelope_filter_hz[0], envelope_filter_hz[1], fs, order=3)
if len(noise_envelope_raw)>0 and not np.all(np.isnan(noise_envelope_raw)):
median_noise_env = np.nanmedian(noise_envelope_raw)
mad_noise_env = np.nanmedian(np.abs(noise_envelope_raw - median_noise_env))
if mad_noise_env == 0 or np.isnan(mad_noise_env): mad_noise_env = 1e-9
noise_envelope_z = 0.6745 * (noise_envelope_raw - median_noise_env) / mad_noise_env
else: noise_envelope_z = None
for start_idx, end_idx in merged_events:
if start_idx >= end_idx: continue
segment_envelope_z = envelope_z[start_idx:end_idx]
if len(segment_envelope_z) == 0: continue
max_power_z_in_segment = np.max(segment_envelope_z)
if max_power_z_in_segment >= high_thresh:
peak_power_idx_in_segment = np.argmax(segment_envelope_z)
peak_sample_abs_power = start_idx + peak_power_idx_in_segment
duration_s = (end_idx - start_idx) / fs
if not (min_dur_ms / 1000 <= duration_s <= max_dur_ms / 1000):
continue
if noise_envelope_z is not None:
if end_idx > len(noise_envelope_z):
pass # print(f"Warning: Ripple end {end_idx} > noise_env len {len(noise_envelope_z)}. Skip noise check.")
elif np.any(noise_envelope_z[start_idx:end_idx] > high_thresh):
continue
peak_val_z_lfp_amp = np.nan
peak_sample_abs_lfp = np.nan
if len(z_filtered_lfp_for_amplitude) > 0 and end_idx <= len(z_filtered_lfp_for_amplitude) and start_idx < len(z_filtered_lfp_for_amplitude):
segment_z_lfp_amp = z_filtered_lfp_for_amplitude[start_idx:end_idx]
if len(segment_z_lfp_amp) > 0:
abs_max_idx_in_segment_amp = np.argmax(np.abs(segment_z_lfp_amp))
peak_val_z_lfp_amp = segment_z_lfp_amp[abs_max_idx_in_segment_amp]
peak_sample_abs_lfp = start_idx + abs_max_idx_in_segment_amp
current_ts_start = timestamps[start_idx] if start_idx < len(timestamps) else np.nan
current_ts_end = timestamps[end_idx-1] if end_idx > 0 and end_idx-1 < len(timestamps) else np.nan
current_ts_peak_power = timestamps[peak_sample_abs_power] if peak_sample_abs_power < len(timestamps) else np.nan
current_ts_peak_lfp = timestamps[int(peak_sample_abs_lfp)] if pd.notna(peak_sample_abs_lfp) and int(peak_sample_abs_lfp) < len(timestamps) else np.nan
final_ripples.append({
'start_sample': start_idx, 'end_sample': end_idx,
'peak_sample_power': peak_sample_abs_power,
'peak_power_zscore': max_power_z_in_segment,
'peak_sample_lfp': peak_sample_abs_lfp,
'peak_lfp_amplitude_zscore': peak_val_z_lfp_amp,
'start_time': current_ts_start, 'end_time': current_ts_end,
'peak_time_power': current_ts_peak_power, 'peak_time_lfp': current_ts_peak_lfp
})
return pd.DataFrame(final_ripples)
def plot_ripple_details(lfp_ca1_avg_for_plot, fs_for_plot, ripple_events_df_for_plot,
window_ms, ripple_band_for_plot, output_dir_plot, session_name_plot):
if ripple_events_df_for_plot.empty:
print("No ripples to plot.")
return
print(f"Generating ripple-triggered plots for {len(ripple_events_df_for_plot)} events...")
window_plot_samples = int(window_ms * fs_for_plot / 1000)
all_ripple_spectrograms, all_ripple_psds = [], []
valid_ripple_count_for_plot = 0
representative_freqs_spec_rip, representative_times_spec_rip_centered, representative_spec_freq_mask_rip = None, None, None
freqs_psd_rip_for_plot = None
for idx, ripple in ripple_events_df_for_plot.iterrows():
peak_idx_abs = ripple['peak_sample_lfp']
if pd.isna(peak_idx_abs): continue
peak_idx_abs = int(peak_idx_abs)
start_plot = peak_idx_abs - window_plot_samples
end_plot = peak_idx_abs + window_plot_samples
if start_plot < 0 or end_plot >= len(lfp_ca1_avg_for_plot): continue
segment_for_plot = lfp_ca1_avg_for_plot[start_plot:end_plot].astype(np.float64)
if len(segment_for_plot) != 2 * window_plot_samples: continue
valid_ripple_count_for_plot += 1
nperseg_spec = min(len(segment_for_plot), max(32, int(fs_for_plot / ripple_band_for_plot[0] * 2.5)))
noverlap_spec = nperseg_spec // 2
if nperseg_spec <= noverlap_spec : noverlap_spec = max(0, nperseg_spec -1)
if nperseg_spec == 0 : continue
try:
current_freqs_spec_rip, current_times_spec_rip, Sxx_rip = scipy.signal.spectrogram(
segment_for_plot, fs=fs_for_plot, window='hann', nperseg=nperseg_spec, noverlap=noverlap_spec,
scaling='density', mode='psd'
)
current_spec_freq_mask_rip = (current_freqs_spec_rip >= ripple_band_for_plot[0]) & (current_freqs_spec_rip <= ripple_band_for_plot[1])
if np.any(current_spec_freq_mask_rip) and Sxx_rip[current_spec_freq_mask_rip, :].size > 0 :
all_ripple_spectrograms.append(10 * np.log10(Sxx_rip[current_spec_freq_mask_rip, :] + 1e-12))
if representative_freqs_spec_rip is None:
representative_freqs_spec_rip = current_freqs_spec_rip
representative_spec_freq_mask_rip = current_spec_freq_mask_rip
representative_times_spec_rip_centered = (current_times_spec_rip - current_times_spec_rip.mean()) * 1000
except ValueError as e:
continue
nperseg_psd_rip = min(len(segment_for_plot), 256)
current_freqs_psd_rip, Pxx_rip = scipy.signal.welch(segment_for_plot, fs=fs_for_plot, nperseg=nperseg_psd_rip, scaling='density')
all_ripple_psds.append(Pxx_rip)
if freqs_psd_rip_for_plot is None :
freqs_psd_rip_for_plot = current_freqs_psd_rip
if not all_ripple_spectrograms or not all_ripple_psds:
print("Not enough valid ripple data for average plots after processing segments.")
return
print(f" Aggregating {len(all_ripple_spectrograms)} spectrograms and {len(all_ripple_psds)} PSDs for averaging.")
min_time_bins_spec = min(s.shape[1] for s in all_ripple_spectrograms) if all_ripple_spectrograms else 0
if min_time_bins_spec == 0:
print(" Spectrograms have inconsistent time bins or are empty. Cannot average for plot.")
return
all_ripple_spectrograms_trimmed = [s[:, :min_time_bins_spec] for s in all_ripple_spectrograms]
avg_ripple_spectrogram = np.mean(all_ripple_spectrograms_trimmed, axis=0)
avg_ripple_psd = np.mean(all_ripple_psds, axis=0)
avg_ripple_psd_db = 10 * np.log10(avg_ripple_psd + 1e-12)
avg_ripple_psd_zscore = zscore(avg_ripple_psd_db, nan_policy='omit') if len(avg_ripple_psd_db) > 1 else avg_ripple_psd_db
avg_ripple_psd_zscore = np.nan_to_num(avg_ripple_psd_zscore)
fig_rip, axes_rip = plt.subplots(2, 1, figsize=(10, 8))
plot_freqs_spec_final = representative_freqs_spec_rip[representative_spec_freq_mask_rip] if representative_freqs_spec_rip is not None and representative_spec_freq_mask_rip is not None and np.any(representative_spec_freq_mask_rip) and len(representative_freqs_spec_rip[representative_spec_freq_mask_rip]) == avg_ripple_spectrogram.shape[0] else np.linspace(ripple_band_for_plot[0],ripple_band_for_plot[1], avg_ripple_spectrogram.shape[0])
plot_times_spec_final = representative_times_spec_rip_centered[:min_time_bins_spec] if representative_times_spec_rip_centered is not None and len(representative_times_spec_rip_centered) >= min_time_bins_spec else np.linspace(-window_ms, window_ms, min_time_bins_spec)
if plot_freqs_spec_final.shape[0] != avg_ripple_spectrogram.shape[0] :
plot_freqs_spec_final = np.linspace(ripple_band_for_plot[0], ripple_band_for_plot[1], avg_ripple_spectrogram.shape[0])
if plot_times_spec_final.shape[0] != avg_ripple_spectrogram.shape[1] :
plot_times_spec_final = np.linspace(-window_ms, window_ms, avg_ripple_spectrogram.shape[1])
im = axes_rip[0].pcolormesh(plot_times_spec_final, plot_freqs_spec_final, avg_ripple_spectrogram,
shading='gouraud', cmap=SPECTROGRAM_CMAP)
axes_rip[0].set_ylabel(f'Frequency ({ripple_band_for_plot[0]}-{ripple_band_for_plot[1]} Hz)')
axes_rip[0].set_xlabel('Time from LFP peak (ms)')
axes_rip[0].set_title(f'Average Ripple Spectrogram (N={valid_ripple_count_for_plot})')
axes_rip[0].axvline(0, color='r', linestyle='--', alpha=0.7)
fig_rip.colorbar(im, ax=axes_rip[0], label='Power (dB/Hz)')
if freqs_psd_rip_for_plot is None and all_ripple_psds:
freqs_psd_rip_for_plot = np.linspace(0, fs_for_plot/2, len(avg_ripple_psd_zscore))
psd_plot_freq_mask = (freqs_psd_rip_for_plot >= ripple_band_for_plot[0]-20) & (freqs_psd_rip_for_plot <= ripple_band_for_plot[1]+20)
axes_rip[1].plot(freqs_psd_rip_for_plot[psd_plot_freq_mask], avg_ripple_psd_zscore[psd_plot_freq_mask])
axes_rip[1].set_xlabel('Frequency (Hz)')
axes_rip[1].set_ylabel('Z-scored Avg Power (dB/Hz)')
axes_rip[1].set_title('Average Ripple PSD (Z-scored)')
axes_rip[1].grid(True, which="both", ls="-", alpha=0.5)
plt.tight_layout()
plot_filename = Path(output_dir_plot) / f"{session_name_plot}_avg_ripple_plots_manual.png"
try:
plt.savefig(plot_filename)
print(f"Saved average ripple plots to {plot_filename}")
except Exception as e_save_plot:
print(f"Error saving ripple plot: {e_save_plot}")
plt.close(fig_rip)
# ------------------------------------------------------------------------------
# Main Execution
# ------------------------------------------------------------------------------
if __name__ == "__main__":
output_dir_specific = Path(OUTPUT_DIR_BASE) / "ripple_analysis_output_py_manual_ds"
output_dir_specific.mkdir(parents=True, exist_ok=True)
session_name_for_files = f"{path_part_for_filename}_{EXPERIMENT_NAME}_{RECORDING_NAME}"
print(f"Starting analysis for: {session_dir}")
print(f"Base for output filenames: {session_name_for_files}")
electrode_map_df = load_electrode_mapping(MAT_FILE_PATH)
if electrode_map_df is None:
sys.exit()
logical_ca1_channels_1based = []
if isinstance(CA1_REGION_IDS, list):
for region_id in CA1_REGION_IDS:
channels_for_region = electrode_map_df[electrode_map_df['RegionID'] == region_id]['Channel'].astype(int).tolist()
logical_ca1_channels_1based.extend(channels_for_region)
elif isinstance(CA1_REGION_IDS, int):
logical_ca1_channels_1based = electrode_map_df[electrode_map_df['RegionID'] == CA1_REGION_IDS]['Channel'].astype(int).tolist()
else:
print(f"Error: CA1_REGION_IDS type. Value: {CA1_REGION_IDS}");
sys.exit()
logical_ca1_channels_1based = sorted(list(set(logical_ca1_channels_1based)))
if not logical_ca1_channels_1based:
print(f"Error: No LOGICAL CA1 channels for RegionIDs {CA1_REGION_IDS}.");
sys.exit()
print(f"Identified LOGICAL CA1 channels (1-16 mapping): {logical_ca1_channels_1based}")
logical_dg_channel_1based_list = []
if DG_REGION_ID is not None:
dg_channels_df = electrode_map_df[electrode_map_df['RegionID'] == DG_REGION_ID]
if not dg_channels_df.empty:
logical_dg_channel_1based_list = [dg_channels_df['Channel'].astype(int).iloc[0]]
print(f"Identified LOGICAL DG noise channel (1-16 mapping): {logical_dg_channel_1based_list[0]}")
else: print(f"Warning: No LOGICAL DG channels for RegionID {DG_REGION_ID}.")
else: print("DG_REGION_ID not set for noise channel.")
physical_ca1_channels_1based = [lc + LOGICAL_TO_PHYSICAL_OFFSET for lc in logical_ca1_channels_1based]
physical_dg_channel_1based_list = [ldc + LOGICAL_TO_PHYSICAL_OFFSET for ldc in logical_dg_channel_1based_list] if logical_dg_channel_1based_list else []
print(f"Translated to PHYSICAL CA1 channels (for CH names in oebin): {physical_ca1_channels_1based}")
if physical_dg_channel_1based_list:
print(f"Translated to PHYSICAL DG noise channel: {physical_dg_channel_1based_list[0]}")
print("\n--- Starting Manual Data Loading via structure.oebin ---")
lfp_data_processed_scaled = None
original_fs_from_oebin = None
lfp_timestamps_loaded = None
loaded_channel_physical_ids_in_order = []
structure_file_path = session_dir / "structure.oebin"
metadata_oebin = None
if not structure_file_path.exists():
print(f"ERROR: structure.oebin not found at {structure_file_path}");
sys.exit()
try:
with open(structure_file_path, 'r') as f: metadata_oebin = json.load(f)
print("Successfully parsed structure.oebin")
except Exception as e:
print(f"ERROR: Could not parse structure.oebin: {e}");
sys.exit()
if metadata_oebin and metadata_oebin.get('continuous') and len(metadata_oebin['continuous']) > 0:
continuous_stream_info = metadata_oebin['continuous'][0]
stream_folder_name = continuous_stream_info.get('folder_name')
original_fs_from_oebin = float(continuous_stream_info.get('sample_rate', 0))
num_channels_total_in_stream = int(continuous_stream_info.get('num_channels', 0))
channels_metadata_list_oebin = continuous_stream_info.get('channels', [])
print(f" Stream Folder: {stream_folder_name}, Original SR: {original_fs_from_oebin} Hz, Total Stream Ch: {num_channels_total_in_stream}")
if not all([stream_folder_name, original_fs_from_oebin > 0, num_channels_total_in_stream > 0]):
print("ERROR: Essential stream info missing from structure.oebin."); sys.exit()
continuous_dat_path = session_dir / "continuous" / stream_folder_name / "continuous.dat"
timestamps_npy_path = session_dir / "continuous" / stream_folder_name / "timestamps.npy"
print(f" Expected continuous.dat: {continuous_dat_path}")
if not continuous_dat_path.exists():
print("ERROR: continuous.dat not found.");
sys.exit()
if not timestamps_npy_path.exists():
print("ERROR: timestamps.npy not found.");
sys.exit()
target_physical_channels_to_extract = sorted(list(set(physical_ca1_channels_1based + (physical_dg_channel_1based_list if physical_dg_channel_1based_list else []))))
print(f" Attempting to extract PHYSICAL 1-based channels: {target_physical_channels_to_extract}")
stream_indices_to_load_0based = []
bit_volts_for_selected_channels = []
num_amplifier_channels_in_stream = sum(1 for ch_meta in channels_metadata_list_oebin if ch_meta.get('channel_name','').upper().startswith("CH"))
for target_ch_physical_num_1based in target_physical_channels_to_extract:
found_in_oebin = False
for oebin_idx_0based, oebin_ch_meta in enumerate(channels_metadata_list_oebin):
oebin_ch_name = oebin_ch_meta.get('channel_name', '').upper()
num_part_str = ''.join(filter(str.isdigit, oebin_ch_name))
if not num_part_str: continue
oebin_ch_num_part = int(num_part_str)
current_ch_matches_target = False
if oebin_ch_name.startswith("CH") and oebin_ch_num_part == target_ch_physical_num_1based:
current_ch_matches_target = True
elif oebin_ch_name.startswith("ADC"):
if (num_amplifier_channels_in_stream + oebin_ch_num_part) == target_ch_physical_num_1based:
current_ch_matches_target = True
if current_ch_matches_target:
stream_indices_to_load_0based.append(oebin_idx_0based)
bit_v = float(oebin_ch_meta.get('bit_volts'))
units = oebin_ch_meta.get('units', '').upper()
if "V" in units and "UV" not in units and abs(bit_v) < 1:
bit_v *= 1e6
bit_volts_for_selected_channels.append(bit_v)
loaded_channel_physical_ids_in_order.append(target_ch_physical_num_1based)
found_in_oebin = True
break
if not found_in_oebin:
print(f"Warning: Target PHYSICAL channel {target_ch_physical_num_1based} not matched in structure.oebin.")
if not stream_indices_to_load_0based:
print("ERROR: No channels mapped for loading.");
sys.exit()
print(f" Mapped 0-based stream indices: {stream_indices_to_load_0based}")
print(f" Corresponding PHYSICAL 1-based IDs loaded: {loaded_channel_physical_ids_in_order}")
try:
raw_data_memmap = np.memmap(continuous_dat_path, dtype='int16', mode='r')
num_total_samples_in_file = len(raw_data_memmap) // num_channels_total_in_stream
valid_length = num_total_samples_in_file * num_channels_total_in_stream
all_channels_data_reshaped = raw_data_memmap[:valid_length].reshape((num_total_samples_in_file,
num_channels_total_in_stream))
lfp_data_processed_scaled = np.zeros((num_total_samples_in_file,
len(stream_indices_to_load_0based)), dtype=np.float32)
for i, (stream_idx, bit_v) in enumerate(zip(stream_indices_to_load_0based, bit_volts_for_selected_channels)):
lfp_data_processed_scaled[:, i] = all_channels_data_reshaped[:, stream_idx].astype(np.float32) * bit_v
del raw_data_memmap, all_channels_data_reshaped
except Exception as e:
print(f"Error reading/processing continuous.dat: {e}");
sys.exit()
try:
lfp_timestamps_loaded = np.load(timestamps_npy_path)
if len(lfp_timestamps_loaded) == 1:
print("Timestamps.npy has one entry; assuming start time and reconstructing.")
start_time_abs = lfp_timestamps_loaded[0]
lfp_timestamps_loaded = start_time_abs + np.arange(num_total_samples_in_file) / original_fs_from_oebin
elif len(lfp_timestamps_loaded) != num_total_samples_in_file:
print(f"Warning: Timestamps ({len(lfp_timestamps_loaded)}) != samples ({num_total_samples_in_file}). Adjusting.")
min_len_ts_data = min(len(lfp_timestamps_loaded), num_total_samples_in_file)
lfp_timestamps_loaded = lfp_timestamps_loaded[:min_len_ts_data]
lfp_data_processed_scaled = lfp_data_processed_scaled[:min_len_ts_data, :]
num_total_samples_in_file = min_len_ts_data
print(f" Adjusted data/timestamps to min length: {min_len_ts_data}")
except Exception as e:
print(f"Error loading timestamps.npy: {e}");
sys.exit()
print("Manual data loading complete.")
else:
print("ERROR: No 'continuous' stream info in structure.oebin.");
sys.exit()
if lfp_data_processed_scaled is None or original_fs_from_oebin is None or lfp_timestamps_loaded is None:
print("Exiting: Manual LFP data loading failed.");
sys.exit()
print(f"LFP data (manual raw). Shape: {lfp_data_processed_scaled.shape}, Original SR: {original_fs_from_oebin} Hz, Timestamps: {len(lfp_timestamps_loaded)}")
ca1_indices_in_loaded = [i for i, ph_id in
enumerate(loaded_channel_physical_ids_in_order)
if ph_id in physical_ca1_channels_1based]
if not ca1_indices_in_loaded:
print("ERROR: CA1 channels not found in manually loaded data array.");
sys.exit()
lfp_ca1_all_raw_sr = lfp_data_processed_scaled[:, ca1_indices_in_loaded]
lfp_ca1_avg_raw_sr = np.mean(lfp_ca1_all_raw_sr, axis=1)
lfp_dg_noise_raw_sr = None
if physical_dg_channel_1based_list:
dg_idx_in_loaded = [i for i, ph_id in
enumerate(loaded_channel_physical_ids_in_order)
if ph_id == physical_dg_channel_1based_list[0]]
if dg_idx_in_loaded: lfp_dg_noise_raw_sr = lfp_data_processed_scaled[:, dg_idx_in_loaded[0]].squeeze()
print(f"Averaged CA1 LFP (raw SR). Shape: {lfp_ca1_avg_raw_sr.shape}")
if lfp_dg_noise_raw_sr is not None:
print(f"DG Noise LFP (raw SR). Shape: {lfp_dg_noise_raw_sr.shape}")
fs_current = original_fs_from_oebin
if fs_current <= TARGET_FS_DOWNSAMPLED:
print(f"Original SR ({fs_current} Hz) is at/below target ({TARGET_FS_DOWNSAMPLED} Hz). No decimation.")
lfp_ca1_avg_final = lfp_ca1_avg_raw_sr.copy()
lfp_dg_noise_final = lfp_dg_noise_raw_sr.copy() if lfp_dg_noise_raw_sr is not None else None
lfp_timestamps_final = lfp_timestamps_loaded.copy()
else:
decimation_factor = int(round(fs_current / TARGET_FS_DOWNSAMPLED))
if decimation_factor < 1: decimation_factor = 1
print(f"\n--- Downsampling LFP from {fs_current} Hz to ~{TARGET_FS_DOWNSAMPLED} Hz (factor: {decimation_factor}) ---")
lfp_ca1_avg_final = scipy.signal.decimate(lfp_ca1_avg_raw_sr, decimation_factor, ftype='fir', zero_phase=True)
print(f"CA1 LFP downsampled. New shape: {lfp_ca1_avg_final.shape}")
if lfp_dg_noise_raw_sr is not None:
lfp_dg_noise_final = scipy.signal.decimate(lfp_dg_noise_raw_sr, decimation_factor, ftype='fir', zero_phase=True)
print(f"DG LFP downsampled. New shape: {lfp_dg_noise_final.shape}")
else: lfp_dg_noise_final = None
lfp_timestamps_final = lfp_timestamps_loaded[::decimation_factor]
fs_current = fs_current / decimation_factor
print(f"Timestamps downsampled. New length: {len(lfp_timestamps_final)}")
print(f"New effective sampling rate (fs_current): {fs_current:.2f} Hz")
print("\n--- Performing Sleep State Scoring on downsampled data ---")
full_spec_sxx, full_spec_freqs, spec_times_centered = compute_spectrogram_custom(
lfp_ca1_avg_final, fs_current, SPECTROGRAM_WINDOW_SEC, SPECTROGRAM_STEP_SEC
)
if full_spec_sxx is None:
print("Spectrogram failed. Exiting.");
sys.exit()
# --- Calculate Sleep States ---
# 1. Existing PCA/Theta Method
freq_mask_for_pca = (full_spec_freqs >= PCA_SPECTROGRAM_FREQ_RANGE[0]) & (full_spec_freqs <= PCA_SPECTROGRAM_FREQ_RANGE[1])
if not np.any(freq_mask_for_pca):
print(f"Error: No freqs in {PCA_SPECTROGRAM_FREQ_RANGE} Hz for PCA. Using full spectrum.");
spec_sxx_for_pca, freqs_for_pca = full_spec_sxx, full_spec_freqs
else:
spec_sxx_for_pca, freqs_for_pca = full_spec_sxx[freq_mask_for_pca, :], full_spec_freqs[freq_mask_for_pca]
spec_abs_times = spec_times_centered + (lfp_timestamps_final[0] if len(lfp_timestamps_final)>0 else 0)
pc1 = compute_pca_custom(spec_sxx_for_pca)
theta_dominance = compute_theta_dominance_custom(full_spec_sxx, full_spec_freqs,
THETA_DOMINANCE_THETA_BAND, THETA_DOMINANCE_TOTAL_BAND)
min_len_metrics = min(len(pc1) if pc1 is not None else 0,
len(theta_dominance) if theta_dominance is not None else 0,
len(spec_abs_times))
# Align arrays
if pc1 is not None: pc1 = pc1[:min_len_metrics]
if theta_dominance is not None: theta_dominance = theta_dominance[:min_len_metrics]
spec_abs_times_aligned = spec_abs_times[:min_len_metrics]
sleep_states, sleep_thresholds = score_sleep_states_custom(pc1, theta_dominance, spec_abs_times_aligned, SPECTROGRAM_STEP_SEC, SCORING_BUFFER_SEC)
# 2. Logic for NREM Detection Selection
is_nrem_state = None
if USE_DT_RATIO_FOR_NREM:
print(f"\n--- [OPTION ENABLED] Estimating NREM using Delta/Theta Ratio > {DT_RATIO_THRESHOLD} ---")
# Calculate Delta Power (Mean power in Delta band)
delta_mask = (full_spec_freqs >= DT_DELTA_BAND[0]) & (full_spec_freqs <= DT_DELTA_BAND[1])
theta_mask = (full_spec_freqs >= DT_THETA_BAND[0]) & (full_spec_freqs <= DT_THETA_BAND[1])
if np.any(delta_mask) and np.any(theta_mask):
dt_delta_power = np.nanmean(full_spec_sxx[delta_mask, :], axis=0)[:min_len_metrics]
dt_theta_power = np.nanmean(full_spec_sxx[theta_mask, :], axis=0)[:min_len_metrics]
# Avoid division by zero
dt_ratio = np.zeros_like(dt_delta_power)
valid_ratio_mask = dt_theta_power > 1e-9
dt_ratio[valid_ratio_mask] = dt_delta_power[valid_ratio_mask] / dt_theta_power[valid_ratio_mask]
is_nrem_state = dt_ratio > DT_RATIO_THRESHOLD
print(f" Identified {np.sum(is_nrem_state)} spectrogram bins as NREM based on Delta/Theta ratio.")
else:
print("Error: Could not calculate Delta or Theta power for ratio. Fallback to scoring.")
if sleep_states is not None: is_nrem_state = sleep_states == 1
else:
# Standard scoring
if sleep_states is not None: is_nrem_state = sleep_states == 1
# Fallback if detection failed completely
if is_nrem_state is None:
print("Warning: NREM state detection failed. No NREM epochs will be processed.")
is_nrem_state = np.zeros(min_len_metrics, dtype=bool)
nrem_diff = np.diff(is_nrem_state.astype(int))
nrem_start_indices_spec = np.where(nrem_diff == 1)[0] + 1
nrem_end_indices_spec = np.where(nrem_diff == -1)[0] + 1
if len(is_nrem_state)>0:
if is_nrem_state[0]: nrem_start_indices_spec = np.insert(nrem_start_indices_spec, 0, 0)
if is_nrem_state[-1]: nrem_end_indices_spec = np.append(nrem_end_indices_spec, len(is_nrem_state))
final_nrem_starts, final_nrem_ends = [], []
if len(nrem_start_indices_spec) > 0 and len(nrem_end_indices_spec) > 0:
# Ensure not empty before potential indexing
# Ensure start_indices are less than end_indices and arrays are of same length for zipping
min_len_nrem_edges = min(len(nrem_start_indices_spec), len(nrem_end_indices_spec))
nrem_start_indices_spec = nrem_start_indices_spec[:min_len_nrem_edges]
nrem_end_indices_spec = nrem_end_indices_spec[:min_len_nrem_edges]
for s, e_idx in zip(nrem_start_indices_spec, nrem_end_indices_spec):
if e_idx > s: final_nrem_starts.append(s); final_nrem_ends.append(e_idx)
nrem_start_indices_spec, nrem_end_indices_spec = np.array(final_nrem_starts), np.array(final_nrem_ends)
nrem_periods_sec = []
lfp_start_time_abs_final = lfp_timestamps_final[0] if len(lfp_timestamps_final) > 0 else 0
lfp_end_time_abs_final = lfp_timestamps_final[-1] if len(lfp_timestamps_final) > 0 else ((len(lfp_ca1_avg_final) / fs_current + lfp_start_time_abs_final) if 'lfp_ca1_avg_final' in locals() and lfp_ca1_avg_final is not None else lfp_start_time_abs_final)
for s_idx, e_idx in zip(nrem_start_indices_spec, nrem_end_indices_spec):
if e_idx > s_idx and s_idx < len(spec_abs_times_aligned) and (e_idx -1) < len(spec_abs_times_aligned): # e_idx-1 for indexing spec_abs_times
start_t = spec_abs_times_aligned[s_idx] - SPECTROGRAM_STEP_SEC / 2
end_t = spec_abs_times_aligned[e_idx - 1] + SPECTROGRAM_STEP_SEC / 2
start_t,end_t = max(start_t,lfp_start_time_abs_final),min(end_t,lfp_end_time_abs_final)
if end_t > start_t: nrem_periods_sec.append((start_t, end_t))
total_nrem_duration_s = sum(e - s for s, e in nrem_periods_sec)
print(f"Total NREM duration: {total_nrem_duration_s:.2f}s from {len(nrem_periods_sec)} episodes.")
all_detected_ripples_df = pd.DataFrame()
if total_nrem_duration_s > 0 and nrem_periods_sec:
print(f"\n--- Detecting ripples in {len(nrem_periods_sec)} NREM epochs ---")
for i, (nrem_s_time, nrem_e_time) in enumerate(nrem_periods_sec):
nrem_s_sample = np.searchsorted(lfp_timestamps_final, nrem_s_time, side='left')
nrem_e_sample = np.searchsorted(lfp_timestamps_final, nrem_e_time, side='right')
if nrem_e_sample <= nrem_s_sample: continue
seg_ca1, seg_ts = lfp_ca1_avg_final[nrem_s_sample:nrem_e_sample], lfp_timestamps_final[nrem_s_sample:nrem_e_sample]
seg_dg = lfp_dg_noise_final[nrem_s_sample:nrem_e_sample] if lfp_dg_noise_final is not None else None
min_len_for_ripple = int(RIPPLE_DURATIONS[2] / 1000 * fs_current * 1.1)
if len(seg_ca1) < min_len_for_ripple: continue
# print(f" Processing NREM epoch {i+1}, duration {len(seg_ca1)/fs_current:.2f}s")
df_rip = find_swr_custom(seg_ca1, seg_ts, fs_current, RIPPLE_THRESHOLDS, RIPPLE_DURATIONS,
RIPPLE_FREQ_RANGE, RIPPLE_ENVELOPE_FILTER_HZ, seg_dg)
if not df_rip.empty:
for col_sample in ['start_sample', 'end_sample', 'peak_sample_power', 'peak_sample_lfp']:
if col_sample in df_rip and pd.api.types.is_numeric_dtype(df_rip[col_sample]):
df_rip[col_sample] = df_rip[col_sample].add(nrem_s_sample, fill_value=0).astype('Int64')
all_detected_ripples_df = pd.concat([all_detected_ripples_df, df_rip], ignore_index=True)
if not all_detected_ripples_df.empty:
all_detected_ripples_df['duration_ms'] = (all_detected_ripples_df['end_time'] - all_detected_ripples_df['start_time']) * 1000.0
all_detected_ripples_df = all_detected_ripples_df.sort_values(by='start_time').reset_index(drop=True)
if len(all_detected_ripples_df) > 0: all_detected_ripples_df['IRI_s'] = np.nan
if len(all_detected_ripples_df) > 1:
iri_values = all_detected_ripples_df['start_time'].iloc[1:].values - all_detected_ripples_df['end_time'].iloc[:-1].values
all_detected_ripples_df.loc[1:, 'IRI_s'] = iri_values
all_detected_ripples_df['num_cycles_approx'] = all_detected_ripples_df.apply(
lambda r: (r['duration_ms']/1000.0)*np.mean(RIPPLE_FREQ_RANGE) if pd.notna(r['duration_ms']) else np.nan, axis=1)
else: print("No NREM / NREM duration zero. Skipping ripple detection.")
if not all_detected_ripples_df.empty:
print(f"\n--- Detected {len(all_detected_ripples_df)} ripples (manual load, downsampled) ---")
ripple_band_filtered_lfp_ca1_avg_final = butter_bandpass_filter(lfp_ca1_avg_final, RIPPLE_FREQ_RANGE[0], RIPPLE_FREQ_RANGE[1], fs_current, order=4)
if len(ripple_band_filtered_lfp_ca1_avg_final) > 0 and not np.all(np.isnan(ripple_band_filtered_lfp_ca1_avg_final)):
z_ripple_band_filtered_lfp_ca1_avg_final = zscore(ripple_band_filtered_lfp_ca1_avg_final, nan_policy='omit')
z_ripple_band_filtered_lfp_ca1_avg_final = np.nan_to_num(z_ripple_band_filtered_lfp_ca1_avg_final)
new_peak_lfp_amplitude_zscore = []
for idx, ripple in all_detected_ripples_df.iterrows():
s_abs_val, e_abs_val = ripple['start_sample'], ripple['end_sample']
if pd.isna(s_abs_val) or pd.isna(e_abs_val): new_peak_lfp_amplitude_zscore.append(np.nan); continue
s_abs, e_abs = int(s_abs_val), int(e_abs_val)
if s_abs>=e_abs or e_abs > len(z_ripple_band_filtered_lfp_ca1_avg_final) or s_abs < 0 :
new_peak_lfp_amplitude_zscore.append(np.nan); continue
segment_z_lfp_abs = z_ripple_band_filtered_lfp_ca1_avg_final[s_abs:e_abs]
if len(segment_z_lfp_abs)>0:
abs_max_idx_in_segment=np.argmax(np.abs(segment_z_lfp_abs))
new_peak_lfp_amplitude_zscore.append(segment_z_lfp_abs[abs_max_idx_in_segment])
else: new_peak_lfp_amplitude_zscore.append(np.nan)
all_detected_ripples_df['peak_lfp_amplitude_zscore'] = new_peak_lfp_amplitude_zscore
else:
all_detected_ripples_df['peak_lfp_amplitude_zscore'] = np.nan
print("Warning: Full CA1 LFP avg (downsampled) problematic for LFP Z-score calculation.")
plot_ripple_details(lfp_ca1_avg_final, fs_current, all_detected_ripples_df,
RIPPLE_PLOT_WINDOW_MS, RIPPLE_FREQ_RANGE,
output_dir_specific, session_name_for_files)
csv_events_path = output_dir_specific / f"{session_name_for_files}_ripple_events_detailed_manual_ds.csv"
all_detected_ripples_df.to_csv(csv_events_path, index=False, float_format='%.4f')
print(f"Saved detailed ripple events to {csv_events_path}")
ripple_rate_hz = len(all_detected_ripples_df)/total_nrem_duration_s if total_nrem_duration_s > 0 else 0
summary_metrics = {
'recording_name': session_name_for_files, 'total_nrem_duration_s': total_nrem_duration_s,
'total_ripples_detected': len(all_detected_ripples_df), 'ripple_rate_nrem_hz': ripple_rate_hz,
'mean_duration_ms': all_detected_ripples_df['duration_ms'].mean(),
'mean_peak_power_zscore': all_detected_ripples_df['peak_power_zscore'].mean(),
'mean_peak_lfp_amplitude_zscore': all_detected_ripples_df['peak_lfp_amplitude_zscore'].mean(),
'mean_IRI_s': all_detected_ripples_df['IRI_s'].mean(skipna=True),
'mean_num_cycles_approx': all_detected_ripples_df['num_cycles_approx'].mean()
}
summary_df = pd.DataFrame([summary_metrics])
csv_summary_path = output_dir_specific / f"{session_name_for_files}_ripple_summary_stats_manual_ds.csv"
summary_df.to_csv(csv_summary_path, index=False, float_format='%.4f')
print(f"Saved ripple summary statistics to {csv_summary_path}")
else: print("No ripples detected. No CSV/plots for ripples.")
if 'root_tk' in locals() and root_tk is not None:
try: root_tk.destroy()
except tk.TclError: pass
print(f"\nAnalysis complete for {session_dir} (using manual data loading and downsampling)")

Ripples_Z-scored_Single Trial Script.py, under CC-BY-4.0 · at the source

Overview

Authors: Yu-Tzu Shih^1,2,3,4,5, Jason Bondoc Alipio^1,2,3,4,5, Zin-Juan Klaft⁶, Nathaniel Green^1,2,3,4, Alok Nath Mohapatra^1,2,3,4, Travis D Goode^1,2,3,4, Muthu Panchanatham^1,2,3,4, Devesh Pathak^1,2,3, Lai Ping Wong⁷, Ruslan Sadreyev⁷, Jung Ho Hyun⁸, Omar Ahmed^9,10,11, Chris Dulla⁶, Amar Sahay^1,2,3,4

ORCID iDs: Yu-Tzu Shih, Zin-Juan Klaft, Alok Nath Mohapatra, Travis D Goode, Jung Ho Hyun, Omar Ahmed, Chris Dulla, Amar Sahay

Center for Regenerative Medicine, Massachusetts General Hospital, Boston, MA, USA
Harvard Stem Cell Institute, Cambridge, MA, USA
Department of Psychiatry, Massachusetts General Hospital, Harvard Medical School, Boston, MA, USA
BROAD Institute of MIT and Harvard, Cambridge, MA, USA
These authors contributed equally: Yu-Tzu Shih, Jason Bondoc Alipio
Department of Neuroscience, Tufts University School of Medicine, Boston, MA, USA
Department of Molecular Biology, Massachusetts General Hospital, Harvard Medical School, Boston, MA, USA
Department of Brain Sciences, Daegu Gyeongbuk Institute of Science and Technology, Daegu, South Korea
Department of Psychology, University of Michigan, Ann Arbor, MI, USA
Neuroscience Graduate Program, University of Michigan, Ann Arbor, MI, USA
Department of Biomedical Engineering, University of Michigan, Ann Arbor, MI, USA

Institutions: Broad Institute (United States); Harvard University (United States); Massachusetts General Hospital (United States); Harvard Stem Cell Institute (United States); Tufts University (United States); Daegu Gyeongbuk Institute of Science and Technology (South Korea); University of Michigan (United States)

Journal: Nature, pages 10.1038/s41586-026-10907-8

Dates: published online 12 August 2026; in print August 2026

Type: Research article · Language: English

License: CC BY-NC-ND

Identifiers: DOI 10.1038/s41586-026-10907-8 · PMID 42587157 · PMCID PMC13531005 · OpenAlex W7202253421

Open access: hybrid, a free copy (OpenAlex)

Status: code verified

Categories: other condition (population), epilepsy (population)

Methods: Spectral & time-frequency, Connectivity, Statistics, Smoothing, state filtering, decompositions, Preprocessing, Evoked potentials

Topic: Neurogenesis and neuroplasticity mechanisms (Developmental Neuroscience, Neuroscience), according to OpenAlex

Funding: NIMH NIH HHS (R01 MH111729, R01 MH131652); NINDS NIH HHS (R01 NS139468, R01 NS113499); NIA NIH HHS (R01 AG076612)

Citations: not cited yet (Europe PMC); 79 references in the paper

Research resources: PV RRID:AB_10000344, GFAP RRID:AB_177521, MEIS2 RRID:AB_2143028, RGS14 RRID:AB_2179918, RGS14 RRID:AB_2179931, TBR1 RRID:AB_2200219, RFP RRID:AB_2209751, GFP RRID:AB_2534023, gephyrin RRID:AB_2619834, FOS RRID:AB_2619946, PV RRID:AB_2721207, RFP RRID:AB_2877097, SYT2 RRID:AB_2916272, BCL11A RRID:AB_444947, SST RRID:AB_831726, IBA1 RRID:AB_839504

Abstract

The abstract is not reproduced here: the paper's license (CC BY-NC-ND) does not allow it. Read it in the paper, at the publisher or on Europe PMC.

Repository

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

Zenodo 18496478

License: CC-BY-4.0

State: the link answers, verified on 27 September 2026

Evidence: files inventoried

Languages: Python (1)

Size: 4 files, 1 script

Software Heritage: not checked

Found in: “Code availability”

Holds: README

Not found: license file, CITATION.cff, environment file, tests, continuous integration, documentation

Tools: Matplotlib (1 file), NumPy (1 file), pandas (1 file), scikit-learn (1 file), SciPy (1 file)

Availability: 1 check, the latest on 27 September 2026: the link answers (HTTP 200)

27 September 2026: the link answers (HTTP 200)

1 file

Ripples_Z-scored_Single Trial Script.py, Python, 938 lines, 3 matches

At the source: doi.org/10.5281/zenodo.18496478

Code availability statement

The paper has a code availability statement. Its license (CC BY-NC-ND) does not allow reproducing it here; in short, from what the harvester recognized in it:

it points to the authors' code: Zenodo 18496478

Read it in the paper: doi.org/10.1038/s41586-026-10907-8.

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:

1 repository of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
1 script, each with its path and the digest of its content;
3 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

Datasets cited

geo:GSE283741, at NCBI GEO; found in “Data availability”

Data availability statement

The paper has a data availability statement. Its license (CC BY-NC-ND) does not allow reproducing it here; in short, from what the harvester recognized in it:

it points to a dataset: NCBI GEO GSE283741

Read it in the paper: doi.org/10.1038/s41586-026-10907-8.

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 3, 28 September 2026

Publisher: n/a → Nature Portfolio

Version 1, 27 September 2026: the first record

Recorded: type, language, journal, pages, dates, 14 authors, 3 funders, 78 references, 16 RRIDs.

Cite

This paper

Shih, Y.-T., Alipio, J. B., Klaft, Z.-J., Green, N., Mohapatra, A. N., Goode, T. D., Panchanatham, M., Pathak, D., Wong, L. P., Sadreyev, R., Hyun, J. H., Ahmed, O., Dulla, C., & Sahay, A. (2026). Procognitive restoration of PV neuron plasticity in neurodevelopmental disorders. Nature, 10.1038/s41586-026-10907-8. https://doi.org/10.1038/s41586-026-10907-8

BibTeX

@article{shih2026procognitive, author = {Shih, Yu-Tzu and Alipio, Jason Bondoc and Klaft, Zin-Juan and Green, Nathaniel and Mohapatra, Alok Nath and Goode, Travis D and Panchanatham, Muthu and Pathak, Devesh and Wong, Lai Ping and Sadreyev, Ruslan and Hyun, Jung Ho and Ahmed, Omar and Dulla, Chris and Sahay, Amar}, title = {{Procognitive restoration of PV neuron plasticity in neurodevelopmental disorders}}, journal = {Nature}, year = {2026}, month = aug, pages = {10.1038/s41586--026--10907--8}, publisher = {Nature Portfolio}, issn = {0028-0836}, doi = {10.1038/s41586-026-10907-8}, url = {https://doi.org/10.1038/s41586-026-10907-8}, pmid = {42587157}, pmcid = {PMC13531005} }

RIS

TY - JOUR AU - Shih, Yu-Tzu AU - Alipio, Jason Bondoc AU - Klaft, Zin-Juan AU - Green, Nathaniel AU - Mohapatra, Alok Nath AU - Goode, Travis D AU - Panchanatham, Muthu AU - Pathak, Devesh AU - Wong, Lai Ping AU - Sadreyev, Ruslan AU - Hyun, Jung Ho AU - Ahmed, Omar AU - Dulla, Chris AU - Sahay, Amar TI - Procognitive restoration of PV neuron plasticity in neurodevelopmental disorders T2 - Nature J2 - Nature PY - 2026 DA - 2026/08/12 SP - 10.1038/s41586 EP - 026-10907-8 SN - 0028-0836 PB - Nature Portfolio DO - 10.1038/s41586-026-10907-8 UR - https://doi.org/10.1038/s41586-026-10907-8 LA - en ER -

CSL-JSON

{ "id": "10.1038/s41586-026-10907-8", "type": "article-journal", "title": "Procognitive restoration of PV neuron plasticity in neurodevelopmental disorders", "container-title": "Nature", "author": [ { "family": "Shih", "given": "Yu-Tzu" }, { "family": "Alipio", "given": "Jason Bondoc" }, { "family": "Klaft", "given": "Zin-Juan" }, { "family": "Green", "given": "Nathaniel" }, { "family": "Mohapatra", "given": "Alok Nath" }, { "family": "Goode", "given": "Travis D" }, { "family": "Panchanatham", "given": "Muthu" }, { "family": "Pathak", "given": "Devesh" }, { "family": "Wong", "given": "Lai Ping" }, { "family": "Sadreyev", "given": "Ruslan" }, { "family": "Hyun", "given": "Jung Ho" }, { "family": "Ahmed", "given": "Omar" }, { "family": "Dulla", "given": "Chris" }, { "family": "Sahay", "given": "Amar" } ], "container-title-short": "Nature", "page": "10.1038/s41586-026-10907-8", "DOI": "10.1038/s41586-026-10907-8", "PMID": "42587157", "PMCID": "PMC13531005", "ISSN": "0028-0836", "publisher": "Nature Portfolio", "URL": "https://doi.org/10.1038/s41586-026-10907-8", "language": "en", "issued": { "date-parts": [ [ 2026, 8, 12 ] ] } }

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

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.