Skip to content

Spike Train Analysis -- 124 Functions (23 Modules)

SC-NeuroCore includes a complete spike train analysis toolkit covering every metric from Elephant, PySpike, SpikeInterface, and NeuroTools. 127 functions across 23 modules (22 spike_stats + 1 explainability). Pure NumPy, zero external dependencies.

Rust acceleration: All 22 spike_stats modules have native Rust implementations with 96 PyO3 bindings (597 tests). See Rust Analysis Engine for the Rust API reference.

Quick Start

Python
from sc_neurocore.analysis import (
    firing_rate, cv_isi, fano_factor, cross_correlation,
    power_spectrum, burst_detection, victor_purpura_distance,
    phase_locking_value, spade_detect, gpfa,
)
import numpy as np

# Generate a Poisson spike train
from sc_neurocore.analysis import homogeneous_poisson
train = homogeneous_poisson(rate_hz=50.0, duration_s=10.0)

# Basic statistics
print(f"Rate: {firing_rate(train):.1f} Hz")
print(f"CV(ISI): {cv_isi(train):.2f}")
print(f"Fano factor: {fano_factor(train):.2f}")

Module Reference

Basic (spike_stats.basic)

Spike extraction, ISI computation, rate, count, binning.

Function Description
spike_times(train, dt) Extract spike times (seconds) from a binary 0/1 array
isi(train, dt) Inter-spike intervals (seconds)
firing_rate(train, dt) Mean firing rate (Hz)
spike_count(train) Total spike count
bin_spike_train(train, bin_size) Bin a binary spike train into spike counts per bin

Variability (spike_stats.variability)

Spike train regularity, complexity, and fractal measures.

Function Description
cv_isi(train, dt) Coefficient of variation of ISI (CV=1 for Poisson, <1 regular)
cv2(train, dt) Local CV2, less sensitive to rate changes (Holt et al. 1996)
local_variation(train, dt) Local variation LV (Shinomoto et al. 2003)
lvr(train, dt, refractoriness_ms) Revised local variation LvR corrected for refractoriness (Shinomoto et al. 2009)
fano_factor(train, window_ms, dt) Fano factor: variance/mean of windowed spike counts
isi_entropy(train, dt, bins) Shannon entropy of ISI distribution (bits)
lempel_ziv_complexity(train) Lempel-Ziv 1976 complexity, normalized by N/log2(N)
approximate_entropy(train, m, r_factor) Approximate entropy ApEn (Pincus 1991)
sample_entropy(train, m, r_factor) Sample entropy SampEn, lower bias (Richman & Moorman 2000)
permutation_entropy(train, order, delay) Bandt-Pompe permutation entropy, normalized to [0,1]
hurst_exponent(train, min_window) Hurst exponent via detrended fluctuation analysis (Peng et al. 1994)
allan_factor(train, dt, n_scales) Allan factor vs window size for fractal clustering detection
rescaled_range(train, min_window) Hurst exponent via R/S analysis (Hurst 1951)
complexity_pdf(train, dt, bins) ISI probability density function (Abeles 1982)
optimal_bin_width(train, dt) Shimazaki-Shinomoto 2007 optimal histogram bin width
optimal_kernel_bandwidth(train, dt) Silverman's rule-of-thumb bandwidth for ISI KDE

Rate Estimation (spike_stats.rate)

Instantaneous and population firing rate estimation.

Function Description
instantaneous_rate(train, dt, kernel, sigma_ms) Kernel-smoothed instantaneous rate (gaussian/exponential/rectangular)
population_rate(trains, dt, sigma_ms) Population-level instantaneous rate across multiple neurons
psth(trials, bin_ms, dt) Peri-stimulus time histogram across trials

Distance Metrics (spike_stats.distance)

Spike train distance and similarity measures.

Function Description
van_rossum_distance(a, b, dt, tau_ms) Van Rossum 2001 exponential-kernel distance
victor_purpura_distance(times_a, times_b, cost_per_s) Victor-Purpura 1996 edit distance
isi_distance(a, b, dt) ISI-distance ratio comparison (Kreuz et al. 2007)
spike_distance(times_a, times_b, t_start, t_end) SPIKE-distance (Kreuz et al. 2013)
spike_sync(times_a, times_b, t_start, t_end) SPIKE-synchronization, normalized to [0,1] (Kreuz et al. 2015)
spike_sync_profile(times_a, times_b, n_bins, ...) Binned SPIKE-synchronization profile
spike_profile(times_a, times_b, n_bins, ...) Binned SPIKE-distance profile
isi_profile(a, b, dt, n_bins) Binned ISI-distance profile
adaptive_spike_distance(times_a, times_b, ...) Adaptive SPIKE-distance with cost interpolation (Kreuz et al. 2013)
schreiber_similarity(a, b, dt, sigma_ms) Smoothed Pearson correlation similarity (Schreiber et al. 2003)
hunter_milton_similarity(times_a, times_b, dt_max) Nearest-neighbour coincidence fraction (Hunter-Milton 2003)
earth_movers_distance(times_a, times_b, ...) EMD between spike time distributions (Rubner et al. 1998)
multi_neuron_victor_purpura(times_list, cost) All-pairs Victor-Purpura distance matrix
generalized_victor_purpura(times_a, times_b, cost_func) Victor-Purpura with arbitrary cost function
spike_distance_matrix(times_list, metric, ...) All-pairs distance matrix (spike_distance/spike_sync/victor_purpura)

Correlation (spike_stats.correlation)

Cross-correlation, synchrony, and covariance measures.

Function Description
cross_correlation(a, b, max_lag_ms, dt) Cross-correlogram between two binary trains
pairwise_correlation(trains, dt) Pairwise Pearson correlation matrix across neurons
event_synchronization(a, b, dt, tau_ms) Quian Quiroga et al. 2002 event synchronization
spike_train_coherence(a, b, dt) Magnitude-squared coherence between spike trains
spike_time_tiling_coefficient(a, b, dt, delta_ms) STTC, corrects for rate bias (Cutts & Eglen 2014)
covariance_matrix(trains, bin_size) Spike count covariance matrix (de la Rocha et al. 2007)
autocorrelation_time(train, dt, max_lag_ms) Autocorrelation time until first zero crossing
noise_correlation(trains, bin_size) Trial-to-trial variability correlation (Averbeck & Lee 2006)
signal_correlation(trains, bin_size) Tuning similarity via Pearson correlation of mean responses
spike_count_covariance(trains, window) Windowed spike count covariance (Kohn & Smith 2005)
joint_psth(a, b, bin_size) Joint PSTH matrix (Aertsen et al. 1989)
coincidence_index(a, b, dt, delta_ms) Rate-corrected coincidence index kappa (Joris et al. 2006)

Spectral (spike_stats.spectral)

Spectral analysis of spike trains.

Function Description
power_spectrum(train, dt) Power spectral density of a binary spike train

Temporal Patterns (spike_stats.temporal)

Burst detection, latency, onset, and change point analysis.

Function Description
burst_detection(train, dt, max_isi_ms, min_spikes) Detect bursts as consecutive short-ISI spikes
first_spike_latency(train, dt) Time to first spike (seconds)
response_onset(train, baseline_steps, dt, threshold_sigma) Detect response onset exceeding baseline + threshold
change_point_detection(train, bin_size, threshold) CUSUM-based firing rate change points (Page 1954)

Stimulus (spike_stats.stimulus)

Spike-triggered analysis and receptive field estimation.

Function Description
spike_triggered_average(stimulus, train, window_steps) STA: average stimulus preceding each spike
spike_triggered_covariance(stimulus, train, window_steps) STC: covariance of pre-spike stimulus (Schwartz et al. 2006)
spatial_information(train, positions, n_bins, dt) Spatial information in bits/spike (Skaggs et al. 1993)
place_field_detection(train, positions, ...) Detect place fields (O'Keefe & Dostrovsky 1971)
tuning_curve(train, stimulus_values, n_bins, dt) Firing rate vs stimulus value (Dayan & Abbott 2001)

LFP Coupling (spike_stats.lfp)

Spike-LFP phase locking, coherence, and phase histograms.

Function Description
phase_locking_value(train, lfp) PLV via Hilbert phase at spike times
spike_field_coherence(train, lfp, dt) Spike-field coherence SFC in frequency domain
spike_phase_histogram(train, lfp, n_bins) Histogram of LFP phase at spike times

Surrogates (spike_stats.surrogates)

Surrogate spike train generation for significance testing.

Function Description
surrogate_isi_shuffle(train, seed) ISI-shuffled surrogate preserving rate + ISI distribution
surrogate_dither(train, dither_ms, dt, seed) Spike jittering surrogate (+/- dither window)
surrogate_trial_shuffle(trains, seed) Trial-order shuffling destroying trial-to-trial correlation
homogeneous_poisson(rate_hz, duration_s, dt, seed) Generate homogeneous Poisson spike train
inhomogeneous_poisson(rate_func, duration_s, dt, seed) Time-varying Poisson via thinning (Lewis & Shedler 1979)
gamma_process(rate_hz, shape, duration_s, dt, seed) Gamma-renewal process (shape=1 Poisson, >1 regular)
compound_poisson_process(rate_hz, burst_mean, ...) Compound Poisson with burst events (Snyder & Miller 1991)
surrogate_joint_isi(train, seed) Joint-ISI surrogate preserving serial correlations (Louis et al. 2010)
surrogate_bin_shuffling(train, bin_size, seed) Within-bin spike shuffling (Hatsopoulos et al. 2003)
surrogate_spike_train_shifting(train, max_shift, seed) Circular shifting surrogate

Information Theory (spike_stats.information)

Information-theoretic measures for spike trains.

Function Description
mutual_information(a, b, bin_size) Mutual information between binned spike trains (bits)
transfer_entropy(source, target, bin_size, lag) Transfer entropy from source to target (bits)
spike_train_entropy(train, bin_size, word_length) Binary word entropy (Strong et al. 1998)
noise_entropy(train, n_trials, bin_size, word_length) Noise entropy via pseudo-trials (de Ruyter van Steveninck et al. 1997)
stimulus_specific_information(counts, stim_ids) SSI in bits (Butts 2003)
kozachenko_leonenko_mi(x, y, k) k-NN mutual information estimator in nats (Kraskov et al. 2004)
time_rescaling_ks_test(times, rate_func, ...) KS goodness-of-fit for point processes (Brown et al. 2002)

Causality (spike_stats.causality)

Granger causality and directed connectivity.

Function Description
pairwise_granger_causality(source, target, ...) Pairwise Granger causality log-likelihood ratio (Granger 1969)
conditional_granger_causality(source, target, cond, ...) Conditional GC controlling for a third signal (Geweke 1984)
spectral_granger_causality(trains, bin_size, order, ...) Frequency-domain GC via VAR model (Geweke 1982)
partial_directed_coherence(trains, ...) PDC: normalized directed connectivity (Baccala & Sameshima 2001)
directed_transfer_function(trains, ...) DTF: normalized transfer function (Kaminski & Blinowska 1991)

Dimensionality Reduction (spike_stats.dimensionality)

Low-dimensional projections of population activity.

Function Description
spike_train_pca(trains, n_components, bin_size) PCA on binned spike count matrix
demixed_pca(trains_by_condition, n_components, ...) Demixed PCA separating condition variance (Kobak et al. 2016)
factor_analysis(trains, n_factors, bin_size, n_iter) Factor analysis via EM (Rubin & Thayer 1982)

Decoding (spike_stats.decoding)

Neural population decoding algorithms.

Function Description
population_vector_decode(trains, directions, window) Georgopoulos population vector decoding
bayesian_decode(counts, tuning_rates, prior) Bayesian MAP decoder (Dayan & Abbott 2001)
maximum_likelihood_decode(counts, tuning_rates) Maximum likelihood Poisson decoder
linear_discriminant_decode(data, labels, test) Fisher LDA decoder (Fisher 1936)
naive_bayes_decode(data, labels, test) Gaussian naive Bayes decoder

Network (spike_stats.network)

Network-level analysis: connectivity, assemblies, synfire chains.

Function Description
functional_connectivity(trains, max_lag_ms, dt) Infer connectivity from peak cross-correlation
unitary_events(trains, bin_size, alpha) Detect significant synchronous bins (Gruen et al. 2002)
cell_assembly_detection(trains, bin_size, threshold) PCA-based cell assembly detection (Lopes-dos-Santos et al. 2013)
synfire_chain_detection(trains, dt, max_delay_ms, ...) Cross-correlation peak ordering (Abeles 1991)

Point Process (spike_stats.point_process)

Point process models and hazard functions.

Function Description
conditional_intensity(train, dt, window_ms) Moving-window MLE conditional intensity (Brown et al. 2004)
isi_hazard_function(train, dt, bins) ISI hazard function h(t) = f(t)/S(t) (Tuckwell 1988)
isi_survivor_function(train, dt, bins) ISI survivor function S(t) = P(ISI > t)
renewal_density(train, dt, bins) Renewal density normalized by mean rate (Cox 1962)

Sorting Quality (spike_stats.sorting_quality)

Spike sorting quality metrics.

Function Description
isolation_distance(cluster, noise) Mahalanobis isolation distance (Harris et al. 2001)
l_ratio(cluster, noise) L-ratio cluster quality (Schmitzer-Torbert et al. 2005)
silhouette_score(features, labels) Mean silhouette score (Rousseeuw 1987)
d_prime(cluster_a, cluster_b) Sensitivity index between two clusters (Green & Swets 1966)
isi_violation_rate(train, dt, refractory_ms) Fraction of ISIs below refractory period (Hill et al. 2011)
presence_ratio(train, n_bins) Fraction of time bins with at least one spike (IBL 2019)
amplitude_cutoff(amplitudes, bins) Estimated missing spike fraction (Hill et al. 2011)
snr(waveforms) Signal-to-noise ratio of spike waveforms (Suner et al. 2005)
nn_hit_rate(cluster, noise, k) k-NN cluster purity (Chung et al. 2017)
drift_metric(waveforms, timestamps, n_bins) Waveform amplitude drift over time (IBL 2019)

Waveform (spike_stats.waveform)

Spike waveform shape analysis.

Function Description
waveform_width(waveform, dt) Trough-to-peak width in seconds (Bartho et al. 2004)
waveform_amplitude(waveform) Peak-to-trough amplitude
waveform_repolarization_slope(waveform, dt) Max dV/dt after trough (Bean 2007)
waveform_recovery_slope(waveform, dt) dV/dt during return to baseline (Bean 2007)
waveform_halfwidth(waveform, dt) Duration at half-minimum amplitude (Bartho et al. 2004)
waveform_pt_ratio(waveform) Peak-to-trough ratio (Bartho et al. 2004)

Statistics (spike_stats.statistics)

Significance testing.

Function Description
significance_bootstrap(func, a, b, n_surrogates, seed) Bootstrap permutation test for pairwise statistics

Patterns (spike_stats.patterns)

Spike directionality, ordering, and higher-order patterns.

Function Description
spike_directionality(times_a, times_b, ...) Asymmetry in [-1,1]: positive means A leads B (Kreuz et al. 2015)
spike_train_order(times_list, ...) Pairwise directionality matrix (Kreuz et al. 2017)
cubic_higher_order(train, dt, max_lag) Third-order cumulant C3(tau1, tau2) (Nikias & Petropulu 1993)

SPADE (spike_stats.spade)

Spike Pattern Detection and Evaluation (Torre et al. 2013).

Function Description
spade_detect(trains, bin_ms, dt, min_support, ...) Detect significant spatiotemporal spike patterns via frequent itemset mining + surrogate testing

GPFA (spike_stats.gpfa)

Gaussian Process Factor Analysis (Yu et al. 2009).

Function Description
gpfa(trains, n_latents, bin_ms, dt, max_iter, ...) Extract smooth latent trajectories via EM with SE-GP priors
gpfa_transform(new_trains, params, bin_ms, dt) Project new data using learned GPFA parameters

Explainability

sc_neurocore.analysis.explainability

SpikeToConceptMapper

XAI Module: Maps spike patterns to semantic concepts.

Source code in src/sc_neurocore/analysis/explainability.py
Python
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
class SpikeToConceptMapper:
    """
    XAI Module: Maps spike patterns to semantic concepts.
    """

    def __init__(self, concept_map: Dict[int, str]):
        self.concept_map = concept_map

    def explain(self, spikes: np.ndarray[Any, Any]) -> str:
        """
        Input: Spike vector (n_neurons,)
        Output: "The agent is thinking about [Concept1, Concept2]"
        """
        active_indices = np.where(spikes > 0)[0]

        concepts = []
        for idx in active_indices:
            if idx in self.concept_map:
                concepts.append(self.concept_map[idx])
            else:
                concepts.append(f"Unknown({idx})")

        if not concepts:
            return "The agent is idle."

        return f"The agent is active on: {', '.join(concepts)}"

explain(spikes)

Input: Spike vector (n_neurons,) Output: "The agent is thinking about [Concept1, Concept2]"

Source code in src/sc_neurocore/analysis/explainability.py
Python
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
def explain(self, spikes: np.ndarray[Any, Any]) -> str:
    """
    Input: Spike vector (n_neurons,)
    Output: "The agent is thinking about [Concept1, Concept2]"
    """
    active_indices = np.where(spikes > 0)[0]

    concepts = []
    for idx in active_indices:
        if idx in self.concept_map:
            concepts.append(self.concept_map[idx])
        else:
            concepts.append(f"Unknown({idx})")

    if not concepts:
        return "The agent is idle."

    return f"The agent is active on: {', '.join(concepts)}"

Integrated Information (Phi*)

Approximate Phi (Barrett & Seth 2011) measures how much a system exceeds the sum of its parts — the consciousness-relevant integrated information.

sc_neurocore.analysis.phi_estimation.phi_star(data, tau=1)

Geometric integrated information (Barrett & Seth 2011).

Phi* = MI(past; future) - max_partition sum MI(past_k; future_k)

Approximated via Gaussian assumption: MI = 0.5 * log(det(Sigma_X) * det(Sigma_Y) / det(Sigma_XY)).

Parameters

data : np.ndarray Shape (n_channels, n_timesteps). Spike counts or continuous signals. tau : int Time lag for past→future mapping.

Returns

float Phi* estimate (bits). Non-negative; 0 = fully reducible.

Source code in src/sc_neurocore/analysis/phi_estimation.py
Python
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
def phi_star(data: np.ndarray, tau: int = 1) -> float:
    """Geometric integrated information (Barrett & Seth 2011).

    Phi* = MI(past; future) - max_partition sum MI(past_k; future_k)

    Approximated via Gaussian assumption: MI = 0.5 * log(det(Sigma_X) * det(Sigma_Y) / det(Sigma_XY)).

    Parameters
    ----------
    data : np.ndarray
        Shape (n_channels, n_timesteps). Spike counts or continuous signals.
    tau : int
        Time lag for past→future mapping.

    Returns
    -------
    float
        Phi* estimate (bits). Non-negative; 0 = fully reducible.
    """
    n, T = data.shape
    if 2 * tau >= T or n < 2:
        return 0.0

    past = data[:, :-tau]
    future = data[:, tau:]

    # Joint mutual information I(past; future)
    mi_whole = _gaussian_mi(past, future)

    # Minimum information partition: try all bipartitions
    # For tractability, only try contiguous splits (first k vs rest)
    mi_parts_min = np.inf
    for k in range(1, n):
        idx_a = list(range(k))
        idx_b = list(range(k, n))
        mi_a = _gaussian_mi(past[idx_a], future[idx_a])
        mi_b = _gaussian_mi(past[idx_b], future[idx_b])
        mi_parts_min = min(mi_parts_min, mi_a + mi_b)

    phi = max(0.0, mi_whole - mi_parts_min)
    return float(phi)

sc_neurocore.analysis.phi_estimation.phi_from_spike_trains(spikes, bin_size=10, tau=1)

Compute Phi* from binary spike trains.

Parameters

spikes : np.ndarray Shape (n_neurons, n_timesteps), binary {0, 1}. bin_size : int Number of timesteps per bin for spike count computation. tau : int Time lag in bins.

Returns

float Phi* in bits.

Source code in src/sc_neurocore/analysis/phi_estimation.py
Python
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
def phi_from_spike_trains(spikes: np.ndarray, bin_size: int = 10, tau: int = 1) -> float:
    """Compute Phi* from binary spike trains.

    Parameters
    ----------
    spikes : np.ndarray
        Shape (n_neurons, n_timesteps), binary {0, 1}.
    bin_size : int
        Number of timesteps per bin for spike count computation.
    tau : int
        Time lag in bins.

    Returns
    -------
    float
        Phi* in bits.
    """
    n_neurons, n_steps = spikes.shape
    n_bins = n_steps // bin_size
    if n_bins < 2 * tau + 2:
        return 0.0

    # Bin spike trains into spike counts
    binned = np.zeros((n_neurons, n_bins))
    for b in range(n_bins):
        binned[:, b] = spikes[:, b * bin_size : (b + 1) * bin_size].sum(axis=1)

    return phi_star(binned, tau=tau)