Code to reproduce analysis of Claar et al (2023) neuropixel data in Burman et al. (2023)

burman_et_al_2023_neuropixel_data_analysis.py

#!/usr/bin/env python
# -*- coding: utf-8 -*

""" Extract data from Claar et al. (2023) & compute for each unit its
- evoked firing rate,
- peri-stimulus histogram entropy, and
- synchrony.

The dataset is available at:

https://dandiarchive.org/dandiset/000458/0.230317.0039

Copyright (C) 2023 by Paul Brodersen.

This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or (at
your option) any later version. This program is distributed in the
hope that it will be useful, but WITHOUT ANY WARRANTY; without even
the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
PURPOSE.  See the GNU General Public License for more details. You
should have received a copy of the GNU General Public License along
with this program. If not, see <https://www.gnu.org/licenses/.

"""

import pathlib
import numpy as np

from pynwb import NWBHDF5IO
from scipy.stats import (
    gaussian_kde,
    entropy as get_entropy,
)


def get_peristimulus_spikes(unit_to_spiketimes, trials, interval):
    peristimulus_spikes = list()
    for unit, spike_times in unit_to_spiketimes.items():
        spike_times = spike_times * 1000 # s -> ms
        for trial, stimulus_time in enumerate(trials["stop_time"].values):
            stimulus_time = stimulus_time * 1000 # s -> ms
            is_within_interval = (spike_times >= stimulus_time + interval[0]) & (spike_times < stimulus_time + interval[1])
            if np.any(is_within_interval):
                spike_times_within_interval = spike_times[is_within_interval] - stimulus_time
                unit_trial_time = np.c_[
                    np.full_like(spike_times_within_interval, unit),
                    np.full_like(spike_times_within_interval, trial),
                    spike_times_within_interval
                ]
                peristimulus_spikes.append(unit_trial_time)
    return np.concatenate(peristimulus_spikes, axis=0)


def get_firing_rate(spike_times, interval, total_trials):
    start, stop = interval
    is_valid = np.logical_and(spike_times >= start, spike_times < stop)
    total_spikes = np.sum(spike_times[is_valid])
    total_time = total_trials * (stop - start) / 1000
    return total_spikes / total_time


def get_peristimulus_histogram(spike_times, interval, bin_width, kde=False):
    if kde:
        evaluate_at = np.arange(interval[0] + bin_width/2, interval[1], bin_width)
        kernel = gaussian_kde(spike_times, bin_width)
        density = kernel(evaluate_at)
        return density / density.sum() * len(spike_times)
    else:
        bins = np.arange(interval[0], interval[1]+bin_width, bin_width)
        indices = np.digitize(spike_times, bins)
        indices -= 1 # index of zero indicates out of bounds samples
        interval_length = interval[1] - interval[0]
        counts = np.bincount(indices[indices>=0], minlength=interval_length)
        counts = counts[:interval_length] # index of len(interval) indicates out of bounds samples
        return counts.astype(float)


def get_peristimulus_histogram_entropy(spike_times, interval, bin_width):
    histogram = get_peristimulus_histogram(spike_times, interval, bin_width)
    probability = histogram / np.sum(histogram)
    return get_entropy(probability)


if __name__ == '__main__':

    data_directory = pathlib.Path("/path/to/claar_et_al_2023/data/")

    # neuropixel data sets
    datasets = [
        "sub-546655_ses-20201023_behavior+ecephys",
        "sub-551397_ses-20210211_behavior+ecephys",
        "sub-551399_ses-20210128_behavior+ecephys",
        "sub-569062_ses-20210218_behavior+ecephys",
        "sub-569064_ses-20210408_behavior+ecephys",
        "sub-569068_ses-20210304_behavior+ecephys",
        "sub-569069_ses-20210312_behavior+ecephys",
        "sub-569072_ses-20210422_behavior+ecephys",
        "sub-569073_ses-20210415_behavior+ecephys",
        "sub-569073_ses-20210416_behavior+ecephys",
        "sub-571619_ses-20210319_behavior+ecephys",
        "sub-571620_ses-20210513_behavior+ecephys",
        "sub-575102_ses-20210603_behavior+ecephys",
        "sub-586466_ses-20210729_behavior+ecephys",
        "sub-586468_ses-20210819_behavior+ecephys",
        "sub-590479_ses-20210921_behavior+ecephys",
        "sub-590480_ses-20211111_behavior+ecephys",
    ]


    for dataset in datasets:
        print()
        print("--------------------------------------------------------------------------------")
        print(dataset)
        print()

        # extract data from nwb files
        filepath = pathlib.Path(data_directory, dataset.split('_')[0], dataset + ".nwb")
        neuropixel_data = NWBHDF5IO(filepath , "r").read()

        # determine stimulated region and subset data to only include RS units from that region
        trial_annotations = neuropixel_data.intervals["trials"].to_dataframe()
        target_regions = trial_annotations[trial_annotations["stimulus_type"] == "electrical"]["estim_target_region"].unique()
        unit_data = neuropixel_data.units.to_dataframe().reset_index()
        is_within_region = unit_data["location"].apply(lambda x : x[:3] in target_regions)
        is_regular_spiking = unit_data["waveform_duration"] > 0.4
        unit_data = unit_data[is_within_region & is_regular_spiking]
        total_units = len(unit_data)

        interval = (2, 12) # ms
        resolution = 1 # ms
        for condition in ["awake", "isoflurane"]:

            # get peristimulus spikes
            trials = trial_annotations[
                (trial_annotations["stimulus_type"] == "electrical").values & \
                (trial_annotations["behavioral_epoch"] == condition).values & \
                trial_annotations["is_valid"].values
            ]
            units, trials, times = get_peristimulus_spikes(unit_data["spike_times"], trials, interval).T

            # compute evoked firing rates
            firing_rate = np.zeros((total_units))
            total_trials = int(np.max(trials) + 1)
            for ii, unit in enumerate(unit_data["id"].values):
                firing_rate[ii] = get_firing_rate(times[units==unit], interval, total_trials)
            unit_data[f"Firing rate [Hz] ({condition})"] = firing_rate

            # compute unit peristimulus histogram entropy
            entropy = np.full((total_units), np.nan)
            for ii, unit in enumerate(unit_data["id"].values):
                unit_times = times[units==unit]
                if len(unit_times) > 1:
                    entropy[ii] = get_peristimulus_histogram_entropy(unit_times, interval, resolution)
            unit_data[f"Entropy [nats] ({condition})"] = entropy

            # compute population coupling
            synchrony = {unit : list() for unit in np.unique(units)}
            for trial in np.unique(trials):
                is_trial = trials == trial
                if np.any(is_trial):
                    population_response = get_peristimulus_histogram(times[is_trial], interval, resolution, kde=False)
                    population_response /= population_response.sum()
                    for unit in np.unique(units[is_trial]):
                        synchrony[unit].extend([population_response[int(time - interval[0])] \
                                                for time in times[np.logical_and(is_trial, units == unit)]])
            synchrony = np.array([np.mean(synchrony[unit]) if unit in synchrony else np.nan for unit in unit_data["id"].values])
            unit_data[f"Synchrony ({condition})"] = synchrony

        # save out unit data
        unit_data["Dataset"] = dataset
        unit_data.to_csv(data_directory / f"{dataset}.csv", index=False)