Spike Analysis

The spikes module provides functions for loading and analyzing spike data from simulations.

bmtool.analysis.spikes.load_spikes_to_df(spike_file, network_name, sort=True, config=None, groupby='pop_name')

Load spike data from an HDF5 file into a pandas DataFrame.

Parameters:

Name	Type	Description	Default
spike_file	str	Path to the HDF5 file containing spike data	required
network_name	str	The name of the network within the HDF5 file from which to load spike data	required
sort	bool	Whether to sort the DataFrame by 'timestamps' (default: True)	True
config	str	Path to configuration file to label the cell type of each spike (default: None)	None
groupby	Union[str, List[str]]	The column(s) to group by (default: 'pop_name')	'pop_name'

Returns:

Type	Description
DataFrame	A pandas DataFrame containing 'node_ids' and 'timestamps' columns from the spike data, with additional columns if a config file is provided

Examples:

>>> df = load_spikes_to_df("spikes.h5", "cortex")
>>> df = load_spikes_to_df("spikes.h5", "cortex", config="config.json", groupby=["pop_name", "model_type"])

Source code in bmtool/analysis/spikes.py

def load_spikes_to_df(
    spike_file: str,
    network_name: str,
    sort: bool = True,
    config: str = None,
    groupby: Union[str, List[str]] = "pop_name",
) -> pd.DataFrame:
    """
    Load spike data from an HDF5 file into a pandas DataFrame.

    Parameters
    ----------
    spike_file : str
        Path to the HDF5 file containing spike data
    network_name : str
        The name of the network within the HDF5 file from which to load spike data
    sort : bool, optional
        Whether to sort the DataFrame by 'timestamps' (default: True)
    config : str, optional
        Path to configuration file to label the cell type of each spike (default: None)
    groupby : Union[str, List[str]], optional
        The column(s) to group by (default: 'pop_name')

    Returns
    -------
    pd.DataFrame
        A pandas DataFrame containing 'node_ids' and 'timestamps' columns from the spike data,
        with additional columns if a config file is provided

    Examples
    --------
    >>> df = load_spikes_to_df("spikes.h5", "cortex")
    >>> df = load_spikes_to_df("spikes.h5", "cortex", config="config.json", groupby=["pop_name", "model_type"])
    """
    with h5py.File(spike_file) as f:
        spikes_df = pd.DataFrame(
            {
                "node_ids": f["spikes"][network_name]["node_ids"],
                "timestamps": f["spikes"][network_name]["timestamps"],
            }
        )

        if sort:
            spikes_df.sort_values(by="timestamps", inplace=True, ignore_index=True)

        if config:
            nodes = load_nodes_from_config(config)
            nodes = nodes[network_name]

            # Convert single string to a list for uniform handling
            if isinstance(groupby, str):
                groupby = [groupby]

            # Ensure all requested columns exist
            missing_cols = [col for col in groupby if col not in nodes.columns]
            if missing_cols:
                raise KeyError(f"Columns {missing_cols} not found in nodes DataFrame.")

            spikes_df = spikes_df.merge(
                nodes[groupby], left_on="node_ids", right_index=True, how="left"
            )

    return spikes_df

bmtool.analysis.spikes.compute_firing_rate_stats(df, groupby='pop_name', start_time=None, stop_time=None)

Computes the firing rates of individual nodes and the mean and standard deviation of firing rates per group.

Args: df (pd.DataFrame): Dataframe containing spike timestamps and node IDs. groupby (str or list of str, optional): Column(s) to group by (e.g., 'pop_name' or ['pop_name', 'layer']). start_time (float, optional): Start time for the analysis window. Defaults to the minimum timestamp in the data. stop_time (float, optional): Stop time for the analysis window. Defaults to the maximum timestamp in the data.

Returns: Tuple[pd.DataFrame, pd.DataFrame]: - The first DataFrame (pop_stats) contains the mean and standard deviation of firing rates per group. - The second DataFrame (individual_stats) contains the firing rate of each individual node.

Source code in bmtool/analysis/spikes.py

def compute_firing_rate_stats(
    df: pd.DataFrame,
    groupby: Union[str, List[str]] = "pop_name",
    start_time: float = None,
    stop_time: float = None,
) -> Tuple[pd.DataFrame, pd.DataFrame]:
    """
    Computes the firing rates of individual nodes and the mean and standard deviation of firing rates per group.

    Args:
        df (pd.DataFrame): Dataframe containing spike timestamps and node IDs.
        groupby (str or list of str, optional): Column(s) to group by (e.g., 'pop_name' or ['pop_name', 'layer']).
        start_time (float, optional): Start time for the analysis window. Defaults to the minimum timestamp in the data.
        stop_time (float, optional): Stop time for the analysis window. Defaults to the maximum timestamp in the data.

    Returns:
        Tuple[pd.DataFrame, pd.DataFrame]:
            - The first DataFrame (`pop_stats`) contains the mean and standard deviation of firing rates per group.
            - The second DataFrame (`individual_stats`) contains the firing rate of each individual node.
    """

    # Ensure groupby is a list
    if isinstance(groupby, str):
        groupby = [groupby]

    # Ensure all columns exist in the dataframe
    for col in groupby:
        if col not in df.columns:
            raise ValueError(f"Column '{col}' not found in dataframe.")

    # Filter dataframe based on start/stop time
    if start_time is not None:
        df = df[df["timestamps"] >= start_time]
    if stop_time is not None:
        df = df[df["timestamps"] <= stop_time]

    # Compute total duration for firing rate calculation
    if start_time is None:
        min_time = df["timestamps"].min()
    else:
        min_time = start_time

    if stop_time is None:
        max_time = df["timestamps"].max()
    else:
        max_time = stop_time

    duration = max_time - min_time  # Avoid division by zero

    if duration <= 0:
        raise ValueError("Invalid time window: Stop time must be greater than start time.")

    # Compute firing rate for each node

    # Compute spike counts per node
    spike_counts = df["node_ids"].value_counts().reset_index()
    spike_counts.columns = ["node_ids", "spike_count"]  # Rename columns

    # Merge with original dataframe to get corresponding labels (e.g., 'pop_name')
    spike_counts = spike_counts.merge(
        df[["node_ids"] + groupby].drop_duplicates(), on="node_ids", how="left"
    )

    # Compute firing rate
    spike_counts["firing_rate"] = spike_counts["spike_count"] / duration * 1000  # scale to Hz
    indivdual_stats = spike_counts

    # Compute mean and standard deviation per group
    pop_stats = spike_counts.groupby(groupby)["firing_rate"].agg(["mean", "std"]).reset_index()

    # Rename columns
    pop_stats.rename(columns={"mean": "firing_rate_mean", "std": "firing_rate_std"}, inplace=True)

    return pop_stats, indivdual_stats

bmtool.analysis.spikes._pop_spike_rate(spike_times, time=None, time_points=None, frequency=False)

Calculate the spike count or frequency histogram over specified time intervals.

Parameters:

Name	Type	Description	Default
spike_times	Union[ndarray, list]	Array or list of spike times in milliseconds	required
time	Optional[Tuple[float, float, float]]	Tuple specifying (start, stop, step) in milliseconds. Used to create evenly spaced time points if time_points is not provided. Default is None.	None
time_points	Optional[Union[ndarray, list]]	Array or list of specific time points for binning. If provided, time is ignored. Default is None.	None
frequency	bool	If True, returns spike frequency in Hz; otherwise, returns spike count. Default is False.	False

Returns:

Type	Description
ndarray	Array of spike counts or frequencies, depending on the frequency flag.

Raises:

Type	Description
ValueError	If both time and time_points are None.

Source code in bmtool/analysis/spikes.py

def _pop_spike_rate(
    spike_times: Union[np.ndarray, list],
    time: Optional[Tuple[float, float, float]] = None,
    time_points: Optional[Union[np.ndarray, list]] = None,
    frequency: bool = False,
) -> np.ndarray:
    """
    Calculate the spike count or frequency histogram over specified time intervals.

    Parameters
    ----------
    spike_times : Union[np.ndarray, list]
        Array or list of spike times in milliseconds
    time : Optional[Tuple[float, float, float]], optional
        Tuple specifying (start, stop, step) in milliseconds. Used to create evenly spaced time points
        if `time_points` is not provided. Default is None.
    time_points : Optional[Union[np.ndarray, list]], optional
        Array or list of specific time points for binning. If provided, `time` is ignored. Default is None.
    frequency : bool, optional
        If True, returns spike frequency in Hz; otherwise, returns spike count. Default is False.

    Returns
    -------
    np.ndarray
        Array of spike counts or frequencies, depending on the `frequency` flag.

    Raises
    ------
    ValueError
        If both `time` and `time_points` are None.
    """
    if time_points is None:
        if time is None:
            raise ValueError("Either `time` or `time_points` must be provided.")
        time_points = np.arange(*time)
        dt = time[2]
    else:
        time_points = np.asarray(time_points).ravel()
        dt = (time_points[-1] - time_points[0]) / (time_points.size - 1)

    bins = np.append(time_points, time_points[-1] + dt)
    spike_rate, _ = np.histogram(np.asarray(spike_times), bins)

    if frequency:
        spike_rate = 1000 / dt * spike_rate

    return spike_rate

bmtool.analysis.spikes.get_population_spike_rate(spike_data, fs=400.0, t_start=0, t_stop=None, config=None, network_name=None, save=False, save_path=None, normalize=False, smooth_window=50, smooth_method='gaussian')

Calculate the population spike rate for each population in the given spike data.

Parameters:

Name	Type	Description	Default
spike_data	DataFrame	A DataFrame containing spike data with columns 'pop_name', 'timestamps', and 'node_ids'	required
fs	float	Sampling frequency in Hz, which determines the time bin size for calculating the spike rate (default: 400.0)	400.0
t_start	float	Start time (in milliseconds) for spike rate calculation (default: 0)	0
t_stop	Optional[float]	Stop time (in milliseconds) for spike rate calculation. If None, defaults to the maximum timestamp in the data	None
config	Optional[str]	Path to a configuration file containing node information, used to determine the correct number of nodes per population. If None, node count is estimated from unique node spikes (default: None)	None
network_name	Optional[str]	Name of the network used in the configuration file, allowing selection of nodes for that network. Required if config is provided (default: None)	None
save	bool	Whether to save the calculated population spike rate to a file (default: False)	False
save_path	Optional[str]	Directory path where the file should be saved if save is True (default: None)	None
normalize	bool	Whether to normalize the spike rates for each population to a range of [0, 1] (default: False)	False
smooth_window	int	Window size for smoothing in number of time bins (default: 50)	50
smooth_method	str	Smoothing method to use: 'gaussian', 'boxcar', or 'exponential' (default: 'gaussian')	'gaussian'

Returns:

Type	Description
DataArray	An xarray DataArray containing the spike rates with dimensions of time, population, and type. The 'type' dimension includes 'raw' and 'smoothed' values. The DataArray includes sampling frequency (fs) as an attribute. If normalize is True, each population's spike rate is scaled to [0, 1].

Raises:

Type	Description
ValueError	If save is True but save_path is not provided. If an invalid smooth_method is specified.

Notes

• If config is None, the function assumes all cells in each population have fired at least once; otherwise, the node count may be inaccurate.
• If normalization is enabled, each population's spike rate is scaled using Min-Max normalization.
• Smoothing is applied using scipy.ndimage's filters based on the specified method.

Source code in bmtool/analysis/spikes.py

def get_population_spike_rate(
    spike_data: pd.DataFrame,
    fs: float = 400.0,
    t_start: float = 0,
    t_stop: Optional[float] = None,
    config: Optional[str] = None,
    network_name: Optional[str] = None,
    save: bool = False,
    save_path: Optional[str] = None,
    normalize: bool = False,
    smooth_window: int = 50,  # Window size for smoothing (in time bins)
    smooth_method: str = "gaussian",  # Smoothing method: 'gaussian', 'boxcar', or 'exponential'
) -> xr.DataArray:
    """
    Calculate the population spike rate for each population in the given spike data.

    Parameters
    ----------
    spike_data : pd.DataFrame
        A DataFrame containing spike data with columns 'pop_name', 'timestamps', and 'node_ids'
    fs : float, optional
        Sampling frequency in Hz, which determines the time bin size for calculating the spike rate (default: 400.0)
    t_start : float, optional
        Start time (in milliseconds) for spike rate calculation (default: 0)
    t_stop : Optional[float], optional
        Stop time (in milliseconds) for spike rate calculation. If None, defaults to the maximum timestamp in the data
    config : Optional[str], optional
        Path to a configuration file containing node information, used to determine the correct number of nodes per population.
        If None, node count is estimated from unique node spikes (default: None)
    network_name : Optional[str], optional
        Name of the network used in the configuration file, allowing selection of nodes for that network.
        Required if `config` is provided (default: None)
    save : bool, optional
        Whether to save the calculated population spike rate to a file (default: False)
    save_path : Optional[str], optional
        Directory path where the file should be saved if `save` is True (default: None)
    normalize : bool, optional
        Whether to normalize the spike rates for each population to a range of [0, 1] (default: False)
    smooth_window : int, optional
        Window size for smoothing in number of time bins (default: 50)
    smooth_method : str, optional
        Smoothing method to use: 'gaussian', 'boxcar', or 'exponential' (default: 'gaussian')

    Returns
    -------
    xr.DataArray
        An xarray DataArray containing the spike rates with dimensions of time, population, and type.
        The 'type' dimension includes 'raw' and 'smoothed' values.
        The DataArray includes sampling frequency (fs) as an attribute.
        If normalize is True, each population's spike rate is scaled to [0, 1].

    Raises
    ------
    ValueError
        If `save` is True but `save_path` is not provided.
        If an invalid smooth_method is specified.

    Notes
    -----
    - If `config` is None, the function assumes all cells in each population have fired at least once;
      otherwise, the node count may be inaccurate.
    - If normalization is enabled, each population's spike rate is scaled using Min-Max normalization.
    - Smoothing is applied using scipy.ndimage's filters based on the specified method.
    """
    import numpy as np
    from scipy import ndimage

    # Validate smoothing method
    if smooth_method not in ["gaussian", "boxcar", "exponential"]:
        raise ValueError(
            f"Invalid smooth_method: {smooth_method}. Choose from 'gaussian', 'boxcar', or 'exponential'."
        )

    pop_spikes = {}
    node_number = {}

    if config is None:
        print(
            "Note: Node number is obtained by counting unique node spikes in the network.\nIf the network did not run for a sufficient duration, or not all cells fired,\nthen this count will not include all nodes so the firing rate will not be of the whole population!"
        )
        print(
            "You can provide a config to calculate the correct amount of nodes! for a true population rate."
        )

    if config:
        if not network_name:
            print(
                "Grabbing first network; specify a network name to ensure correct node population is selected."
            )

    # Get t_stop if not provided
    if t_stop is None:
        t_stop = spike_data["timestamps"].max()

    # Get population names and prepare data
    populations = spike_data["pop_name"].unique()
    for pop_name in populations:
        ps = spike_data[spike_data["pop_name"] == pop_name]

        if config:
            nodes = load_nodes_from_config(config)
            if network_name:
                nodes = nodes[network_name]
            else:
                nodes = list(nodes.values())[0] if nodes else {}
            nodes = nodes[nodes["pop_name"] == pop_name]
            node_number[pop_name] = nodes.index.nunique()

        else:
            node_number[pop_name] = ps["node_ids"].nunique()

        filtered_spikes = spike_data[
            (spike_data["pop_name"] == pop_name)
            & (spike_data["timestamps"] > t_start)
            & (spike_data["timestamps"] < t_stop)
        ]
        pop_spikes[pop_name] = filtered_spikes

    # Calculate time points
    time = np.arange(t_start, t_stop, 1000 / fs)  # Convert sampling frequency to time steps

    # Calculate spike rates for each population
    spike_rates = []
    for p in populations:
        raw_rate = _pop_spike_rate(pop_spikes[p]["timestamps"], (t_start, t_stop, 1000 / fs))
        rate = fs / node_number[p] * raw_rate
        spike_rates.append(rate)

    spike_rates_array = np.array(spike_rates).T  # Transpose to have time as first dimension

    # Calculate smoothed version for each population
    smoothed_rates = []

    for i in range(spike_rates_array.shape[1]):
        pop_rate = spike_rates_array[:, i]

        if smooth_method == "gaussian":
            # Gaussian smoothing (sigma is approximately window/6 for a Gaussian filter)
            sigma = smooth_window / 6
            smoothed_pop_rate = ndimage.gaussian_filter1d(pop_rate, sigma=sigma)
        elif smooth_method == "boxcar":
            # Boxcar/uniform smoothing
            kernel = np.ones(smooth_window) / smooth_window
            smoothed_pop_rate = ndimage.convolve1d(pop_rate, kernel, mode="nearest")
        elif smooth_method == "exponential":
            # Exponential smoothing
            alpha = 2 / (smooth_window + 1)  # Equivalent to window size in exponential smoothing
            smoothed_pop_rate = np.zeros_like(pop_rate)
            smoothed_pop_rate[0] = pop_rate[0]
            for t in range(1, len(pop_rate)):
                smoothed_pop_rate[t] = alpha * pop_rate[t] + (1 - alpha) * smoothed_pop_rate[t - 1]

        smoothed_rates.append(smoothed_pop_rate)

    smoothed_rates_array = np.array(smoothed_rates).T  # Transpose to have time as first dimension

    # Stack raw and smoothed data
    combined_data = np.stack([spike_rates_array, smoothed_rates_array], axis=2)

    # Create DataArray with the additional 'type' dimension
    spike_rate_array = xr.DataArray(
        combined_data,
        coords={"time": time, "population": populations, "type": ["raw", "smoothed"]},
        dims=["time", "population", "type"],
        attrs={
            "fs": fs,
            "normalized": False,
            "smooth_method": smooth_method,
            "smooth_window": smooth_window,
        },
    )

    # Normalize if requested
    if normalize:
        # Apply normalization for each population and each type (raw/smoothed)
        for pop_idx in range(len(populations)):
            for type_idx, type_name in enumerate(["raw", "smoothed"]):
                pop_data = spike_rate_array.sel(population=populations[pop_idx], type=type_name)
                min_val = pop_data.min(dim="time")
                max_val = pop_data.max(dim="time")

                # Handle case where min == max (constant signal)
                if max_val != min_val:
                    spike_rate_array.loc[:, populations[pop_idx], type_name] = (
                        pop_data - min_val
                    ) / (max_val - min_val)

        spike_rate_array.attrs["normalized"] = True

    # Save if requested
    if save:
        if save_path is None:
            raise ValueError("save_path must be provided if save is True.")

        os.makedirs(save_path, exist_ok=True)
        save_file = os.path.join(save_path, "spike_rate.h5")
        spike_rate_array.to_netcdf(save_file)

    return spike_rate_array