Joint Sharp Wave Ripple Detection¶

%reload_ext autoreload
%autoreload 2
%matplotlib inline

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from scipy import ndimage, signal

import neuro_py as npy
from neuro_py.detectors.sharp_wave_ripple import detect_sharp_wave_ripples

npy.plotting.set_plotting_defaults()

RIPPLE_FREQUENCY = 165.0


def normalize(signal_in):
    signal_in = np.asarray(signal_in, dtype=float)
    return (signal_in - signal_in.mean()) / signal_in.std()


def brown_noise(n_samples, seed=None):
    """Generate normalized brown noise with power decreasing at low frequencies."""
    rng = np.random.default_rng(seed)
    uneven = n_samples % 2
    spectrum = rng.normal(size=n_samples // 2 + 1 + uneven) + 1j * rng.normal(
        size=n_samples // 2 + 1 + uneven
    )
    scale = np.arange(1, spectrum.size + 1)
    noise = np.fft.irfft(spectrum / scale).real
    if uneven:
        noise = noise[:-1]
    return normalize(noise)


def simulate_swr_session(
    fs=1250.0,
    duration=1.5,
    swr_times=(0.45, 1.02),
    ripple_amplitude=2.0,
    ripple_duration=0.100,
    sharp_wave_amplitude=-2.8,
    sharp_wave_duration=0.140,
    ripple_noise_amplitude=1.1,
    sharp_wave_noise_amplitude=0.8,
    seed=7,
):
    timestamps = np.arange(0.0, duration, 1.0 / fs)
    ripple_channel = ripple_noise_amplitude * brown_noise(timestamps.size, seed=seed)
    sharp_wave_channel = sharp_wave_noise_amplitude * brown_noise(
        timestamps.size, seed=seed + 1
    )

    ripple_carrier = np.sin(2 * np.pi * RIPPLE_FREQUENCY * timestamps)
    for center in swr_times:
        ripple_envelope = np.exp(
            -((timestamps - center) ** 2) / (2 * (ripple_duration / 6.0) ** 2)
        )
        sharp_wave_envelope = np.exp(
            -((timestamps - center) ** 2) / (2 * (sharp_wave_duration / 6.0) ** 2)
        )
        ripple_channel += ripple_amplitude * ripple_carrier * ripple_envelope
        sharp_wave_channel += sharp_wave_amplitude * sharp_wave_envelope

    truth = pd.DataFrame(
        {
            "event_type": "swr",
            "true_time": np.asarray(swr_times, dtype=float),
            "true_start": np.asarray(swr_times, dtype=float) - ripple_duration / 2.0,
            "true_stop": np.asarray(swr_times, dtype=float) + ripple_duration / 2.0,
        }
    )
    return timestamps, ripple_channel, sharp_wave_channel, truth


sim_timestamps, sim_ripple, sim_sharp_wave, sim_truth = simulate_swr_session()

sim_events = detect_sharp_wave_ripples(
    ripple_signal=sim_ripple,
    sharp_wave_signal=sim_sharp_wave,
    fs=1250.0,
    timestamps=sim_timestamps,
    ripple_channel=0,
    low_threshold=0.6,
    high_threshold=1.8,
    sharp_wave_low_threshold=0.4,
    sharp_wave_high_threshold=1.4,
    smooth_sigma=0.002,
    min_duration=0.020,
    max_duration=0.120,
    sharp_wave_min_duration=0.020,
    sharp_wave_max_duration=0.220,
    peak_window=0.060,
    merge_gap=0.010,
    boundary_mode="union",
    save_mat=False,
)

sim_events[["start", "stop", "peaks", "peakNormedPower", "sharp_wave_peakNormedPower"]]

def nearest_detection_time(event_time, detected_peaks):
    if len(detected_peaks) == 0:
        return np.nan, np.nan
    distances = np.abs(detected_peaks - event_time)
    nearest = int(np.argmin(distances))
    return float(detected_peaks[nearest]), float(distances[nearest])


detected_peaks = sim_events["peaks"].to_numpy(dtype=float)
rows = []
for event in sim_truth.itertuples(index=False):
    nearest_peak, error = nearest_detection_time(event.true_time, detected_peaks)
    rows.append(
        {
            "event_type": event.event_type,
            "true_time": event.true_time,
            "nearest_detected_peak": nearest_peak,
            "error_ms": error * 1000.0,
            "detected": bool(error <= 0.025) if not np.isnan(error) else False,
        }
    )

sim_comparison = pd.DataFrame(rows)
display(sim_comparison)

def plot_simulated_swr_result(events, truth, timestamps, ripple_raw, sharp_wave_raw):
    fig, axes = plt.subplots(
        2,
        1,
        figsize=npy.plotting.set_size("nature_double", fraction=1.3, ratio=0.3),
        sharex=True,
        dpi=180,
    )
    ripple_ax, sharp_wave_ax = axes

    for row in truth.itertuples(index=False):
        for ax in axes:
            ax.axvspan(
                row.true_start,
                row.true_stop,
                color="#16a34a",
                alpha=0.16,
                lw=0,
                zorder=-1000,
            )
    for start, stop in events[["start", "stop"]].to_numpy():
        for ax in axes:
            ax.axvspan(start, stop, color="#f59e0b", alpha=0.30, lw=0, zorder=-900)

    ripple_ax.plot(timestamps, ripple_raw, color="#111827", lw=1.1)
    ripple_ax.set_ylabel("Ripple\nchannel")

    sharp_wave_ax.plot(timestamps, sharp_wave_raw, color="#047857", lw=1.2)
    sharp_wave_ax.set_ylabel("Sharp-wave\nchannel")
    sharp_wave_ax.set_xlabel("Time (s)")

    for ax in axes:
        ax.margins(x=0)
        ax.spines[["top", "right"]].set_visible(False)

    ripple_ax.text(
        0.01,
        1,
        "green = simulated SWR, orange = detected SWR",
        transform=ripple_ax.transAxes,
        ha="left",
        va="top",
        fontsize=10,
        bbox={
            "boxstyle": "round,pad=0.3",
            "facecolor": "white",
            "alpha": 0.5,
            "edgecolor": "none",
        },
    )
    return fig


fig = plot_simulated_swr_result(
    sim_events,
    sim_truth,
    sim_timestamps,
    sim_ripple,
    sim_sharp_wave,
)
fig.suptitle("Simulated joint sharp wave ripple detection", y=1.04, fontsize=15);

basepath = r"S:\data\HMC\HMC1\day8" # 12691.7

ripples = detect_sharp_wave_ripples(
    basepath=basepath,
    low_threshold=1.0,
    high_threshold=2.5,
    sharp_wave_low_threshold=0.25,
    sharp_wave_high_threshold=1.5,
    min_duration=0.020,
    max_duration=0.250,
    sharp_wave_min_duration=0.005,
    sharp_wave_max_duration=0.500,
    merge_gap=0.001,
    noise_threshold=10,
    min_inter_event_interval=0.025,
    peak_window=0.150,
    boundary_mode="union",
    save_mat=False, # CHANGE TO TRUE TO SAVE MATLAB FILE
)

ripples

summary_columns = [
    "start",
    "stop",
    "peaks",
    "duration",
    "ripple_duration",
    "sharp_wave_duration",
    "frequency",
    "peakNormedPower",
    "sharp_wave_peakNormedPower",
    "noise_peakNormedPower",
]

print(f"Detected {len(ripples):,} joint SWR events")
display(ripples[summary_columns].head(10))

Detected 11,376 joint SWR events

n_max_time = 20 * 60
print(
    f"Number of detected events in first {n_max_time} seconds: {(ripples.stop < n_max_time).sum()}"
)
channel_tags = npy.io.load_channel_tags(basepath)


def first_tagged_channel(tags, *keys):
    for key in keys:
        if key in tags:
            return int(np.atleast_1d(tags[key]["channels"])[0] - 1)
    return None


def zscore(values):
    values = np.asarray(values, dtype=float)
    scale = np.nanstd(values)
    if scale == 0 or np.isnan(scale):
        return np.zeros_like(values, dtype=float)
    return (values - np.nanmean(values)) / scale


ripple_channel = (
    int(ripples["ripple_channel"].iloc[0])
    if (
        not ripples.empty
        and "ripple_channel" in ripples.columns
        and not np.isnan(ripples["ripple_channel"].iloc[0])
    )
    else first_tagged_channel(channel_tags, "ripple", "Ripple", "CA1sp", "ca1sp")
)
sharp_wave_channel = first_tagged_channel(
    channel_tags, "SharpWave", "sharpwave", "sharp_wave", "CA1sr", "ca1sr"
)

if ripple_channel is None or sharp_wave_channel is None:
    raise ValueError(
        "This tutorial plot needs both ripple and sharp-wave channel tags."
    )

events_for_plot = (
    ripples.iloc[ripples.stop.values < n_max_time, :]
    .sort_values(["peakNormedPower", "sharp_wave_peakNormedPower"], ascending=False)
    .head(6)
)
plot_margin = 0.25
plot_start = max(0.0, float(events_for_plot["start"].min() - plot_margin))
plot_stop = float(events_for_plot["stop"].max() + plot_margin)
plot_epoch = np.array([[plot_start, plot_stop]])

lfp = npy.io.LFPLoader(
    basepath,
    channels=[ripple_channel, sharp_wave_channel],
    ext="lfp",
    epoch=plot_epoch,
)
fs = float(lfp.fs)
timestamps = np.asarray(lfp.abscissa_vals)
lfp_data = np.asarray(lfp.data)
raw = np.asarray(lfp_data[0]).squeeze()
sharp_wave_raw = np.asarray(lfp_data[1]).squeeze()

ripple_sos = signal.butter(4, (120.0, 250.0), btype="bandpass", fs=fs, output="sos")
ripple_filtered = signal.sosfiltfilt(ripple_sos, raw)
ripple_envelope = np.abs(signal.hilbert(ripple_filtered))
ripple_envelope = ndimage.gaussian_filter1d(
    ripple_envelope, sigma=0.004 * fs, mode="nearest"
)
ripple_power = zscore(ripple_envelope)

sharp_sos = signal.butter(4, (2.0, 50.0), btype="bandpass", fs=fs, output="sos")
ripple_low = signal.sosfiltfilt(sharp_sos, raw)
sharp_low = signal.sosfiltfilt(sharp_sos, sharp_wave_raw)
sharp_wave_diff = ripple_low - sharp_low
sharp_wave_power = zscore(sharp_wave_diff)

print(f"Ripple channel: {ripple_channel}")
print(f"Sharp-wave channel: {sharp_wave_channel}")
print(f"Loaded plotting window: {plot_start:.2f} to {plot_stop:.2f} s")
print(f"Loaded LFP shape: {lfp_data.shape}")
print(f"Sampling rate: {fs:.1f} Hz")

Number of detected events in first 1200 seconds: 459
Ripple channel: 261
Sharp-wave channel: 282
Loaded plotting window: 227.08 to 1029.31 s
Loaded LFP shape: (2, 1002791)
Sampling rate: 1250.0 Hz

def plot_joint_swr_summary(
    events,
    events_for_plot,
    timestamps,
    ripple_raw,
    sharp_wave_raw,
    ripple_power,
    sharp_wave_power,
    fs,
    ripple_high_threshold=2.5,
    sharp_wave_high_threshold=1.5,
    example_window=0.150,
):
    fig = plt.figure(
        figsize=npy.plotting.set_size("nature_double", fraction=1.65),
        constrained_layout=True,
        dpi=200,
    )
    outer = fig.add_gridspec(2, 1, height_ratios=[0.9, 2.2], hspace=0.08)

    ax_overview = fig.add_subplot(outer[0, 0])
    ax_overview.plot(
        timestamps, ripple_raw, color="#111827", lw=0.8, label="Ripple channel raw"
    )
    ax_overview.plot(
        timestamps,
        sharp_wave_raw,
        color="#047857",
        lw=0.8,
        alpha=0.85,
        label="Sharp-wave channel raw",
    )
    for start, stop in events.loc[
        (events["stop"] >= timestamps[0]) & (events["start"] <= timestamps[-1]),
        ["start", "stop"],
    ].to_numpy():
        ax_overview.axvspan(start, stop, color="#f59e0b", alpha=0.16, lw=0)
    ax_overview.set_title("Raw LFP in the loaded plotting window")
    ax_overview.set_xlabel("Time (s)")
    ax_overview.set_ylabel("LFP")
    ax_overview.legend(loc="upper right", frameon=False)

    event_grid = outer[1, 0].subgridspec(2, 3, hspace=0.34, wspace=0.22)
    event_axes = event_grid.subplots().ravel()
    half_window = int(round(example_window * fs))

    for panel_index, (container_ax, (_, event)) in enumerate(
        zip(event_axes, events_for_plot.sort_values("peaks").iterrows())
    ):
        container_ax.axis("off")
        subgrid = container_ax.get_subplotspec().subgridspec(
            2, 1, height_ratios=[1.2, 0.8], hspace=0.05
        )
        raw_ax = fig.add_subplot(subgrid[0, 0])
        feature_ax = fig.add_subplot(subgrid[1, 0], sharex=raw_ax)

        peak_idx = int(np.argmin(np.abs(timestamps - event["peaks"])))
        start_idx = max(0, peak_idx - half_window)
        stop_idx = min(len(timestamps), peak_idx + half_window + 1)
        local_time = (timestamps[start_idx:stop_idx] - event["peaks"]) * 1000.0

        raw_ax.plot(
            local_time,
            ripple_raw[start_idx:stop_idx],
            color="#111827",
            lw=1.0,
            label="Ripple raw",
        )
        raw_ax.plot(
            local_time,
            sharp_wave_raw[start_idx:stop_idx],
            color="#047857",
            lw=1.0,
            alpha=0.85,
            label="Sharp-wave raw",
        )
        raw_ax.axvspan(
            (event["start"] - event["peaks"]) * 1000.0,
            (event["stop"] - event["peaks"]) * 1000.0,
            color="#f59e0b",
            alpha=0.16,
            lw=0,
        )
        raw_ax.axvline(0.0, color="#2563eb", ls="--", lw=0.9)
        raw_ax.set_title(
            f"Ripple z={event['peakNormedPower']:.1f} | SW z={event['sharp_wave_peakNormedPower']:.1f}"
        )
        raw_ax.set_ylabel("Raw LFP")
        raw_ax.tick_params(labelbottom=False)

        feature_ax.plot(
            local_time,
            ripple_power[start_idx:stop_idx],
            color="#991b1b",
            lw=1.1,
            label="Ripple feature",
        )
        feature_ax.plot(
            local_time,
            sharp_wave_power[start_idx:stop_idx],
            color="#047857",
            lw=1.1,
            label="Sharp-wave feature",
        )
        feature_ax.axhline(
            ripple_high_threshold, color="#991b1b", ls="--", lw=0.8, alpha=0.75
        )
        feature_ax.axhline(
            sharp_wave_high_threshold, color="#047857", ls=":", lw=0.9, alpha=0.85
        )
        feature_ax.axvline(0.0, color="#2563eb", ls="--", lw=0.8, alpha=0.8)
        feature_ax.set_xlabel("Time from ripple trough (ms)")
        feature_ax.set_ylabel("z")

        if panel_index == 0:
            raw_ax.legend(loc="upper left", frameon=False, fontsize=8)
            feature_ax.legend(loc="upper left", frameon=False, fontsize=8)

    for ax in event_axes[len(events_for_plot) :]:
        ax.axis("off")

    return fig


fig = plot_joint_swr_summary(
    ripples,
    events_for_plot,
    timestamps,
    raw,
    sharp_wave_raw,
    ripple_power,
    sharp_wave_power,
    fs,
    ripple_high_threshold=2.5,
    sharp_wave_high_threshold=1.5,
)
fig.suptitle("Joint sharp wave ripple detection summary", y=1.02, fontsize=15);

detect_sharp_wave_ripples(
    basepath=basepath,
    low_threshold=1.0,
    high_threshold=2.5,
    sharp_wave_low_threshold=0.25,
    sharp_wave_high_threshold=1.5,
    min_duration=0.020,
    max_duration=0.250,
    sharp_wave_min_duration=0.005,
    sharp_wave_max_duration=0.500,
    noise_threshold=10,
    merge_gap=0.001,
    min_inter_event_interval=0.025,
    peak_window=0.150,
    boundary_mode="union",
    save_mat=True,
    overwrite=True,
)