harding

GMM-based OFF-period detection for multi-channel extracellular recordings, after Harding et al. (2023), BMC Neuroscience.

Harding et al. extract low-amplitude segments (LAS) from a single MUA amplitude channel: fit a Gaussian mixture to a 2-D feature space of multi-scale smoothed amplitude, classify each sample, then refine against negative half-waves of a wake-baseline-subtracted signal. This package generalizes that to probes — fit per channel, combine the per-channel detections into a spatial mask, and label connected components as discrete OFF periods.

It depends only on packages from PyPI. There is no project layout, no subject registry, and no hypnogram library: you pass a recording, a table of time intervals, and paths of your choosing.

Install

uv add harding

Use

Fitting is expensive and detection is not, so they are separate phases.

import numpy as np
import spikeinterface as si

from harding import baseline, detect, gmm, preprocess

# 0. MUA amplitude envelope. The chain is lazy; save it once and fit against
#    the saved copy rather than re-evaluating it for every phase.
raw = si.read_spikeglx("/path/to/recording")
rec = preprocess.build_preprocessing_chain(raw, resample_rate=500)
rec = rec.save(folder="mua_traces.zarr", format="zarr", n_jobs=16)

channel_ids = np.asarray(rec.get_channel_ids())
y_coords = rec.get_channel_locations()[:, 1]        # depth, in microns
fs = rec.get_sampling_frequency()

# Time intervals in the recording's own time base: any (n, 2) array, or a
# DataFrame with start_time / end_time columns.
wake_bouts = [[120.0, 900.0], [3600.0, 4200.0]]
nrem_bouts = [[1000.0, 2400.0], [5000.0, 7200.0]]

# 1. Fit, once per recording. Both steps are parallel over channels.
bl = baseline.compute_wake_baseline(rec, wake_bouts, channel_ids, n_jobs=8)
baseline.save_wake_baseline(bl, "fits/wake_baseline.nc")

fits = gmm.fit_gmms(rec, nrem_bouts, channel_ids, fs, n_jobs=8)
gmm.save_gmm_fits(fits, "fits/gmm")

# 2. Detect, per window of interest.
offs, label_ixs = detect.detect_offs_spatial(
    rec,
    condition_bouts=nrem_bouts,
    channel_ids=channel_ids,
    y_coords=y_coords,
    gmm_fits=fits,
    wake_baseline=bl,
    n_jobs=8,
)

offs is a DataFrame, one row per OFF period, whose columns are documented in harding.morphology.Off: timing, depth span, area, convexity, trace-value summaries, centre of mass, and onset/offset edge-synchrony measures. label_ixs maps each row's label to the (time_indices, channel_indices) it occupies, for plotting overlays.

What each module does

Module	Purpose
preprocess	Lazy SpikeInterface chain → rectified, resampled MUA envelope
bouts	The time-interval table every entry point accepts
baseline	Per-channel wake baseline (mean + median)
gmm	Per-channel GMM fitting, K selection, persistence
detect	Per-channel classification, half-wave refinement, spatial combination
morphology	Mask cleaning, connected components, per-event properties
sampling	Reproducible row subsampling

Parameters

Defaults follow the paper except where noted.

Parameter	Default	Source
Bandpass	300–5000 Hz	Paper
Resample rate	500 Hz	Paper (~498 Hz)
Heavy smoothing	σ = 12.4 ms	Paper (62 ms window, wf = 2.5)
Light smoothing	σ = 4.4 ms	Paper (22 ms window, wf = 2.5)
GMM components	K = 1..8, Calinski–Harabasz	Paper
NREM sample fraction	5%	Conservative; paper uses 1%
Wake baseline	Median	More robust than the mean

Scope

Estimating probe drift is out of scope. build_preprocessing_chain takes an optional motion_correction callable applied before rectification, so you can supply your own — for example from spikeinterface.sortingcomponents.motion.

morphology.py reimplements generic spatial morphology (mask cleaning, connected components, event properties) rather than importing it from a larger analysis package, so that this package depends on nothing outside PyPI.

Harding's original single-channel MATLAB implementation is OFFAD.