BMTK Stimulus Generation Tutorial¶

This notebook demonstrates how to use the StimulusBuilder to generate neural stimulus inputs for BMTK simulations, including baseline activity, population-specific shell input, and assembly-based patterned stimuli with analysis and visualization.

Import Required Libraries¶

In [ ]:

import os
import sys
import numpy as np
import pandas as pd
import importlib
import h5py

from bmtool.stimulus.core import StimulusBuilder
from bmtool.util import util
from bmtk.utils.reports.spike_trains import PoissonSpikeGenerator

from bmtool.analysis.spikes import load_spikes_to_df, compute_firing_rate_stats, get_population_spike_rate
from bmtool.bmplot.spikes import plot_firing_rate_histogram, plot_firing_rate_distribution, raster
import matplotlib.pyplot as plt
import json

import os import sys import numpy as np import pandas as pd import importlib import h5py from bmtool.stimulus.core import StimulusBuilder from bmtool.util import util from bmtk.utils.reports.spike_trains import PoissonSpikeGenerator from bmtool.analysis.spikes import load_spikes_to_df, compute_firing_rate_stats, get_population_spike_rate from bmtool.bmplot.spikes import plot_firing_rate_histogram, plot_firing_rate_distribution, raster import matplotlib.pyplot as plt import json

Configure Paths and Seeds¶

Set up the configuration file path, output directory, and random seeds for reproducibility.

In [ ]:

# IMPORTANT: Update this path to point to your BMTK configuration file
config_path = 'path/to/your/simulation_config.json'  # Adjust this path
output_dir = 'stimulus_output'
os.makedirs(output_dir, exist_ok=True)

# Set seeds for reproducible data generation
np.random.seed(42)
net_seed = 123  # Seed for assembly generation
psg_seed = 2    # Seed for Poisson spike generation

# IMPORTANT: Update this path to point to your BMTK configuration file config_path = 'path/to/your/simulation_config.json' # Adjust this path output_dir = 'stimulus_output' os.makedirs(output_dir, exist_ok=True) # Set seeds for reproducible data generation np.random.seed(42) net_seed = 123 # Seed for assembly generation psg_seed = 2 # Seed for Poisson spike generation

Initialize StimulusBuilder¶

Create a StimulusBuilder instance with the configuration file and random seeds.

In [ ]:

sb = StimulusBuilder(config=config_path, net_seed=net_seed, psg_seed=psg_seed)

sb = StimulusBuilder(config=config_path, net_seed=net_seed, psg_seed=psg_seed)

Generate Baseline Activity¶

Create baseline spiking activity across all nodes with a lognormal firing rate distribution.

In [ ]:

baseline_file = os.path.join(output_dir, 'baseline.h5')

baseline_params = {
    'base': {'mean_firing_rate': 20.0, 'stdev': 2.0}
}

sb.generate_background(
    output_path=baseline_file,
    network_name='baseline',
    population_params=baseline_params,
    t_start=0.0,
    t_stop=15.0
)

baseline_file = os.path.join(output_dir, 'baseline.h5') baseline_params = { 'base': {'mean_firing_rate': 20.0, 'stdev': 2.0} } sb.generate_background( output_path=baseline_file, network_name='baseline', population_params=baseline_params, t_start=0.0, t_stop=15.0 )

Verify Shell Input Statistics¶

Load shell spikes and visualize population-specific firing rate statistics. Each population (PN, PV, SST) can have its own firing rate distribution.

In [ ]:

shell_file = os.path.join(output_dir, 'shell.h5')

shell_params = {
    'PN': {'mean_firing_rate': 1.9, 'stdev': 1.8},
    'PV': {'mean_firing_rate': 7.5, 'stdev': 6.4},
    'SST': {'mean_firing_rate': 5.0, 'stdev': 6.0}
}

sb.generate_background(
    output_path=shell_file,
    network_name='shell',
    population_params=shell_params,
    t_start=0.0,
    t_stop=15.0
)

shell_file = os.path.join(output_dir, 'shell.h5') shell_params = { 'PN': {'mean_firing_rate': 1.9, 'stdev': 1.8}, 'PV': {'mean_firing_rate': 7.5, 'stdev': 6.4}, 'SST': {'mean_firing_rate': 5.0, 'stdev': 6.0} } sb.generate_background( output_path=shell_file, network_name='shell', population_params=shell_params, t_start=0.0, t_stop=15.0 )

Generate Shell Input¶

Create population-specific background activity with different firing rates for PN, PV, and SST cell types.

In [ ]:

shell_file = os.path.join(output_dir, 'shell.h5')

shell_params = {
    'PN': (1.9, 1.8),
    'PV': (7.5, 6.4),
    'SST': (5.0, 6.0)
}

sb.generate_shell_input(
    output_path=shell_file,
    network_name='shell', 
    shell_params=shell_params,
    t_start=0.0,
    t_stop=15.0,
    distribution='lognormal'
)

shell_file = os.path.join(output_dir, 'shell.h5') shell_params = { 'PN': (1.9, 1.8), 'PV': (7.5, 6.4), 'SST': (5.0, 6.0) } sb.generate_shell_input( output_path=shell_file, network_name='shell', shell_params=shell_params, t_start=0.0, t_stop=15.0, distribution='lognormal' )

Verify Shell Input Statistics¶

Load shell input spikes and visualize population-specific firing rate statistics.

In [ ]:

shell = os.path.join(output_dir, 'shell.h5')
spikes_df = load_spikes_to_df(shell, network_name='shell', config=config_path)
pop_stats, individual_stats = compute_firing_rate_stats(spikes_df)
plot_firing_rate_distribution(individual_stats, groupby='pop_name', plot_type='box', logscale=False)
plt.show()

shell = os.path.join(output_dir, 'shell.h5') spikes_df = load_spikes_to_df(shell, network_name='shell', config=config_path) pop_stats, individual_stats = compute_firing_rate_stats(spikes_df) plot_firing_rate_distribution(individual_stats, groupby='pop_name', plot_type='box', logscale=False) plt.show()

Create Thalamic Assemblies¶

Partition thalamic neurons into assemblies based on the pulse_group_id property, creating separate groups for each pulse pattern.

In [ ]:

sb.create_assemblies(
    name='thalamus_groups',
    network_name='thalamus',
    method='property',
    property_name='pulse_group_id'
)

assembly_list = sb.assemblies['thalamus_groups']
n_assemblies = len(assembly_list)
print(f"Created {n_assemblies} assemblies from pulse_group_id")

sb.create_assemblies( name='thalamus_groups', network_name='thalamus', method='property', property_name='pulse_group_id' ) assembly_list = sb.assemblies['thalamus_groups'] n_assemblies = len(assembly_list) print(f"Created {n_assemblies} assemblies from pulse_group_id")

Generate Long Stimulus Pattern¶

Create a stimulus with the 'long' pattern, where each assembly receives a sustained burst during its designated cycle.

In [ ]:

long_file = os.path.join(output_dir, 'long.h5')

# First, create assemblies from node properties
sb.create_assemblies(
    name='thalamus_groups',
    network_name='thalamus',
    method='property',
    property_name='pulse_group_id'
)

assembly_list = sb.assemblies['thalamus_groups']
n_assemblies = len(assembly_list)
print(f"Created {n_assemblies} assemblies from pulse_group_id")

# Then generate stimulus using the assembly name
firing_rate_params = (0.0, 50.0, 0.0)  # (off_rate, burst_rate, silent_rate)

sb.generate_stimulus(
    output_path=long_file,
    pattern_type='long',
    assembly_name='thalamus_groups',
    population='thalamus',
    firing_rate=firing_rate_params,
    t_start=1.0,
    t_stop=15.0,
    on_time=1.0,
    off_time=0.5
)

long_file = os.path.join(output_dir, 'long.h5') # First, create assemblies from node properties sb.create_assemblies( name='thalamus_groups', network_name='thalamus', method='property', property_name='pulse_group_id' ) assembly_list = sb.assemblies['thalamus_groups'] n_assemblies = len(assembly_list) print(f"Created {n_assemblies} assemblies from pulse_group_id") # Then generate stimulus using the assembly name firing_rate_params = (0.0, 50.0, 0.0) # (off_rate, burst_rate, silent_rate) sb.generate_stimulus( output_path=long_file, pattern_type='long', assembly_name='thalamus_groups', population='thalamus', firing_rate=firing_rate_params, t_start=1.0, t_stop=15.0, on_time=1.0, off_time=0.5 )

Visualize Long Stimulus Response¶

Create a helper function to analyze and visualize stimulus-evoked responses with raster plots and firing rate time series for each pulse group.

In [ ]:

def visualize_stimulus(h5_path, network_name='thalamus', pulse_group_key='pulse_group_id'):
    """
    Load stimulus spikes and create firing rate and raster visualizations.
    
    Args:
        h5_path: Path to the stimulus HDF5 file
        network_name: Name of the network in the file
        pulse_group_key: Column name for grouping (e.g., 'pulse_group_id')
    """
    spikes_df = load_spikes_to_df(h5_path, network_name=network_name, config=config_path, 
                                   groupby=['pop_name', pulse_group_key])

    # Get unique pulse group IDs (skip NaNs)
    pulse_groups = sorted([x for x in spikes_df[pulse_group_key].unique() if not pd.isna(x)])
    if not pulse_groups:
        pulse_groups = [None]
        
    populations = sorted(spikes_df['pop_name'].unique())

    # Create grid of subplots
    n_groups = len(pulse_groups)
    cols = min(3, n_groups)
    rows = (n_groups + cols - 1) // cols if n_groups > 0 else 1
    
    fig, axes = plt.subplots(rows, cols, figsize=(15, 4 * rows), squeeze=False)
    axes = axes.flatten()

    # Plot firing rate time series for each pulse group
    for pulse_idx, pulse_group in enumerate(pulse_groups):
        ax = axes[pulse_idx]
        
        # Filter data by pulse group
        if pulse_group is not None:
            mask = (spikes_df[pulse_group_key] == pulse_group)
            label = f'Pulse Group {pulse_group}'
        else:
            mask = slice(None)
            label = 'All nodes'
            
        filtered_df = spikes_df[mask]
        if filtered_df.empty: 
            continue
        
        # Compute and plot spike rate
        spike_rate_subset = get_population_spike_rate(filtered_df, fs=400, network_name=network_name)
        spike_rate_subset.sel(type='smoothed').plot(ax=ax, label=label)
        ax.set_title(f'{network_name} - {label}')
        ax.set_ylabel('Firing Rate (Hz)')
        ax.set_xlabel('Time (ms)')
        ax.grid(True, alpha=0.3)

    # Hide unused subplots
    for idx in range(n_groups, len(axes)):
        axes[idx].set_visible(False)

    plt.tight_layout()
    plt.show()

    # Create spike raster plot
    fig, ax = plt.subplots(figsize=(15, 8))
    
    # Sort nodes for better visualization
    sort_cols = []
    if pulse_group_key in spikes_df.columns: 
        sort_cols.append(pulse_group_key)
    sort_cols.append('node_ids')
    
    sorted_nodes = spikes_df.sort_values(sort_cols)['node_ids'].unique()
    node_to_y = {node_id: i for i, node_id in enumerate(sorted_nodes)}
    
    spikes_df_plot = spikes_df.copy()
    spikes_df_plot['node_y'] = spikes_df_plot['node_ids'].map(node_to_y)
    
    # Color code by population
    unique_pops = spikes_df_plot['pop_name'].unique()
    cmap = plt.get_cmap("tab10")
    color_map = {pop: cmap(i / len(unique_pops)) for i, pop in enumerate(unique_pops)}
    
    for pop_name, group in spikes_df_plot.groupby('pop_name'):
        ax.scatter(group["timestamps"], group['node_y'], color=color_map[pop_name], s=0.5, alpha=0.6)
    
    handles = [ax.scatter([], [], color=color_map[pop], label=pop, s=20) for pop in unique_pops]
    ax.legend(handles=handles, title="Population", loc="upper right")
    
    ax.set_xlabel('Time (ms)')
    ax.set_ylabel('Node ID (sorted)')
    ax.set_title(f'{network_name} Spike Raster')
    
    plt.tight_layout()
    plt.show()


# Visualize the long pattern stimulus
visualize_stimulus(long_file, network_name='thalamus')

def visualize_stimulus(h5_path, network_name='thalamus', pulse_group_key='pulse_group_id'): """ Load stimulus spikes and create firing rate and raster visualizations. Args: h5_path: Path to the stimulus HDF5 file network_name: Name of the network in the file pulse_group_key: Column name for grouping (e.g., 'pulse_group_id') """ spikes_df = load_spikes_to_df(h5_path, network_name=network_name, config=config_path, groupby=['pop_name', pulse_group_key]) # Get unique pulse group IDs (skip NaNs) pulse_groups = sorted([x for x in spikes_df[pulse_group_key].unique() if not pd.isna(x)]) if not pulse_groups: pulse_groups = [None] populations = sorted(spikes_df['pop_name'].unique()) # Create grid of subplots n_groups = len(pulse_groups) cols = min(3, n_groups) rows = (n_groups + cols - 1) // cols if n_groups > 0 else 1 fig, axes = plt.subplots(rows, cols, figsize=(15, 4 * rows), squeeze=False) axes = axes.flatten() # Plot firing rate time series for each pulse group for pulse_idx, pulse_group in enumerate(pulse_groups): ax = axes[pulse_idx] # Filter data by pulse group if pulse_group is not None: mask = (spikes_df[pulse_group_key] == pulse_group) label = f'Pulse Group {pulse_group}' else: mask = slice(None) label = 'All nodes' filtered_df = spikes_df[mask] if filtered_df.empty: continue # Compute and plot spike rate spike_rate_subset = get_population_spike_rate(filtered_df, fs=400, network_name=network_name) spike_rate_subset.sel(type='smoothed').plot(ax=ax, label=label) ax.set_title(f'{network_name} - {label}') ax.set_ylabel('Firing Rate (Hz)') ax.set_xlabel('Time (ms)') ax.grid(True, alpha=0.3) # Hide unused subplots for idx in range(n_groups, len(axes)): axes[idx].set_visible(False) plt.tight_layout() plt.show() # Create spike raster plot fig, ax = plt.subplots(figsize=(15, 8)) # Sort nodes for better visualization sort_cols = [] if pulse_group_key in spikes_df.columns: sort_cols.append(pulse_group_key) sort_cols.append('node_ids') sorted_nodes = spikes_df.sort_values(sort_cols)['node_ids'].unique() node_to_y = {node_id: i for i, node_id in enumerate(sorted_nodes)} spikes_df_plot = spikes_df.copy() spikes_df_plot['node_y'] = spikes_df_plot['node_ids'].map(node_to_y) # Color code by population unique_pops = spikes_df_plot['pop_name'].unique() cmap = plt.get_cmap("tab10") color_map = {pop: cmap(i / len(unique_pops)) for i, pop in enumerate(unique_pops)} for pop_name, group in spikes_df_plot.groupby('pop_name'): ax.scatter(group["timestamps"], group['node_y'], color=color_map[pop_name], s=0.5, alpha=0.6) handles = [ax.scatter([], [], color=color_map[pop], label=pop, s=20) for pop in unique_pops] ax.legend(handles=handles, title="Population", loc="upper right") ax.set_xlabel('Time (ms)') ax.set_ylabel('Node ID (sorted)') ax.set_title(f'{network_name} Spike Raster') plt.tight_layout() plt.show() # Visualize the long pattern stimulus visualize_stimulus(long_file, network_name='thalamus')

Generate Short Stimulus Pattern¶

Create a stimulus with the 'short' pattern, where multiple brief bursts are delivered sequentially to each assembly within each cycle.

In [ ]:

short_file = os.path.join(output_dir, 'short.h5')
firing_rate_params = (0.0, 50.0, 0.0)  # (off_rate, burst_rate, silent_rate)

sb.generate_stimulus(
    output_path=short_file,
    pattern_type='short',
    assembly_name='thalamus_groups',  # Use the same assembly created earlier
    population='thalamus',
    firing_rate=firing_rate_params,
    t_start=1.0,
    t_stop=15.0,
    on_time=1.0,
    off_time=0.5
)

short_file = os.path.join(output_dir, 'short.h5') firing_rate_params = (0.0, 50.0, 0.0) # (off_rate, burst_rate, silent_rate) sb.generate_stimulus( output_path=short_file, pattern_type='short', assembly_name='thalamus_groups', # Use the same assembly created earlier population='thalamus', firing_rate=firing_rate_params, t_start=1.0, t_stop=15.0, on_time=1.0, off_time=0.5 )

Visualize Short Stimulus Response¶

Analyze and compare the short stimulus pattern response across different pulse groups.

In [ ]:

visualize_stimulus(short_file, network_name='thalamus')

visualize_stimulus(short_file, network_name='thalamus')

Summary¶

This tutorial demonstrated the key features of the StimulusBuilder module:

1. Unified Background Generation: Use generate_background() with population_params to create activity grouped by any node property (pop_name, layer, etc).

2. Flexible Distribution Types: Each group can independently use either constant rates or statistical distributions (lognormal, normal). Omit stdev for constant-rate groups.

3. Property-Based Grouping: Control how nodes are grouped using the groupby parameter (default: 'pop_name'). This allows grouping by cell type, layer, region, or custom properties.

4. Assembly-Based Stimulus Generation: Partition neurons into groups using multiple strategies (random, property-based, grid-based) and deliver coordinated stimulus patterns.

5. Firing Rate Patterns: Use different temporal patterns ('long' vs 'short') to deliver precisely timed stimuli to neural assemblies.

6. Validation and Analysis: Verify generated stimuli using firing rate statistics, histograms, and raster plots to ensure they meet your experimental design requirements.

Explore other firing patterns ('ramp', 'join', 'fade', 'loop') in the Stimulus Module Overview to further customize stimulus delivery for your simulations.