Spike Analysis Tutorial#

This tutorial covers working with spike-sorted neural data and integrating it with EEG features in NeuRodent.

Overview#

Spike Analysis Results (SAR) integrates spike-sorted single-unit data with continuous EEG analysis:

  1. Load spike-sorted data from MountainSort or other spike sorters

  2. Compute peri-event time histograms (PETH)

  3. Analyze spike-LFP relationships

  4. Integrate with WAR for combined analysis

import sys
from pathlib import Path
import logging

import numpy as np
import matplotlib.pyplot as plt

from neurodent import core, visualization

logging.basicConfig(level=logging.INFO)
logger = logging.getLogger()

1. Loading Spike-Sorted Data#

NeuRodent works with spike-sorted data from various spike sorters via SpikeInterface:

# Load data
data_path = Path("/path/to/spike/data")
animal_id = "animal_001"

lro = core.LongRecordingOrganizer(
    base_folder=data_path,
    animal_id=animal_id,
    mode="bin"
)

ao = visualization.AnimalOrganizer(lro)

2. MountainSort Analysis#

If using MountainSort for spike sorting:

# Create MountainSort analyzer
msa = core.MountainSortAnalyzer(
    lro=lro,
    output_folder=Path("./mountainsort_output")
)

# Run spike sorting (if not already done)
# msa.run_mountainsort()

# Load sorting results
sorting = msa.load_sorting()
print(f"Found {len(sorting.get_unit_ids())} units")

3. Computing Spike Analysis#

Compute spike-related features:

# Compute spike analysis
sar = ao.compute_spike_analysis(
    multiprocess_mode='serial'
)

print(f"Spike analysis computed for {len(sar)} sessions")

4. Converting to MNE Format#

SAR can be converted to MNE format for advanced analysis:

# Convert each SAR to MNE format
for sar_session in sar:
    # Convert to MNE Raw object
    sar_session.convert_to_mne(
        chunk_len=1440  # Length in minutes
    )
    
    # Save as FIF format
    output_path = Path("./mne_output") / sar_session.animal_day
    sar_session.save_fif_and_json(output_path)
    
print("Converted to MNE format")

5. Loading Saved SAR#

Load previously saved spike analysis:

# Load SAR from MNE format
sar_path = Path("./mne_output/animal_001_2023-12-15")
sar_loaded = visualization.SpikeAnalysisResult.load_fif_and_json(sar_path)

# Access MNE object
mne_obj = sar_loaded.result_mne
print(f"Loaded MNE object: {mne_obj}")
print(f"Channels: {mne_obj.ch_names}")
print(f"Sample rate: {mne_obj.info['sfreq']} Hz")

6. Peri-Event Analysis with MNE#

Analyze EEG around spike times:

import mne

# Extract events from annotations (spike times)
events, event_dict = mne.events_from_annotations(mne_obj)

print(f"Found {len(events)} events")
print(f"Event types: {list(event_dict.keys())}")

# Create epochs around spike times
epochs = mne.Epochs(
    mne_obj,
    events,
    event_id=event_dict,
    tmin=-0.5,  # 500 ms before spike
    tmax=1.0,   # 1 s after spike
    baseline=(-0.5, -0.1),
    preload=True
)

print(f"Created {len(epochs)} epochs")

7. Spike-Triggered Averages#

Compute spike-triggered average of LFP:

# Compute evoked response (spike-triggered average)
evoked = epochs.average()

# Plot spike-triggered average
fig = evoked.plot(spatial_colors=True, gfp=True)
plt.title("Spike-Triggered Average")
plt.show()

8. Time-Frequency Analysis#

Analyze oscillations around spike times:

# Compute time-frequency representation
freqs = np.arange(1, 50, 1)  # 1-50 Hz
n_cycles = freqs / 2  # Number of cycles increases with frequency

# Compute TFR using Morlet wavelets
power = epochs.compute_tfr(
    method='morlet',
    freqs=freqs,
    n_cycles=n_cycles,
    use_fft=True,
    average=True
)

# Plot TFR
power.plot(
    picks='eeg',
    baseline=(-0.5, -0.1),
    mode='percent',
    title="Peri-Spike Time-Frequency"
)
plt.show()

9. Spike-LFP Phase Locking#

Analyze phase relationships between spikes and LFP oscillations:

# Filter for specific frequency band (e.g., theta: 4-8 Hz)
epochs_filtered = epochs.copy().filter(
    l_freq=4,
    h_freq=8,
    method='iir'
)

# Apply Hilbert transform to get instantaneous phase
epochs_hilbert = epochs_filtered.apply_hilbert()

# Extract phase at spike time (t=0)
spike_phases = np.angle(epochs_hilbert.get_data()[:, :, epochs.time_as_index(0)[0]])

# Plot phase distribution
plt.figure(figsize=(10, 5))
plt.subplot(1, 2, 1)
plt.hist(spike_phases.flatten(), bins=50, density=True)
plt.xlabel("Phase (radians)")
plt.ylabel("Density")
plt.title("Spike Phase Distribution")

# Polar plot
plt.subplot(1, 2, 2, projection='polar')
plt.hist(spike_phases.flatten(), bins=50, density=True)
plt.title("Spike Phase (Polar)")
plt.tight_layout()
plt.show()

10. Integrating with WAR#

Combine spike analysis with windowed EEG analysis:

# Compute both WAR and SAR
war = ao.compute_windowed_analysis(
    features=['rms', 'psdband', 'nspike'],  # Include spike count
    multiprocess_mode='serial'
)

# Access spike counts from WAR
spike_counts = war.nspike
print(f"Spike count data shape: {spike_counts.shape}")

# Plot spike counts over time
ap = visualization.AnimalPlotter(war)
fig = ap.plot_feature_over_time('nspike')
plt.title("Spike Counts Over Time")
plt.show()

Summary#

This tutorial covered:

  1. Loading spike-sorted data

  2. Computing spike analysis (SAR)

  3. Converting to MNE format

  4. Peri-event analysis

  5. Spike-triggered averages

  6. Time-frequency analysis around spikes

  7. Spike-LFP phase locking

  8. Integrating spikes with windowed analysis

Next Steps#