Detecting Spikes in Continuous Recordings

Import Libraries

import numpy as np
import pandas as pd
from matplotlib import pyplot as plt
from scipy import signal

Utility Functions

Define the utility functions required for this session.

class utils:
    @staticmethod
    def plot_frequency_response(b, a, cutoff):
        w, h = signal.freqz(b, a, fs=30000)
        plt.semilogx(w, np.abs(h), label="Filter")
        if not isinstance(cutoff, list):
            cutoff = [cutoff]
        for c in cutoff:
            plt.axvline(c, color="red", label="Cutoff")
        plt.xlabel("Frequency [Hz]")
        plt.ylabel("Gain [a.u.]")
        plt.legend()

Download Data

Download the data required for this session

import requests
url = "https://uni-bonn.sciebo.de/public.php/webdav/oLE33JLS7GcERgG"
print("Downloading Data ...")
response = requests.get(("/").join(url.split("/")[:-1]), auth=(url.split("/")[-1], "ibots"))
response.raise_for_status()
with open("flash_spikes.parquet", "wb") as file:
    file.write(response.content)
print("Done!")

Downloading Data ...
Done!

In this notebook, our aim is to extract spike waveforms from raw electrophysiological recordings from a tetrode (a probe with four separate electrodes). In this section we are going to load the data into a pandas data frame and explore its contents.

Exercise: Load the file "raw_tetrode_recording.parquet" into a pandas data frame and print the first 10 rows

df = pd.read_parquet("raw_tetrode_recording.parquet")
df.head(10)

Exercise: what is the duration of the recording and how many samples (i.e. rows) are there? What is the sampling rate of the recording?

sampling_rate = len(df) / df.index.max(0)
sampling_rate

30000.015104174978

Example: Plot the data recorded at "Ch1" between tmin=0 and tmax=3.

tmin, tmax = 0, 3
data = df["Ch1"].loc[tmin:tmax]
plt.plot(data.index, data)

Exercise: Plot the data recorded at "Ch2" between tmin=10 and tmax=20.

Solution

tmin, tmax = 10, 20
data = df["Ch2"].loc[tmin:tmax]
plt.plot(data.index, data)

Exercise: Compute the difference between the data recorded at channels "Ch1" and "Ch2" between tmin=10 and tmax=20, and plot it.

Solution

tmin, tmax = 10, 20
data = df["Ch1"].loc[tmin:tmax] - df["Ch2"].loc[tmin:tmax]
plt.plot(data.index, data)
plt.ylabel("Voltage [muV]")

Exercise: Plot the data recorded at channels "Ch1" and "Ch3" between tmin=0.05 and tmax=0.06 in the same plot.

Solution

tmin, tmax = 0.05, 0.06
for ch in ["Ch1", "Ch3"]:
    data = df[ch].loc[tmin:tmax]
    plt.plot(data.index, data, label=ch)
plt.legend()
plt.xlabel("Time [s]")
plt.ylabel("Voltage [muV]")

Extracellular electrophysiological recordings contain raw voltage traces from every electrode. These recordings contain the spikes we want to analyze but they also contain other signals like local field potentials (LFPs) and artifacts. To extract the spikes we have to filter out the part of the signal that contains the spikes. In this section, we will explore how to construct infinite-impulse response (IIR) filter and apply them to the data to extract specific parts of the frequency spectrum.

Exercises

Example: Create a Butterworth "lowpass" filter with order=4, cutoff=500 and fs=30000 and plot its frequency response.

cutoff = 500
order = 4
b, a = signal.butter(order, cutoff, btype="lowpass", fs=30000)
utils.plot_frequency_response(b, a, cutoff=cutoff)

Example: Use signal.filtfilt() to apply the filter coefficients b and a to the data recorded at "Ch1" between tmin=0 and tmax=3 and plot the raw and filtered signal in a 2-by-1 subplot grid.

tmin, tmax = 0, 3
data = df["Ch1"].loc[tmin:tmax]
data_filt = signal.filtfilt(b, a, data)

plt.subplot(2, 1, 1)
plt.plot(data.index, data)
plt.subplot(2, 1, 2)
plt.plot(data.index, data_filt)

Exercise: Create a Butterworth "lowpass" filter with order=4, cutoff=300 and fs=30000 and plot its frequency response.

Solution

cutoff = 300
order = 4
b, a = signal.butter(order, cutoff, btype="lowpass", fs=30000)
utils.plot_frequency_response(b, a, cutoff=cutoff)

Exercise: Use signal.filtfilt() to apply the filter coefficients b and a to the data recorded at "Ch2" between tmin=0 and tmax=3 and plot the raw and filtered signal in a 2-by-1 subplot grid.

Solution

tmin, tmax = 0, 3
data = df["Ch2"].loc[tmin:tmax]
data_filt = signal.filtfilt(b, a, data)

plt.subplot(2, 1, 1)
plt.plot(data.index, data)
plt.subplot(2, 1, 2)
plt.plot(data.index, data_filt)

Exercise: Create a Butterworth "highpass" filter with order=6, cutoff=1000 and fs=30000 and plot its frequency response.

Solution

cutoff = 1000
order = 6
b, a = signal.butter(order, cutoff, btype="highpass", fs=30000)
utils.plot_frequency_response(b, a, cutoff=cutoff)

Exercise: Use signal.filtfilt() to apply the filter coefficients b and a to the data recorded at "Ch3" between tmin=15 and tmax=20 and plot the raw and filtered signal in a 2-by-1 subplot grid.

Solution

tmin, tmax = 15, 20
data = df["Ch3"].loc[tmin:tmax]
data_filt = signal.filtfilt(b, a, data)

plt.subplot(2, 1, 1)
plt.plot(data.index, data)
plt.subplot(2, 1, 2)
plt.plot(data.index, data_filt)
plt.tight_layout()

Exercise: Create a Butterworth "bandpass" filter with order=6, cutoff=[500, 5000] and and fs=30000 and plot its frequency response.

Solution

cutoff = [500, 5000]
order = 6
b, a = signal.butter(order, cutoff, btype="bandpass", fs=30000)
utils.plot_frequency_response(b, a, cutoff=cutoff)

Exercise: Use signal.filtfilt() to apply the filter coefficients b and a to the data recorded at "Ch4" between tmin=50 and tmax=60 and plot the raw and filtered signal in a 2-by-1 subplot grid.

Solution

tmin, tmax = 50, 60
data = df["Ch4"].loc[tmin:tmax]
data_filt = signal.filtfilt(b, a, data)

plt.subplot(2, 1, 1)
plt.plot(data.index, data)
plt.subplot(2, 1, 2)
plt.plot(data.index, data_filt)

Exercise: Use signal.filtfilt() to apply the filter coefficients b and a to the whole data frame df across axis=0 and return the filtered data. Then, use df.loc[:]=data to overwrite df with the filtered data.

data = signal.filtfilt(b, a, df, axis=0)
df.loc[:] = data

To detect spikes, we need to estimate the level of noise in the data and define a threshold (i.e. 5 times the noise level) to categorize signals that cross this threshold as spikes. The noise level is calculated as the median absolute deviation (MAD) of the signal:

This is a useful measure because the median is more robust to large deviations compared to the mean. Because an action potential is a depolarization (i.e. a positive spike) in the cell, it is reflected as a negative spike in the extracellular recording. In this section, we are going to explore how to compare the signal against a noise threshold to identify the times at which spikes occur.

Exercises

Example: Compute the robust noise level sigma for first channel as the median absolute deviation (MAD) of the data recorded at "Ch1".

x = df["Ch1"]
sigma = np.median(np.abs(x - x.median(axis=0)), axis=0)
sigma

23.840796979425427

Exercise: Compute the robust noise level sigma for the data recorded at "Ch2". Does this mean the signal recorded at "Ch2" is more or less noisy than the signal at "Ch1"?

x = df["Ch2"]
sigma = np.median(np.abs(x - x.median(axis=0)), axis=0)
sigma

25.290102706075356

Example: Find which samples in df are below the threshold of b * sigma for b=2 and assign them to a variable is_below_thresh. Convert that variable to integers and print the total number of samples below the threshold.

b = 2
threshold = -1 * (np.mean(df, axis=0) + b * sigma)
is_below_thresh = (df < threshold).any(axis=1)
is_below_thresh = is_below_thresh.astype(int)
sum(is_below_thresh)

320593

Example: Get the time points where difference between the values of is_below_threshold is 1, i.e. where the signal crosses from above to below the threshold and print the number of crossing_times.

is_crossing = is_below_thresh.diff() == 1
crossing_times = df[is_crossing].index
len(crossing_times)

85050

Exercise: Set b=4 and find which samples in df are below the threshold of b * sigma for b=4 and assign them to a variable is_below_thresh. Convert that variable to integers and print the total number of samples below the threshold.

b = 4
threshold = -1 * (np.mean(df, axis=0) + b * sigma)
is_below_thresh = (df < threshold).any(axis=1)
is_below_thresh = is_below_thresh.astype(int)
sum(is_below_thresh)

43714

Exercise: Get the time points where difference between the values of is_below_threshold is 1, i.e. where the signal crosses from above to below the threshold and print the number of crossing_times.

is_crossing = is_below_thresh.diff() == 1
crossing_times = df[is_crossing].index
len(crossing_times)

12315

Exercise: Set b=5, get the times where the signal crosses the b*sigma threshold and print the total number of crossing_times.

b = 5
threshold = -1 * (np.mean(df, axis=0) + b * sigma)
is_below_thresh = (df < threshold).any(axis=1)
is_below_thresh = is_below_thresh.astype(int)
is_crossing = is_below_thresh.diff() == 1
crossing_times = df[is_crossing].index
len(crossing_times)

5062

Example: For the first crossing time, get the lowest point (i.e. the negative peak) within a window=2*1e-3 seconds and get the channel voltages at this peak_time.

window = 2 * 1e-3
tmin = crossing_times[0]
tmax = tmin + window
peak_time = df.loc[tmin:tmax].idxmin().min()
df.loc[peak_time]

Ch1    -80.497902
Ch2    -82.152938
Ch3   -143.714896
Ch4    -85.106740
Name: 0.011866667284722287, dtype: float64

Exercise: Get the peak time for the 10th crossing time and get the channel voltages at the peak time. Which channel recorded the largest signal for this spike?

window = 2 * 1e-3
tmin = crossing_times[9]
tmax = tmin + window
peak_time = df.loc[tmin:tmax].idxmin().min()
df.loc[peak_time]

Ch1   -117.819185
Ch2    -88.697246
Ch3   -182.508835
Ch4    -93.611148
Name: 0.06706667015972258, dtype: float64

Exercise: Complete the loop below to get the peak_time for every crossing time c_t and append it to the list spike_times.

window = 2 * 1e-3
spike_times = []
for c_t in crossing_times:
    tmin = # your code here
    tmax = # your code here
    peak_time= # your code here
    spike_times.append(peak_time)
len(crossing_times), len(spike_times)

window = 2 * 1e-3
spike_times = []
for c_t in crossing_times:
    tmin = c_t
    tmax = tmin + window
    peak_time = df.loc[tmin:tmax].idxmin().min()
    spike_times.append(peak_time)
len(crossing_times), len(spike_times)

(5062, 5062)

Example: Plot the recording at "Ch1" and mark the spike_times with red dots. Then, set the plt.xlim() so the plot shows the time between 0.03 and 0.04 seconds.

ch = "Ch1"
plt.plot(df[ch].index, df[ch])
plt.plot(df[ch].loc[spike_times].index, df[ch].loc[spike_times], "r.")
plt.xlim(0.03, 0.04)

(0.03, 0.04)

Exercise: Plot the recording at "Ch4" and mark the spike_times with red dots. Then, set the plt.xlim() so the plot shows the time between 0.03 and 0.04 seconds.

Solution

ch = "Ch4"
plt.plot(df[ch].index, df[ch])
plt.plot(df[ch].loc[spike_times].index, df[ch].loc[spike_times], "r.")
plt.xlim(0.03, 0.04)

(0.03, 0.04)

Exercise: Fill the plt.plot() commands in the loop below to plot the recording at each channel ch and mark the spike_times with red dots.

for i, ch in enumerate(df.columns):
    plt.subplot(4, 1, i + 1)
    plt.plot() # your code here
    plt.plot() # your code here
    plt.xlim(0.03, 0.04)

Solution

for i, ch in enumerate(df.columns):
    plt.subplot(4, 1, i + 1)
    plt.plot(df[ch].index, df[ch])
    plt.plot(df[ch].loc[spike_times].index, df[ch].loc[spike_times], "r.")
    plt.xlim(0.03, 0.04)

Exercise: Plot the recording at "Ch3" and mark the spike_times with red dots. Then, set the plt.xlim() so the plot shows the time between 1 and 1.5 seconds.

Solution

ch = "Ch3"
plt.plot(df[ch].index, df[ch])
plt.plot(df[ch].loc[spike_times].index, df[ch].loc[spike_times], "r.")
plt.xlim(1, 1.5)

(1.0, 1.5)

Now that we have the times at which spikes occur, we can extract the spike waveforms that we will use for spike sorting from the signal. Those spikes will be stored in a 3D matrix where the first dimension represents the spikes, the second represents the points in time and the third dimension represents the channels. By visualizing the waveforms across channels, we will see the differences in spike shape, amplitude and arrival times that we can use to sort the spikes later on.

Exercises

Example: Get the waveform of the 1st spike by extracting the signal recorded within +/- 2*-1e3 seconds around the spike time.

window = 2 * 1e-3
tmin = spike_times[0] - window
tmax = spike_times[0] + window
waveform = df.loc[tmin:tmax]

Example: Plot the extracted waveform for every channel and add a plt.legend().

plt.plot(waveform.index, waveform, label=waveform.columns)
plt.legend()

Exercise: Get the waveform of the 1st spike by extracting the signal recorded within +/- 0.5*-1e3 seconds around the spike time and plot the extracted spike for every channel.

Solution

window = 0.5 * 1e-3
tmin = spike_times[0] - window
tmax = spike_times[0] + window
waveform = df.loc[tmin:tmax]
plt.plot(waveform.index, waveform, label=waveform.columns)
plt.legend()

Exercise: Get the waveform of the 2nd spike by extracting the signal recorded within +/- 0.5*-1e3 seconds around the spike time and plot the extracted spike for every channel.

Solution

window = 0.5 * 1e-3
tmin = spike_times[1] - window
tmax = spike_times[1] + window
waveform = df.loc[tmin:tmax]
plt.plot(waveform.index, waveform, label=waveform.columns)
plt.legend()

Exercise: Complete the loop below by extracting the waveform for every item s in the list of spike_times.

window = 0.5 * 1e-3
waveforms = []
for s in spike_times:
    waveform = # your code here
    waveforms.append(waveform.to_numpy())
waveforms = np.stack(waveforms)
waveforms.shape

window = 0.5 * 1e-3
waveforms = []
for s in spike_times:
    waveform = df.loc[s - window : s + window]
    waveforms.append(waveform.to_numpy())
waveforms = np.stack(waveforms)
waveforms.shape

(5062, 29, 4)

Example: Plot the first 5 waveforms for each channel.

for i in range(4):
    plt.subplot(2, 2, i + 1)
    plt.plot(waveforms[:5, :, i].T, linewidth=0.5)
plt.tight_layout()

Exercise: Plot the last 5 waveforms for each channel.

Solution

for i in range(4):
    plt.subplot(2, 2, i + 1)
    plt.plot(waveforms[-5:, :, i].T, linewidth=0.5)
plt.tight_layout()

Exercise: Compute the mean waveform for every channel (i.e. waveforms.mean(axis=0)) and plot them in the same graph.

Solution

plt.plot(waveforms.mean(axis=0), linewidth=0.5, label=df.columns)
plt.xlabel("Time [ms]")
plt.ylabel("Voltage [muV]")
plt.legend()