braintools.analysis module

Neuroscience data analysis.

braintools.analysis groups the measures used to analyse recorded or simulated neural activity. Functions live in six subpackages and are always accessed through them (they are not re-exported at this level):

Subpackage	Contents
braintools.analysis.spike	ISI statistics (CV, CV2, LV, Fano), PSTH, correlograms, spike-triggered average, firing rates, spike-train distances, synchrony
braintools.analysis.spectral	Welch and multitaper spectra, spectrograms, Morlet wavelets, band-pass filtering, band power, PLV/PLI/wPLI, coherence, cross-frequency coupling
braintools.analysis.field	LFP forward model, current source density
braintools.analysis.population	PCA, CCA, RSA, decoding, mutual information, Kuramoto order, Granger causality, avalanches, branching ratio, functional connectivity
braintools.analysis.stats	correlation helpers, spike jitter, ISI shuffling, phase randomisation, bootstrap confidence intervals
braintools.analysis.graph	structural graph measures of connectivity: degree distribution, clustering, path length, modularity, motifs, spectral radius, small-world index

Apart from graph, these functions were previously part of braintools.metric. The old names (for example braintools.metric.firing_rate) still work in 0.4 but emit a BraintoolsDeprecationWarning when called and will be removed in 1.0.

Time conventions. Time arguments accept brainunit.Quantity values, which are converted automatically. Plain numbers follow each function’s documented convention: spike-matrix measures take dt in milliseconds (the brainstate convention), whereas firing_rate(), phase_locking_value() and every function in spectral take seconds. Passing quantities avoids having to remember which.

Examples

Spike-train analysis on a spike matrix of shape (n_time, n_neurons):

>>> import brainunit as u
>>> import jax.numpy as jnp
>>> from braintools.analysis import spike
>>> spikes = jnp.zeros((200, 5)).at[50:53, :].set(1.0)
>>> rate = spike.firing_rate(spikes, width=10 * u.ms, dt=1 * u.ms)
>>> rate.shape
(200,)
>>> sync = spike.spike_train_synchrony(spikes, window_size=5 * u.ms, dt=1 * u.ms)
>>> round(float(sync), 3)
1.0

Spectral analysis of an LFP sampled every millisecond:

>>> from braintools.analysis import spectral
>>> t = jnp.arange(4000) * 0.001  # seconds
>>> lfp = jnp.sin(2 * jnp.pi * 10 * t)
>>> freqs, psd = spectral.power_spectral_density(lfp, dt=1 * u.ms, nperseg=1000)
>>> float(freqs[jnp.argmax(psd)])
10.0

Population and statistics helpers:

>>> from braintools.analysis import population, stats
>>> activities = jnp.column_stack([jnp.sin(t), jnp.sin(t), jnp.cos(t)])
>>> population.functional_connectivity(activities).shape
(3, 3)
>>> fc = population.functional_connectivity(activities)
>>> round(float(stats.matrix_correlation(fc, fc)), 3)
1.0

See also

For runnable, narrative walkthroughs, see the tutorials:

• Firing Rates and Population Synchrony

• Spike-Train Distances and Synchrony

• Local Field Potential Analysis

• The analysis toolbox: spikes, spectra and populations

Spike Trains

Spike-train analysis: statistics, PSTHs, correlograms, distances and synchrony.

Two input conventions are used:

• event-based measures take a braintools.spike.SpikeTrains (build one from a dense raster with from_dense());

• the rate and population-synchrony measures take a spike matrix of shape (n_time_steps, n_neurons) (time on axis 0) with the sampling step dt.

• Interval statistics: isi(), cv(), cv2(), lv(), fano_factor().

• Event-aligned measures: psth(), spike_triggered_average().

• Correlograms: cross_correlogram(), auto_correlogram().

• Rates: firing_rate().

• Distances between two trains: isi_distance(), spike_distance(), victor_purpura_distance(), van_rossum_distance().

• Population synchrony: cross_correlation(), spike_train_synchrony(), burst_synchrony_index(), spike_time_tiling_coefficient(), correlation_index().

• Locking to a rhythm: phase_locking_value().

Examples

>>> import brainunit as u
>>> import jax.numpy as jnp
>>> from braintools.analysis import spike
>>> spikes = jnp.zeros((100, 4)).at[20:25, :].set(1.0)
>>> rate = spike.firing_rate(spikes, width=5 * u.ms, dt=1 * u.ms)
>>> rate.shape
(100,)
>>> float(spike.victor_purpura_distance(jnp.array([1.0, 2.0]), jnp.array([1.0, 2.0])))
0.0

Interval statistics from a SpikeTrains:

>>> from braintools.spike import SpikeTrains
>>> raster = jnp.zeros((100, 2)).at[::10, 0].set(1.0).at[::25, 1].set(1.0)
>>> trains = SpikeTrains.from_dense(raster, dt=1 * u.ms)
>>> spike.cv(trains)
Array([0., 0.], dtype=float32)

firing_rate	Calculate the smoothed population firing rate from spike data.
victor_purpura_distance	Calculate Victor-Purpura distance between two spike trains.
van_rossum_distance	Calculate van Rossum distance between two spike trains.
spike_train_synchrony	Calculate spike train synchrony using the SPIKE-synchronization measure.
burst_synchrony_index	Calculate burst synchrony index based on co-occurring burst events.
phase_locking_value	Calculate phase-locking value (PLV) for spike synchronization.
spike_time_tiling_coefficient	Calculate Spike Time Tiling Coefficient (STTC).
correlation_index	Calculate correlation index for spike train synchrony.
cross_correlation	Calculate cross-correlation index between neurons.

Event-based measures on braintools.spike.SpikeTrains: interval statistics, event-aligned averages, correlograms and spike-train distances:

isi	Compute the inter-spike intervals of every neuron.
cv	Compute the coefficient of variation of the inter-spike intervals.
cv2	Compute the local coefficient of variation CV2 of Holt et al. (1996).
lv	Compute the local variation LV of Shinomoto et al. (2003).
fano_factor	Compute the Fano factor, the variance-to-mean ratio of spike counts.
psth	Compute the event-aligned peri-stimulus time histogram (PSTH).
spike_triggered_average	Average a sampled signal around the spikes of every neuron.
cross_correlogram	Compute the spike-time cross-correlogram of paired neurons.
auto_correlogram	Compute the spike-time auto-correlogram of every neuron.
isi_distance	Compute the ISI-distance between two populations, neuron by neuron.
spike_distance	Compute the SPIKE-distance between two populations, neuron by neuron.

Spectral Analysis

Spectral analysis of continuous signals such as local field potentials.

All functions take the sampling interval dt (a plain number is in seconds; a brainunit.Quantity is converted) and return frequencies in Hz.

• Spectra: power_spectral_density() (Welch), multitaper_psd(), spectrogram(), band_power(), spectral_entropy().

• Time-frequency and filtering: morlet_transform(), bandpass(), analytic_signal().

• Coupling between signals: plv(), pli(), wpli(), coherence_analysis(), lfp_phase_coherence().

• Cross-frequency coupling: phase_amplitude_coupling(), theta_gamma_coupling().

Examples

>>> import brainunit as u
>>> import jax.numpy as jnp
>>> from braintools.analysis import spectral
>>> t = jnp.arange(4000) * 0.001
>>> lfp = jnp.sin(2 * jnp.pi * 10 * t)
>>> freqs, psd = spectral.power_spectral_density(lfp, dt=1 * u.ms, nperseg=1000)
>>> float(freqs[jnp.argmax(psd)])
10.0

power_spectral_density	Estimate the one-sided power spectral density (PSD) using Welch's method.
coherence_analysis	Compute the magnitude-squared coherence between two LFP signals (Welch).
phase_amplitude_coupling	Compute phase-amplitude coupling (PAC) via the Tort modulation index.
theta_gamma_coupling	Compute theta-gamma coupling strength using standard frequency bands.
spectral_entropy	Compute the normalized spectral entropy of an LFP signal.
lfp_phase_coherence	Compute pairwise phase coherence between LFP channels in a frequency band.

Multitaper spectra, time-frequency analysis, filtering and phase synchrony between continuous signals:

spectrogram	Compute the short-time Fourier power spectrogram along time axis 0.
multitaper_psd	Estimate the power spectral density with DPSS (Slepian) multitapers.
morlet_transform	Compute the complex Morlet wavelet transform along time axis 0.
bandpass	Zero-phase Butterworth band-pass filter along time axis 0.
analytic_signal	Compute the analytic signal by the FFT Hilbert transform.
band_power	Integrate the Welch power spectral density over a frequency band.
plv	Compute the phase-locking value between two continuous signals.
pli	Compute the phase-lag index between two continuous signals.
wpli	Compute the weighted phase-lag index between two continuous signals.

Field Potentials

Extracellular field potentials: forward models and laminar analysis.

• unitary_LFP() builds the local field potential generated by a spiking population with the kernel method of Telenczuk et al. (2020).

• current_source_density() estimates current sources and sinks from a laminar LFP recording (second spatial difference).

Examples

>>> import jax.numpy as jnp
>>> from braintools.analysis import field
>>> lfp_laminar = jnp.tile(jnp.sin(jnp.linspace(0, jnp.pi, 16)), (50, 1))
>>> csd = field.current_source_density(lfp_laminar, electrode_spacing=0.1)
>>> csd.shape
(50, 14)

unitary_LFP	Calculate unitary local field potentials (uLFP) from spike train data.
current_source_density	Compute current source density (CSD) from laminar LFP recordings.

Population Activity

Population-level analysis of continuous multichannel activity.

The functions take activity of shape (n_time, n_signals), such as membrane potentials, firing rates or regional signals.

• voltage_fluctuation(): Golomb’s \(\chi^2\) synchrony index.

• functional_connectivity(): pairwise Pearson correlation matrix.

• functional_connectivity_dynamics(): correlation between sliding-window connectivity matrices.

• Dimensionality and alignment: pca(), participation_ratio(), cca(), rdm(), rsa().

• Decoding and information: decode(), mutual_information().

• Collective dynamics: kuramoto_order_parameter(), granger_causality(), avalanches(), branching_ratio().

Examples

>>> import jax.numpy as jnp
>>> from braintools.analysis import population
>>> t = jnp.linspace(0, 10, 200)
>>> activities = jnp.column_stack([jnp.sin(t), jnp.sin(t), -jnp.sin(t)])
>>> fc = population.functional_connectivity(activities)
>>> fc.shape
(3, 3)
>>> round(float(fc[0, 2]), 3)
-1.0

voltage_fluctuation	Calculate neuronal synchronization via voltage variance analysis.
functional_connectivity	Compute functional connectivity matrix from time series data.
functional_connectivity_dynamics	Compute functional connectivity dynamics (FCD) matrix.

Dimensionality, representational alignment, decoding and information:

pca	Principal component analysis by mean-centred SVD.
participation_ratio	Participation ratio (effective dimensionality) of population activity.
cca	Canonical correlation analysis via SVD of the whitened cross-covariance.
rdm	Representational dissimilarity matrix between conditions.
rsa	Representational similarity: correlation between two RDMs.
decode	Cross-validated linear decoding of integer labels from population activity.
mutual_information	Plug-in mutual information (in nats) between two discrete variables.

Collective dynamics:

kuramoto_order_parameter	Compute the Kuramoto order parameter of a population of oscillators.
granger_causality	Compute pairwise conditional Granger causality from a VAR model.
avalanches	Detect neuronal avalanches in a population count series.
branching_ratio	Estimate the branching ratio of a population count series.

Result containers (named tuples, documented by their Attributes):

PCAResult	Result of pca() (a pytree, so it passes through jax.jit).
CCAResult	Result of cca() (a pytree, so it passes through jax.jit).
AvalancheResult	Neuronal avalanches detected in a population count series.

Statistics

General statistics used across neural data analysis.

• matrix_correlation(): Pearson correlation between the upper triangles of two square matrices (for example, two connectivity matrices).

• weighted_correlation(): weighted Pearson correlation of two series.

• Resampling (building blocks for significance tests, often called “surrogate data” in the literature): jitter_spikes(), shuffle_isi(), phase_randomize().

• bootstrap_ci(): bootstrap confidence intervals of any statistic.

Examples

>>> import jax.numpy as jnp
>>> from braintools.analysis import stats
>>> x = jnp.array([1.0, 2.0, 3.0, 4.0, 5.0])
>>> w = jnp.array([1.0, 1.0, 2.0, 2.0, 1.0])
>>> round(float(stats.weighted_correlation(x, 2 * x, w)), 3)
1.0

matrix_correlation	Compute Pearson correlation of upper triangular elements of two matrices.
weighted_correlation	Compute weighted Pearson correlation between two data series.

Resampling for significance testing, and bootstrap confidence intervals:

jitter_spikes	Randomly displace every spike by an independent offset.
shuffle_isi	Randomly permute the inter-spike intervals of every neuron.
phase_randomize	Randomise the Fourier phases of a real signal, keeping its power spectrum.
bootstrap_ci	Bootstrap confidence interval of a statistic.

Graph Measures

Structural graph measures of connectivity.

The functions take a braintools.conn.ConnectionResult, a SciPy sparse matrix, or a dense adjacency matrix (NumPy, JAX or a nested list), and work on sparse matrices throughout, so they scale to large networks without building an N x N array. They run on the host and are not jit/vmap/grad compatible. Self-connections are handled as in networkx: they are counted by degree_distribution() (once in the in- and out-degree, twice in the total degree), modularity() and the reciprocity of motif_counts(), and kept by spectral_radius(); they are ignored by clustering_coefficient(), characteristic_path_length(), the triad census of motif_counts() and small_world_index().

• degree_distribution(): in-, out- or total degrees and their histogram.

• clustering_coefficient(): directed (Fagiolo) or undirected clustering, optionally weighted.

• characteristic_path_length(): mean shortest path length and global efficiency, exact or from sampled sources.

• modularity(): Newman (undirected) or Leicht-Newman (directed) modularity of a partition.

• motif_counts(): the 16-class triad census, dyad counts and reciprocity.

• spectral_radius(): largest eigenvalue magnitude of the weight matrix.

• small_world_index(): small-world coefficient sigma against Erdos-Renyi graphs with the same numbers of nodes and edges.

Examples

>>> import brainunit as u
>>> from braintools.conn import FixedProb
>>> from braintools.analysis import graph
>>> conn = FixedProb(0.1, weight=0.5 * u.nS, seed=0)(pre_size=200, post_size=200)
>>> dd = graph.degree_distribution(conn, direction='in')
>>> dd.degrees.shape
(200,)
>>> 0.0 < graph.clustering_coefficient(conn) < 0.2
True

degree_distribution	Compute node degrees and their histogram.
clustering_coefficient	Compute the local clustering coefficient, or its average over nodes.
characteristic_path_length	Compute the characteristic path length and the global efficiency.
modularity	Compute the modularity of a given partition.
motif_counts	Count the directed triad classes (triad census) and dyads.
spectral_radius	Compute the largest absolute eigenvalue of the connectivity matrix.
small_world_index	Compute the small-world index \(\sigma\) against random graphs.

Result containers (named tuples, documented by their Attributes):

DegreeDistribution	Result of degree_distribution().
PathLengthResult	Result of characteristic_path_length().
MotifCounts	Result of motif_counts().
SmallWorldResult	Result of small_world_index().