Analysis Module
The Analysis module provides a comprehensive suite of tools for processing and analyzing simulation results from BMTK models. It is organized into several submodules, each focusing on specific aspects of neural data analysis.
Overview
Spike Analysis
- Loading and processing spike data
- Computing firing rate statistics
- Population spike rate analysis
- Spike train analysis tools
LFP/ECP Analysis
- Loading and processing LFP/ECP data
- Spectral analysis and filtering
- Time-frequency analysis
- Signal quality metrics
Entrainment Analysis
- Phase-locking analysis
- Spike-field coherence
- Population entrainment metrics
- Cross-frequency coupling
Network Connectivity
- Connection statistics
- Network topology analysis
- Synaptic weight distributions
- Connectivity visualization
Spike Analysis
Load and analyze spike data from simulation output:
import pandas as pd
import matplotlib.pyplot as plt
from bmtool.analysis.spikes import load_spikes_to_df, compute_firing_rate_stats, get_population_spike_rate
from bmtool.bmplot import raster, plot_firing_rate_pop_stats, plot_firing_rate_distribution
# Load spike data from a simulation
spikes_df = load_spikes_to_df(
spike_file='output/spikes.h5',
network_name='network',
config='config.json', # Optional: to label cell types
groupby='pop_name'
)
# Get basic spike statistics by population
pop_stats, individual_stats = compute_firing_rate_stats(
df=spikes_df,
groupby='pop_name',
start_time=500,
stop_time=1500
)
print("Population firing rate statistics:")
print(pop_stats)
# Calculate population spike rates over time
population_rates = get_population_spike_rate(
spikes=spikes_df,
fs=400.0, # Sampling frequency in Hz
t_start=0,
t_stop=2000,
config='config.json', # Optional
network_name='network' # Optional
)
# Plot population rates
for pop_name, rates in population_rates.items():
plt.plot(rates, label=pop_name)
plt.xlabel('Time (ms)')
plt.ylabel('Firing Rate (Hz)')
plt.legend()
plt.title('Population Firing Rates')
plt.show()
Raster Plots
Create raster plots to visualize spike patterns using the BMPlot module:
import matplotlib.pyplot as plt
from bmtool.analysis.spikes import load_spikes_to_df
from bmtool.bmplot import raster
# Load spike data
spikes_df = load_spikes_to_df(
spike_file='output/spikes.h5',
network_name='network',
config='config.json'
)
# Create a basic raster plot
fig, ax = plt.subplots(figsize=(10, 6))
raster(
spikes_df=spikes_df,
groupby='pop_name',
time_range=(0, 2000),
ax=ax
)
plt.show()
# Plot firing rate statistics
fig, ax = plt.subplots(figsize=(10, 6))
plot_firing_rate_pop_stats(
firing_stats=pop_stats,
groupby='pop_name',
ax=ax
)
plt.show()
# Plot firing rate distributions
fig, ax = plt.subplots(figsize=(10, 6))
plot_firing_rate_distribution(
individual_stats=individual_stats,
groupby='pop_name',
ax=ax
)
plt.show()
LFP/ECP Analysis
Analyze LFP (Local Field Potential) and ECP (Extracellular Potential) data:
import numpy as np
import matplotlib.pyplot as plt
import xarray as xr
from bmtool.analysis.lfp import (
load_ecp_to_xarray,
ecp_to_lfp,
slice_time_series,
cwt_spectrogram_xarray,
plot_spectrogram,
butter_bandpass_filter,
fit_fooof
)
# Load ECP data
ecp_data = load_ecp_to_xarray('output/ecp.h5', demean=True)
# Convert ECP to LFP with filtering
lfp_data = ecp_to_lfp(
ecp_data=ecp_data,
cutoff=250, # Cutoff frequency in Hz
fs=10000 # Sampling frequency in Hz
)
# Slice data to specific time range
lfp_slice = slice_time_series(lfp_data, time_ranges=(500, 1500))
# Calculate spectrogram
spectrogram = cwt_spectrogram_xarray(
x=lfp_slice.sel(channel=0).data,
fs=10000,
freq_range=(1, 100),
nNotes=8
)
# Plot spectrogram
fig, ax = plt.subplots(figsize=(10, 6))
plot_spectrogram(
sxx_xarray=spectrogram,
log_power=True,
ax=ax
)
plt.show()
Frequency Analysis and Phase Locking
Analyze frequency content and phase locking between signals:
import numpy as np
import matplotlib.pyplot as plt
from bmtool.analysis.lfp import (
butter_bandpass_filter,
calculate_spike_lfp_plv,
calculate_signal_signal_plv,
calculate_ppc
)
from bmtool.analysis.spikes import load_spikes_to_df
# Load spike data and LFP data
spikes_df = load_spikes_to_df('output/spikes.h5', network_name='network')
lfp_data = load_ecp_to_xarray('output/ecp.h5', demean=True)
# Filter LFP to specific frequency band (e.g., theta: 4-8 Hz)
lfp_signal = lfp_data.sel(channel=0).data
fs = 10000 # Hz
filtered_lfp = butter_bandpass_filter(
data=lfp_signal,
lowcut=4,
highcut=8,
fs=fs
)
# Extract spike times for a specific population
population_spikes = spikes_df[spikes_df['pop_name'] == 'Pyramidal']
spike_times = population_spikes['timestamps'].to_numpy()
# Calculate phase-locking value between spikes and LFP
plv = calculate_spike_lfp_plv(
spike_times=spike_times,
lfp_signal=filtered_lfp,
spike_fs=1000, # Spike time unit in milliseconds
lfp_fs=fs,
fmin=4,
fmax=8
)
print(f"Phase-locking value: {plv}")
# Calculate pairwise phase consistency
ppc = calculate_ppc(
spike_times=spike_times,
lfp_signal=filtered_lfp,
spike_fs=1000,
lfp_fs=fs,
fmin=4,
fmax=8
)
print(f"Pairwise phase consistency: {ppc}")
Signal Processing
Apply filters and transformations to time series data:
import numpy as np
import matplotlib.pyplot as plt
from bmtool.analysis.lfp import (
butter_bandpass_filter,
wavelet_filter,
fit_fooof,
generate_resd_from_fooof,
calculate_SNR
)
# Apply band-pass filter
filtered_signal = butter_bandpass_filter(
data=lfp_signal,
lowcut=30,
highcut=80,
fs=10000
)
# Apply wavelet filter centered at a specific frequency
gamma_filtered = wavelet_filter(
x=lfp_signal,
freq=40, # Center frequency in Hz
fs=10000, # Sampling rate
bandwidth=10 # Bandwidth in Hz
)
# Calculate power spectrum and fit FOOOF model
from scipy import signal
# Calculate power spectrum
freqs, pxx = signal.welch(lfp_signal, fs=10000, nperseg=4096)
# Fit FOOOF model to extract oscillatory and aperiodic components
fooof_model = fit_fooof(
f=freqs,
pxx=pxx,
freq_range=[1, 100],
peak_width_limits=[1, 8],
max_n_peaks=6
)
# Get the residuals between the original spectrum and the aperiodic fit
resid_spectra, idx_freqs = generate_resd_from_fooof(fooof_model)
# Calculate signal-to-noise ratio in a specific frequency band
snr = calculate_SNR(fooof_model, freq_band=(30, 80))
print(f"Signal-to-noise ratio in gamma band: {snr}")